From 3db9929af442baf0b57d9c5bf1c8cba0a8272e03 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 6 Dec 2024 14:11:56 +0100 Subject: [PATCH 01/28] =?UTF-8?q?add=20reflecting=20boundary=20conditions?= =?UTF-8?q?=20for=20Sv=C3=A4rd-Kalisch=20equations?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../svaerd_kalisch_1d_basic_reflecting.jl | 49 +++++ src/DispersiveShallowWater.jl | 2 +- src/equations/bbm_bbm_1d.jl | 12 +- src/equations/equations.jl | 7 + src/equations/svaerd_kalisch_1d.jl | 179 ++++++++++++++++-- test/test_svaerd_kalisch_1d.jl | 32 ++++ 6 files changed, 257 insertions(+), 24 deletions(-) create mode 100644 examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl new file mode 100644 index 00000000..8d08b9fd --- /dev/null +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl @@ -0,0 +1,49 @@ +using OrdinaryDiffEqTsit5 +using DispersiveShallowWater +using SummationByPartsOperators: MattssonNordström2004, derivative_operator + +############################################################################### +# Semidiscretization of the Svärd-Kalisch equations + +equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0, + alpha = 0.0, beta = 1 / 3, gamma = 0.0) + +initial_condition = initial_condition_manufactured_reflecting +source_terms = source_terms_manufactured_reflecting +boundary_conditions = boundary_condition_reflecting + +# create homogeneous mesh +coordinates_min = 0.0 +coordinates_max = 1.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with SBP operators of accuracy order 4 +accuracy_order = 4 +D1 = derivative_operator(MattssonNordström2004(), + derivative_order = 1, accuracy_order = accuracy_order, + xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N) +D2 = derivative_operator(MattssonNordström2004(), + derivative_order = 2, accuracy_order = accuracy_order, + xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N) +solver = Solver(D1, D2) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions, + source_terms = source_terms) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi; interval = 100, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_callback, summary_callback) + +saveat = range(tspan..., length = 100) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 54633bf0..12237564 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -19,7 +19,7 @@ using BandedMatrices: BandedMatrix using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation -using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, cholesky, cholesky!, +using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, lu!, cholesky, cholesky!, issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 8f2f8610..e7836e0b 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -338,8 +338,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D{BathymetryFlat}, # homogeneous Neumann boundary conditions if solver.D1 isa DerivativeOperator || - solver.D1 isa UniformCoupledOperator || - solver.D1 isa FourierDerivativeOperator + solver.D1 isa UniformCoupledOperator D1_b = BandedMatrix(solver.D1) invImD2n = lu(I + 1 / 6 * D^2 * inv(M) * D1_b' * PdM * D1_b) elseif solver.D1 isa UpwindOperators @@ -382,8 +381,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D, # homogeneous Neumann boundary conditions if solver.D1 isa DerivativeOperator || - solver.D1 isa UniformCoupledOperator || - solver.D1 isa FourierDerivativeOperator + solver.D1 isa UniformCoupledOperator D1_b = BandedMatrix(solver.D1) invImD2n = lu(I + 1 / 6 * inv(M) * D1_b' * PdM * K * D1_b) elseif solver.D1 isa UpwindOperators @@ -434,8 +432,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition, end # energy and mass conservative semidiscretization if solver.D1 isa PeriodicDerivativeOperator || - solver.D1 isa UniformPeriodicCoupledOperator || - solver.D1 isa FourierDerivativeOperator + solver.D1 isa UniformPeriodicCoupledOperator @trixi_timeit timer() "deta hyperbolic" begin mul!(deta, solver.D1, tmp1) end @@ -501,8 +498,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition, end # energy and mass conservative semidiscretization if solver.D1 isa DerivativeOperator || - solver.D1 isa UniformCoupledOperator || - solver.D1 isa FourierDerivativeOperator + solver.D1 isa UniformCoupledOperator @trixi_timeit timer() "deta hyperbolic" begin mul!(deta, solver.D1, tmp1) end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8cb78ed0..a8c677e2 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -581,6 +581,13 @@ function solve_system_matrix!(dv, system_matrix, rhs, return nothing end +# Svärd-Kalisch equations for reflecting boundary conditions uses LU factorization +function solve_system_matrix!(dv, system_matrix, ::SvaerdKalischEquations1D, cache) + (; factorization) = cache + lu!(factorization, system_matrix) + ldiv!(factorization, dv) +end + function solve_system_matrix!(dv, system_matrix, ::Union{BBMEquation1D, BBMBBMEquations1D}) ldiv!(system_matrix, dv) end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 964a509f..d295cf08 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -140,6 +140,79 @@ function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D) return SVector(dq1, dq2, zero(dq1)) end +""" + initial_condition_manufactured_reflecting(x, t, equations::SvaerdKalischEquations1D, mesh) + +A smooth manufactured solution for reflecting boundary conditions in combination +with [`source_terms_manufactured_reflecting`](@ref). +""" +function initial_condition_manufactured_reflecting(x, t, + equations::SvaerdKalischEquations1D, + mesh) + eta = exp(2 * t) * cospi(x) + v = exp(t) * x * sinpi(x) + b = -5 - 2 * cospi(2 * x) + D = equations.eta0 - b + return SVector(eta, v, D) +end + +""" + source_terms_manufactured_reflecting(q, x, t, equations::SvaerdKalischEquations1D, mesh) + +A smooth manufactured solution for reflecting boundary conditions in combination +with [`initial_condition_manufactured_reflecting`](@ref). +""" +function source_terms_manufactured_reflecting(q, x, t, equations::SvaerdKalischEquations1D) + g = gravity_constant(equations) + eta0 = still_water_surface(q, equations) + beta = equations.beta + a1 = sinpi(2 * x) + a2 = cospi(2 * x) + a9 = exp(t) + a11 = eta0 + 2.0 * a2 + 5.0 + a13 = 0.2 * eta0 + 0.4 * a2 + 1 + a22 = sinpi(x) + a23 = cospi(x) + a24 = exp(2 * t) + a25 = a24 * a23 + a26 = a25 + 2.0 * a2 + 5.0 + a27 = a9 * a22 + a28 = pi * x * a9 * a23 + a27 + a29 = -pi * a24 * a22 - 4.0 * pi * a1 + a30 = a26 * a24 + a31 = a26 * a27 + dq1 = x * a29 * a27 + pi * x * a26 * a9 * a23 + a31 + 2 * a25 + dq2 = 100.0 * pi * beta * a28 * a13^2 * a1 + 40.0 * pi * beta * a28 * a13 * a11 * a1 - + 25.0 * beta * (-pi^2 * x * a27 + 2 * pi * a9 * a23) * a13^2 * a11 - + pi * g * a30 * a22 + x^2 * a29 * a24 * a22^2 / 2 + pi * x^2 * a30 * a22 * a23 + + x * a28 * a31 / 2 + x * a30 * a22^2 + x * a31 - + x * (x * a29 * a27 + pi * x * a26 * a9 * a23 + a31) * a27 / 2 + + return SVector(dq1, dq2, zero(dq1)) +end + +# For periodic boundary conditions +function assemble_system_matrix!(cache, h, D1, + ::SvaerdKalischEquations1D) + (; M_h, minus_MD1betaD1) = cache + + @.. M_h = h + scale_by_mass_matrix!(M_h, D1) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + return Symmetric(Diagonal(M_h) + minus_MD1betaD1) +end + +# For reflecting boundary conditions +function assemble_system_matrix!(cache, h, + ::SvaerdKalischEquations1D) + (; D1betaD1d) = cache + return Diagonal((@view h[2:(end - 1)])) - D1betaD1d +end + function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, ::BoundaryConditionPeriodic, @@ -185,7 +258,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat system_matrix = let cache = (; M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, - D1, D1mat, equations) + D1, equations) end elseif D1 isa PeriodicUpwindOperators D1_central = D1.central @@ -193,7 +266,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, minus_MD1betaD1 = D1mat_minus' * Diagonal(M_beta) * D1mat_minus system_matrix = let cache = (; M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, - D1, D1mat_minus, equations) + D1, equations) end else throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) @@ -211,18 +284,48 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, return cache end -function assemble_system_matrix!(cache, h, D1, D1mat, - ::SvaerdKalischEquations1D) - (; M_h, minus_MD1betaD1) = cache - - @.. M_h = h - scale_by_mass_matrix!(M_h, D1) - - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - return Symmetric(Diagonal(M_h) + minus_MD1betaD1) +# Reflecting boundary conditions assume alpha = gamma = 0 +function create_cache(mesh, equations::SvaerdKalischEquations1D, + solver, initial_condition, + ::BoundaryConditionReflecting, + RealT, uEltype) + if !isapprox(equations.alpha, 0, atol = 1e-14) || + !isapprox(equations.gamma, 0, atol = 1e-14) + throw(ArgumentError("Reflecting boundary conditions for Svärd-Kalisch equations only implemented for alpha = gamma = 0")) + end + D1 = solver.D1 + N = nnodes(mesh) + # Assume D is independent of time and compute D evaluated at mesh points once. + D = Array{RealT}(undef, N) + x = grid(solver) + for i in eachnode(solver) + D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations) + end + h = ones(RealT, N) + hv = zero(h) + b = zero(h) + eta_x = zero(h) + v_x = zero(h) + h_v_x = zero(h) + hv2_x = zero(h) + beta_hat = equations.beta * D .^ 3 + if D1 isa DerivativeOperator || + D1 isa UniformCoupledOperator + D1_central = D1 + D1mat = sparse(D1_central) + Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2)) + D1betaD1d = sparse(D1mat * Diagonal(beta_hat) * D1mat * Pd)[2:(end - 1), :] + system_matrix = let cache = (; D1betaD1d) + assemble_system_matrix!(cache, h, equations) + end + factorization = lu(system_matrix) + else + throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) + end + factorization = lu(system_matrix) + cache = (; factorization, D1betaD1d, D, h, hv, b, eta_x, v_x, + h_v_x, hv2_x, beta_hat, D1_central, D1) + return cache end # Discretization that conserves @@ -315,7 +418,7 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver) @trixi_timeit timer() "assemble system matrix" begin - system_matrix = assemble_system_matrix!(cache, h, D1, D1_central, equations) + system_matrix = assemble_system_matrix!(cache, h, D1, equations) end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv @@ -326,6 +429,52 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, return nothing end +function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, + initial_condition, ::BoundaryConditionReflecting, source_terms, + solver, cache) + (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central, D1) = cache + + g = gravity_constant(equations) + eta, v = q.x + deta, dv, dD = dq.x + fill!(dD, zero(eltype(dD))) + + @trixi_timeit timer() "deta hyperbolic" begin + @.. b = equations.eta0 - D + @.. h = eta - b + @.. hv = h * v + + mul!(deta, D1_central, -hv) + end + + # split form + @trixi_timeit timer() "dv hyperbolic" begin + mul!(eta_x, D1_central, eta) + mul!(v_x, D1_central, v) + mul!(h_v_x, D1_central, hv) + @.. tmp1 = hv * v + mul!(hv2_x, D1_central, tmp1) + @.. dv = -(0.5 * (hv2_x + hv * v_x - v * h_v_x) + + g * h * eta_x) + end + + # no split form + # dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ + # g * h .* eta_x) + + @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, + solver) + @trixi_timeit timer() "assemble system matrix" begin + system_matrix = assemble_system_matrix!(cache, h, equations) + end + @trixi_timeit timer() "solving elliptic system" begin + solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, equations, cache) + dv[1] = dv[end] = zero(eltype(dv)) + end + + return nothing +end + # The modified entropy/energy takes the whole `q` for every point in space """ DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D, cache) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 5b919b9a..00306113 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -17,6 +17,38 @@ end @test_allocations(semi, sol, allocs=90_000) end +@testitem "svaerd_kalisch_1d_basic_reflecting" setup=[Setup, SvaerdKalischEquations1D, AdditionalImports] begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), + tspan=(0.0, 1.0), + l2=[4.067039325687744e-5 6.705583873243185e-8 0.0], + linf=[0.00028373482419530305 1.4602668701318988e-7 0.0], + cons_error=[1.0484869654379814e-9 0.5469460930247622 0.0], + change_waterheight=1.0484869654379814e-9, + change_velocity=0.5469460930247622, + change_entropy=14.059432069002527) + + @test_allocations(semi, sol, allocs=650_000) + + # # test upwind operators + # D1 = upwind_operators(Mattsson2017; derivative_order = 1, + # accuracy_order = accuracy_order, xmin = mesh.xmin, + # xmax = mesh.xmax, + # N = mesh.N) + # D2 = sparse(D1.plus) * sparse(D1.minus) + # solver = Solver(D1, D2) + # @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), + # tspan=(0.0, 1.0), + # solver=solver, + # l2=[0.002345799818043513 3.254313503127441e-8 0.0], + # linf=[0.05625950533062252 6.815531732318192e-7 0.0], + # cons_error=[1.6607871307809518e-9 0.5469460993745239 0.0], + # change_waterheight=-1.6607871307809518e-9, + # change_velocity=0.5469460993745239, + # change_entropy=132.10938489083918) + + # @test_allocations(semi, sol, allocs=2_000) +end + @testitem "svaerd_kalisch_1d_dingemans" setup=[ Setup, SvaerdKalischEquations1D, From 94dad028ec0f8f09560a92a19f666e4b819525da Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Fri, 6 Dec 2024 14:28:31 +0100 Subject: [PATCH 02/28] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/DispersiveShallowWater.jl | 3 ++- test/test_svaerd_kalisch_1d.jl | 6 +++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 12237564..2d70b67a 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -19,7 +19,8 @@ using BandedMatrices: BandedMatrix using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation -using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, lu!, cholesky, cholesky!, +using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, lu!, cholesky, + cholesky!, issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 00306113..e1185eae 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -17,7 +17,11 @@ end @test_allocations(semi, sol, allocs=90_000) end -@testitem "svaerd_kalisch_1d_basic_reflecting" setup=[Setup, SvaerdKalischEquations1D, AdditionalImports] begin +@testitem "svaerd_kalisch_1d_basic_reflecting" setup=[ + Setup, + SvaerdKalischEquations1D, + AdditionalImports +] begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), l2=[4.067039325687744e-5 6.705583873243185e-8 0.0], From 0d5ca48c717e2833d4e300c32cce43a9a027dc8e Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 6 Dec 2024 14:39:32 +0100 Subject: [PATCH 03/28] add back accidently removed FourierDerivativeOperator --- src/equations/bbm_bbm_1d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index e7836e0b..99cfd788 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -432,7 +432,8 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition, end # energy and mass conservative semidiscretization if solver.D1 isa PeriodicDerivativeOperator || - solver.D1 isa UniformPeriodicCoupledOperator + solver.D1 isa UniformPeriodicCoupledOperator || + solver.D1 isa FourierDerivativeOperator @trixi_timeit timer() "deta hyperbolic" begin mul!(deta, solver.D1, tmp1) end From c7f3d4e3bc9c68ff187379d0fed8d3a353a54d5a Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Sun, 8 Dec 2024 15:38:31 +0100 Subject: [PATCH 04/28] also implement upwind operators --- .../svaerd_kalisch_1d_basic_reflecting.jl | 9 ++--- src/equations/svaerd_kalisch_1d.jl | 22 +++++++---- test/test_svaerd_kalisch_1d.jl | 37 +++++++++---------- 3 files changed, 34 insertions(+), 34 deletions(-) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl index 8d08b9fd..0590c419 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl @@ -4,7 +4,7 @@ using SummationByPartsOperators: MattssonNordström2004, derivative_operator ############################################################################### # Semidiscretization of the Svärd-Kalisch equations - +# For reflecting boundary conditions, alpha and gamma need to be 0 equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0, alpha = 0.0, beta = 1 / 3, gamma = 0.0) @@ -23,10 +23,7 @@ accuracy_order = 4 D1 = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = accuracy_order, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N) -D2 = derivative_operator(MattssonNordström2004(), - derivative_order = 2, accuracy_order = accuracy_order, - xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N) -solver = Solver(D1, D2) +solver = Solver(D1, nothing) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, @@ -41,7 +38,7 @@ summary_callback = SummaryCallback() analysis_callback = AnalysisCallback(semi; interval = 100, extra_analysis_errors = (:conservation_error,), extra_analysis_integrals = (waterheight_total, - velocity, entropy)) + entropy_modified)) callbacks = CallbackSet(analysis_callback, summary_callback) saveat = range(tspan..., length = 100) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index d295cf08..b681f23f 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -233,7 +233,6 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, alpha_eta_x_x = zero(h) y_x = zero(h) v_y_x = zero(h) - yv_x = zero(h) vy = zero(h) vy_x = zero(h) y_v_x = zero(h) @@ -273,7 +272,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, end factorization = cholesky(system_matrix) cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x, - alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx, + alpha_eta_x_x, y_x, v_y_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, alpha_hat, beta_hat, gamma_hat, tmp2, D1_central, M, D1, M_h) if D1 isa PeriodicUpwindOperators @@ -309,16 +308,22 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, h_v_x = zero(h) hv2_x = zero(h) beta_hat = equations.beta * D .^ 3 + Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2)) if D1 isa DerivativeOperator || D1 isa UniformCoupledOperator D1_central = D1 D1mat = sparse(D1_central) - Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2)) D1betaD1d = sparse(D1mat * Diagonal(beta_hat) * D1mat * Pd)[2:(end - 1), :] system_matrix = let cache = (; D1betaD1d) assemble_system_matrix!(cache, h, equations) end - factorization = lu(system_matrix) + elseif D1 isa UpwindOperators + D1_central = D1.central + D1betaD1d = sparse(sparse(D1.plus) * Diagonal(beta_hat) * + sparse(D1.minus) * Pd)[2:(end - 1), :] + system_matrix = let cache = (; D1betaD1d) + assemble_system_matrix!(cache, h, equations) + end else throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) end @@ -340,9 +345,9 @@ end function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, initial_condition, ::BoundaryConditionPeriodic, source_terms, solver, cache) - (; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x, + (; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, alpha_hat, gamma_hat, - tmp1, tmp2, D1_central, M, D1) = cache + tmp1, tmp2, D1_central, D1) = cache g = gravity_constant(equations) eta, v = q.x @@ -432,7 +437,7 @@ end function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, initial_condition, ::BoundaryConditionReflecting, source_terms, solver, cache) - (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central, D1) = cache + (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central) = cache g = gravity_constant(equations) eta, v = q.x @@ -507,7 +512,8 @@ See also [`energy_total_modified`](@ref). @.. b = equations.eta0 - D @.. h = eta - b - if D1 isa PeriodicUpwindOperators + if D1 isa PeriodicUpwindOperators || + D1 isa UpwindOperators mul!(v_x, D1.minus, v) else mul!(v_x, D1, v) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index e1185eae..169d629a 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -28,29 +28,26 @@ end linf=[0.00028373482419530305 1.4602668701318988e-7 0.0], cons_error=[1.0484869654379814e-9 0.5469460930247622 0.0], change_waterheight=1.0484869654379814e-9, - change_velocity=0.5469460930247622, - change_entropy=14.059432069002527) + change_entropy_modified=459.90372362418947) @test_allocations(semi, sol, allocs=650_000) - # # test upwind operators - # D1 = upwind_operators(Mattsson2017; derivative_order = 1, - # accuracy_order = accuracy_order, xmin = mesh.xmin, - # xmax = mesh.xmax, - # N = mesh.N) - # D2 = sparse(D1.plus) * sparse(D1.minus) - # solver = Solver(D1, D2) - # @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), - # tspan=(0.0, 1.0), - # solver=solver, - # l2=[0.002345799818043513 3.254313503127441e-8 0.0], - # linf=[0.05625950533062252 6.815531732318192e-7 0.0], - # cons_error=[1.6607871307809518e-9 0.5469460993745239 0.0], - # change_waterheight=-1.6607871307809518e-9, - # change_velocity=0.5469460993745239, - # change_entropy=132.10938489083918) - - # @test_allocations(semi, sol, allocs=2_000) + # test upwind operators + D1 = upwind_operators(Mattsson2017; derivative_order = 1, + accuracy_order = accuracy_order, xmin = mesh.xmin, + xmax = mesh.xmax, + N = mesh.N) + solver = Solver(D1, nothing) + @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), + tspan=(0.0, 1.0), + solver=solver, + l2=[5.1845375611538653e-5 4.111955051600758e-9 0.0], + linf=[0.0003710888676375923 8.797788351305735e-8 0.0], + cons_error=[1.7004244096716813e-9 0.5469460935005923 0.0], + change_waterheight=-1.7004244096716813e-9, + change_entropy_modified=459.9037221456176) + + @test_allocations(semi, sol, allocs=650_000) end @testitem "svaerd_kalisch_1d_dingemans" setup=[ From 9e9d472ad3726d15c7ccf6ae8fcd8177bef1e8db Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Sun, 8 Dec 2024 15:40:27 +0100 Subject: [PATCH 05/28] reformat --- src/DispersiveShallowWater.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 2d70b67a..a8c1f41d 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -19,9 +19,8 @@ using BandedMatrices: BandedMatrix using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation -using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, lu!, cholesky, - cholesky!, - issuccess +using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, + lu, lu!, cholesky, cholesky!, issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @recipe, @series From 874028b7346fef9013e3a292ead30ae688e4cba6 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Tue, 10 Dec 2024 12:38:31 +0100 Subject: [PATCH 06/28] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/equations/svaerd_kalisch_1d.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index b681f23f..1184ca09 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -288,8 +288,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, ::BoundaryConditionReflecting, RealT, uEltype) - if !isapprox(equations.alpha, 0, atol = 1e-14) || - !isapprox(equations.gamma, 0, atol = 1e-14) + if !iszero(equations.alpha) || !iszero(equations.gamma,) throw(ArgumentError("Reflecting boundary conditions for Svärd-Kalisch equations only implemented for alpha = gamma = 0")) end D1 = solver.D1 @@ -307,7 +306,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, v_x = zero(h) h_v_x = zero(h) hv2_x = zero(h) - beta_hat = equations.beta * D .^ 3 + beta_hat = equations.beta .* D .^ 3 Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2)) if D1 isa DerivativeOperator || D1 isa UniformCoupledOperator @@ -437,6 +436,8 @@ end function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, initial_condition, ::BoundaryConditionReflecting, source_terms, solver, cache) + # When constructing the `cache`, we assert that alpha and gamma are zero. + # We use this explicitly in the code below. (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central) = cache g = gravity_constant(equations) @@ -463,10 +464,6 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, g * h * eta_x) end - # no split form - # dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ - # g * h .* eta_x) - @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver) @trixi_timeit timer() "assemble system matrix" begin From 25ad669685329f0acf0eedbef99ed816f9646647 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Tue, 10 Dec 2024 12:39:37 +0100 Subject: [PATCH 07/28] Update src/equations/svaerd_kalisch_1d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/svaerd_kalisch_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 1184ca09..01ed7bcc 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -288,7 +288,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, ::BoundaryConditionReflecting, RealT, uEltype) - if !iszero(equations.alpha) || !iszero(equations.gamma,) + if !iszero(equations.alpha) || !iszero(equations.gamma) throw(ArgumentError("Reflecting boundary conditions for Svärd-Kalisch equations only implemented for alpha = gamma = 0")) end D1 = solver.D1 From 8fe7401ce3600a243fe764d604f7874ffdd52646 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 10 Dec 2024 16:40:26 +0100 Subject: [PATCH 08/28] set stricter time integration tolerances --- .../svaerd_kalisch_1d_basic_reflecting.jl | 2 +- test/test_svaerd_kalisch_1d.jl | 20 +++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl index 0590c419..a6edb3dd 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl @@ -42,5 +42,5 @@ analysis_callback = AnalysisCallback(semi; interval = 100, callbacks = CallbackSet(analysis_callback, summary_callback) saveat = range(tspan..., length = 100) -sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, +sol = solve(ode, Tsit5(), abstol = 1e-12, reltol = 1e-12, save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 169d629a..c24d762d 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -24,11 +24,11 @@ end ] begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), - l2=[4.067039325687744e-5 6.705583873243185e-8 0.0], - linf=[0.00028373482419530305 1.4602668701318988e-7 0.0], - cons_error=[1.0484869654379814e-9 0.5469460930247622 0.0], - change_waterheight=1.0484869654379814e-9, - change_entropy_modified=459.90372362418947) + l2=[5.402315467757905e-6 6.705583190754449e-8 0.0], + linf=[9.787584609100008e-5 1.4602668577112787e-7 0.0], + cons_error=[1.0484865187415942e-9 0.5469460930247753 0.0], + change_waterheight=1.0484865187415942e-9, + change_entropy_modified=459.9037236233692) @test_allocations(semi, sol, allocs=650_000) @@ -41,11 +41,11 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), solver=solver, - l2=[5.1845375611538653e-5 4.111955051600758e-9 0.0], - linf=[0.0003710888676375923 8.797788351305735e-8 0.0], - cons_error=[1.7004244096716813e-9 0.5469460935005923 0.0], - change_waterheight=-1.7004244096716813e-9, - change_entropy_modified=459.9037221456176) + l2=[5.862278093002509e-6 4.111933117101271e-9 0.0], + linf=[3.135229084794133e-5 8.797788287467911e-8 0.0], + cons_error=[1.700425282207037e-9 0.5469460935004404 0.0], + change_waterheight=-1.700425282207037e-9, + change_entropy_modified=459.9037221442477) @test_allocations(semi, sol, allocs=650_000) end From 75ec37f9d563931b963dd32eb1e678d5b885055e Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 10 Dec 2024 17:04:20 +0100 Subject: [PATCH 09/28] lower tolerance --- test/test_svaerd_kalisch_1d.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index c24d762d..9a919506 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -28,7 +28,8 @@ end linf=[9.787584609100008e-5 1.4602668577112787e-7 0.0], cons_error=[1.0484865187415942e-9 0.5469460930247753 0.0], change_waterheight=1.0484865187415942e-9, - change_entropy_modified=459.9037236233692) + change_entropy_modified=459.9037236233692, + atol=1e-11) # to make CI pass @test_allocations(semi, sol, allocs=650_000) @@ -45,7 +46,8 @@ end linf=[3.135229084794133e-5 8.797788287467911e-8 0.0], cons_error=[1.700425282207037e-9 0.5469460935004404 0.0], change_waterheight=-1.700425282207037e-9, - change_entropy_modified=459.9037221442477) + change_entropy_modified=459.9037221442477, + atol=1e-11) # to make CI pass @test_allocations(semi, sol, allocs=650_000) end From 8859d06e8b57a35e86fc4c6c3aec2d583fa99ccd Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 10:29:08 +0100 Subject: [PATCH 10/28] use Cholesky decomposition --- src/DispersiveShallowWater.jl | 2 +- src/equations/equations.jl | 7 --- src/equations/svaerd_kalisch_1d.jl | 84 ++++++++++++++++++------------ 3 files changed, 52 insertions(+), 41 deletions(-) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index a8c1f41d..83801d31 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -20,7 +20,7 @@ using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, - lu, lu!, cholesky, cholesky!, issuccess + lu, cholesky, cholesky!, issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @recipe, @series diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a8c677e2..8cb78ed0 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -581,13 +581,6 @@ function solve_system_matrix!(dv, system_matrix, rhs, return nothing end -# Svärd-Kalisch equations for reflecting boundary conditions uses LU factorization -function solve_system_matrix!(dv, system_matrix, ::SvaerdKalischEquations1D, cache) - (; factorization) = cache - lu!(factorization, system_matrix) - ldiv!(factorization, dv) -end - function solve_system_matrix!(dv, system_matrix, ::Union{BBMEquation1D, BBMBBMEquations1D}) ldiv!(system_matrix, dv) end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 82b620a5..e746743c 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -192,9 +192,10 @@ function source_terms_manufactured_reflecting(q, x, t, equations::SvaerdKalischE end # For periodic boundary conditions -function assemble_system_matrix!(cache, h, D1, - ::SvaerdKalischEquations1D) - (; M_h, minus_MD1betaD1) = cache +function assemble_system_matrix!(cache, h, + ::SvaerdKalischEquations1D, + ::BoundaryConditionPeriodic) + (; D1, M_h, minus_MD1betaD1) = cache @.. M_h = h scale_by_mass_matrix!(M_h, D1) @@ -208,14 +209,23 @@ end # For reflecting boundary conditions function assemble_system_matrix!(cache, h, - ::SvaerdKalischEquations1D) - (; D1betaD1d) = cache - return Diagonal((@view h[2:(end - 1)])) - D1betaD1d + ::SvaerdKalischEquations1D, + ::BoundaryConditionReflecting) + (; D1, M_h, minus_MD1betaD1) = cache + + @.. M_h = h + scale_by_mass_matrix!(M_h, D1) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + return Symmetric(Diagonal((@view M_h[2:(end - 1)])) + minus_MD1betaD1) end function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, - ::BoundaryConditionPeriodic, + boundary_conditions::BoundaryConditionPeriodic, RealT, uEltype) D1 = solver.D1 # Assume D is independent of time and compute D evaluated at mesh points once. @@ -245,8 +255,8 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, beta_hat = equations.beta * D .^ 3 gamma_hat = equations.gamma * sqrt.(g * D) .* D .^ 3 tmp2 = zero(h) - M = mass_matrix(D1) - M_h = zero(h) + M_h = copy(h) + scale_by_mass_matrix!(M_h, D1) M_beta = copy(beta_hat) scale_by_mass_matrix!(M_beta, D1) if D1 isa PeriodicDerivativeOperator || @@ -255,17 +265,17 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, D1_central = D1 D1mat = sparse(D1_central) minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat - system_matrix = let cache = (; M_h, minus_MD1betaD1) + system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, - D1, equations) + equations, boundary_conditions) end elseif D1 isa PeriodicUpwindOperators D1_central = D1.central D1mat_minus = sparse(D1.minus) minus_MD1betaD1 = D1mat_minus' * Diagonal(M_beta) * D1mat_minus - system_matrix = let cache = (; M_h, minus_MD1betaD1) + system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, - D1, equations) + equations, boundary_conditions) end else throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) @@ -274,7 +284,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, - alpha_hat, beta_hat, gamma_hat, tmp2, D1_central, M, D1, M_h) + alpha_hat, beta_hat, gamma_hat, tmp2, D1_central, D1, M_h) if D1 isa PeriodicUpwindOperators eta_x_upwind = zero(h) v_x_upwind = zero(h) @@ -286,7 +296,7 @@ end # Reflecting boundary conditions assume alpha = gamma = 0 function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, - ::BoundaryConditionReflecting, + boundary_conditions::BoundaryConditionReflecting, RealT, uEltype) if !iszero(equations.alpha) || !iszero(equations.gamma) throw(ArgumentError("Reflecting boundary conditions for Svärd-Kalisch equations only implemented for alpha = gamma = 0")) @@ -307,28 +317,33 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, h_v_x = zero(h) hv2_x = zero(h) beta_hat = equations.beta .* D .^ 3 + M_h = copy(h) + scale_by_mass_matrix!(M_h, D1) + M_beta = copy(beta_hat) + scale_by_mass_matrix!(M_beta, D1) Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2)) if D1 isa DerivativeOperator || D1 isa UniformCoupledOperator D1_central = D1 D1mat = sparse(D1_central) - D1betaD1d = sparse(D1mat * Diagonal(beta_hat) * D1mat * Pd)[2:(end - 1), :] - system_matrix = let cache = (; D1betaD1d) - assemble_system_matrix!(cache, h, equations) + minus_MD1betaD1 = sparse(D1mat' * Diagonal(M_beta) * D1mat * Pd)[2:(end - 1), :] + system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) + assemble_system_matrix!(cache, h, equations, boundary_conditions) end elseif D1 isa UpwindOperators D1_central = D1.central - D1betaD1d = sparse(sparse(D1.plus) * Diagonal(beta_hat) * - sparse(D1.minus) * Pd)[2:(end - 1), :] - system_matrix = let cache = (; D1betaD1d) - assemble_system_matrix!(cache, h, equations) + D1mat_minus = sparse(D1.minus) + minus_MD1betaD1 = sparse(D1mat_minus' * Diagonal(M_beta) * D1mat_minus * Pd)[2:(end - 1), + :] + system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) + assemble_system_matrix!(cache, h, equations, boundary_conditions) end else throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) end - factorization = lu(system_matrix) - cache = (; factorization, D1betaD1d, D, h, hv, b, eta_x, v_x, - h_v_x, hv2_x, beta_hat, D1_central, D1) + factorization = cholesky(system_matrix) + cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x, + h_v_x, hv2_x, beta_hat, D1_central, D1, M_h) return cache end @@ -342,8 +357,8 @@ end # [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669) # TODO: Simplify for the case of flat bathymetry and use higher-order operators function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, - initial_condition, ::BoundaryConditionPeriodic, source_terms, - solver, cache) + initial_condition, boundary_conditions::BoundaryConditionPeriodic, + source_terms, solver, cache) (; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, alpha_hat, gamma_hat, tmp1, tmp2, D1_central, D1) = cache @@ -422,7 +437,8 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver) @trixi_timeit timer() "assemble system matrix" begin - system_matrix = assemble_system_matrix!(cache, h, D1, equations) + system_matrix = assemble_system_matrix!(cache, h, equations, + boundary_conditions) end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv @@ -434,11 +450,11 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, - initial_condition, ::BoundaryConditionReflecting, source_terms, - solver, cache) + initial_condition, boundary_conditions::BoundaryConditionReflecting, + source_terms, solver, cache) # When constructing the `cache`, we assert that alpha and gamma are zero. # We use this explicitly in the code below. - (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central) = cache + (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central, D1) = cache g = gravity_constant(equations) eta, v = q.x @@ -467,10 +483,12 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver) @trixi_timeit timer() "assemble system matrix" begin - system_matrix = assemble_system_matrix!(cache, h, equations) + system_matrix = assemble_system_matrix!(cache, h, equations, boundary_conditions) end @trixi_timeit timer() "solving elliptic system" begin - solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, equations, cache) + tmp1 .= dv + solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, (@view tmp1[2:(end - 1)]), equations, D1, + cache) dv[1] = dv[end] = zero(eltype(dv)) end From 40468f3650fb5436f90a3abfbad23ab248c44ffd Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 11 Dec 2024 10:29:55 +0100 Subject: [PATCH 11/28] Update src/equations/svaerd_kalisch_1d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/svaerd_kalisch_1d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index e746743c..b2a65da0 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -487,7 +487,8 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv - solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, (@view tmp1[2:(end - 1)]), equations, D1, + solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, + (@view tmp1[2:(end - 1)]), equations, D1, cache) dv[1] = dv[end] = zero(eltype(dv)) end From 40f52cb9e6eb8e2d7f53f9fd2ef91ea0fb38ec18 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 12:11:43 +0100 Subject: [PATCH 12/28] temporarily add lu decomposition again --- src/DispersiveShallowWater.jl | 2 +- src/equations/equations.jl | 7 ++++ src/equations/svaerd_kalisch_1d.jl | 60 +++++++++++++++++++++++++++++- 3 files changed, 67 insertions(+), 2 deletions(-) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 83801d31..a8c1f41d 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -20,7 +20,7 @@ using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, - lu, cholesky, cholesky!, issuccess + lu, lu!, cholesky, cholesky!, issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @recipe, @series diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8cb78ed0..85e098ef 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -581,6 +581,13 @@ function solve_system_matrix!(dv, system_matrix, rhs, return nothing end +# Svärd-Kalisch equations for reflecting boundary conditions uses LU factorization +function solve_system_matrix_lu!(dv, system_matrix_lu, ::SvaerdKalischEquations1D, cache) + (; factorization_lu) = cache + lu!(factorization_lu, system_matrix_lu) + ldiv!(factorization_lu, dv) +end + function solve_system_matrix!(dv, system_matrix, ::Union{BBMEquation1D, BBMBBMEquations1D}) ldiv!(system_matrix, dv) end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index b2a65da0..a8e2ecfd 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -223,6 +223,13 @@ function assemble_system_matrix!(cache, h, return Symmetric(Diagonal((@view M_h[2:(end - 1)])) + minus_MD1betaD1) end +function assemble_system_matrix_lu!(cache, h, + ::SvaerdKalischEquations1D, + ::BoundaryConditionReflecting) + (; D1betaD1d) = cache + return Diagonal((@view h[2:(end - 1)])) - D1betaD1d +end + function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, boundary_conditions::BoundaryConditionPeriodic, @@ -326,10 +333,15 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, D1 isa UniformCoupledOperator D1_central = D1 D1mat = sparse(D1_central) + D1betaD1d = sparse(D1mat * Diagonal(beta_hat) * D1mat * Pd)[2:(end - 1), :] minus_MD1betaD1 = sparse(D1mat' * Diagonal(M_beta) * D1mat * Pd)[2:(end - 1), :] system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, equations, boundary_conditions) end + system_matrix_lu = let cache = (; D1betaD1d) + assemble_system_matrix_lu!(cache, h, + equations, boundary_conditions) + end elseif D1 isa UpwindOperators D1_central = D1.central D1mat_minus = sparse(D1.minus) @@ -342,7 +354,9 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) end factorization = cholesky(system_matrix) - cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x, + factorization_lu = lu(system_matrix_lu) + cache = (; factorization, factorization_lu, minus_MD1betaD1, D1betaD1d, D, h, hv, b, + eta_x, v_x, h_v_x, hv2_x, beta_hat, D1_central, D1, M_h) return cache end @@ -496,6 +510,50 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, return nothing end +function rhs_lu!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, + initial_condition, boundary_conditions::BoundaryConditionReflecting, + source_terms, solver, cache) + # When constructing the `cache`, we assert that alpha and gamma are zero. + # We use this explicitly in the code below. + (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central, D1) = cache + + g = gravity_constant(equations) + eta, v = q.x + deta, dv, dD = dq.x + fill!(dD, zero(eltype(dD))) + + @trixi_timeit timer() "deta hyperbolic" begin + @.. b = equations.eta0 - D + @.. h = eta - b + @.. hv = h * v + + mul!(deta, D1_central, -hv) + end + + # split form + @trixi_timeit timer() "dv hyperbolic" begin + mul!(eta_x, D1_central, eta) + mul!(v_x, D1_central, v) + mul!(h_v_x, D1_central, hv) + @.. tmp1 = hv * v + mul!(hv2_x, D1_central, tmp1) + @.. dv = -(0.5 * (hv2_x + hv * v_x - v * h_v_x) + + g * h * eta_x) + end + + @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, + solver) + @trixi_timeit timer() "assemble system matrix" begin + system_matrix_lu = assemble_system_matrix_lu!(cache, h, equations, boundary_conditions) + end + @trixi_timeit timer() "solving elliptic system" begin + solve_system_matrix_lu!((@view dv[2:(end - 1)]), system_matrix_lu, equations, cache) + dv[1] = dv[end] = zero(eltype(dv)) + end + + return nothing +end + # The modified entropy/energy takes the whole `q` for every point in space """ DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D, cache) From f01761985db1fee95551c9a83655e24d79e6c405 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 15:31:49 +0100 Subject: [PATCH 13/28] fix scale_by_mass_matrix! --- src/equations/equations.jl | 2 +- src/equations/svaerd_kalisch_1d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 85e098ef..1fe16d5f 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -567,7 +567,7 @@ function solve_system_matrix!(dv, system_matrix, rhs, if issuccess(factorization) scale_by_mass_matrix!(rhs, D1) # see https://github.com/JoshuaLampert/DispersiveShallowWater.jl/issues/122 - dv .= factorization \ rhs + dv .= system_matrix \ rhs[2:(end - 1)] else # The factorization may fail if the time step is too large # and h becomes negative. diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index a8e2ecfd..054b60fb 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -502,7 +502,7 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, - (@view tmp1[2:(end - 1)]), equations, D1, + tmp1, equations, D1, cache) dv[1] = dv[end] = zero(eltype(dv)) end From ca345a5f9eadab8306e6673df4ca7e54d393c3e4 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 15:37:23 +0100 Subject: [PATCH 14/28] remove LU decomposition again --- src/DispersiveShallowWater.jl | 2 +- src/equations/equations.jl | 9 +---- src/equations/svaerd_kalisch_1d.jl | 60 +----------------------------- 3 files changed, 3 insertions(+), 68 deletions(-) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index a8c1f41d..83801d31 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -20,7 +20,7 @@ using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, - lu, lu!, cholesky, cholesky!, issuccess + lu, cholesky, cholesky!, issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @recipe, @series diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 1fe16d5f..6a527573 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -567,7 +567,7 @@ function solve_system_matrix!(dv, system_matrix, rhs, if issuccess(factorization) scale_by_mass_matrix!(rhs, D1) # see https://github.com/JoshuaLampert/DispersiveShallowWater.jl/issues/122 - dv .= system_matrix \ rhs[2:(end - 1)] + dv .= factorization \ rhs[2:(end - 1)] else # The factorization may fail if the time step is too large # and h becomes negative. @@ -581,13 +581,6 @@ function solve_system_matrix!(dv, system_matrix, rhs, return nothing end -# Svärd-Kalisch equations for reflecting boundary conditions uses LU factorization -function solve_system_matrix_lu!(dv, system_matrix_lu, ::SvaerdKalischEquations1D, cache) - (; factorization_lu) = cache - lu!(factorization_lu, system_matrix_lu) - ldiv!(factorization_lu, dv) -end - function solve_system_matrix!(dv, system_matrix, ::Union{BBMEquation1D, BBMBBMEquations1D}) ldiv!(system_matrix, dv) end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 054b60fb..791e0b53 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -223,13 +223,6 @@ function assemble_system_matrix!(cache, h, return Symmetric(Diagonal((@view M_h[2:(end - 1)])) + minus_MD1betaD1) end -function assemble_system_matrix_lu!(cache, h, - ::SvaerdKalischEquations1D, - ::BoundaryConditionReflecting) - (; D1betaD1d) = cache - return Diagonal((@view h[2:(end - 1)])) - D1betaD1d -end - function create_cache(mesh, equations::SvaerdKalischEquations1D, solver, initial_condition, boundary_conditions::BoundaryConditionPeriodic, @@ -333,15 +326,10 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, D1 isa UniformCoupledOperator D1_central = D1 D1mat = sparse(D1_central) - D1betaD1d = sparse(D1mat * Diagonal(beta_hat) * D1mat * Pd)[2:(end - 1), :] minus_MD1betaD1 = sparse(D1mat' * Diagonal(M_beta) * D1mat * Pd)[2:(end - 1), :] system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, equations, boundary_conditions) end - system_matrix_lu = let cache = (; D1betaD1d) - assemble_system_matrix_lu!(cache, h, - equations, boundary_conditions) - end elseif D1 isa UpwindOperators D1_central = D1.central D1mat_minus = sparse(D1.minus) @@ -354,9 +342,7 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, throw(ArgumentError("unknown type of first-derivative operator: $(typeof(D1))")) end factorization = cholesky(system_matrix) - factorization_lu = lu(system_matrix_lu) - cache = (; factorization, factorization_lu, minus_MD1betaD1, D1betaD1d, D, h, hv, b, - eta_x, v_x, + cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, beta_hat, D1_central, D1, M_h) return cache end @@ -510,50 +496,6 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, return nothing end -function rhs_lu!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, - initial_condition, boundary_conditions::BoundaryConditionReflecting, - source_terms, solver, cache) - # When constructing the `cache`, we assert that alpha and gamma are zero. - # We use this explicitly in the code below. - (; D, h, hv, b, eta_x, v_x, h_v_x, hv2_x, tmp1, D1_central, D1) = cache - - g = gravity_constant(equations) - eta, v = q.x - deta, dv, dD = dq.x - fill!(dD, zero(eltype(dD))) - - @trixi_timeit timer() "deta hyperbolic" begin - @.. b = equations.eta0 - D - @.. h = eta - b - @.. hv = h * v - - mul!(deta, D1_central, -hv) - end - - # split form - @trixi_timeit timer() "dv hyperbolic" begin - mul!(eta_x, D1_central, eta) - mul!(v_x, D1_central, v) - mul!(h_v_x, D1_central, hv) - @.. tmp1 = hv * v - mul!(hv2_x, D1_central, tmp1) - @.. dv = -(0.5 * (hv2_x + hv * v_x - v * h_v_x) + - g * h * eta_x) - end - - @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, - solver) - @trixi_timeit timer() "assemble system matrix" begin - system_matrix_lu = assemble_system_matrix_lu!(cache, h, equations, boundary_conditions) - end - @trixi_timeit timer() "solving elliptic system" begin - solve_system_matrix_lu!((@view dv[2:(end - 1)]), system_matrix_lu, equations, cache) - dv[1] = dv[end] = zero(eltype(dv)) - end - - return nothing -end - # The modified entropy/energy takes the whole `q` for every point in space """ DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D, cache) From 3fa969f64749247e54be6d9d313953250fb32afe Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 15:48:34 +0100 Subject: [PATCH 15/28] scale_by_mass_matrix! outside of solve_system_matrix! --- src/equations/equations.jl | 4 +--- src/equations/serre_green_naghdi_1d.jl | 4 ++++ src/equations/svaerd_kalisch_1d.jl | 4 +++- test/test_svaerd_kalisch_1d.jl | 20 ++++++++++---------- 4 files changed, 18 insertions(+), 14 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 6a527573..1c19c37f 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -565,9 +565,8 @@ function solve_system_matrix!(dv, system_matrix, rhs, (; factorization) = cache cholesky!(factorization, system_matrix; check = false) if issuccess(factorization) - scale_by_mass_matrix!(rhs, D1) # see https://github.com/JoshuaLampert/DispersiveShallowWater.jl/issues/122 - dv .= factorization \ rhs[2:(end - 1)] + dv .= factorization \ rhs else # The factorization may fail if the time step is too large # and h becomes negative. @@ -575,7 +574,6 @@ function solve_system_matrix!(dv, system_matrix, rhs, end else factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(rhs, D1) ldiv!(dv, factorization, rhs) end return nothing diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 7c4d1984..6beaf4d6 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -399,6 +399,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla end @trixi_timeit timer() "solving elliptic system" begin + scale_by_mass_matrix!(tmp, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end @@ -484,6 +485,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat end @trixi_timeit timer() "solving elliptic system" begin + scale_by_mass_matrix!(tmp1, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end @@ -589,6 +591,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin + scale_by_mass_matrix!(tmp1, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end @@ -702,6 +705,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin + scale_by_mass_matrix!(tmp1, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 791e0b53..01f1bb09 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -442,6 +442,7 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv + scale_by_mass_matrix!(tmp1, D1) solve_system_matrix!(dv, system_matrix, tmp1, equations, D1, cache) end @@ -487,8 +488,9 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv + scale_by_mass_matrix!(tmp1, D1) solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, - tmp1, equations, D1, + (@view tmp1[2:(end - 1)]), equations, D1, cache) dv[1] = dv[end] = zero(eltype(dv)) end diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 9a919506..9dddbf7d 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -24,11 +24,11 @@ end ] begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), - l2=[5.402315467757905e-6 6.705583190754449e-8 0.0], - linf=[9.787584609100008e-5 1.4602668577112787e-7 0.0], - cons_error=[1.0484865187415942e-9 0.5469460930247753 0.0], - change_waterheight=1.0484865187415942e-9, - change_entropy_modified=459.9037236233692, + l2=[5.402315494643131e-6 6.70558305594986e-8 0.0], + linf=[9.787585531206844e-5 1.460266869784954e-7 0.0], + cons_error=[1.0484871583691058e-9 0.5469460930247998 0.0], + change_waterheight=1.0484871583691058e-9, + change_entropy_modified=459.90372362340514, atol=1e-11) # to make CI pass @test_allocations(semi, sol, allocs=650_000) @@ -42,11 +42,11 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), solver=solver, - l2=[5.862278093002509e-6 4.111933117101271e-9 0.0], - linf=[3.135229084794133e-5 8.797788287467911e-8 0.0], - cons_error=[1.700425282207037e-9 0.5469460935004404 0.0], - change_waterheight=-1.700425282207037e-9, - change_entropy_modified=459.9037221442477, + l2=[5.862278175937948e-6 4.11195454078554e-9 0.0], + linf=[3.135228725170691e-5 8.797787950237668e-8 0.0], + cons_error=[1.700425028441774e-9 0.5469460935005555 0.0], + change_waterheight=-1.700425028441774e-9, + change_entropy_modified=459.9037221442321, atol=1e-11) # to make CI pass @test_allocations(semi, sol, allocs=650_000) From 82becee1561a2cecd2418e0d2cd59c3edd4cf128 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 15:49:16 +0100 Subject: [PATCH 16/28] try CI without tolerances --- test/test_svaerd_kalisch_1d.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 9dddbf7d..f1823ba1 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -28,8 +28,7 @@ end linf=[9.787585531206844e-5 1.460266869784954e-7 0.0], cons_error=[1.0484871583691058e-9 0.5469460930247998 0.0], change_waterheight=1.0484871583691058e-9, - change_entropy_modified=459.90372362340514, - atol=1e-11) # to make CI pass + change_entropy_modified=459.90372362340514) @test_allocations(semi, sol, allocs=650_000) @@ -46,8 +45,7 @@ end linf=[3.135228725170691e-5 8.797787950237668e-8 0.0], cons_error=[1.700425028441774e-9 0.5469460935005555 0.0], change_waterheight=-1.700425028441774e-9, - change_entropy_modified=459.9037221442321, - atol=1e-11) # to make CI pass + change_entropy_modified=459.9037221442321) @test_allocations(semi, sol, allocs=650_000) end From 60e800ea3638b4d918e9dedd705827e875d8fe96 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 15:50:48 +0100 Subject: [PATCH 17/28] reorder imports --- src/DispersiveShallowWater.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 83801d31..54633bf0 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -19,8 +19,8 @@ using BandedMatrices: BandedMatrix using DiffEqBase: DiffEqBase, SciMLBase, terminate! using FastBroadcast: @.. using Interpolations: Interpolations, linear_interpolation -using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, - lu, cholesky, cholesky!, issuccess +using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, cholesky, cholesky!, + issuccess using PolynomialBases: PolynomialBases using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @recipe, @series From 08c04e694f1e6dd1c425174be269a3547b3c189b Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 15:54:11 +0100 Subject: [PATCH 18/28] fix variable name --- src/equations/serre_green_naghdi_1d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 6beaf4d6..5af96671 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -485,7 +485,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp1, D1) + scale_by_mass_matrix!(tmp, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end @@ -591,7 +591,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp1, D1) + scale_by_mass_matrix!(tmp, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end @@ -705,7 +705,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp1, D1) + scale_by_mass_matrix!(tmp, D1) solve_system_matrix!(dv, system_matrix, tmp, equations, D1, cache) end From 94991f34a3f6f35b23d23734aab88646edbaa75c Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 16:02:19 +0100 Subject: [PATCH 19/28] refactor to reduce code redundancy --- src/equations/equations.jl | 16 ++++++++++++++++ src/equations/serre_green_naghdi_1d.jl | 20 ++++++++------------ src/equations/svaerd_kalisch_1d.jl | 9 +++------ 3 files changed, 27 insertions(+), 18 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 1c19c37f..ef81e73d 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -557,6 +557,22 @@ abstract type AbstractSerreGreenNaghdiEquations{NDIMS, NVARS} <: include("serre_green_naghdi_1d.jl") include("hyperbolic_serre_green_naghdi_1d.jl") +function solve_system_matrix!(dv, system_matrix, rhs, + ::Union{SvaerdKalischEquations1D, + SerreGreenNaghdiEquations1D}, + D1, cache, ::BoundaryConditionPeriodic) + scale_by_mass_matrix!(rhs, D1) + solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache) +end + +function solve_system_matrix!(dv, system_matrix, rhs, + ::Union{SvaerdKalischEquations1D, + SerreGreenNaghdiEquations1D}, + D1, cache, ::BoundaryConditionReflecting) + scale_by_mass_matrix!(rhs, D1) + solve_system_matrix!(dv, system_matrix, (@view rhs[2:(end - 1)]), equations, D1, cache) +end + function solve_system_matrix!(dv, system_matrix, rhs, ::Union{SvaerdKalischEquations1D, SerreGreenNaghdiEquations1D}, diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 5af96671..4dcc91a8 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -399,9 +399,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp, D1) - solve_system_matrix!(dv, system_matrix, tmp, - equations, D1, cache) + solve_system_matrix!(dv, system_matrix, + tmp, equations, D1, cache, boundary_conditions) end return nothing @@ -485,9 +484,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp, D1) - solve_system_matrix!(dv, system_matrix, tmp, - equations, D1, cache) + solve_system_matrix!(dv, system_matrix, + tmp, equations, D1, cache, boundary_conditions) end return nothing @@ -591,9 +589,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp, D1) - solve_system_matrix!(dv, system_matrix, tmp, - equations, D1, cache) + solve_system_matrix!(dv, system_matrix, + tmp, equations, D1, cache, boundary_conditions) end return nothing @@ -705,9 +702,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - scale_by_mass_matrix!(tmp, D1) - solve_system_matrix!(dv, system_matrix, tmp, - equations, D1, cache) + solve_system_matrix!(dv, system_matrix, + tmp, equations, D1, cache, boundary_conditions) end return nothing diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 01f1bb09..057b49b1 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -442,9 +442,8 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv - scale_by_mass_matrix!(tmp1, D1) - solve_system_matrix!(dv, system_matrix, tmp1, - equations, D1, cache) + solve_system_matrix!(dv, system_matrix, + tmp1, equations, D1, cache, boundary_conditions) end return nothing @@ -488,10 +487,8 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv - scale_by_mass_matrix!(tmp1, D1) solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, - (@view tmp1[2:(end - 1)]), equations, D1, - cache) + tmp1, equations, D1, cache, boundary_conditions) dv[1] = dv[end] = zero(eltype(dv)) end From 320d8830282950f6f0dbf41c05d1c07dd8335ec2 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 16:03:34 +0100 Subject: [PATCH 20/28] fix --- src/equations/equations.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ef81e73d..4c8f7ec0 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -558,16 +558,16 @@ include("serre_green_naghdi_1d.jl") include("hyperbolic_serre_green_naghdi_1d.jl") function solve_system_matrix!(dv, system_matrix, rhs, - ::Union{SvaerdKalischEquations1D, - SerreGreenNaghdiEquations1D}, + equations::Union{SvaerdKalischEquations1D, + SerreGreenNaghdiEquations1D}, D1, cache, ::BoundaryConditionPeriodic) scale_by_mass_matrix!(rhs, D1) solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache) end function solve_system_matrix!(dv, system_matrix, rhs, - ::Union{SvaerdKalischEquations1D, - SerreGreenNaghdiEquations1D}, + equations::Union{SvaerdKalischEquations1D, + SerreGreenNaghdiEquations1D}, D1, cache, ::BoundaryConditionReflecting) scale_by_mass_matrix!(rhs, D1) solve_system_matrix!(dv, system_matrix, (@view rhs[2:(end - 1)]), equations, D1, cache) From a2d84975399d7165a9236420bdcf81d33a036314 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 16:09:18 +0100 Subject: [PATCH 21/28] pass boundary_conditions --- src/equations/serre_green_naghdi_1d.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 4dcc91a8..4f35a0bd 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -314,19 +314,20 @@ end function rhs!(dq, q, t, mesh, equations::SerreGreenNaghdiEquations1D, initial_condition, - ::BoundaryConditionPeriodic, + boundary_conditions::BoundaryConditionPeriodic, source_terms::Nothing, solver, cache) if cache.D1 isa PeriodicUpwindOperators - rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type) + rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type, boundary_conditions) else - rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type) + rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type, boundary_conditions) end return nothing end -function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat) +function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat, + boundary_conditions::BoundaryConditionPeriodic) # Unpack physical parameters and SBP operator `D1` as well as the # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) @@ -406,7 +407,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla return nothing end -function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat) +function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat, + boundary_conditions::BoundaryConditionPeriodic) # Unpack physical parameters and SBP operator `D1` as well as the # SBP upwind operator in sparse matrix form `D1mat_minus` g = gravity_constant(equations) @@ -492,7 +494,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat end function rhs_sgn_central!(dq, q, equations, source_terms, cache, - ::Union{BathymetryMildSlope, BathymetryVariable}) + ::Union{BathymetryMildSlope, BathymetryVariable}, + boundary_conditions::BoundaryConditionPeriodic) # Unpack physical parameters and SBP operator `D1` as well as the # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) @@ -597,7 +600,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, - ::Union{BathymetryMildSlope, BathymetryVariable}) + ::Union{BathymetryMildSlope, BathymetryVariable}, + boundary_conditions::BoundaryConditionPeriodic) # Unpack physical parameters and SBP operator `D1` as well as the # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) From 9bfebe360a99fc5cc8e9ac101ab4f8e57b985928 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 16:10:53 +0100 Subject: [PATCH 22/28] format --- src/equations/serre_green_naghdi_1d.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 4f35a0bd..d9e162a3 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -318,9 +318,11 @@ function rhs!(dq, q, t, mesh, source_terms::Nothing, solver, cache) if cache.D1 isa PeriodicUpwindOperators - rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type, boundary_conditions) + rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type, + boundary_conditions) else - rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type, boundary_conditions) + rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type, + boundary_conditions) end return nothing From 2d5083310507058115a7065846ebde89d0626b3b Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 16:32:47 +0100 Subject: [PATCH 23/28] adjust tolerances --- test/test_svaerd_kalisch_1d.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index f1823ba1..9085eebc 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -28,7 +28,8 @@ end linf=[9.787585531206844e-5 1.460266869784954e-7 0.0], cons_error=[1.0484871583691058e-9 0.5469460930247998 0.0], change_waterheight=1.0484871583691058e-9, - change_entropy_modified=459.90372362340514) + change_entropy_modified=459.90372362340514, + atol=1e-11) # to make CI pass @test_allocations(semi, sol, allocs=650_000) @@ -45,7 +46,8 @@ end linf=[3.135228725170691e-5 8.797787950237668e-8 0.0], cons_error=[1.700425028441774e-9 0.5469460935005555 0.0], change_waterheight=-1.700425028441774e-9, - change_entropy_modified=459.9037221442321) + change_entropy_modified=459.9037221442321, + atol=1e-10) # to make CI pass @test_allocations(semi, sol, allocs=650_000) end From 209cc76080ff015c2e8cd5618eeae652346de1ed Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 11 Dec 2024 17:09:58 +0100 Subject: [PATCH 24/28] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/equations/equations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 4c8f7ec0..d003e40a 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -570,7 +570,7 @@ function solve_system_matrix!(dv, system_matrix, rhs, SerreGreenNaghdiEquations1D}, D1, cache, ::BoundaryConditionReflecting) scale_by_mass_matrix!(rhs, D1) - solve_system_matrix!(dv, system_matrix, (@view rhs[2:(end - 1)]), equations, D1, cache) + solve_system_matrix!(dv, system_matrix, (@view rhs[(begin + 1):(end - 1)]), equations, D1, cache) end function solve_system_matrix!(dv, system_matrix, rhs, From 230812189db4659558cfd0a7831f3a78039900fa Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 11 Dec 2024 17:11:34 +0100 Subject: [PATCH 25/28] Update src/equations/equations.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/equations.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index d003e40a..7f9d2e54 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -570,7 +570,8 @@ function solve_system_matrix!(dv, system_matrix, rhs, SerreGreenNaghdiEquations1D}, D1, cache, ::BoundaryConditionReflecting) scale_by_mass_matrix!(rhs, D1) - solve_system_matrix!(dv, system_matrix, (@view rhs[(begin + 1):(end - 1)]), equations, D1, cache) + solve_system_matrix!(dv, system_matrix, (@view rhs[(begin + 1):(end - 1)]), equations, + D1, cache) end function solve_system_matrix!(dv, system_matrix, rhs, From 3302cd0c21320d9a42bb46ea583c3ca1f87e20a3 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 17:15:27 +0100 Subject: [PATCH 26/28] begin + 1instead of 2 --- src/equations/svaerd_kalisch_1d.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 057b49b1..71e633c6 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -220,7 +220,7 @@ function assemble_system_matrix!(cache, h, # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. - return Symmetric(Diagonal((@view M_h[2:(end - 1)])) + minus_MD1betaD1) + return Symmetric(Diagonal((@view M_h[(begin + 1):(end - 1)])) + minus_MD1betaD1) end function create_cache(mesh, equations::SvaerdKalischEquations1D, @@ -326,15 +326,16 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D, D1 isa UniformCoupledOperator D1_central = D1 D1mat = sparse(D1_central) - minus_MD1betaD1 = sparse(D1mat' * Diagonal(M_beta) * D1mat * Pd)[2:(end - 1), :] + minus_MD1betaD1 = sparse(D1mat' * Diagonal(M_beta) * + D1mat * Pd)[(begin + 1):(end - 1), :] system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, equations, boundary_conditions) end elseif D1 isa UpwindOperators D1_central = D1.central D1mat_minus = sparse(D1.minus) - minus_MD1betaD1 = sparse(D1mat_minus' * Diagonal(M_beta) * D1mat_minus * Pd)[2:(end - 1), - :] + minus_MD1betaD1 = sparse(D1mat_minus' * Diagonal(M_beta) * + D1mat_minus * Pd)[(begin + 1):(end - 1), :] system_matrix = let cache = (; D1, M_h, minus_MD1betaD1) assemble_system_matrix!(cache, h, equations, boundary_conditions) end @@ -487,7 +488,7 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D, end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv - solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix, + solve_system_matrix!((@view dv[(begin + 1):(end - 1)]), system_matrix, tmp1, equations, D1, cache, boundary_conditions) dv[1] = dv[end] = zero(eltype(dv)) end From 7358093cff699397ac2583ffaad4cde2a9021f0f Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 17:17:17 +0100 Subject: [PATCH 27/28] set tolerances in tests --- .../svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl | 2 +- test/test_svaerd_kalisch_1d.jl | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl index a6edb3dd..0590c419 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl @@ -42,5 +42,5 @@ analysis_callback = AnalysisCallback(semi; interval = 100, callbacks = CallbackSet(analysis_callback, summary_callback) saveat = range(tspan..., length = 100) -sol = solve(ode, Tsit5(), abstol = 1e-12, reltol = 1e-12, +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 9085eebc..ce64b389 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -24,6 +24,8 @@ end ] begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), + abstol=1e-12, + reltol=1e-12, # this example is relatively unstable with higher tolerances l2=[5.402315494643131e-6 6.70558305594986e-8 0.0], linf=[9.787585531206844e-5 1.460266869784954e-7 0.0], cons_error=[1.0484871583691058e-9 0.5469460930247998 0.0], @@ -42,6 +44,8 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"), tspan=(0.0, 1.0), solver=solver, + abstol=1e-12, + reltol=1e-12, # this example is relatively unstable with higher tolerances l2=[5.862278175937948e-6 4.11195454078554e-9 0.0], linf=[3.135228725170691e-5 8.797787950237668e-8 0.0], cons_error=[1.700425028441774e-9 0.5469460935005555 0.0], From 7092265877d86b41c3cffbb346a03b269510c7e9 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Wed, 11 Dec 2024 17:20:47 +0100 Subject: [PATCH 28/28] add test for ArgumentError --- test/test_unit.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/test_unit.jl b/test/test_unit.jl index 6952ea6b..8e26df8e 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -250,6 +250,15 @@ end U_modified_total = @inferred DispersiveShallowWater.integrate(U_modified, semi) @test isapprox(U_modified_total, e_modified_total) end + + @testset "reflecting boundary conditions" begin + initial_condition = initial_condition_manufactured_reflecting + boundary_conditions = boundary_condition_reflecting + mesh = Mesh1D(-1.0, 1.0, 10) + solver = Solver(mesh, 4) + @test_throws ArgumentError Semidiscretization(mesh, equations, initial_condition, + solver; boundary_conditions) + end end @testitem "SerreGreenNaghdiEquations1D" setup=[Setup] begin