Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/laplace/secondderivative.jl @ 436:cffdac9c612d bugfix/tensor_application_multiplication

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Close branch before merge
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 19 Oct 2020 20:54:23 +0200
parents dacbcba33d7d
children 011ca1639153
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
1 """
293
f63232aeb1c6 Move laplace.jl from DiffOps to SbpOperators. Rename constandlaplace to secondderivative
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 292
diff changeset
2 SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1}
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
3 Implements the Laplace tensor operator `L` with constant grid spacing and coefficients
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
4 in 1D dimension
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
5 """
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
6
348
7fe43d902a27 Start trying to change LazyTensors
Jonatan Werpers <jonatan@werpers.com>
parents: 333
diff changeset
7 struct SecondDerivative{T,N,M,K} <: TensorMapping{T,1,1}
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
8 h_inv::T # The grid spacing could be included in the stencil already. Preferable?
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
9 innerStencil::Stencil{T,N}
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
10 closureStencils::NTuple{M,Stencil{T,K}}
356
0844069ab5ff Reinclude SbpOperators and fix most of the code and tests there.
Jonatan Werpers <jonatan@werpers.com>
parents: 348
diff changeset
11 size::NTuple{1,Int}
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
12 end
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
13 export SecondDerivative
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
14
356
0844069ab5ff Reinclude SbpOperators and fix most of the code and tests there.
Jonatan Werpers <jonatan@werpers.com>
parents: 348
diff changeset
15 function SecondDerivative(grid::EquidistantGrid{1}, innerStencil, closureStencils)
381
dacbcba33d7d Refactor EquidistantGrid to not store spacing or inverse spacing
Jonatan Werpers <jonatan@werpers.com>
parents: 380
diff changeset
16 h_inv = inverse_spacing(grid)[1]
380
81053b1992b6 Remove parity field from Secondderivative
Jonatan Werpers <jonatan@werpers.com>
parents: 361
diff changeset
17 return SecondDerivative(h_inv, innerStencil, closureStencils, size(grid))
356
0844069ab5ff Reinclude SbpOperators and fix most of the code and tests there.
Jonatan Werpers <jonatan@werpers.com>
parents: 348
diff changeset
18 end
0844069ab5ff Reinclude SbpOperators and fix most of the code and tests there.
Jonatan Werpers <jonatan@werpers.com>
parents: 348
diff changeset
19
0844069ab5ff Reinclude SbpOperators and fix most of the code and tests there.
Jonatan Werpers <jonatan@werpers.com>
parents: 348
diff changeset
20 LazyTensors.range_size(D2::SecondDerivative) = D2.size
0844069ab5ff Reinclude SbpOperators and fix most of the code and tests there.
Jonatan Werpers <jonatan@werpers.com>
parents: 348
diff changeset
21 LazyTensors.domain_size(D2::SecondDerivative) = D2.size
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
22
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
23 #TODO: The 1D tensor mappings should not have to dispatch on 1D tuples if we write LazyTensor.apply for vararg right?!?!
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
24 # Currently have to index the Tuple{Index} in each method in order to call the stencil methods which is ugly.
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
25 # I thought I::Vararg{Index,R} fell back to just Index for R = 1
299
27a0bca5e1f2 Add apply_transpose and fix minor issues
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 293
diff changeset
26
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
27 # Apply for different regions Lower/Interior/Upper or Unknown region
328
9cc5d1498b2d Refactor 1D diagonal inner product in quadrature.jl to separate file. Write tests for quadratures. Clean up laplace and secondderivative
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 311
diff changeset
28 function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Lower}) where T
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
29 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(I)], v, Int(I))
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
30 end
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
31
328
9cc5d1498b2d Refactor 1D diagonal inner product in quadrature.jl to separate file. Write tests for quadratures. Clean up laplace and secondderivative
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 311
diff changeset
32 function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Interior}) where T
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
33 return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(I))
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
34 end
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
35
328
9cc5d1498b2d Refactor 1D diagonal inner product in quadrature.jl to separate file. Write tests for quadratures. Clean up laplace and secondderivative
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 311
diff changeset
36 function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
37 N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v))
380
81053b1992b6 Remove parity field from Secondderivative
Jonatan Werpers <jonatan@werpers.com>
parents: 361
diff changeset
38 return @inbounds D2.h_inv*D2.h_inv*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I))
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
39 end
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
40
328
9cc5d1498b2d Refactor 1D diagonal inner product in quadrature.jl to separate file. Write tests for quadratures. Clean up laplace and secondderivative
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 311
diff changeset
41 function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T
287
dd621017b695 Change apply_2nd_derivative to Lazy Tensor apply
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 286
diff changeset
42 N = length(v) # TODO: Use domain_size here instead?
311
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
43 r = getregion(Int(index), closuresize(D2), N)
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
44 I = Index(Int(index), r)
f2d6ec89dfc5 Make apply dispatch on Index instead of Tuples of Index
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 299
diff changeset
45 return LazyTensors.apply(D2, v, I)
286
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
46 end
7247e85dc1e8 Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
47
328
9cc5d1498b2d Refactor 1D diagonal inner product in quadrature.jl to separate file. Write tests for quadratures. Clean up laplace and secondderivative
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 311
diff changeset
48 closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} = M