Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/volume_operator.jl @ 875:067a322e4f73 laplace_benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/laplace_opset
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jan 2022 10:55:08 +0100
parents ae28f1d7ef5e
children b41180efb6c2 469ed954b493 0a856fb96db4
comparison
equal deleted inserted replaced
874:7e9ebd572deb 875:067a322e4f73
1 """ 1 """
2 volume_operator(grid,inner_stencil,closure_stencils,parity,direction) 2 volume_operator(grid, inner_stencil, closure_stencils, parity, direction)
3
3 Creates a volume operator on a `Dim`-dimensional grid acting along the 4 Creates a volume operator on a `Dim`-dimensional grid acting along the
4 specified coordinate `direction`. The action of the operator is determined by the 5 specified coordinate `direction`. The action of the operator is determined by
5 stencils `inner_stencil` and `closure_stencils`. 6 the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the
6 When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned. 7 corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the
7 When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product 8 returned operator is the appropriate outer product of a one-dimensional
8 of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`, 9 operators and `IdentityMapping`s, e.g for `Dim=3` the volume operator in the
9 the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`. 10 y-direction is `I⊗op⊗I`.
10 """ 11 """
11 function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M} 12 function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction)
12 #TODO: Check that direction <= Dim? 13 #TODO: Check that direction <= Dim?
13 14
14 # Create 1D volume operator in along coordinate direction 15 # Create 1D volume operator in along coordinate direction
15 op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) 16 op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity)
16 # Create 1D IdentityMappings for each coordinate direction 17 # Create 1D IdentityMappings for each coordinate direction
17 one_d_grids = restrict.(Ref(grid), Tuple(1:Dim)) 18 one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
18 Is = IdentityMapping{T}.(size.(one_d_grids)) 19 Is = IdentityMapping{eltype(grid)}.(size.(one_d_grids))
19 # Formulate the correct outer product sequence of the identity mappings and 20 # Formulate the correct outer product sequence of the identity mappings and
20 # the volume operator 21 # the volume operator
21 parts = Base.setindex(Is, op, direction) 22 parts = Base.setindex(Is, op, direction)
22 return foldl(⊗, parts) 23 return foldl(⊗, parts)
23 end 24 end
32 size::NTuple{1,Int} 33 size::NTuple{1,Int}
33 parity::Parity 34 parity::Parity
34 end 35 end
35 36
36 function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity) 37 function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity)
37 return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity) 38 return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity)
38 end 39 end
39 40
40 closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M 41 closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M
41 42
42 LazyTensors.range_size(op::VolumeOperator) = op.size 43 LazyTensors.range_size(op::VolumeOperator) = op.size