Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 873:9929c99754fb laplace_benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add some benchmarks using the Laplace Operator Set
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 19 Jan 2022 13:15:45 +0100
parents 1970ebceabe4
children 067a322e4f73
comparison
equal deleted inserted replaced
865:545a6c1a0a0e 873:9929c99754fb
20 H::TensorMapping # Inner product operator 20 H::TensorMapping # Inner product operator
21 H_inv::TensorMapping # Inverse inner product operator 21 H_inv::TensorMapping # Inverse inner product operator
22 e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators. 22 e::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary restriction operators.
23 d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators 23 d::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Normal derivative operators
24 H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product? 24 H_boundary::StaticDict{<:BoundaryIdentifier,<:TensorMapping} # Boundary quadrature operators # TODO: Boundary inner product?
25 order::Int
25 end 26 end
26 export Laplace 27 export Laplace
27 28
28 function Laplace(grid::AbstractGrid, fn; order) 29 function Laplace(grid::AbstractGrid, fn; order)
29 # TODO: Removed once we can construct the volume and 30 # TODO: Removed once we can construct the volume and
46 n_ids = length(ids) 47 n_ids = length(ids)
47 e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid,e_closure_stencil,ids[i]),n_ids) 48 e_pairs = ntuple(i -> ids[i] => boundary_restriction(grid,e_closure_stencil,ids[i]),n_ids)
48 d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids) 49 d_pairs = ntuple(i -> ids[i] => normal_derivative(grid,d_closure_stencil,ids[i]),n_ids)
49 Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids) 50 Hᵧ_pairs = ntuple(i -> ids[i] => inner_product(boundary_grid(grid,ids[i]),H_closure_stencils),n_ids)
50 51
51 return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs)) 52 return Laplace(Δ, H, H⁻¹, StaticDict(e_pairs), StaticDict(d_pairs), StaticDict(Hᵧ_pairs), order)
52 end 53 end
53 54
54 # TODO: Consider pretty printing of the following form 55 # TODO: Consider pretty printing of the following form
55 # Base.show(io::IO, L::Laplace{T, Dim}) where {T,Dim,TM} = print(io, "Laplace{$T, $Dim, $TM}(", L.D, L.H, L.H_inv, L.e, L.d, L.H_boundary, ")") 56 # Base.show(io::IO, L::Laplace{T, Dim}) where {T,Dim,TM} = print(io, "Laplace{$T, $Dim, $TM}(", L.D, L.H, L.H_inv, L.e, L.d, L.H_boundary, ")")
56 57
57 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D) 58 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
58 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D) 59 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
59 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...) 60 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
60 61
62 """
63 closure_size(L::Lapalace)
64
65 Returns the inner product operator associated with `L`
66
67 """
68 closure_size(L::Laplace) = closure_size(L.D)
69 export closure_size
61 70
62 """ 71 """
63 inner_product(L::Lapalace) 72 inner_product(L::Lapalace)
64 73
65 Returns the inner product operator associated with `L` 74 Returns the inner product operator associated with `L`
122 boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id] 131 boundary_quadrature(L::Laplace,id::BoundaryIdentifier) = L.H_boundary[id]
123 boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N) 132 boundary_quadrature(L::Laplace,ids::NTuple{N,BoundaryIdentifier}) where N = ntuple(i->L.H_boundary[ids[i]],N)
124 boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N) 133 boundary_quadrature(L::Laplace,ids::Vararg{BoundaryIdentifier,N}) where N = ntuple(i->L.H_boundary[ids[i]],N)
125 export boundary_quadrature 134 export boundary_quadrature
126 135
136 abstract type BoundaryConditionType end
137 struct NeumannBC <: BoundaryConditionType end
138 struct DirichletBC <: BoundaryConditionType end
139 export NeumannBC
140
141 boundary_condition(L, id, ::NeumannBC) = neumann_bc(L, id)
142 boundary_condition(L, id, ::DirichletBC) = dirichlet_bc(L, id)
143 export boundary_condition
144
145 function neumann_bc(L::Laplace, id::BoundaryIdentifier)
146 H_inv = inverse_inner_product(L)
147 e = boundary_restriction(L,id)
148 d = normal_derivative(L,id)
149 H_b = boundary_quadrature(L,id)
150 tau = H_inv∘e'∘H_b
151 closure = tau∘d
152 # TODO: Return penalty once we have implemented scalar scaling of the operators.
153 return closure
154 end
155
156 function dirichlet_bc(L::Laplace, id::BoundaryIdentifier)
157 error("Not implemented")
158 # H_inv = inverse_inner_product(L)
159 # e = boundary_restriction(L,id)
160 # d = normal_derivative(L,id)
161 # H_b = boundary_quadrature(L,id)
162 # gamma = borrowing_parameter(L)
163 # tuning = 1.2
164 # S = ScalingOperator(tuning * -1.0/gamma)
165 # tau = H_inv∘(S∘e' + d')∘H_b
166 # closure = tau∘e
167 # penalty = ScalingOperator(-1)∘tau
168 # return (closure, penalty)
169 end
170
171 # function borrowing_parameter(L)
172 # if L.order == 2
173 # return 0.4
174 # elseif L.order == 4
175 # return 0.2508
176 # elseif L.order == 6
177 # return 0.1878
178 # elseif L.order == 8
179 # return 0.0015
180 # elseif L.order == 10
181 # return 0.0351
182 # end
183 # end
127 184
128 """ 185 """
129 laplace(grid, inner_stencil, closure_stencils) 186 laplace(grid, inner_stencil, closure_stencils)
130 187
131 Creates the Laplace operator operator `Δ` as a `TensorMapping` 188 Creates the Laplace operator operator `Δ` as a `TensorMapping`