Mercurial > repos > public > sbplib_julia
changeset 1384:851d1e4ab3de
Merge feature/variable_derivatives
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 08 Jun 2023 15:52:22 +0200 |
| parents | 8e59c85a25c0 (current diff) a2ebd7fa011a (diff) |
| children | 896d978725fa 4d628c83987e 3d6425c36d32 7694b35d137d bdcdbd4ea9cd 86026367a9ff 9c689a627244 |
| files | |
| diffstat | 10 files changed, 826 insertions(+), 57 deletions(-) [+] |
line wrap: on
line diff
--- a/Notes.md Thu Jun 08 09:00:48 2023 +0200 +++ b/Notes.md Thu Jun 08 15:52:22 2023 +0200 @@ -71,59 +71,6 @@ dictionary-structure containing stencils, tuples, scalars and other types ready for input to the methods creating the operators. -## Variable second derivative - -2020-12-08 after discussion with Vidar: -We will have to handle the variable second derivative in a new variant of -VolumeOperator, "SecondDerivativeVariable?". Somehow it needs to know about -the coefficients. They should be provided as an AbstractVector. Where they are -provided is another question. It could be that you provide a reference to the -array to the constructor of SecondDerivativeVariable. If that array is mutable -you are free to change it whenever and the changes should propagate -accordingly. Another option is that the counter part to "Laplace" for this -variable second derivate returns a function or acts like a functions that -takes an Abstract array and returns a SecondDerivativeVariable with the -appropriate array. This would allow syntax like `D2(a)*v`. Can this be made -performant? - -For the 1d case we can have a constructor -`SecondDerivativeVariable(D2::SecondDerivativeVariable, a)` that just creates -a copy with a different `a`. - -Apart from just the second derivative in 1D we need operators for higher -dimensions. What happens if a=a(x,y)? Maybe this can be solved orthogonally to -the `D2(a)*v` issue, meaning that if a constant nD version of -SecondDerivativeVariable is available then maybe it can be wrapped to support -function like syntax. We might have to implement `SecondDerivativeVariable` -for N dimensions which takes a N dimensional a. If this could be easily -closured to allow D(a) syntax we would have come a long way. - -For `Laplace` which might use a variable D2 if it is on a curvilinear grid we -might want to choose how to calculate the metric coefficients. They could be -known on closed form, they could be calculated from the grid coordinates or -they could be provided as a vector. Which way you want to do it might change -depending on for example if you are memory bound or compute bound. This choice -cannot be done on the grid since the grid shouldn't care about the computer -architecture. The most sensible option seems to be to have an argument to the -`Laplace` function which controls how the coefficients are gotten from the -grid. The argument could for example be a function which is to be applied to -the grid. - -What happens if the grid or the varible coefficient is dependent on time? -Maybe it becomes important to support `D(a)` or even `D(t,a)` syntax in a more -general way. - -``` -g = TimeDependentGrid() -L = Laplace(g) -function Laplace(g::TimeDependentGrid) - g_logical = logical(g) # g_logical is time independent - ... Build a L(a) assuming we can do that ... - a(t) = metric_coeffs(g,t) - return t->L(a(t)) -end -``` - ## Known size of range and domain? Is there any reason to use a trait to differentiate between fixed size and unknown size?
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/benchmark/benchmark_laplace.jl Thu Jun 08 15:52:22 2023 +0200 @@ -0,0 +1,220 @@ +using Sbplib +using Sbplib.SbpOperators +using Sbplib.Grids +using Sbplib.RegionIndices +using BenchmarkTools + +# TODO: Move the below benchmarks into the benchmark suite + +const operator_path = sbp_operators_path()*"standard_diagonal.toml" + +function benchmark_const_coeff_1d(;N = 100, order = 4) + stencil_set = read_stencil_set(operator_path; order=order) + g = equidistant_grid(N, 0., 1.) + D = second_derivative(g, stencil_set) + u = rand(size(g)...) + u_xx = rand(size(g)...) + + b_naive = @benchmark $u_xx .= $D*$u + b_reg = @benchmark $apply_region_1d!($u_xx,$u,$D) + b_thrd = @benchmark $apply_region_threaded_1d!($u_xx,$u,$D) + print_benchmark_result("Naive apply",b_naive) + print_benchmark_result("Region apply",b_reg) + print_benchmark_result("Threaded region apply",b_thrd) +end + +function benchmark_var_coeff_1d(;N = 100, order = 4) + stencil_set = read_stencil_set(operator_path; order=order) + g = equidistant_grid(N, 0., 1.) + c = rand(size(g)...) + c_lz = eval_on(g, x -> 0.5) + D = second_derivative_variable(g, c, stencil_set) + D_lz = second_derivative_variable(g, c_lz, stencil_set) + u = rand(size(g)...) + u_xx = rand(size(g)...) + + b_naive = @benchmark $u_xx .= $D*$u + b_naive_lz = @benchmark $u_xx .= $D_lz*$u + b_reg = @benchmark $apply_region_1d!($u_xx,$u,$D) + b_reg_lz = @benchmark $apply_region_1d!($u_xx,$u,$D_lz) + b_thrd = @benchmark $apply_region_threaded_1d!($u_xx,$u,$D) + b_thrd_lz = @benchmark $apply_region_threaded_1d!($u_xx,$u,$D_lz) + print_benchmark_result("Naive apply",b_naive) + print_benchmark_result("Naive apply lazy coeff",b_naive_lz) + print_benchmark_result("Region apply",b_reg) + print_benchmark_result("Region apply lazy coeff",b_reg_lz) + print_benchmark_result("Threaded region apply",b_thrd) + print_benchmark_result("Threaded region apply lazy coeff",b_thrd_lz) +end + +function benchmark_const_coeff_2d(;N = 100, order = 4) + stencil_set = read_stencil_set(operator_path; order=order) + g = equidistant_grid((N,N), (0.,0.,),(1.,1.)) + D = Laplace(g, stencil_set) + u = rand(size(g)...) + u_xx = rand(size(g)...) + if order == 2 + clz_sz = 1 + elseif order == 4 + clz_sz = 4 + else + error() + end + + b_naive = @benchmark $u_xx .= $D*$u + b_reg = @benchmark $apply_region_2d!($u_xx,$u,$D,$clz_sz) + b_thrd = @benchmark $apply_region_threaded_2d!($u_xx,$u,$D,$clz_sz) + print_benchmark_result("Naive apply",b_naive) + print_benchmark_result("Region apply",b_reg) + print_benchmark_result("Threaded region apply",b_thrd) +end + +function benchmark_var_coeff_2d(;N = 100, order = 4) + stencil_set = read_stencil_set(operator_path; order=order) + g = equidistant_grid((N,N), (0.,0.,),(1.,1.)) + c = rand(size(g)...) + c_lz = eval_on(g, x-> 0.5) + D = second_derivative_variable(g, c, stencil_set, 1) + second_derivative_variable(g, c, stencil_set, 2) + D_lz = second_derivative_variable(g, c_lz, stencil_set, 1) + second_derivative_variable(g, c_lz, stencil_set, 2) + u = rand(size(g)...) + u_xx = rand(size(g)...) + + if order == 2 + clz_sz = 1 + elseif order == 4 + clz_sz = 6 + else + error() + end + + # Check correctnesss + # u_xx .= D*u + # u_xx_tmp = zeros(size(u_xx)...) + # u_xx_tmp .= u_xx + # apply_region_threaded_2d!(u_xx, u, D, clz_sz) + + # @show sum(abs.(u_xx_tmp .- u_xx)) + # @show pointer(u_xx_tmp) == pointer(u_xxs + + + b_naive = @benchmark $u_xx .= $D*$u + b_naive_lz = @benchmark $u_xx .= $D_lz*$u + b_reg = @benchmark $apply_region_2d!($u_xx,$u,$D, $clz_sz) + b_reg_lz = @benchmark $apply_region_2d!($u_xx,$u,$D_lz, $clz_sz) + b_thrd = @benchmark $apply_region_threaded_2d!($u_xx,$u,$D, $clz_sz) + b_thrd_lz = @benchmark $apply_region_threaded_2d!($u_xx,$u,$D_lz, $clz_sz) + print_benchmark_result("Naive apply",b_naive) + print_benchmark_result("Naive apply lazy coeff",b_naive_lz) + print_benchmark_result("Region apply",b_reg) + print_benchmark_result("Region apply lazy coeff",b_reg_lz) + print_benchmark_result("Threaded region apply",b_thrd) + print_benchmark_result("Threaded region apply lazy coeff",b_thrd_lz) +end + +function print_benchmark_result(title_str,res) + if title_str[1] != ' ' + title_str = lpad(title_str,length(title_str)+1, " ") + end + if title_str[end] != ' ' + title_str = rpad(title_str,length(title_str)+1, " ") + end + tot_len = 76 + pad_len = Int(tot_len/2) + header = lpad(title_str,pad_len,"*") + header = rpad(header,tot_len,"*") + bottom = repeat("*",tot_len) + println(header) + display(res) + println(bottom) + return +end + +function apply_region_1d!(u_xx, u, D) + clz_sz = SbpOperators.closure_size(D) + tm = D*u + for i ∈ @view eachindex(u)[1:clz_sz] + u_xx[i] = tm[Index{Lower}(i)] + end + for i ∈ @view eachindex(u)[clz_sz+1:end-clz_sz] + u_xx[i] = tm[Index{Interior}(i)] + end + for i ∈ @view eachindex(u)[end-clz_sz+1:end] + u_xx[i] = tm[Index{Upper}(i)] + end +end + +function apply_region_threaded_1d!(u_xx, u, D) + clz_sz = SbpOperators.closure_size(D) + tm = D*u + for i ∈ @view eachindex(u)[1:clz_sz] + u_xx[i] = tm[Index{Lower}(i)] + end + Threads.@threads for i ∈ @view eachindex(u)[clz_sz+1:end-clz_sz] + u_xx[i] = tm[Index{Interior}(i)] + end + for i ∈ @view eachindex(u)[end-clz_sz+1:end] + u_xx[i] = tm[Index{Upper}(i)] + end +end + +function apply_region_2d!(u_xx, u, D, clz_sz) + tm = D*u + for I ∈ @view CartesianIndices(u)[1:clz_sz,1:clz_sz] + u_xx[I] = tm[Index{Lower}(I[1]),Index{Lower}(I[2])] + end + for I ∈ @view CartesianIndices(u)[1:clz_sz,clz_sz+1:end-clz_sz] + u_xx[I] = tm[Index{Lower}(I[1]),Index{Interior}(I[2])] + end + for I ∈ @view CartesianIndices(u)[1:clz_sz,end-clz_sz+1:end] + u_xx[I] = tm[Index{Lower}(I[1]),Index{Upper}(I[2])] + end + for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,1:clz_sz] + u_xx[I] = tm[Index{Interior}(I[1]),Index{Lower}(I[2])] + end + for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,clz_sz+1:end-clz_sz] + u_xx[I] = tm[Index{Interior}(I[1]),Index{Interior}(I[2])] + end + for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,end-clz_sz+1:end] + u_xx[I] = tm[Index{Interior}(I[1]),Index{Upper}(I[2])] + end + for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,1:clz_sz] + u_xx[I] = tm[Index{Upper}(I[1]),Index{Lower}(I[2])] + end + for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,clz_sz+1:end-clz_sz] + u_xx[I] = tm[Index{Upper}(I[1]),Index{Interior}(I[2])] + end + for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,end-clz_sz+1:end] + u_xx[I] = tm[Index{Upper}(I[1]),Index{Upper}(I[2])] + end +end + +function apply_region_threaded_2d!(u_xx, u, D, clz_sz) + tm = D*u + for I ∈ @view CartesianIndices(u)[1:clz_sz,1:clz_sz] + u_xx[I] = tm[Index{Lower}(I[1]),Index{Lower}(I[2])] + end + for I ∈ @view CartesianIndices(u)[1:clz_sz,clz_sz+1:end-clz_sz] + u_xx[I] = tm[Index{Lower}(I[1]),Index{Interior}(I[2])] + end + for I ∈ @view CartesianIndices(u)[1:clz_sz,end-clz_sz+1:end] + u_xx[I] = tm[Index{Lower}(I[1]),Index{Upper}(I[2])] + end + for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,1:clz_sz] + u_xx[I] = tm[Index{Interior}(I[1]),Index{Lower}(I[2])] + end + Threads.@threads for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,clz_sz+1:end-clz_sz] + u_xx[I] = tm[Index{Interior}(I[1]),Index{Interior}(I[2])] + end + for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,end-clz_sz+1:end] + u_xx[I] = tm[Index{Interior}(I[1]),Index{Upper}(I[2])] + end + for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,1:clz_sz] + u_xx[I] = tm[Index{Upper}(I[1]),Index{Lower}(I[2])] + end + for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,clz_sz+1:end-clz_sz] + u_xx[I] = tm[Index{Upper}(I[1]),Index{Interior}(I[2])] + end + for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,end-clz_sz+1:end] + u_xx[I] = tm[Index{Upper}(I[1]),Index{Upper}(I[2])] + end +end \ No newline at end of file
--- a/benchmark/benchmarks.jl Thu Jun 08 09:00:48 2023 +0200 +++ b/benchmark/benchmarks.jl Thu Jun 08 15:52:22 2023 +0200 @@ -1,10 +1,13 @@ using BenchmarkTools + using Sbplib using Sbplib.Grids using Sbplib.SbpOperators using Sbplib.RegionIndices using Sbplib.LazyTensors +using LinearAlgebra + const SUITE = BenchmarkGroup() @@ -69,6 +72,30 @@ SUITE["derivatives"]["second_derivative"]["3D"]["z"] = @benchmarkable $u3 .= $Dz*$v3 +SUITE["derivatives"]["second_derivative_variable"] = BenchmarkGroup() + +c1 = map(x->sin(x)+2, g1) +D₂ = second_derivative_variable(g1, c1, stencil_set) +SUITE["derivatives"]["second_derivative_variable"]["1D"] = @benchmarkable $u1 .= $D₂*$v1 + +c2 = map(x->sin(x[1] + x[2])+2, g2) +Dx = second_derivative_variable(g2, c2, stencil_set, 1) +Dy = second_derivative_variable(g2, c2, stencil_set, 2) +SUITE["derivatives"]["second_derivative_variable"]["2D"] = BenchmarkGroup() +SUITE["derivatives"]["second_derivative_variable"]["2D"]["x"] = @benchmarkable $u2 .= $Dx*$v2 +SUITE["derivatives"]["second_derivative_variable"]["2D"]["y"] = @benchmarkable $u2 .= $Dy*$v2 + +c3 = map(x->sin(norm(x))+2, g3) +Dx = second_derivative_variable(g3, c3, stencil_set, 1) +Dy = second_derivative_variable(g3, c3, stencil_set, 2) +Dz = second_derivative_variable(g3, c3, stencil_set, 3) +SUITE["derivatives"]["second_derivative_variable"]["3D"] = BenchmarkGroup() +SUITE["derivatives"]["second_derivative_variable"]["3D"]["x"] = @benchmarkable $u3 .= $Dx*$v3 +SUITE["derivatives"]["second_derivative_variable"]["3D"]["y"] = @benchmarkable $u3 .= $Dy*$v3 +SUITE["derivatives"]["second_derivative_variable"]["3D"]["z"] = @benchmarkable $u3 .= $Dz*$v3 + + + SUITE["derivatives"]["addition"] = BenchmarkGroup()
--- a/src/SbpOperators/SbpOperators.jl Thu Jun 08 09:00:48 2023 +0200 +++ b/src/SbpOperators/SbpOperators.jl Thu Jun 08 15:52:22 2023 +0200 @@ -5,6 +5,7 @@ export read_stencil_set export get_stencil_set export parse_stencil +export parse_nested_stencil export parse_scalar export parse_tuple export sbp_operators_path @@ -19,6 +20,7 @@ export normal_derivative export first_derivative export second_derivative +export second_derivative_variable export undivided_skewed04 using Sbplib.RegionIndices @@ -30,6 +32,8 @@ even = 1 end +export closure_size + include("stencil.jl") include("stencil_set.jl") include("volumeops/volume_operator.jl") @@ -37,6 +41,7 @@ include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/first_derivative.jl") include("volumeops/derivatives/second_derivative.jl") +include("volumeops/derivatives/second_derivative_variable.jl") include("volumeops/derivatives/dissipation.jl") include("volumeops/laplace/laplace.jl") include("volumeops/inner_products/inner_product.jl")
--- a/src/SbpOperators/operators/standard_diagonal.toml Thu Jun 08 09:00:48 2023 +0200 +++ b/src/SbpOperators/operators/standard_diagonal.toml Thu Jun 08 15:52:22 2023 +0200 @@ -20,6 +20,10 @@ H.inner = "1" H.closure = ["1/2"] +e.closure = ["1"] +d1.closure = {s = ["3/2", "-2", "1/2"], c = 1} + + D1.inner_stencil = ["-1/2", "0", "1/2"] D1.closure_stencils = [ {s = ["-1", "1"], c = 1}, @@ -30,8 +34,10 @@ {s = ["1", "-2", "1"], c = 1}, ] -e.closure = ["1"] -d1.closure = {s = ["3/2", "-2", "1/2"], c = 1} +D2variable.inner_stencil = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]] +D2variable.closure_stencils = [ + {s = [["2", "-1", "0"],["-3", "1", "0"],["1","0","0"]], c = 1}, +] [[stencil_set]] @@ -40,6 +46,9 @@ H.inner = "1" H.closure = ["17/48", "59/48", "43/48", "49/48"] +e.closure = ["1"] +d1.closure = {s = ["11/6", "-3", "3/2", "-1/3"], c = 1} + D1.inner_stencil = ["1/12","-2/3","0","2/3","-1/12"] D1.closure_stencils = [ {s = [ "-24/17", "59/34", "-4/17", "-3/34", "0", "0"], c = 1}, @@ -56,5 +65,88 @@ {s = [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"], c = 4}, ] +D2variable.inner_stencil = [ + ["-1/8", "1/6", "-1/8", "0", "0" ], + [ "1/6", "1/2", "1/2", "1/6", "0" ], + ["-1/24", "-5/6", "-3/4", "-5/6", "-1/24"], + [ "0", "1/6", "1/2", "1/2", "1/6" ], + [ "0", "0", "-1/8", "1/6", "-1/8" ], +] +D2variable.closure_stencils = [ + {c = 1, s = [ + [ "920/289", "-59/68", "-81031200387/366633756146", "-69462376031/733267512292", "0", "0", "0", "0" ], + ["-1740/289", "0", "6025413881/7482321554", "1612249989/7482321554", "0", "0", "0", "0" ], + [ "1128/289", "59/68", "-6251815797/8526366422", "-639954015/17052732844", "0", "0", "0", "0" ], + [ "-308/289", "0", "1244724001/7482321554", "-752806667/7482321554", "0", "0", "0", "0" ], + [ "0", "0", "-148737261/10783345769", "148737261/10783345769", "0", "0", "0", "0" ], + [ "0", "0", "-3/833", "3/833", "0", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + ]}, + {c = 2, s = [ + [ "12/17", "0", "102125659/440136562", "27326271/440136562", "0", "0", "0", "0" ], + [ "-59/68", "0", "-156920047993625/159775733917868", "-12001237118451/79887866958934", "0", "0", "0", "0" ], + [ "2/17", "0", "1489556735319/1857857371138", "149729180391/1857857371138", "0", "0", "0", "0" ], + [ "3/68", "0", "-13235456910147/159775733917868", "3093263736297/79887866958934", "0", "0", "0", "0" ], + [ "0", "0", "67535018271/2349643145851", "-67535018271/2349643145851", "0", "0", "0", "0" ], + [ "0", "0", "441/181507", "-441/181507", "0", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + ]}, + {c = 3, s = [ + [ "-96/731", "59/172", "-6251815797/21566691538", "-639954015/43133383076", "0", "0", "0", "0" ], + [ "118/731", "0", "87883847383821/79887866958934", "8834021643069/79887866958934", "0", "0", "0", "0" ], + [ "-16/731", "-59/172", "-1134866646907639536627/727679167377258785038", "-13777050223300597/23487032885926596", "-26254/557679", "0", "0", "0" ], + [ "-6/731", "0", "14509020271326561681/14850595252597118062", "17220493277981/79887866958934", "1500708/7993399", "0", "0", "0" ], + [ "0", "0", "-4841930283098652915/21402328452272317207", "31597236232005/115132514146699", "-26254/185893", "0", "0", "0" ], + [ "0", "0", "-2318724711/1653303156799", "960119/1147305747", "13564/23980197", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + ]}, + {c = 4, s = [ + [ "-36/833", "0", "1244724001/21566691538", "-752806667/21566691538", "0", "0", "0", "0" ], + [ "177/3332", "0", "-780891957698673/7829010961975532", "3724542049827/79887866958934", "0", "0", "0", "0" ], + [ "-6/833", "0", "14509020271326561681/16922771334354855466", "2460070468283/13005001597966", "1500708/9108757", "0", "0", "0" ], + [ "-9/3332", "0", "-217407431400324796377/207908333536359652868", "-1950062198436997/3914505480987766", "-7476412/9108757", "-2/49", "0", "0" ], + [ "0", "0", "4959271814984644613/21402328452272317207", "47996144728947/115132514146699", "4502124/9108757", "8/49", "0", "0" ], + [ "0", "0", "-2258420001/1653303156799", "-1063649/8893843", "1473580/9108757", "-6/49", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + ]}, + {c = 5, s = [ + [ "0", "0", "-49579087/10149031312", "49579087/10149031312", "0", "0", "0", "0" ], + [ "0", "0", "1328188692663/37594290333616", "-1328188692663/37594290333616", "0", "0", "0", "0" ], + [ "0", "0", "-1613976761032884305/7963657098519931984", "10532412077335/42840005263888", "-564461/4461432", "0", "0", "0" ], + [ "0", "0", "4959271814984644613/20965546238960637264", "15998714909649/37594290333616", "375177/743572", "1/6", "0", "0" ], + [ "0", "0", "-8386761355510099813/128413970713633903242", "-2224717261773437/2763180339520776", "-280535/371786", "-5/6", "-1/24", "0" ], + [ "0", "0", "13091810925/13226425254392", "35039615/213452232", "1118749/2230716", "1/2", "1/6", "0" ], + [ "0", "0", "0", "0", "-1/8", "1/6", "-1/8", "0" ], + [ "0", "0", "0", "0", "0", "0", "0", "0" ], + ]}, + {c = 6, s = [ + [ "0", "0", "-1/784", "1/784", "0", "0", "0", "0" ], + [ "0", "0", "8673/2904112", "-8673/2904112", "0", "0", "0", "0" ], + [ "0", "0", "-33235054191/26452850508784", "960119/1280713392", "3391/6692148", "0", "0", "0" ], + [ "0", "0", "-752806667/539854092016", "-1063649/8712336", "368395/2230716", "-1/8", "0", "0" ], + [ "0", "0", "13091810925/13226425254392", "35039615/213452232", "1118749/2230716", "1/2", "1/6", "0" ], + [ "0", "0", "-660204843/13226425254392", "-3290636/80044587", "-5580181/6692148", "-3/4", "-5/6", "-1/24"], + [ "0", "0", "0", "0", "1/6", "1/2", "1/2", "1/6" ], + [ "0", "0", "0", "0", "0", "-1/8", "1/6", "-1/8" ], + ]} +] + + + +[[stencil_set]] + +order = 6 + +H.inner = "1" +H.closure = ["13649/43200", "12013/8640", "2711/4320", "5359/4320", "7877/8640", "43801/43200"] + + + + + e.closure = ["1"] -d1.closure = {s = ["11/6", "-3", "3/2", "-1/3"], c = 1} +d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"]
--- a/src/SbpOperators/stencil_set.jl Thu Jun 08 09:00:48 2023 +0200 +++ b/src/SbpOperators/stencil_set.jl Thu Jun 08 15:52:22 2023 +0200 @@ -110,6 +110,33 @@ end end + +""" + parse_nested_stencil(parsed_toml) + +Accept parsed TOML and read it as a nested tuple. + +See also [`read_stencil_set`](@ref), [`parse_stencil`](@ref). +""" +function parse_nested_stencil(parsed_toml) + if parsed_toml isa Array + weights = parse_stencil.(parsed_toml) + return CenteredNestedStencil(weights...) + end + + center = parsed_toml["c"] + weights = parse_tuple.(parsed_toml["s"]) + return NestedStencil(weights...; center) +end + +""" + parse_nested_stencil(T, parsed_toml) + +Parse the input as a nested stencil with element type `T`. +""" +parse_nested_stencil(T, parsed_toml) = NestedStencil{T}(parse_nested_stencil(parsed_toml)) + + """ parse_scalar(parsed_toml)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Thu Jun 08 15:52:22 2023 +0200 @@ -0,0 +1,188 @@ +""" + 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, 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, stencil_set) + return second_derivative_variable(TensorGrid(g), coeff, stencil_set, 1) +end + +function second_derivative_variable(g::TensorGrid, coeff, 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(coeff, scaled_inner_stencil, scaled_closure_stencils, dir) +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 + + +""" + SecondDerivativeVariable{Dir,T,D,...} <: LazyTensor{T,D,D} + +A second derivative operator in direction `Dir` with a variable coefficient. +""" +struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: LazyTensor{T,D,D} + inner_stencil::IStencil + closure_stencils::NTuple{M,CStencil} + coefficient::TArray + + function SecondDerivativeVariable(coefficient::AbstractArray, inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, dir) where {T,M} + D = ndims(coefficient) + IStencil = typeof(inner_stencil) + CStencil = eltype(closure_stencils) + TArray = typeof(coefficient) + return new{dir,T,D,M,IStencil,CStencil,TArray}(inner_stencil, closure_stencils, coefficient) + end +end + +derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir + +closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils) + +LazyTensors.range_size(op::SecondDerivativeVariable) = size(op.coefficient) +LazyTensors.domain_size(op::SecondDerivativeVariable) = size(op.coefficient) + + +function derivative_view(op, a, I) + d = derivative_direction(op) + + Iview = Base.setindex(I,:,d) + return @view a[Iview...] +end + +function apply_lower(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i) +end + +function apply_interior(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + return apply_stencil(op.inner_stencil, c̃, ṽ, i) +end + +function apply_upper(op::SecondDerivativeVariable, v, I...) + ṽ = derivative_view(op, v, I) + c̃ = derivative_view(op, op.coefficient, I) + + i = I[derivative_direction(op)] + sz = domain_size(op)[derivative_direction(op)] + stencil = op.closure_stencils[sz-i+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I::Vararg{Index}) + if I[derivative_direction(op)] isa Index{Lower} + return apply_lower(op, v, Int.(I)...) + elseif I[derivative_direction(op)] isa Index{Upper} + return apply_upper(op, v, Int.(I)...) + elseif I[derivative_direction(op)] isa Index{Interior} + return apply_interior(op, v, Int.(I)...) + else + error("Invalid region") + end +end + +function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...) + dir = derivative_direction(op) + sz = domain_size(op)[dir] + + i = I[dir] + + I = map(i->Index(i, Interior), I) + if 0 < i <= closure_size(op) + I = Base.setindex(I, Index(i, Lower), dir) + return LazyTensors.apply(op, v, I...) + elseif closure_size(op) < i <= sz-closure_size(op) + I = Base.setindex(I, Index(i, Interior), dir) + return LazyTensors.apply(op, v, I...) + elseif sz-closure_size(op) < i <= sz + I = Base.setindex(I, Index(i, Upper), dir) + return LazyTensors.apply(op, v, I...) + else + error("Bounds error") # This should be `throw(BoundsError())` but the type inference is so fragile that it doesn't work. Needs investigation. / Jonatan 2023-06-08 + end +end + + +# 2D Specific implementations to avoid type instability +# TBD: Can this be solved by fixing the general methods instead? + + +## x-direction +function apply_lower(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i) +end + +function apply_interior(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, i) +end + +function apply_upper(op::SecondDerivativeVariable{1}, v, i, j) + ṽ = @view v[:,j] + c̃ = @view op.coefficient[:,j] + + sz = domain_size(op)[derivative_direction(op)] + stencil = op.closure_stencils[sz-i+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i) +end + + +## y-direction +function apply_lower(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + return @inbounds apply_stencil(op.closure_stencils[j], c̃, ṽ, j) +end + +function apply_interior(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, j) +end + +function apply_upper(op::SecondDerivativeVariable{2}, v, i, j) + ṽ = @view v[i,:] + c̃ = @view op.coefficient[i,:] + + sz = domain_size(op)[derivative_direction(op)] + stencil = op.closure_stencils[sz-j+1] + return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, j) +end
--- a/src/SbpOperators/volumeops/volume_operator.jl Thu Jun 08 09:00:48 2023 +0200 +++ b/src/SbpOperators/volumeops/volume_operator.jl Thu Jun 08 15:52:22 2023 +0200 @@ -12,7 +12,7 @@ function VolumeOperator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity) return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity) -end +end # TBD: Remove this function? closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M
--- a/test/SbpOperators/stencil_set_test.jl Thu Jun 08 09:00:48 2023 +0200 +++ b/test/SbpOperators/stencil_set_test.jl Thu Jun 08 15:52:22 2023 +0200 @@ -4,6 +4,7 @@ using Sbplib.SbpOperators import Sbplib.SbpOperators.Stencil +import Sbplib.SbpOperators.NestedStencil @testset "readoperator" begin toml_str = """ @@ -170,3 +171,18 @@ @test SbpOperators.parse_rational(2) isa Rational @test SbpOperators.parse_rational(2) == 2//1 end + +@testset "parse_nested_stencil" begin + toml = TOML.parse(""" + s1 = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]] + s2 = {s = [[ "2", "-1", "0"],[ "-3", "1", "0"],["1", "0", "0"]], c = 1} + s3 = {s = [[ "2", "-1", "0"],[ "-3", "1", "0"],["1", "0", "0"]], c = 2} + """) + + @test parse_nested_stencil(toml["s1"]) == CenteredNestedStencil((1//2, 1//2, 0//1),( -1//2, -1//1, -1//2),(0//1, 1//2, 1//2)) + @test parse_nested_stencil(toml["s2"]) == NestedStencil((2//1, -1//1, 0//1),( -3//1, 1//1, 0//1),(1//1, 0//1, 0//1), center = 1) + @test parse_nested_stencil(toml["s3"]) == NestedStencil((2//1, -1//1, 0//1),( -3//1, 1//1, 0//1),(1//1, 0//1, 0//1), center = 2) + + @test parse_nested_stencil(Float64, toml["s1"]) == CenteredNestedStencil((1/2, 1/2, 0.),( -1/2, -1., -1/2),(0., 1/2, 1/2)) + @test parse_nested_stencil(Int, toml["s2"]) == NestedStencil((2, -1, 0),( -3, 1, 0),(1, 0, 0), center = 1) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl Thu Jun 08 15:52:22 2023 +0200 @@ -0,0 +1,247 @@ +using Test + +using Sbplib.Grids +using Sbplib.LazyTensors +using Sbplib.SbpOperators +using Sbplib.RegionIndices +using Sbplib.SbpOperators: NestedStencil, CenteredNestedStencil, SecondDerivativeVariable + +using LinearAlgebra + +@testset "second_derivative_variable" begin + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) + + @testset "1D" begin + g = equidistant_grid(11, 0., 1.) + c = [ 1., 3., 6., 10., 15., 21., 28., 36., 45., 55., 66.] + + @testset "checking c" begin + c_short = rand(5) + c_long = rand(16) + c_higher_dimension = rand(11,11) + + @test_throws DimensionMismatch("the size (5,) of the coefficient does not match the size (11,) of the grid") second_derivative_variable(g, c_short, stencil_set) + @test_throws DimensionMismatch("the size (16,) of the coefficient does not match the size (11,) of the grid") second_derivative_variable(g, c_long, stencil_set) + @test_throws ArgumentError("The coefficient has dimension 2 while the grid is dimension 1") second_derivative_variable(g, c_higher_dimension, stencil_set) + end + + @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̄, stencil_set) + 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̄, stencil_set, 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 = ( + NestedStencil(( 2., -1., 0.),(-3., 1., 0.), (1., 0., 0.), center = 1), + ) + + @testset "1D" begin + c = [ 1., 3., 6., 10., 15., 21., 28., 36., 45., 55., 66.] + @testset "Constructors" begin + @test SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) isa LazyTensor + + D₂ᶜ = SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) + @test range_dim(D₂ᶜ) == 1 + @test domain_dim(D₂ᶜ) == 1 + + end + + @testset "sizes" begin + D₂ᶜ = SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) + @test closure_size(D₂ᶜ) == 1 + @test range_size(D₂ᶜ) == (11,) + @test domain_size(D₂ᶜ) == (11,) + end + + @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₂ᶜ = SecondDerivativeVariable(c̄, interior_stencil, closure_stencils, 1) + 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 + + @testset "type stability" begin + g = equidistant_grid(11, 0., 10.) # h = 1 + c̄ = eval_on(g,x-> -1) + v̄ = eval_on(g,x->1.) + + D₂ᶜ = SecondDerivativeVariable(c̄, interior_stencil, closure_stencils, 1) + + @inferred SbpOperators.apply_lower(D₂ᶜ, v̄, 1) + @inferred SbpOperators.apply_interior(D₂ᶜ, v̄, 5) + @inferred SbpOperators.apply_upper(D₂ᶜ, v̄, 11) + @inferred (D₂ᶜ*v̄)[Index(1,Lower)] + 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 "Constructors" begin + @test SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) isa LazyTensor + @test SecondDerivativeVariable(c, interior_stencil, closure_stencils, 2) isa LazyTensor + + D₂ᶜ = SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) + @test range_dim(D₂ᶜ) == 2 + @test domain_dim(D₂ᶜ) == 2 + end + + @testset "sizes" begin + D₂ᶜ = SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) + @test range_size(D₂ᶜ) == (11,9) + @test domain_size(D₂ᶜ) == (11,9) + @test closure_size(D₂ᶜ) == 1 + + D₂ᶜ = SecondDerivativeVariable(c, interior_stencil, closure_stencils, 2) + @test range_size(D₂ᶜ) == (11,9) + @test domain_size(D₂ᶜ) == (11,9) + @test closure_size(D₂ᶜ) == 1 + end + + @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₂ᶜ = SecondDerivativeVariable(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) + end + end +end +