Mercurial > repos > public > sbplib_julia
changeset 1635:b62770cec7d0 update/julia_1.10.3
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 25 Jun 2024 15:32:19 +0200 |
| parents | 7f5e4d2a5112 (current diff) 29c7a787e67d (diff) |
| children | d3ae8df957c1 |
| files | benchmark/benchmarks.jl |
| diffstat | 14 files changed, 390 insertions(+), 10 deletions(-) [+] |
line wrap: on
line diff
diff -r 7f5e4d2a5112 -r b62770cec7d0 Notes.md --- a/Notes.md Sat May 25 16:02:55 2024 -0700 +++ b/Notes.md Tue Jun 25 15:32:19 2024 +0200 @@ -1,5 +1,26 @@ # Notes +## How to dispatch for different operators +We have a problem in how dispatch for different operators work. + * We want to keep the types simple and flat (Awkward to forward `apply`) + * We want to dispatch SATs on the parameters of the continuous operator. (a * div for example) + * We want to allow keeping the same stencil_set across different calls. (maybe not so bad for the user to be responsible) + +Could remove the current opset idea and introduce a description of continuous operators + ```julia +abstract type DifferentialOperator end + +struct Laplace <: DifferentialOperator end +struct Advection <: DifferentialOperator + v +end + +difference_operator(::Laplace, grid, stencil_set) = ... # Returns a plain LazyTensor. Replaces the current `laplace()` function. +sat_tensors(::Laplace, grid, stencil_set, bc) = ... + +sat(::DifferentialOperator, grid, stencil_set, bc) = ... + ``` + ## Reading operators Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and
diff -r 7f5e4d2a5112 -r b62770cec7d0 Project.toml --- a/Project.toml Sat May 25 16:02:55 2024 -0700 +++ b/Project.toml Tue Jun 25 15:32:19 2024 +0200 @@ -1,7 +1,7 @@ name = "Sbplib" uuid = "5a373a26-915f-4769-bcab-bf03835de17b" authors = ["Jonatan Werpers <jonatan@werpers.com>", "Vidar Stiernström <vidar.stiernstrom@it.uu.se>, and contributors"] -version = "0.1.0" +version = "0.1.1" [deps] StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff -r 7f5e4d2a5112 -r b62770cec7d0 TODO.md --- a/TODO.md Sat May 25 16:02:55 2024 -0700 +++ b/TODO.md Tue Jun 25 15:32:19 2024 +0200 @@ -21,6 +21,10 @@ See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing) - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe - [ ] Move export statements to top of each module + - [ ] Implement apply_transpose for + - [ ] ElementwiseTensorOperation + - [ ] VolumeOperator + - [ ] Laplace - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative. Kan vi använda promote istället?
diff -r 7f5e4d2a5112 -r b62770cec7d0 benchmark/benchmarks.jl --- a/benchmark/benchmarks.jl Sat May 25 16:02:55 2024 -0700 +++ b/benchmark/benchmarks.jl Tue Jun 25 15:32:19 2024 +0200 @@ -15,7 +15,7 @@ ll(d) = ntuple(i->0., d) lu(d) = ntuple(i->1., d) -g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)[1]) +g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)...) g2 = equidistant_grid(ll(2), lu(2), sz(2)...) g3 = equidistant_grid(ll(3), lu(3), sz(3)...)
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/Grids/grid.jl --- a/src/Grids/grid.jl Sat May 25 16:02:55 2024 -0700 +++ b/src/Grids/grid.jl Tue Jun 25 15:32:19 2024 +0200 @@ -74,7 +74,7 @@ end end -Base.size(cv) = size(cv.v) +Base.size(cv::ArrayComponentView) = size(cv.v) Base.getindex(cv::ArrayComponentView, i::Int) = cv.v[i][cv.component_index...] Base.getindex(cv::ArrayComponentView, I::Vararg{Int}) = cv.v[I...][cv.component_index...] IndexStyle(::Type{<:ArrayComponentView{<:Any,<:Any,AT}}) where AT = IndexStyle(AT)
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/LazyTensors/lazy_tensor_operations.jl --- a/src/LazyTensors/lazy_tensor_operations.jl Sat May 25 16:02:55 2024 -0700 +++ b/src/LazyTensors/lazy_tensor_operations.jl Tue Jun 25 15:32:19 2024 +0200 @@ -121,6 +121,7 @@ Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm) Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm +Base.:-(tm::LazyTensor) = (-one(eltype(tm)))*tm """ InflatedTensor{T,R,D} <: LazyTensor{T,R,D}
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/SbpOperators/SbpOperators.jl --- a/src/SbpOperators/SbpOperators.jl Sat May 25 16:02:55 2024 -0700 +++ b/src/SbpOperators/SbpOperators.jl Tue Jun 25 15:32:19 2024 +0200 @@ -11,7 +11,6 @@ export sbp_operators_path # Operators -export boundary_quadrature export boundary_restriction export inner_product export inverse_inner_product @@ -22,20 +21,34 @@ export second_derivative export second_derivative_variable export undivided_skewed04 - -using Sbplib.RegionIndices -using Sbplib.LazyTensors -using Sbplib.Grids +export closure_size @enum Parity begin odd = -1 even = 1 end -export closure_size +# Boundary conditions +export BoundaryCondition +export NeumannCondition +export DirichletCondition +export discretize_data +export boundary_data +export boundary +export sat +export sat_tensors + +# Using +using Sbplib.RegionIndices +using Sbplib.LazyTensors +using Sbplib.Grids + +# Includes include("stencil.jl") include("stencil_set.jl") +include("boundary_conditions/boundary_condition.jl") +include("boundary_conditions/sat.jl") include("volumeops/volume_operator.jl") include("volumeops/stencil_operator_distinct_closures.jl") include("volumeops/constant_interior_scaling_operator.jl")
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/SbpOperators/boundary_conditions/boundary_condition.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/boundary_conditions/boundary_condition.jl Tue Jun 25 15:32:19 2024 +0200 @@ -0,0 +1,62 @@ +""" + BoundaryCondition + +Description of a boundary condition. Implementations describe the kind of +boundary condition, what boundary the condition applies to, and any associated +data. Should implement [`boundary`](@ref) and may implement +[`boundary_data`](@ref) if applicable. + +For examples see [`DirichletCondition`](@ref) and [`NeumannCondition`](@ref) +""" +abstract type BoundaryCondition end + +""" + boundary(::BoundaryCondition) + +The boundary identifier of the BoundaryCondition. +""" +function boundary end + +""" + boundary_data(::BoundaryCondition) + +If implemented, the data associated with the BoundaryCondition. +""" +function boundary_data end + +""" + discretize_data(grid, bc::BoundaryCondition) + +The data of `bc` as a lazily evaluated grid function on the boundary grid +specified by `boundary(bc)`. +""" +function discretize_data(grid, bc::BoundaryCondition) + return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc)) +end + +""" + DirichletCondition{DT,BID} + +A Dirichlet condition with `data::DT` on the boundary +specified by the boundary identifier `BID`. +""" +struct DirichletCondition{DT,BID} <: BoundaryCondition + data::DT + boundary::BID +end +boundary_data(bc::DirichletCondition) = bc.data +boundary(bc::DirichletCondition) = bc.boundary + +""" + NeumannCondition{DT,BID} + +A Neumann condition with `data::DT` on the boundary +specified by the boundary identifier `BID`. +""" +struct NeumannCondition{DT,BID} <: BoundaryCondition + data::DT + boundary::BID +end +boundary_data(bc::NeumannCondition) = bc.data +boundary(bc::NeumannCondition) = bc.boundary +
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/SbpOperators/boundary_conditions/sat.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/boundary_conditions/sat.jl Tue Jun 25 15:32:19 2024 +0200 @@ -0,0 +1,25 @@ +""" + sat_tensors(op, grid, bc::BoundaryCondition; kwargs...) + +The penalty tensor and boundary operator used to construct a +simultaneous-approximation-term for imposing `bc` related to `op`. + +For `penalty_tensor, L = sat_tensors(...)` then `SAT(u,g) = +penalty_tensor*(L*u - g)` where `g` is the boundary data. +""" +function sat_tensors end + + +""" + sat(op, grid, bc::BoundaryCondition; kwargs...) + +Simultaneous-Approximation-Term for a general `bc` to `op`. Returns a function +`SAT(u,g)` weakly imposing `bc` when added to `op*u`. + +Internally `sat_tensors(op, grid, bc; ...)` is called to construct the +necessary parts for the SAT. +""" +function sat(op, grid, bc::BoundaryCondition; kwargs...) + penalty_tensor, L = sat_tensors(op, grid, bc; kwargs...) + return SAT(u, g) = penalty_tensor*(L*u - g) +end
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/SbpOperators/operators/standard_diagonal.toml --- a/src/SbpOperators/operators/standard_diagonal.toml Sat May 25 16:02:55 2024 -0700 +++ b/src/SbpOperators/operators/standard_diagonal.toml Tue Jun 25 15:32:19 2024 +0200 @@ -34,11 +34,14 @@ {s = ["1", "-2", "1"], c = 1}, ] +D2.positivity = {theta_M = "0.3636363636", theta_R = "1.000000538455350", m_b = "2"} + 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]] order = 4 @@ -65,6 +68,8 @@ {s = [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"], c = 4}, ] +D2.positivity = {theta_M = "0.2505765857", theta_R = "0.577587500088313", m_b = "4"} + D2variable.inner_stencil = [ ["-1/8", "1/6", "-1/8", "0", "0" ], [ "1/6", "1/2", "1/2", "1/6", "0" ], @@ -136,7 +141,6 @@ ] - [[stencil_set]] order = 6 @@ -150,3 +154,5 @@ e.closure = ["1"] d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"] + +D2.positivity = {theta_M = "0.1878687080", theta_R = "0.3697", m_b = "7"}
diff -r 7f5e4d2a5112 -r b62770cec7d0 src/SbpOperators/volumeops/laplace/laplace.jl --- a/src/SbpOperators/volumeops/laplace/laplace.jl Sat May 25 16:02:55 2024 -0700 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Tue Jun 25 15:32:19 2024 +0200 @@ -52,3 +52,78 @@ return Δ end laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set) + +""" + sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) + +The operators required to construct the SAT for imposing a Dirichlet +condition. `H_tuning` and `R_tuning` are used to specify the strength of the +penalty. + +See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref). +""" +function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.) + id = boundary(bc) + set = Δ.stencil_set + H⁻¹ = inverse_inner_product(g,set) + Hᵧ = inner_product(boundary_grid(g, id), set) + e = boundary_restriction(g, set, id) + d = normal_derivative(g, set, id) + B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning) + penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ + return penalty_tensor, e +end + +""" + sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) + +The operators required to construct the SAT for imposing a Neumann condition. + +See also: [`sat`](@ref), [`NeumannCondition`](@ref). +""" +function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition) + id = boundary(bc) + set = Δ.stencil_set + H⁻¹ = inverse_inner_product(g,set) + Hᵧ = inner_product(boundary_grid(g, id), set) + e = boundary_restriction(g, set, id) + d = normal_derivative(g, set, id) + + penalty_tensor = -H⁻¹∘e'∘Hᵧ + return penalty_tensor, d +end + +""" + positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) + +Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive +definite with respect to the boundary quadrature. Here `d` is the normal +derivative and `e` is the boundary restriction operator. `B` can then be used +to form a symmetric and energy stable penalty for a Dirichlet condition. The +parameters `H_tuning` and `R_tuning` are used to specify the strength of the +penalty and must be greater than 1. For details we refer to +https://doi.org/10.1016/j.jcp.2020.109294 +""" +function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) + @assert(H_tuning ≥ 1.) + @assert(R_tuning ≥ 1.) + Nτ_H, τ_R = positivity_limits(Δ,g,bc) + return H_tuning*Nτ_H + R_tuning*τ_R +end + +# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then +# change bc::BoundaryCondition to id::BoundaryIdentifier +function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition) + h = spacing(g) + θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1]) + θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"]) + + τ_H = 1/(h*θ_H) + τ_R = 1/(h*θ_R) + return τ_H, τ_R +end + +function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition) + τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc) + return τ_H*ndims(g), τ_R +end
diff -r 7f5e4d2a5112 -r b62770cec7d0 test/SbpOperators/boundary_conditions/boundary_condition_test.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl Tue Jun 25 15:32:19 2024 +0200 @@ -0,0 +1,39 @@ +using Test + +using Sbplib.Grids +using Sbplib.RegionIndices +using Sbplib.SbpOperators + +@testset "BoundaryCondition" begin + grid_1d = equidistant_grid(0.0, 1.0, 11) + grid_2d = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 15) + grid_3d = equidistant_grid((0.0, 0.0, 0.0), (1.0,1.0, 1.0), 11, 15, 13) + (id_l,_) = boundary_identifiers(grid_1d) + (_,_,_,id_n) = boundary_identifiers(grid_2d) + (_,_,_,_,id_b,_) = boundary_identifiers(grid_3d) + + g = 3.14 + f(x,y,z) = x^2+y^2+z^2 + @testset "Constructors" begin + @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower} + @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,Lower}} + end + + @testset "boundary" begin + @test boundary(DirichletCondition(g,id_l)) == id_l + @test boundary(NeumannCondition(f,id_b)) == id_b + end + + @testset "boundary_data" begin + @test boundary_data(DirichletCondition(g,id_l)) == g + @test boundary_data(NeumannCondition(f,id_b)) == f + end + + @testset "discretize_data" begin + @test fill(g) ≈ discretize_data(grid_1d,DirichletCondition(g,id_l)) + @test g*ones(11,1) ≈ discretize_data(grid_2d,DirichletCondition(g,id_n)) + X = repeat(0:1/10:1, inner = (1,15)) + Y = repeat(0:1/14:1, outer = (1,11)) + @test map((x,y)->f(x,y,0), X,Y') ≈ discretize_data(grid_3d,NeumannCondition(f,id_b)) + end +end
diff -r 7f5e4d2a5112 -r b62770cec7d0 test/SbpOperators/boundary_conditions/sat_test.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SbpOperators/boundary_conditions/sat_test.jl Tue Jun 25 15:32:19 2024 +0200 @@ -0,0 +1,71 @@ +using Test + +using Sbplib.Grids +using Sbplib.LazyTensors +using Sbplib.SbpOperators + +stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) + +struct MockOp end + +function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition; a = 1.) + e = boundary_restriction(g, stencil_set, boundary(bc)) + L = a*e + sat_op = e' + return sat_op, L +end + +function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::NeumannCondition) + e = boundary_restriction(g, stencil_set, boundary(bc)) + d = normal_derivative(g, stencil_set, boundary(bc)) + L = d + sat_op = e' + return sat_op, L +end + +@testset "sat" begin + op = MockOp() + @testset "1D" begin + grid = equidistant_grid(0., 1., 11) + l, r = boundary_identifiers(grid) + u = eval_on(grid, x-> 1. + 2x^2) + dc = DirichletCondition(1.0, l) + g_l = discretize_data(grid, dc) + SAT_l = sat(op, grid, dc) + @test SAT_l(u, g_l) ≈ zeros((size(grid))) atol = 1e-13 + + nc = NeumannCondition(4.0, r) + g_r = discretize_data(grid, nc) + SAT_r = sat(op, grid, nc) + @test SAT_r(u, g_r) ≈ zeros((size(grid))) atol = 1e-13 + end + @testset "2D" begin + grid = equidistant_grid((0.,0.), (1.,1.), 11, 13) + W, E, S, N = boundary_identifiers(grid) + u = eval_on(grid, (x,y) -> x+y^2) + + dc_W = DirichletCondition(1.0, W) + SAT_W = sat(op, grid, dc_W) + g_W = discretize_data(grid, dc_W) + r_W = zeros(size(grid)) + r_W[1,:] .= map(y -> (y^2-1.), range(0., 1., length=13)) + @test SAT_W(u, g_W) ≈ r_W atol = 1e-13 + + dc_E = DirichletCondition(2, E) + SAT_E = sat(op, grid, dc_E; a = 2.) + g_E = discretize_data(grid, dc_E) + r_E = zeros(size(grid)) + r_E[end,:] .= map(y -> (2*(1. + y^2)-2.), range(0., 1., length=13)) + @test SAT_E(u, g_E) ≈ r_E atol = 1e-13 + + nc_S = NeumannCondition(.0, S) + SAT_S = sat(op, grid, nc_S) + g_S = discretize_data(grid, nc_S) + @test SAT_S(u, g_S) ≈ zeros(size(grid)) atol = 1e-13 + + nc_N = NeumannCondition(2.0, N) + SAT_N = sat(op, grid, nc_N) + g_N = discretize_data(grid, nc_N) + @test SAT_N(u, g_N) ≈ zeros(size(grid)) atol = 1e-13 + end +end
diff -r 7f5e4d2a5112 -r b62770cec7d0 test/SbpOperators/volumeops/laplace/laplace_test.jl --- a/test/SbpOperators/volumeops/laplace/laplace_test.jl Sat May 25 16:02:55 2024 -0700 +++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl Tue Jun 25 15:32:19 2024 +0200 @@ -88,3 +88,66 @@ end end +@testset "sat_tensors" begin + # TODO: The following tests should be implemented + # 1. Symmetry D'H == H'D (test_broken below) + # 2. Test eigenvalues of and/or solution to Poisson + # 3. Test tuning of Dirichlet conditions + # + # These tests are likely easiest to implement once + # we have support for generating matrices from tensors. + + operator_path = sbp_operators_path()*"standard_diagonal.toml" + orders = (2,4) + tols = (5e-2,5e-4) + sz = (201,401) + g = equidistant_grid((0.,0.), (1.,1.), sz...) + + # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO) + # of the operator D = Δ + SAT, where SAT is the tensor composition of the + # operators from sat_tensor. Note that SAT*u should approximate 0 for the + # conditions chosen. + + @testset "Dirichlet" begin + for (o, tol) ∈ zip(orders,tols) + stencil_set = read_stencil_set(operator_path; order=o) + Δ = Laplace(g, stencil_set) + H = inner_product(g, stencil_set) + u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y))) + Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y))) + D = Δ + for id ∈ boundary_identifiers(g) + D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id))) + end + e = D*u .- Δu + # Accuracy + @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol + # Symmetry + r = randn(size(u)) + @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D. + end + end + + @testset "Neumann" begin + @testset "Dirichlet" begin + for (o, tol) ∈ zip(orders,tols) + stencil_set = read_stencil_set(operator_path; order=o) + Δ = Laplace(g, stencil_set) + H = inner_product(g, stencil_set) + u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y))) + Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y))) + D = Δ + for id ∈ boundary_identifiers(g) + D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id))) + end + e = D*u .- Δu + # Accuracy + @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol + # Symmetry + r = randn(size(u)) + @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D. + end + end + end +end +