Mercurial > repos > public > sbplib_julia
changeset 1370:4ef8fb75d144 feature/variable_derivatives
Start splitting out a second_derivative_variable function
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Fri, 26 May 2023 21:39:08 +0200 |
| parents | ff64acfc1ec9 |
| children | d0e48c2e6aad |
| files | src/SbpOperators/SbpOperators.jl src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl |
| diffstat | 3 files changed, 164 insertions(+), 97 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl Fri May 26 14:59:37 2023 +0200 +++ b/src/SbpOperators/SbpOperators.jl Fri May 26 21:39:08 2023 +0200 @@ -20,6 +20,7 @@ export normal_derivative export first_derivative export second_derivative +export second_derivative_variable export undivided_skewed04 using Sbplib.RegionIndices
--- a/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Fri May 26 14:59:37 2023 +0200 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Fri May 26 21:39:08 2023 +0200 @@ -1,3 +1,51 @@ +""" + second_derivative_variable(g, coeff ..., [direction]) + +The variable second derivative operator as a `LazyTensor` on the given grid. +`coeff` is a grid function of the variable coefficient. + +Approximates the d/dξ c d/dξ on `g` along the coordinate dimension specified +by `direction`. +""" +function second_derivative_variable end + + +function second_derivative_variable(g::TensorGrid, coeff::AbstractArray, inner_stencil::NestedStencil, closure_stencils, dir) + check_coefficient(g, coeff) + + Δxᵢ = spacing(g.grids[dir]) + scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2) + scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2) + return SecondDerivativeVariable{dir, ndims(g)}(scaled_inner_stencil, scaled_closure_stencils, coeff) +end + +function second_derivative_variable(g::EquidistantGrid, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils) + return second_derivative_variable(TensorGrid(g), coeff, inner_stencil, closure_stencils, 1) +end + +function second_derivative_variable(g::TensorGrid, coeff::AbstractArray, stencil_set, dir::Int) + inner_stencil = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"]) + closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"]) + + return second_derivative_variable(g, coeff, inner_stencil, closure_stencils, dir) +end + +function second_derivative_variable(g::EquidistantGrid, coeff::AbstractArray, stencil_set) + return second_derivative_variable(TensorGrid(g), coeff, stencil_set, 1) +end + + +function check_coefficient(g, coeff) + if ndims(g) != ndims(coeff) + throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(ndims(g))")) + end + + if size(g) != size(coeff) + throw(DimensionMismatch("the size $(size(coeff)) of the coefficient does not match the size $(size(g)) of the grid")) + end +end + + # REVIEW: Fixa docs """ SecondDerivativeVariable{Dir,T,D,...} <: LazyTensor{T,D,D} @@ -17,60 +65,6 @@ end end -function SecondDerivativeVariable(g::TensorGrid, coeff::AbstractArray, inner_stencil::NestedStencil, closure_stencils, dir) - check_coefficient(g, coeff) - - Δxᵢ = spacing(g.grids[dir]) - scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2) - scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2) - return SecondDerivativeVariable{dir, ndims(g)}(scaled_inner_stencil, scaled_closure_stencils, coeff) -end - -function SecondDerivativeVariable(g::EquidistantGrid, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils) - return SecondDerivativeVariable(TensorGrid(g), coeff, inner_stencil, closure_stencils, 1) -end - -@doc raw""" - SecondDerivativeVariable(g::EquidistantGrid, coeff::AbstractArray, stencil_set, dir) - -Create a `LazyTensor` for the second derivative with a variable coefficient -`coeff` on `grid` from the stencils in `stencil_set`. The direction is -determined by `dir`. - -`coeff` is a grid function on `grid`. - -# Example -With -``` -D = SecondDerivativeVariable(g, c, stencil_set, 2) -``` -then `D*u` approximates -```math -\frac{\partial}{\partial y} c(x,y) \frac{\partial u}{\partial y}, -``` -on ``(0,1)⨯(0,1)`` represented by `g`. -""" -function SecondDerivativeVariable(g::TensorGrid, coeff::AbstractArray, stencil_set, dir::Int) - inner_stencil = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"]) - closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"]) - - return SecondDerivativeVariable(g, coeff, inner_stencil, closure_stencils, dir) -end - -function SecondDerivativeVariable(g::EquidistantGrid, coeff::AbstractArray, stencil_set) - return SecondDerivativeVariable(TensorGrid(g), coeff, stencil_set, 1) -end - -function check_coefficient(g, coeff) - if ndims(g) != ndims(coeff) - throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(ndims(g))")) - end - - if size(g) != size(coeff) - throw(DimensionMismatch("the size $(size(coeff)) of the coefficient does not match the size $(size(g)) of the grid")) - end -end - derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl Fri May 26 14:59:37 2023 +0200 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl Fri May 26 21:39:08 2023 +0200 @@ -8,6 +8,121 @@ using LinearAlgebra +@testset "second_derivative_variable" begin + interior_stencil = CenteredNestedStencil((1/2, 1/2, 0.),(-1/2, -1., -1/2),( 0., 1/2, 1/2)) + closure_stencils = [ + NestedStencil(( 2., -1., 0.),(-3., 1., 0.), (1., 0., 0.), center = 1), + ] + + @testset "1D" begin + g = equidistant_grid(11, 0., 1.) + c = [ 1., 3., 6., 10., 15., 21., 28., 36., 45., 55., 66.] + + + @testset "application" begin + + function apply_to_functions(; v, c) + g = equidistant_grid(11, 0., 10.) # h = 1 + c̄ = eval_on(g,c) + v̄ = eval_on(g,v) + + D₂ᶜ = second_derivative_variable(g, c̄, interior_stencil, closure_stencils) + return D₂ᶜ*v̄ + end + + @test apply_to_functions(v=x->1., c=x-> -1.) == zeros(11) + @test apply_to_functions(v=x->1., c=x-> -x ) == zeros(11) + @test apply_to_functions(v=x->x, c=x-> 1.) == zeros(11) + @test apply_to_functions(v=x->x, c=x-> -x ) == -ones(11) + @test apply_to_functions(v=x->x^2, c=x-> 1.) == 2ones(11) + end + end + + @testset "2D" begin + g = equidistant_grid((11,9), (0.,0.), (10.,8.)) # h = 1 + c = eval_on(g, (x,y)->x+y) + + @testset "application" begin + function apply_to_functions(dir; v, c) + g = equidistant_grid((11,9), (0.,0.), (10.,8.)) # h = 1 + c̄ = eval_on(g,c) + v̄ = eval_on(g,v) + + D₂ᶜ = second_derivative_variable(g, c̄, interior_stencil, closure_stencils,dir) + return D₂ᶜ*v̄ + end + + # x-direction + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)-> -1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)->- x ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->x, c=(x,y)-> -x ) == -ones(11,9) + @test apply_to_functions(1,v=(x,y)->x^2, c=(x,y)-> 1.) == 2ones(11,9) + + @test apply_to_functions(1,v=(x,y)->1., c=(x,y)->- y ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y, c=(x,y)-> -y ) == zeros(11,9) + @test apply_to_functions(1,v=(x,y)->y^2, c=(x,y)-> 1.) == zeros(11,9) + + # y-direction + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)-> -1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)->- y ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->y, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->y, c=(x,y)-> -y ) == -ones(11,9) + @test apply_to_functions(2,v=(x,y)->y^2, c=(x,y)-> 1.) == 2ones(11,9) + + @test apply_to_functions(2,v=(x,y)->1., c=(x,y)->- x ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> -x ) == zeros(11,9) + @test apply_to_functions(2,v=(x,y)->x^2, c=(x,y)-> 1.) == zeros(11,9) + + + @testset "standard diagonal operators" begin + c(x,y) = exp(x) + exp(1.5(1-y)) + v(x,y) = sin(x) + cos(1.5(1-y)) + + Dxv(x,y) = cos(x)*exp(x) - (exp(x) + exp(1.5 - 1.5y))*sin(x) + Dyv(x,y) = -1.5(1.5exp(x) + 1.5exp(1.5 - 1.5y))*cos(1.5 - 1.5y) - 2.25exp(1.5 - 1.5y)*sin(1.5 - 1.5y) + + g₁ = equidistant_grid((60,67), (0.,0.), (1.,2.)) + g₂ = refine(g₁,2) + + c̄₁ = eval_on(g₁, c) + c̄₂ = eval_on(g₂, c) + + v̄₁ = eval_on(g₁, v) + v̄₂ = eval_on(g₂, v) + + + function convergence_rate_estimate(stencil_set, dir, Dv_true) + D₁ = second_derivative_variable(g₁, c̄₁, stencil_set, dir) + D₂ = second_derivative_variable(g₂, c̄₂, stencil_set, dir) + + Dv̄₁ = D₁*v̄₁ + Dv̄₂ = D₂*v̄₂ + + Dv₁ = eval_on(g₁,Dv_true) + Dv₂ = eval_on(g₂,Dv_true) + + e₁ = norm(Dv̄₁ - Dv₁)/norm(Dv₁) + e₂ = norm(Dv̄₂ - Dv₂)/norm(Dv₂) + + return log2(e₁/e₂) + end + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) + @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 1.5 rtol = 1e-1 + @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 1.5 rtol = 1e-1 + + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) + @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 2.5 rtol = 1e-1 + @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 2.5 rtol = 2e-1 + end + end + end +end + + @testset "SecondDerivativeVariable" begin interior_stencil = CenteredNestedStencil((1/2, 1/2, 0.),(-1/2, -1., -1/2),( 0., 1/2, 1/2)) closure_stencils = [ @@ -139,49 +254,6 @@ @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> 1.) == zeros(11,9) @test apply_to_functions(2,v=(x,y)->x, c=(x,y)-> -x ) == zeros(11,9) @test apply_to_functions(2,v=(x,y)->x^2, c=(x,y)-> 1.) == zeros(11,9) - - - @testset "standard diagonal operators" begin - c(x,y) = exp(x) + exp(1.5(1-y)) - v(x,y) = sin(x) + cos(1.5(1-y)) - - Dxv(x,y) = cos(x)*exp(x) - (exp(x) + exp(1.5 - 1.5y))*sin(x) - Dyv(x,y) = -1.5(1.5exp(x) + 1.5exp(1.5 - 1.5y))*cos(1.5 - 1.5y) - 2.25exp(1.5 - 1.5y)*sin(1.5 - 1.5y) - - g₁ = equidistant_grid((60,67), (0.,0.), (1.,2.)) - g₂ = refine(g₁,2) - - c̄₁ = eval_on(g₁, c) - c̄₂ = eval_on(g₂, c) - - v̄₁ = eval_on(g₁, v) - v̄₂ = eval_on(g₂, v) - - - function convergence_rate_estimate(stencil_set, dir, Dv_true) - D₁ = SecondDerivativeVariable(g₁, c̄₁, stencil_set, dir) - D₂ = SecondDerivativeVariable(g₂, c̄₂, stencil_set, dir) - - Dv̄₁ = D₁*v̄₁ - Dv̄₂ = D₂*v̄₂ - - Dv₁ = eval_on(g₁,Dv_true) - Dv₂ = eval_on(g₂,Dv_true) - - e₁ = norm(Dv̄₁ - Dv₁)/norm(Dv₁) - e₂ = norm(Dv̄₂ - Dv₂)/norm(Dv₂) - - return log2(e₁/e₂) - end - - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 2) - @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 1.5 rtol = 1e-1 - @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 1.5 rtol = 1e-1 - - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) - @test convergence_rate_estimate(stencil_set, 1, Dxv) ≈ 2.5 rtol = 1e-1 - @test convergence_rate_estimate(stencil_set, 2, Dyv) ≈ 2.5 rtol = 2e-1 - end end end end