Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 1395:bdcdbd4ea9cd feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default. Comment out broken tests for boundary_conditions at sat
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 26 Jul 2023 21:35:50 +0200
parents 470a70a6c1e6 08f06bfacd5c
children 35840a0681d1
comparison
equal deleted inserted replaced
1217:ea2e8254820a 1395:bdcdbd4ea9cd
1 """ 1 """
2 Laplace{T, Dim, TM} <: LazyTensor{T, Dim, Dim} 2 Laplace{T, Dim, TM} <: LazyTensor{T, Dim, Dim}
3 3
4 Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a 4 The Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a
5 `LazyTensor`. Additionally `Laplace` stores the `StencilSet` 5 `LazyTensor`.
6 used to construct the `LazyTensor `.
7 """ 6 """
8 struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim} 7 struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim}
9 D::TM # Difference operator 8 D::TM # Difference operator
10 stencil_set::StencilSet # Stencil set of the operator 9 stencil_set::StencilSet # Stencil set of the operator
11 end 10 end
12 11
13 """ 12 """
14 Laplace(grid::Equidistant, stencil_set) 13 Laplace(g::Grid, stencil_set::StencilSet)
15 14
16 Creates the `Laplace` operator `Δ` on `grid` given a `stencil_set`. 15 Creates the `Laplace` operator `Δ` on `g` given `stencil_set`.
17 16
18 See also [`laplace`](@ref). 17 See also [`laplace`](@ref).
19 """ 18 """
20 function Laplace(grid::EquidistantGrid, stencil_set::StencilSet) 19 function Laplace(g::Grid, stencil_set::StencilSet)
21 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 20 Δ = laplace(g, stencil_set)
22 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) 21 return Laplace(Δ, stencil_set)
23 Δ = laplace(grid, inner_stencil,closure_stencils)
24 return Laplace(Δ,stencil_set)
25 end 22 end
26 23
27 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D) 24 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
28 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D) 25 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)
29 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...) 26 LazyTensors.apply(L::Laplace, v::AbstractArray, I...) = LazyTensors.apply(L.D,v,I...)
30 27
31 # TODO: Implement pretty printing of Laplace once pretty printing of LazyTensors is implemented. 28 # TODO: Implement pretty printing of Laplace once pretty printing of LazyTensors is implemented.
32 # Base.show(io::IO, L::Laplace) = ... 29 # Base.show(io::IO, L::Laplace) = ...
33 30
34 """ 31 """
35 laplace(grid::EquidistantGrid, inner_stencil, closure_stencils) 32 laplace(g::Grid, stencil_set)
36 33
37 Creates the Laplace operator operator `Δ` as a `LazyTensor` 34 Creates the Laplace operator operator `Δ` as a `LazyTensor` on `g`.
38 35
39 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using 36 `Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `g`. The
40 the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` 37 approximation depends on the type of grid and the stencil set.
41 for the points in the closure regions.
42
43 On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a
44 multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s
45 where the sum is carried out lazily.
46 38
47 See also: [`second_derivative`](@ref). 39 See also: [`second_derivative`](@ref).
48 """ 40 """
49 function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils) 41 function laplace end
50 Δ = second_derivative(grid, inner_stencil, closure_stencils, 1) 42 function laplace(g::TensorGrid, stencil_set)
51 for d = 2:ndims(grid) 43 # return mapreduce(+, enumerate(g.grids)) do (i, gᵢ)
52 Δ += second_derivative(grid, inner_stencil, closure_stencils, d) 44 # Δᵢ = laplace(gᵢ, stencil_set)
45 # LazyTensors.inflate(Δᵢ, size(g), i)
46 # end
47
48 Δ = LazyTensors.inflate(laplace(g.grids[1], stencil_set), size(g), 1)
49 for d = 2:ndims(g)
50 Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d)
53 end 51 end
54 return Δ 52 return Δ
55 end 53 end
54 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
55
56 56
57 """ 57 """
58 sat_tensors(Δ::Laplace, g::EquidistantGrid, bc::NeumannCondition) 58 sat_tensors(Δ::Laplace, g::TensorGrid, bc::NeumannCondition)
59 59
60 Returns anonymous functions for construction the `LazyTensorApplication`s 60 Returns anonymous functions for construction the `LazyTensorApplication`s
61 recuired in order to impose a Neumann boundary condition. 61 recuired in order to impose a Neumann boundary condition.
62 62
63 See also: [`sat`,`NeumannCondition`](@ref). 63 See also: [`sat`,`NeumannCondition`](@ref).
64 """ 64 """
65 function BoundaryConditions.sat_tensors(Δ::Laplace, g::EquidistantGrid, bc::NeumannCondition) 65 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
66 id = bc.id 66 id = bc.id
67 set = Δ.stencil_set 67 set = Δ.stencil_set
68 H⁻¹ = inverse_inner_product(g,set) 68 H⁻¹ = inverse_inner_product(g,set)
69 Hᵧ = inner_product(boundary_grid(g, id), set) 69 Hᵧ = inner_product(boundary_grid(g, id), set)
70 e = boundary_restriction(g, set, id) 70 e = boundary_restriction(g, set, id)
71 d = normal_derivative(g, set, id) 71 d = normal_derivative(g, set, id)
72 72
73 closure(u) = H⁻¹*e'*Hᵧ*d*u 73 closure(u) = H⁻¹*e'*Hᵧ*d*u
74 penalty(g) = -H⁻¹*e'*Hᵧ*g 74 penalty(g) = -H⁻¹*e'*Hᵧ*g
75 return closure, penalty 75 return closure, penalty
76 end 76 end