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..0590c419 --- /dev/null +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl @@ -0,0 +1,46 @@ +using OrdinaryDiffEqTsit5 +using DispersiveShallowWater +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) + +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) +solver = Solver(D1, nothing) + +# 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, + entropy_modified)) +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/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 8f2f8610..99cfd788 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 @@ -501,8 +499,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..7f9d2e54 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -557,6 +557,23 @@ 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, + 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, + equations::Union{SvaerdKalischEquations1D, + 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) +end + function solve_system_matrix!(dv, system_matrix, rhs, ::Union{SvaerdKalischEquations1D, SerreGreenNaghdiEquations1D}, @@ -565,7 +582,6 @@ 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 else @@ -575,7 +591,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..d9e162a3 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -314,19 +314,22 @@ 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) @@ -399,14 +402,15 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla end @trixi_timeit timer() "solving elliptic system" begin - 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 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) @@ -484,15 +488,16 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat end @trixi_timeit timer() "solving elliptic system" begin - 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 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) @@ -589,15 +594,16 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - 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 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) @@ -702,8 +708,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - 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 bb594838..71e633c6 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -140,9 +140,92 @@ 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, + ::SvaerdKalischEquations1D, + ::BoundaryConditionPeriodic) + (; 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(M_h) + minus_MD1betaD1) +end + +# For reflecting boundary conditions +function assemble_system_matrix!(cache, h, + ::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[(begin + 1):(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. @@ -160,7 +243,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) @@ -173,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 || @@ -183,26 +265,26 @@ 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, D1mat, 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, D1mat_minus, equations) + equations, boundary_conditions) end else 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, - 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) + 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) @@ -211,18 +293,59 @@ 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 +# Reflecting boundary conditions assume alpha = gamma = 0 +function create_cache(mesh, equations::SvaerdKalischEquations1D, + solver, initial_condition, + 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")) + 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 + M_h = copy(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) + 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) + 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)[(begin + 1):(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 = 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 # Discretization that conserves @@ -235,11 +358,11 @@ 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) - (; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x, + 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, M, D1) = cache + tmp1, tmp2, D1_central, D1) = cache g = gravity_constant(equations) eta, v = q.x @@ -315,12 +438,59 @@ 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, equations, + boundary_conditions) + end + @trixi_timeit timer() "solving elliptic system" begin + tmp1 .= dv + solve_system_matrix!(dv, system_matrix, + tmp1, equations, D1, cache, boundary_conditions) + end + + return nothing +end + +function rhs!(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 = assemble_system_matrix!(cache, h, equations, boundary_conditions) end @trixi_timeit timer() "solving elliptic system" begin tmp1 .= dv - solve_system_matrix!(dv, system_matrix, tmp1, - equations, D1, cache) + 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 return nothing @@ -358,7 +528,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 5b919b9a..ce64b389 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -17,6 +17,45 @@ 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), + 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], + change_waterheight=1.0484871583691058e-9, + change_entropy_modified=459.90372362340514, + atol=1e-11) # to make CI pass + + @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) + solver = Solver(D1, nothing) + @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], + change_waterheight=-1.700425028441774e-9, + change_entropy_modified=459.9037221442321, + atol=1e-10) # to make CI pass + + @test_allocations(semi, sol, allocs=650_000) +end + @testitem "svaerd_kalisch_1d_dingemans" setup=[ Setup, SvaerdKalischEquations1D, 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