Mercurial > repos > public > sbplib_julia
comparison src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 881:aa4875f9a530 feature/variable_derivatives
Start implementing the variable second derivative
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 20 Jan 2022 15:18:14 +0100 |
| parents | |
| children | 9098fc936776 |
comparison
equal
deleted
inserted
replaced
| 880:f74189c954d6 | 881:aa4875f9a530 |
|---|---|
| 1 export SecondDerivativeVariable | |
| 2 | |
| 3 # """ | |
| 4 # SecondDerivativeVariable(grid, inner_stencil, closure_stencils, parity, direction) | |
| 5 | |
| 6 # Creates a volume operator on a `Dim`-dimensional grid acting along the | |
| 7 # specified coordinate `direction`. The action of the operator is determined by | |
| 8 # the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the | |
| 9 # corresponding `SecondDerivativeVariable` tensor mapping is returned. When `Dim>1`, the | |
| 10 # returned operator is the appropriate outer product of a one-dimensional | |
| 11 # operators and `IdentityMapping`s, e.g for `Dim=3` the volume operator in the | |
| 12 # y-direction is `I⊗op⊗I`. | |
| 13 # """ | |
| 14 # function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction) | |
| 15 # #TODO: Check that direction <= Dim? | |
| 16 | |
| 17 # # Create 1D volume operator in along coordinate direction | |
| 18 # op = SecondDerivativeVariable(restrict(grid, direction), inner_stencil, closure_stencils, parity) | |
| 19 # # Create 1D IdentityMappings for each coordinate direction | |
| 20 # one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) | |
| 21 # Is = IdentityMapping{eltype(grid)}.(size.(one_d_grids)) | |
| 22 # # Formulate the correct outer product sequence of the identity mappings and | |
| 23 # # the volume operator | |
| 24 # parts = Base.setindex(Is, op, direction) | |
| 25 # return foldl(⊗, parts) | |
| 26 # end | |
| 27 | |
| 28 """ | |
| 29 SecondDerivativeVariable{T,N,M,K} <: TensorOperator{T,1} | |
| 30 Implements a one-dimensional constant coefficients volume operator | |
| 31 """ | |
| 32 struct SecondDerivativeVariable{T,N,M,K} <: TensorMapping{T,1,1} | |
| 33 inner_stencil::NestedStencil{T,N} | |
| 34 closure_stencils::NTuple{M,NestedStencil{T,K}} | |
| 35 size::NTuple{1,Int} | |
| 36 end | |
| 37 | |
| 38 function SecondDerivativeVariable(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) | |
| 39 return SecondDerivativeVariable(inner_stencil, Tuple(closure_stencils), size(grid)) | |
| 40 end | |
| 41 | |
| 42 closure_size(::SecondDerivativeVariable{T,N,M}) where {T,N,M} = M | |
| 43 | |
| 44 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size | |
| 45 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size | |
| 46 | |
| 47 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Lower}) where T | |
| 48 return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i)) | |
| 49 end | |
| 50 | |
| 51 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Interior}) where T | |
| 52 return apply_stencil(op.inner_stencil, v, Int(i)) | |
| 53 end | |
| 54 | |
| 55 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Upper}) where T | |
| 56 return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i)) | |
| 57 end | |
| 58 | |
| 59 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i) where T | |
| 60 r = getregion(i, closure_size(op), op.size[1]) | |
| 61 return LazyTensors.apply(op, v, Index(i, r)) | |
| 62 end |