diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic_reflecting.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic_reflecting.jl new file mode 100644 index 00000000..456d4357 --- /dev/null +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic_reflecting.jl @@ -0,0 +1,48 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: MattssonNordström2004, derivative_operator + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) + +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 = 10, + 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/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index ba810735..831219e1 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -101,6 +101,47 @@ function source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D return SVector(dq1, dq2, zero(dq1)) end +""" + initial_condition_manufactured_reflecting(x, t, equations::BBMBBMVariableEquations1D, 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::BBMBBMVariableEquations1D, + mesh) + eta = exp(2 * t) * cospi(x) + v = exp(t) * x * sinpi(x) + D = 5 + 2 * cospi(2 * x) + return SVector(eta, v, D) +end + +""" + source_terms_manufactured_reflecting(q, x, t, equations::BBMBBMVariableEquations1D, 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::BBMBBMVariableEquations1D) + g = equations.gravity + a1 = cospi(2 * x) + a2 = sinpi(2 * x) + a8 = cospi(x) + a9 = sinpi(x) + a10 = exp(t) + a11 = exp(2 * t) + a12 = cospi(3 * x) + a13 = sinpi(3 * x) + dq1 = (pi * x * a11 * a1 + 4 * pi * x * a8 + 3 * pi * x * a12 + + 20 * pi^2 * (1 - a1)^2 * a10 * a8 / 3 + a11 * a2 / 2 + 2 * a10 * a8 + + 7 * pi^2 * a10 * a8 / 3 + 14 * pi^2 * a10 * a12 + 4 * a9 + a13) * a10 + dq2 = (-pi * g * a10 * a9 + pi * x^2 * a10 * a2 / 2 + x * a10 * a9^2 + x * a9 + + pi * (400 * pi * x * a9^5 - 824 * pi * x * a9^3 + 385 * pi * x * a9 - + 160 * a9^4 * a8 + 336 * a9^2 * a8 - 98 * a8) / 6) * a10 + + return SVector(dq1, dq2, zero(dq1)) +end + """ initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh) @@ -171,6 +212,43 @@ function create_cache(mesh, equations::BBMBBMVariableEquations1D, return (invImDKD = invImDKD, invImD2K = invImD2K, D = D) end +function create_cache(mesh, equations::BBMBBMVariableEquations1D, + solver, initial_condition, + ::BoundaryConditionReflecting, + RealT, uEltype) + # Assume D is independent of time and compute D evaluated at mesh points once. + D = Array{RealT}(undef, nnodes(mesh)) + x = grid(solver) + for i in eachnode(solver) + D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations) + end + K = Diagonal(D .^ 2) + K_i = Diagonal(D[2:(end - 1)] .^ 2) + N = nnodes(mesh) + M = mass_matrix(solver.D1) + Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2)) + D2d = (sparse(solver.D2) * Pd)[2:(end - 1), :] + # homogeneous Dirichlet boundary conditions + invImD2Kd = inv(I - 1 / 6 * D2d * K_i) + m = diag(M) + m[1] = 0 + m[end] = 0 + PdM = Diagonal(m) + + # homogeneous Neumann boundary conditions + if solver.D1 isa DerivativeOperator || + solver.D1 isa UniformCoupledOperator + D1_b = BandedMatrix(solver.D1) + invImDKDn = inv(I + 1 / 6 * inv(M) * D1_b' * PdM * K * D1_b) + elseif solver.D1 isa UpwindOperators + D1plus_b = BandedMatrix(solver.D1.plus) + invImDKDn = inv(I + 1 / 6 * inv(M) * D1plus_b' * PdM * K * D1plus_b) + else + @error "unknown type of first-derivative operator: $(typeof(solver.D1))" + end + return (invImD2Kd = invImD2Kd, invImDKDn = invImDKDn, D = D) +end + # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see # - Joshua Lampert and Hendrik Ranocha (2024) # Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations @@ -211,6 +289,46 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D, return nothing end +function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D, + initial_condition, ::BoundaryConditionReflecting, source_terms, + solver, cache) + @unpack invImDKDn, invImD2Kd, D = cache + + q = wrap_array(u_ode, mesh, equations, solver) + dq = wrap_array(du_ode, mesh, equations, solver) + + eta = view(q, 1, :) + v = view(q, 2, :) + deta = view(dq, 1, :) + dv = view(dq, 2, :) + dD = view(dq, 3, :) + fill!(dD, zero(eltype(dD))) + + # energy and mass conservative semidiscretization + if solver.D1 isa DerivativeOperator || + solver.D1 isa UniformCoupledOperator + @timeit timer() "deta hyperbolic" deta[:]=-solver.D1 * (D .* v + eta .* v) + @timeit timer() "dv hyperbolic" dv[:]=-solver.D1 * + (equations.gravity * eta + 0.5 * v .^ 2) + elseif solver.D1 isa UpwindOperators + @timeit timer() "deta hyperbolic" deta[:]=-solver.D1.minus * (D .* v + eta .* v) + @timeit timer() "dv hyperbolic" dv[:]=-solver.D1.plus * + (equations.gravity * eta + 0.5 * v .^ 2) + else + @error "unknown type of first-derivative operator: $(typeof(solver.D1))" + end + + @timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver) + + @timeit timer() "deta elliptic" deta[:]=invImDKDn * deta + @timeit timer() "dv elliptic" begin + dv[2:(end - 1)] = invImD2Kd * dv[2:(end - 1)] + dv[1] = dv[end] = zero(eltype(dv)) + end + + return nothing +end + @inline function prim2cons(q, equations::BBMBBMVariableEquations1D) eta, v, D = q diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl index 03f2abf2..e5c20fb2 100644 --- a/test/test_bbm_bbm_variable_bathymetry_1d.jl +++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl @@ -97,6 +97,50 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") change_velocity=-3.1478183774857893e-15, change_entropy=3.417533493699221e-7) end + + @trixi_testset "bbm_bbm_variable_bathymetry_1d_basic_reflecting" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "bbm_bbm_variable_bathymetry_1d_basic_reflecting.jl"), + tspan=(0.0, 1.0), + l2=[0.00022223640634220205 7.62845117195942e-9 0.0], + linf=[0.00029474379052363275 1.7386610817737846e-8 0.0], + cons_error=[9.711471991876855e-10 0.5469460962477994 0.0], + change_waterheight=9.711471991876855e-10, + change_velocity=0.5469460962477994, + change_entropy=132.10935771964688, + atol=1e-10, + atol_ints=1e-10) # in order to make CI pass + + # test upwind operators + using SummationByPartsOperators: upwind_operators, Mattsson2017 + using SparseArrays: sparse + using OrdinaryDiffEq: solve + 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) + semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions, + source_terms = source_terms) + ode = semidiscretize(semi, (0.0, 1.0)) + sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) + atol = 1e-7 # in order to make CI pass + rtol = 1e-12 + errs = errors(analysis_callback) + l2 = [0.0023458014985457366 3.261497256675908e-8] + l2_measured = errs.l2_error[:, end] + for (l2_expected, l2_actual) in zip(l2, l2_measured) + @test isapprox(l2_expected, l2_actual, atol = atol, rtol = rtol) + end + linf = [0.05625954556488111 6.815520172120948e-7] + linf_measured = errs.linf_error[:, end] + for (linf_expected, linf_actual) in zip(linf, linf_measured) + @test isapprox(linf_expected, linf_actual, atol = atol, rtol = rtol) + end + end end end # module