From a7dcc23e5bd03ab0fb453f1b257cc6603c865a55 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Mon, 20 Nov 2023 12:45:58 +0100 Subject: [PATCH] Add support for source terms (#65) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * add support for source terms * fix calculation of source terms * add test for manufactured solution of variable BBM-BBM * fix test of convergence order * adjust tolerance * source terms for BBMBBMEquations * format * WIP: add source terms for Svärd-Kalisch equations * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * add references * add support of source terms for Svärd-Kalisch equations with constant bathymetry * lower tolerances * fix tolerances * fix tolerances --------- Co-authored-by: Hendrik Ranocha --- .../bbm_bbm_1d/bbm_bbm_1d_manufactured.jl | 40 +++++++++ ...bbm_variable_bathymetry_1d_manufactured.jl | 40 +++++++++ .../svaerd_kalisch_1d_manufactured.jl | 43 ++++++++++ src/DispersiveShallowWater.jl | 6 +- src/equations/bbm_bbm_1d.jl | 66 +++++++++++--- .../bbm_bbm_variable_bathymetry_1d.jl | 86 +++++++++++++++---- src/equations/svaerd_kalisch_1d.jl | 81 +++++++++++++---- src/semidiscretization.jl | 36 ++++---- src/solver.jl | 17 ++++ test/test_bbm_bbm_1d.jl | 13 +++ test/test_bbm_bbm_variable_bathymetry_1d.jl | 17 +++- test/test_svaerd_kalisch_1d.jl | 12 +++ test/test_unit.jl | 3 +- 13 files changed, 395 insertions(+), 65 deletions(-) create mode 100644 examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl create mode 100644 examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl create mode 100644 examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl diff --git a/examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl new file mode 100644 index 00000000..c051cb3e --- /dev/null +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl @@ -0,0 +1,40 @@ +using OrdinaryDiffEq +using DispersiveShallowWater + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 4.0) + +initial_condition = initial_condition_manufactured +source_terms = source_terms_manufactured +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = 0.0 +coordinates_max = 1.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) + +# 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) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_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/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl new file mode 100644 index 00000000..12204888 --- /dev/null +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl @@ -0,0 +1,40 @@ +using OrdinaryDiffEq +using DispersiveShallowWater + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) + +initial_condition = initial_condition_manufactured +source_terms = source_terms_manufactured +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = 0.0 +coordinates_max = 1.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) + +# 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) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_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/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl new file mode 100644 index 00000000..781d4bf3 --- /dev/null +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl @@ -0,0 +1,43 @@ +using OrdinaryDiffEq +using DispersiveShallowWater + +############################################################################### +# Semidiscretization of the Svärd-Kalisch equations + +equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0, + alpha = 0.0004040404040404049, + beta = 0.49292929292929294, + gamma = 0.15707070707070708) + +initial_condition = initial_condition_manufactured +source_terms = source_terms_manufactured +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = 0.0 +coordinates_max = 1.0 +N = 128 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) + +# 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) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_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 d378b2b0..d866f52f 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -23,8 +23,8 @@ import SummationByPartsOperators: grid, xmin, xmax include("boundary_conditions.jl") include("mesh.jl") -include("solver.jl") include("equations/equations.jl") +include("solver.jl") include("semidiscretization.jl") include("callbacks_step/callbacks_step.jl") include("visualization.jl") @@ -47,7 +47,9 @@ export Semidiscretization, semidiscretize, grid export boundary_condition_periodic -export initial_condition_convergence_test, initial_condition_dingemans +export initial_condition_convergence_test, + initial_condition_manufactured, source_terms_manufactured, + initial_condition_dingemans export AnalysisCallback, RelaxationCallback export tstops, errors, integrals diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 8e8046b5..fa5844aa 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -53,6 +53,41 @@ function initial_condition_convergence_test(x, t, equations::BBMBBMEquations1D, return SVector(eta, v) end +""" + initial_condition_manufactured(x, t, equations::BBMBBMEquations1D, mesh) + +A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref). +""" +function initial_condition_manufactured(x, t, + equations::BBMBBMEquations1D, + mesh) + eta = exp(t) * cospi(2 * (x - 2 * t)) + v = exp(t / 2) * sinpi(2 * (x - t / 2)) + return SVector(eta, v) +end + +""" + source_terms_manufactured(q, x, t, equations::BBMBBMEquations1D, mesh) + +A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref). +""" +function source_terms_manufactured(q, x, t, equations::BBMBBMEquations1D) + g = equations.gravity + D = equations.D + a3 = cospi(t - 2 * x) + a4 = sinpi(t - 2 * x) + a5 = sinpi(2 * t - 4 * x) + a6 = sinpi(4 * t - 2 * x) + a7 = cospi(4 * t - 2 * x) + dq1 = -2 * pi^2 * D^2 * (4 * pi * a6 - a7) * exp(t) / 3 + 2 * pi * D * exp(t / 2) * a3 - + 2 * pi * exp(3 * t / 2) * a4 * a6 + 2 * pi * exp(3 * t / 2) * a3 * a7 - + 4 * pi * exp(t) * a6 + exp(t) * a7 + dq2 = -pi^2 * D^2 * (a4 + 2 * pi * a3) * exp(t / 2) / 3 + 2 * pi * g * exp(t) * a6 - + exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5 + + return SVector(dq1, dq2) +end +# function create_cache(mesh, equations::BBMBBMEquations1D, solver, @@ -61,15 +96,14 @@ function create_cache(mesh, uEltype) if solver.D1 isa PeriodicDerivativeOperator || solver.D1 isa UniformPeriodicCoupledOperator - invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ Matrix(solver.D1) + invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2)) elseif solver.D1 isa PeriodicUpwindOperators - invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ - Matrix(solver.D1.central) + invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2)) else @error "unknown type of first-derivative operator" end - tmp1 = Array{RealT}(undef, nnodes(mesh)) - return (invImD2_D = invImD2_D, tmp1 = tmp1) + tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback` + return (invImD2 = invImD2, tmp1 = tmp1) end # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see @@ -77,8 +111,8 @@ end # A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations # [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119) function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_condition, - ::BoundaryConditionPeriodic, solver, cache) - @unpack invImD2_D, tmp1 = cache + ::BoundaryConditionPeriodic, source_terms, solver, cache) + @unpack invImD2 = cache q = wrap_array(u_ode, mesh, equations, solver) dq = wrap_array(du_ode, mesh, equations, solver) @@ -89,11 +123,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_cond dv = view(dq, 2, :) # energy and mass conservative semidiscretization - @. tmp1 = -(equations.D * v + eta * v) - mul!(deta, invImD2_D, tmp1) + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator + deta[:] = -solver.D1 * (equations.D * v + eta .* v) + dv[:] = -solver.D1 * (equations.gravity * eta + 0.5 * v .^ 2) + elseif solver.D1 isa PeriodicUpwindOperators + deta[:] = -solver.D1.central * (equations.D * v + eta .* v) + dv[:] = -solver.D1.central * (equations.gravity * eta + 0.5 * v .^ 2) + else + @error "unknown type of first-derivative operator" + end + + calc_sources!(dq, q, t, source_terms, equations, solver) - @. tmp1 = -(equations.gravity * eta + 0.5 * v^2) - mul!(dv, invImD2_D, tmp1) + deta[:] = invImD2 * deta + dv[:] = invImD2 * dv return nothing end diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index ee4729e4..ad63e6a0 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -41,8 +41,7 @@ For details see Example 5 in Section 3 from (here adapted for dimensional equati Exact Traveling-Wave Solutions to Bidirectional Wave Equations [DOI: 10.1023/A:1026667903256](https://doi.org/10.1023/A:1026667903256) """ -function initial_condition_convergence_test(x, - t, +function initial_condition_convergence_test(x, t, equations::BBMBBMVariableEquations1D, mesh) g = equations.gravity @@ -58,6 +57,49 @@ function initial_condition_convergence_test(x, return SVector(eta, v, D) end +""" + initial_condition_manufactured(x, t, equations::BBMBBMVariableEquations1D, mesh) + +A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref). +""" +function initial_condition_manufactured(x, t, + equations::BBMBBMVariableEquations1D, + mesh) + eta = exp(t) * cospi(2 * (x - 2 * t)) + v = exp(t / 2) * sinpi(2 * (x - t / 2)) + D = 5.0 + 2.0 * cospi(2 * x) + return SVector(eta, v, D) +end + +""" + source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D, mesh) + +A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref). +""" +function source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D) + g = equations.gravity + a1 = cospi(2 * x) + a2 = sinpi(2 * x) + a3 = cospi(t - 2 * x) + a4 = sinpi(t - 2 * x) + a5 = sinpi(2 * t - 4 * x) + a6 = sinpi(4 * t - 2 * x) + a7 = cospi(4 * t - 2 * x) + dq1 = -2 * pi^2 * (4 * pi * a6 - a7) * (2 * a1 + 5)^2 * exp(t) / 3 + + 8 * pi^2 * (a6 + 4 * pi * a7) * (2 * a1 + 5) * exp(t) * a2 / 3 + + 2 * pi * (2 * a1 + 5) * exp(t / 2) * a3 - 2 * pi * exp(3 * t / 2) * a4 * a6 + + 2 * pi * exp(3 * t / 2) * a3 * a7 + 4 * pi * exp(t / 2) * a2 * a4 - + 4 * pi * exp(t) * a6 + exp(t) * a7 + dq2 = 2 * pi * g * exp(t) * a6 - + pi^2 * + (8 * (2 * pi * a4 - a3) * (2 * a1 + 5) * a2 + + (a4 + 2 * pi * a3) * (2 * a1 + 5)^2 + + 4 * (a4 + 2 * pi * a3) * (16 * sinpi(x)^4 - 26 * sinpi(x)^2 + 7)) * exp(t / 2) / + 3 - exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5 + + return SVector(dq1, dq2, 0.0) +end + """ initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh) @@ -100,7 +142,7 @@ end function create_cache(mesh, equations::BBMBBMVariableEquations1D, - solver::Solver, + solver, initial_condition, RealT, uEltype) @@ -113,18 +155,16 @@ function create_cache(mesh, K = Diagonal(D .^ 2) if solver.D1 isa PeriodicDerivativeOperator || solver.D1 isa UniformPeriodicCoupledOperator - invImDKD_D = (I - 1 / 6 * sparse(solver.D1) * K * sparse(solver.D1)) \ - Matrix(solver.D1) - invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1) + invImDKD = inv(I - 1 / 6 * Matrix(solver.D1) * K * Matrix(solver.D1)) + invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K) elseif solver.D1 isa PeriodicUpwindOperators - invImDKD_D = (I - 1 / 6 * sparse(solver.D1.minus) * K * sparse(solver.D1.plus)) \ - Matrix(solver.D1.minus) - invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1.plus) + invImDKD = inv(I - 1 / 6 * Matrix(solver.D1.minus) * K * Matrix(solver.D1.plus)) + invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K) else @error "unknown type of first-derivative operator" end - tmp1 = Array{RealT}(undef, nnodes(mesh)) - return (invImDKD_D = invImDKD_D, invImD2K_D = invImD2K_D, tmp1 = tmp1) + tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback` + return (invImDKD = invImDKD, invImD2K = invImD2K, tmp1 = tmp1) end # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see @@ -133,9 +173,9 @@ end # [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119) # Here, adapted for spatially varying bathymetry. function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D, - initial_condition, - ::BoundaryConditionPeriodic, solver, cache) - @unpack invImDKD_D, invImD2K_D, tmp1 = cache + initial_condition, ::BoundaryConditionPeriodic, source_terms, + solver, cache) + @unpack invImDKD, invImD2K = cache q = wrap_array(u_ode, mesh, equations, solver) dq = wrap_array(du_ode, mesh, equations, solver) @@ -148,11 +188,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D, dD = view(dq, 3, :) fill!(dD, zero(eltype(dD))) - @. tmp1 = -(D * v + eta * v) - mul!(deta, invImDKD_D, tmp1) + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator + deta[:] = -solver.D1 * (D .* v + eta .* v) + dv[:] = -solver.D1 * (equations.gravity * eta + 0.5 * v .^ 2) + elseif solver.D1 isa PeriodicUpwindOperators + deta[:] = -solver.D1.minus * (D .* v + eta .* v) + dv[:] = -solver.D1.plus * (equations.gravity * eta + 0.5 * v .^ 2) + else + @error "unknown type of first-derivative operator" + end + + calc_sources!(dq, q, t, source_terms, equations, solver) - @. tmp1 = -(equations.gravity * eta + 0.5 * v^2) - mul!(dv, invImD2K_D, tmp1) + deta[:] = invImDKD * deta + dv[:] = invImD2K * dv return nothing end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 633637ba..cf15d74e 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -84,9 +84,55 @@ function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, return SVector(eta, v, D) end +""" + initial_condition_manufactured(x, t, equations::SvaerdKalischEquations1D, mesh) + +A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref). +""" +function initial_condition_manufactured(x, t, + equations::SvaerdKalischEquations1D, + mesh) + eta = exp(t) * cospi(2 * (x - 2 * t)) + v = exp(t / 2) * sinpi(2 * (x - t / 2)) + D = 3.0 + return SVector(eta, v, D) +end + +""" + source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D, mesh) + +A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref). +""" +function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D) + g = equations.gravity + D = q[3, 1] # D is constant, thus simply take the first entry + eta0 = equations.eta0 + alpha = equations.alpha + beta = equations.beta + gamma = equations.gamma + a1 = cospi(t - 2 * x) + a2 = sinpi(t - 2 * x) + a3 = cospi(4 * t - 2 * x) + a4 = sinpi(4 * t - 2 * x) + + dq1 = 8 * pi^3 * alpha * sqrt(g * (D + eta0)) * (D + eta0)^2 * exp(t) * a4 + + 2 * pi * (D + exp(t) * a3) * exp(t / 2) * a1 - 2 * pi * exp(3 * t / 2) * a2 * a4 - + 4 * pi * exp(t) * a4 + exp(t) * a3 + dq2 = 2 * pi * D * g * exp(t) * a4 - D * exp(t / 2) * a2 / 2 - + pi * D * exp(t / 2) * a1 - 2 * pi * D * exp(t) * a2 * a1 + + 8 * pi^3 * alpha * (D + eta0)^2 * sqrt(D * g + eta0 * g) * exp(3 * t / 2) * a1 * + a3 - 2 * pi^2 * beta * (D + eta0)^3 * exp(t / 2) * a2 - + 4 * pi^3 * beta * (D + eta0)^3 * exp(t / 2) * a1 + + 2 * pi * g * exp(2 * t) * a4 * a3 + + 8.0 * pi^3 * gamma * (D + eta0)^3 * sqrt(D * g + eta0 * g) * exp(t / 2) * a1 - + exp(3 * t / 2) * a2 * a3 / 2 - pi * exp(3 * t / 2) * a1 * a3 - + 2 * pi * exp(2 * t) * a2 * a1 * a3 + return SVector(dq1, dq2, 0.0) +end +# function create_cache(mesh, equations::SvaerdKalischEquations1D, - solver::Solver, + solver, initial_condition, RealT, uEltype) @@ -123,8 +169,8 @@ end # Discretization that conserves the mass (for eta and v) and is energy-bounded for periodic boundary conditions function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, - initial_condition, - ::BoundaryConditionPeriodic, solver, cache) + initial_condition, ::BoundaryConditionPeriodic, source_terms, + solver, cache) @unpack hmD1betaD1, D1betaD1, d, h, hv, alpha_hat, beta_hat, gamma_hat, tmp1, tmp2, D1_central = cache q = wrap_array(u_ode, mesh, equations, solver) dq = wrap_array(du_ode, mesh, equations, solver) @@ -137,7 +183,7 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, dD = view(dq, 3, :) fill!(dD, zero(eltype(dD))) - @. h = eta + D + h = eta .+ D hv = h .* v if solver.D1 isa PeriodicDerivativeOperator || @@ -164,18 +210,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, hmD1betaD1 = Diagonal(h) - D1betaD1 # split form - tmp2 = -(0.5 * (D1_central * (hv .* v) + hv .* D1v - v .* (D1_central * hv)) + - equations.gravity * h .* D1eta + - 0.5 * (vD1y - D1vy - yD1v) - - 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - - 0.5 * solver.D2 * (gamma_hat .* D1v)) - # not split form - # tmp2 = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ - # equations.gravity * h .* D1eta + - # vD1y - D1vy - - # 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - - # 0.5 * solver.D2 * (gamma_hat .* D1v)) - dv[:] = hmD1betaD1 \ tmp2 + dv[:] = -(0.5 * (D1_central * (hv .* v) + hv .* D1v - v .* (D1_central * hv)) + + equations.gravity * h .* D1eta + + 0.5 * (vD1y - D1vy - yD1v) - + 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - + 0.5 * solver.D2 * (gamma_hat .* D1v)) + + # no split form + # dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ + # equations.gravity * h .* D1eta + + # vD1y - D1vy - + # 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - + # 0.5 * solver.D2 * (gamma_hat .* D1v)) + + calc_sources!(dq, q, t, source_terms, equations, solver) + dv[:] = hmD1betaD1 \ dv return nothing end diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl index 05aadf7f..eda6bdb5 100644 --- a/src/semidiscretization.jl +++ b/src/semidiscretization.jl @@ -6,7 +6,7 @@ A `struct` containing everything needed to describe a spatial semidiscretization of an equation. """ struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions, - Solver, Cache} + SourceTerms, Solver, Cache} mesh::Mesh equations::Equations @@ -15,30 +15,34 @@ struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions, initial_condition::InitialCondition boundary_conditions::BoundaryConditions + source_terms::SourceTerms solver::Solver cache::Cache function Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions, - Solver, + SourceTerms, Solver, Cache}(mesh::Mesh, equations::Equations, initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, + source_terms::SourceTerms, solver::Solver, cache::Cache) where {Mesh, Equations, InitialCondition, - BoundaryConditions, Solver, - Cache} + BoundaryConditions, SourceTerms, + Solver, Cache} @assert ndims(mesh) == ndims(equations) @assert xmin(mesh) == xmin(solver.D1) @assert xmax(mesh) == xmax(solver.D1) @assert nnodes(mesh) == length(grid(solver)) - new(mesh, equations, initial_condition, boundary_conditions, solver, cache) + new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, + cache) end end """ Semidiscretization(mesh, equations, initial_condition, solver; + source_terms=nothing, boundary_conditions=boundary_condition_periodic, RealT=real(solver), uEltype=RealT, @@ -47,6 +51,7 @@ end Construct a semidiscretization of a PDE. """ function Semidiscretization(mesh, equations, initial_condition, solver; + source_terms = nothing, boundary_conditions = boundary_condition_periodic, # `RealT` is used as real type for node locations etc. # while `uEltype` is used as element type of solutions etc. @@ -57,12 +62,10 @@ function Semidiscretization(mesh, equations, initial_condition, solver; initial_cache...) Semidiscretization{typeof(mesh), typeof(equations), typeof(initial_condition), - typeof(boundary_conditions), typeof(solver), typeof(cache)}(mesh, - equations, - initial_condition, - boundary_conditions, - solver, - cache) + typeof(boundary_conditions), typeof(source_terms), + typeof(solver), typeof(cache)}(mesh, equations, initial_condition, + boundary_conditions, source_terms, + solver, cache) end function Base.show(io::IO, semi::Semidiscretization) @@ -73,6 +76,7 @@ function Base.show(io::IO, semi::Semidiscretization) print(io, ", ", semi.equations) print(io, ", ", semi.initial_condition) print(io, ", ", semi.boundary_conditions) + print(io, ", ", semi.source_terms) print(io, ", ", semi.solver) print(io, ", cache(") for (idx, key) in enumerate(keys(semi.cache)) @@ -93,7 +97,8 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) println(io, " mesh: ", semi.mesh) println(io, " equations: ", get_name(semi.equations)) println(io, " initial condition: ", semi.initial_condition) - print(io, " boundary condition: ", semi.boundary_conditions) + println(io, " boundary condition: ", semi.boundary_conditions) + print(io, " source terms: ", semi.source_terms) end end @@ -196,10 +201,11 @@ function wrap_array(u_ode, semi::Semidiscretization) end function rhs!(du_ode, u_ode, semi::Semidiscretization, t) - @unpack mesh, equations, initial_condition, boundary_conditions, solver, cache = semi + @unpack mesh, equations, initial_condition, boundary_conditions, solver, source_terms, cache = semi - rhs!(du_ode, u_ode, t, mesh, equations, initial_condition, boundary_conditions, solver, - cache) + rhs!(du_ode, u_ode, t, mesh, equations, initial_condition, boundary_conditions, + source_terms, + solver, cache) return nothing end diff --git a/src/solver.jl b/src/solver.jl index 06e8a6d8..1bc2d645 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -114,3 +114,20 @@ function calc_error_norms(u_ode, t, initial_condition, mesh::Mesh1D, equations, end return l2_error, linf_error end + +function calc_sources!(dq, q, t, source_terms::Nothing, + equations::AbstractEquations{1}, solver::Solver) + return nothing +end + +function calc_sources!(dq, q, t, source_terms, + equations::AbstractEquations{1}, solver::Solver) + x = grid(solver) + for i in eachnode(solver) + local_source = source_terms(view(q, :, i), x[i], t, equations) + for v in eachvariable(equations) + dq[v, i] += local_source[v] + end + end + return nothing +end diff --git a/test/test_bbm_bbm_1d.jl b/test/test_bbm_bbm_1d.jl index 22236e2c..f7a72458 100644 --- a/test/test_bbm_bbm_1d.jl +++ b/test/test_bbm_bbm_1d.jl @@ -41,6 +41,19 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d") change_velocity=-4.547473508864641e-13, change_entropy=0.0) end + + @trixi_testset "bbm_bbm_1d_manufactured" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_manufactured.jl"), + tspan=(0.0, 1.0), + l2=[4.365176233405813e-9 6.7151849982388e-10], + linf=[6.226559268185383e-9 9.698699621196738e-10], + cons_error=[3.873483998828586e-12 2.2986355942039745e-11], + change_waterheight=3.873483998828586e-12, + change_velocity=2.2986355942039745e-11, + change_entropy=17.387441847193436, + atol=1e-10, + atol_ints=1e-10) # in order to make CI pass + end end end # module diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl index 8bcea168..7910ab09 100644 --- a/test/test_bbm_bbm_variable_bathymetry_1d.jl +++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl @@ -54,7 +54,22 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") cons_error=[7.194245199571014e-14 1.2388735870505485e-12 0.0], change_waterheight=-1.5765166949677223e-13, change_velocity=1.8791999206735355e-13, - change_entropy=2.1316282072803006e-14) + change_entropy=2.1316282072803006e-14, + atol=1e-11) # in order to make CI pass) + end + + @trixi_testset "bbm_bbm_variable_bathymetry_1d_manufactured" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "bbm_bbm_variable_bathymetry_1d_manufactured.jl"), + tspan=(0.0, 1.0), + l2=[6.867287425380051e-9 3.446178245128195e-9 0.0], + linf=[1.1667709465257303e-8 5.917459633408839e-9 0.0], + cons_error=[1.359075785939412e-11 3.8711139735371144e-13 0.0], + change_waterheight=-1.359075785939412e-11, + change_velocity=-3.8711139735371144e-13, + change_entropy=17.81701226932122, + atol=1e-10, + atol_ints=1e-11) # in order to make CI pass end @trixi_testset "bbm_bbm_variable_bathymetry_1d_well_balanced" begin diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 935e22b4..96bc6969 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -8,6 +8,18 @@ include("test_util.jl") EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d") @testset "SvaerdKalisch1D" begin + @trixi_testset "svaerd_kalisch_1d_manufactured" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "svaerd_kalisch_1d_manufactured.jl"), + tspan=(0.0, 0.1), + l2=[7.467887060263923e-5 2.7796353838948894e-8 0.0], + linf=[0.0001613144395267163 4.344495230235168e-8 0.0], + cons_error=[2.3635607360183997e-16 8.084235123776567e-10 0.0], + change_waterheight=-2.3635607360183997e-16, + change_entropy=0.1342289500320556, + atol=1e-4) # in order to make CI pass + end + @trixi_testset "svaerd_kalisch_1d_dingemans" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_dingemans.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index ebe64afc..bee4c3ac 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -221,10 +221,9 @@ using SparseArrays: sparse, SparseMatrixCSC accuracy_orders = [2, 4, 6] for accuracy_order in accuracy_orders - eoc_mean_values, _ = convergence_test(default_example(), 2, initial_N = 128, + eoc_mean_values, _ = convergence_test(default_example(), 2, N = 512, tspan = (0.0, 1.0), accuracy_order = accuracy_order) - show(eoc_mean_values) @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5) @test isapprox(eoc_mean_values[:linf][2], accuracy_order, atol = 0.5) @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5)