From 99583ab70a616ea5f409a5e0f4c980ef052b2b08 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Aug 2024 10:15:52 +0200 Subject: [PATCH] Use FastBroadcast.jl (#143) * add hyperbolic_serre_green_naghdi_conservation.jl * use FastBroadcast.jl for hyperbolic SGN * minor performance tuning of hyperbolic SGN * use FastBroadcast.jl for standard SGN * relax FastBroadcast.jl compat to 0.3 * allow FastBroadcast 0.2.8 Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --------- Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- Project.toml | 2 + ...erbolic_serre_green_naghdi_conservation.jl | 75 +++++ .../serre_green_naghdi_conservation.jl | 1 - src/DispersiveShallowWater.jl | 1 + .../hyperbolic_serre_green_naghdi_1d.jl | 57 ++-- src/equations/serre_green_naghdi_1d.jl | 256 +++++++++--------- test/test_hyperbolic_serre_green_naghdi_1d.jl | 30 ++ 7 files changed, 261 insertions(+), 161 deletions(-) create mode 100644 examples/hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_conservation.jl diff --git a/Project.toml b/Project.toml index 43ba140d..f4342347 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.4.4-pre" [deps] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PolynomialBases = "c74db56a-226d-5e98-8bb0-a6049094aeea" @@ -25,6 +26,7 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" [compat] BandedMatrices = "1" DiffEqBase = "6.143" +FastBroadcast = "0.2.8, 0.3" Interpolations = "0.14.2, 0.15" LinearAlgebra = "1" PolynomialBases = "0.4.15" diff --git a/examples/hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_conservation.jl b/examples/hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_conservation.jl new file mode 100644 index 00000000..79a41211 --- /dev/null +++ b/examples/hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_conservation.jl @@ -0,0 +1,75 @@ +# This elixir contains an artificial setup that can be used to check the +# conservation properties of the equations and numerical methods as well as +# a possible directional bias (if the velocity is set to zero). See +# - Hendrik Ranocha and Mario Ricchiuto (2024) +# Structure-preserving approximations of the Serre-Green-Naghdi +# equations in standard and hyperbolic form +# [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) + +using OrdinaryDiffEq +using DispersiveShallowWater + +############################################################################### +# Semidiscretization of the hyperbolic Serre-Green-Naghdi equations + +bathymetry_type = bathymetry_mild_slope +equations = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type; + lambda = 500.0, + gravity_constant = 9.81, + eta0 = 1.0) + +function initial_condition_conservation_test(x, t, + equations::HyperbolicSerreGreenNaghdiEquations1D, + mesh) + eta = 1 + exp(-x^2) + v = 1.0e-2 # set this to zero to test a directional bias + b = 0.25 * cospi(x / 75) + + # We use the feature that we can only return the physical variables + # used by the `hyperbolic_approximation_limit`, i.e., the + # `SerreGreenNaghdiEquations1D.` + D = equations.eta0 - b + return SVector(eta, v, D) +end + +# create homogeneous mesh +coordinates_min = -150.0 +coordinates_max = +150.0 +N = 1_000 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 2 +accuracy_order = 2 +solver = Solver(mesh, accuracy_order) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, + initial_condition_conservation_test, solver; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 35.0) +ode = semidiscretize(semi, tspan) + +# The callbacks support an additional `io` argument to write output to a file +# or any other IO stream. The default is stdout. We use this here to enable +# setting it to `devnull` to benchmark the full simulation including the time +# to compute the errors etc. but without the time to write the output to the +# terminal. +io = stdout +summary_callback = SummaryCallback(io) +analysis_callback = AnalysisCallback(semi; interval = 10, io, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy, + entropy_modified)) + +callbacks = CallbackSet(analysis_callback, summary_callback) + +# optimized time integration methods like this one are much more efficient +# for stiff problems (λ big) than standard methods like Tsit5() +alg = RDPK3SpFSAL35() +sol = solve(ode, alg; + save_everystep = false, callback = callbacks) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl index 709f3958..5ddbcc9d 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl @@ -8,7 +8,6 @@ using OrdinaryDiffEq using DispersiveShallowWater -using SummationByPartsOperators: upwind_operators, periodic_derivative_operator ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 54fad696..67967ff4 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -17,6 +17,7 @@ module DispersiveShallowWater 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 diff --git a/src/equations/hyperbolic_serre_green_naghdi_1d.jl b/src/equations/hyperbolic_serre_green_naghdi_1d.jl index eff0810a..ec8b1d3e 100644 --- a/src/equations/hyperbolic_serre_green_naghdi_1d.jl +++ b/src/equations/hyperbolic_serre_green_naghdi_1d.jl @@ -137,11 +137,11 @@ function set_approximation_variables!(q, mesh, eta, v, D, w, H = q.x # H ≈ h - @. H = eta + D - equations.eta0 # h = (h + b) + (eta0 - b) - eta0 + @.. H = eta + D - equations.eta0 # h = (h + b) + (eta0 - b) - eta0 # w ≈ -h v_x mul!(w, D1, v) - @. w = -H * w + @.. w = -H * w return nothing end @@ -322,8 +322,8 @@ function rhs!(dq, q, t, mesh, # Compute all derivatives required below (; h, b, b_x, H_over_h, h_x, v_x, hv_x, v2_x, h_hpb_x, H_x, H2_h_x, w_x, hvw_x, tmp) = cache - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b if !(bathymetry_type isa BathymetryFlat) mul!(b_x, D1, b) end @@ -335,37 +335,37 @@ function rhs!(dq, q, t, mesh, mul!(v_x, D1, v) # hv2_x = D1 * (h * v) - @. tmp = h * v + @.. tmp = h * v mul!(hv_x, D1, tmp) # v2_x = D1 * (v.^2) - @. tmp = v^2 + @.. tmp = v^2 mul!(v2_x, D1, tmp) # h_hpb_x = D1 * (h .* eta) - @. tmp = h * eta + @.. tmp = h * eta mul!(h_hpb_x, D1, tmp) # H_x = D1 * H mul!(H_x, D1, H) # H2_h_x = D1 * (H^2 / h) - @. H_over_h = H / h - @. tmp = H * H_over_h + @.. H_over_h = H / h + @.. tmp = H * H_over_h mul!(H2_h_x, D1, tmp) # w_x = D1 * w mul!(w_x, D1, w) # hvw_x = D1 * (h * v * w) - @. tmp = h * v * w + @.. tmp = h * v * w mul!(hvw_x, D1, tmp) # Plain: h_t + (h v)_x = 0 # # Split form for energy conservation: # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) + @.. dh = -(h_x * v + h * v_x) # Plain: h v_t + h v v_x + g (h + b) h_x # + ... = 0 @@ -377,16 +377,14 @@ function rhs!(dq, q, t, mesh, # + λ/2 b_x - λ/2 H/h b_x = 0 lambda_6 = lambda / 6 lambda_3 = lambda / 3 - @. dv = -(g * h_hpb_x - g * (h + b) * h_x - + - 0.5 * h * v2_x - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - 0.5 * h * v * v_x - + lambda_6 * (H_over_h * H_over_h * h_x - H2_h_x) - + lambda_3 * (1 - H_over_h) * H_x) / h + @.. dv = -(g * h_hpb_x - g * eta * h_x + + 0.5 * (h * v2_x - v^2 * h_x) + + 0.5 * v * (hv_x - h * v_x) + + lambda_6 * (H_over_h^2 * h_x - H2_h_x) + + lambda_3 * (1 - H_over_h) * H_x) / h if !(bathymetry_type isa BathymetryFlat) lambda_2 = lambda / 2 - @. dv -= lambda_2 * (1 - H_over_h) * b_x / h + @.. dv -= lambda_2 * (1 - H_over_h) * b_x / h end # Plain: h w_t + h v w_x = λ - λ H / h @@ -394,19 +392,14 @@ function rhs!(dq, q, t, mesh, # Split form for energy conservation: # h w_t + 1/2 (h v w)_x + 1/2 h v w_x # - 1/2 h_x v w - 1/2 h w v_x = λ - λ H / h - @. dw = (-(0.5 * hvw_x - + - 0.5 * h * v * w_x - - - 0.5 * h_x * v * w - - - 0.5 * h * w * v_x) + lambda * (1 - H_over_h)) / h + @.. dw = (-0.5 * (hvw_x + h * v * w_x + dh * w) + + lambda * (1 - H_over_h)) / h # No special split form for energy conservation required: # H_t + v H_x + 3/2 v b_x = w - @. dH = -v * H_x + w + @.. dH = w - v * H_x if !(bathymetry_type isa BathymetryFlat) - @. dH -= 1.5 * v * b_x + @.. dH -= 1.5 * v * b_x end end @@ -468,14 +461,14 @@ function energy_total_modified(q_global, # `q_global` is an `ArrayPartition`. It collects the individual arrays for # the total water height `eta = h + b` and the velocity `v`. eta, v, D, w, H = q_global.x - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b e = zero(h) # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2 + λ/6 h (1 - H/h)^2 - @. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * w^2 + - lambda / 6 * h * (1 - H / h)^2 + @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * w^2 + + lambda / 6 * h * (1 - H / h)^2 return e end diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index af9cc82f..a9db959d 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -258,7 +258,7 @@ function create_cache(mesh, # had access to the initial condition and the exact value of the # bathymetry here. let x = grid(D1) - @. b_x = x^3 + @.. b_x = x^3 end if D1 isa PeriodicUpwindOperators @@ -313,9 +313,9 @@ function assemble_system_matrix!(cache, h, D1, D1mat, equations::SerreGreenNaghdiEquations1D{BathymetryFlat}) (; M_h, M_h3_3) = cache - @. M_h = h + @.. M_h = h scale_by_mass_matrix!(M_h, D1) - @. M_h3_3 = (1 / 3) * h^3 + @.. M_h3_3 = (1 / 3) * h^3 scale_by_mass_matrix!(M_h3_3, D1) # Floating point errors accumulate a bit and the system matrix @@ -335,12 +335,12 @@ function assemble_system_matrix!(cache, h, b_x, D1, D1mat, elseif equations.bathymetry_type isa BathymetryVariable factor = 1.0 end - @. M_h_p_h_bx2 = h + factor * h * b_x^2 + @.. M_h_p_h_bx2 = h + factor * h * b_x^2 scale_by_mass_matrix!(M_h_p_h_bx2, D1) inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 + @.. M_h3_3 = inv3 * h^3 scale_by_mass_matrix!(M_h3_3, D1) - @. M_h2_bx = 0.5 * h^2 * b_x + @.. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to @@ -419,50 +419,50 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b mul!(h_x, D1, h) mul!(v_x, D1, v) - @. tmp = h^2 + @.. tmp = h^2 mul!(h2_x, D1, tmp) - @. tmp = h * v + @.. tmp = h * v mul!(hv_x, D1, tmp) - @. tmp = v^2 + @.. tmp = v^2 mul!(v2_x, D1, tmp) - @. tmp = h^2 * v * v_x + @.. tmp = h^2 * v * v_x mul!(h2_v_vx_x, D1, tmp) - @. tmp = h * v_x + @.. tmp = h * v_x mul!(h_vx_x, D1, tmp) inv6 = 1 / 6 - @. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x - - - inv6 * h * h2_v_vx_x - - - inv6 * h^2 * v * h_vx_x) + @.. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x + - + inv6 * h * h2_v_vx_x + - + inv6 * h^2 * v * h_vx_x) mul!(p_x, D1, tmp) # Plain: h_t + (h v)_x = 0 # # Split form for energy conservation: # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) + @.. dh = -(h_x * v + h * v_x) # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -(g * h2_x - g * h * h_x - + - 0.5 * h * v2_x - - - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - - - 0.5 * h * v * v_x - + - p_x) + @.. tmp = -(g * h2_x - g * h * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + + p_x) end # The code below is equivalent to @@ -502,53 +502,53 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b mul!(h_x, D1, h) mul!(v_x, D1, v) mul!(v_x_upwind, D1_upwind.minus, v) - @. tmp = h^2 + @.. tmp = h^2 mul!(h2_x, D1, tmp) - @. tmp = h * v + @.. tmp = h * v mul!(hv_x, D1, tmp) - @. tmp = v^2 + @.. tmp = v^2 mul!(v2_x, D1, tmp) - @. tmp = h^2 * v * v_x + @.. tmp = h^2 * v * v_x mul!(h2_v_vx_x, D1, tmp) - @. tmp = h * v_x + @.. tmp = h * v_x mul!(h_vx_x, D1, tmp) # p_+ - @. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind + @.. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind mul!(p_x, D1_upwind.plus, tmp) # p_0 minv6 = -1 / 6 - @. tmp = minv6 * (h * h2_v_vx_x - + - h^2 * v * h_vx_x) + @.. tmp = minv6 * (h * h2_v_vx_x + + + h^2 * v * h_vx_x) mul!(p_x, D1, tmp, 1.0, 1.0) # Plain: h_t + (h v)_x = 0 # # Split form for energy conservation: # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) + @.. dh = -(h_x * v + h * v_x) # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -(g * h2_x - g * h * h_x - + - 0.5 * h * v2_x - - - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - - - 0.5 * h * v * v_x - + - p_x) + @.. tmp = -(g * h2_x - g * h * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + + p_x) end # The code below is equivalent to @@ -590,70 +590,70 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, (; psi) = cache end - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b mul!(b_x, D1, b) mul!(h_x, D1, h) mul!(v_x, D1, v) - @. tmp = h * eta + @.. tmp = h * eta mul!(h_hpb_x, D1, tmp) - @. tmp = h * v + @.. tmp = h * v mul!(hv_x, D1, tmp) - @. tmp = v^2 + @.. tmp = v^2 mul!(v2_x, D1, tmp) - @. tmp = h^2 * v * v_x + @.. tmp = h^2 * v * v_x mul!(h2_v_vx_x, D1, tmp) - @. tmp = h * v_x + @.. tmp = h * v_x mul!(h_vx_x, D1, tmp) inv6 = 1 / 6 - @. p_h = (0.5 * h * (h * v_x + h_x * v) * v_x - - - inv6 * h2_v_vx_x - - - inv6 * h * v * h_vx_x) - @. tmp = h * b_x * v^2 + @.. p_h = (0.5 * h * (h * v_x + h_x * v) * v_x + - + inv6 * h2_v_vx_x + - + inv6 * h * v * h_vx_x) + @.. tmp = h * b_x * v^2 mul!(p_x, D1, tmp) - @. p_h += 0.25 * p_x + @.. p_h += 0.25 * p_x if equations.bathymetry_type isa BathymetryVariable - @. psi = 0.125 * p_x + @.. psi = 0.125 * p_x end - @. tmp = b_x * v + @.. tmp = b_x * v mul!(p_x, D1, tmp) - @. p_h += 0.25 * h * v * p_x + @.. p_h += 0.25 * h * v * p_x if equations.bathymetry_type isa BathymetryVariable - @. psi += 0.125 * h * v * p_x + @.. psi += 0.125 * h * v * p_x end - @. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v + @.. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v if equations.bathymetry_type isa BathymetryVariable - @. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v + @.. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v end - @. tmp = p_h * h + @.. tmp = p_h * h mul!(p_x, D1, tmp) # Plain: h_t + (h v)_x = 0 # # Split form for energy conservation: # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) + @.. dh = -(h_x * v + h * v_x) # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -(g * h_hpb_x - g * eta * h_x - + - 0.5 * h * v2_x - - - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - - - 0.5 * h * v * v_x - + p_x - + 1.5 * p_h * b_x) + @.. tmp = -(g * h_hpb_x - g * eta * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) if equations.bathymetry_type isa BathymetryVariable - @. tmp = tmp - psi * b_x + @.. tmp = tmp - psi * b_x end end @@ -698,76 +698,76 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, (; psi) = cache end - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b mul!(b_x, D1, b) mul!(h_x, D1, h) mul!(v_x, D1, v) mul!(v_x_upwind, D1_upwind.minus, v) - @. tmp = h * eta + @.. tmp = h * eta mul!(h_hpb_x, D1, tmp) - @. tmp = h * v + @.. tmp = h * v mul!(hv_x, D1, tmp) - @. tmp = v^2 + @.. tmp = v^2 mul!(v2_x, D1, tmp) - @. tmp = h^2 * v * v_x + @.. tmp = h^2 * v * v_x mul!(h2_v_vx_x, D1, tmp) - @. tmp = h * v_x + @.. tmp = h * v_x mul!(h_vx_x, D1, tmp) # p_0 minv6 = -1 / 6 - @. p_h = minv6 * (h2_v_vx_x - + - h * v * h_vx_x) - @. tmp = h * b_x * v^2 + @.. p_h = minv6 * (h2_v_vx_x + + + h * v * h_vx_x) + @.. tmp = h * b_x * v^2 mul!(p_x, D1, tmp) - @. p_h += 0.25 * p_x + @.. p_h += 0.25 * p_x if equations.bathymetry_type isa BathymetryVariable - @. psi = 0.125 * p_x + @.. psi = 0.125 * p_x end - @. tmp = b_x * v + @.. tmp = b_x * v mul!(p_x, D1, tmp) - @. p_h += 0.25 * h * v * p_x + @.. p_h += 0.25 * h * v * p_x if equations.bathymetry_type isa BathymetryVariable - @. psi += 0.125 * h * v * p_x + @.. psi += 0.125 * h * v * p_x end - @. p_0 = p_h * h + @.. p_0 = p_h * h mul!(p_x, D1, p_0) # p_+ - @. tmp = (0.5 * h * (h * v_x + h_x * v) * v_x_upwind - - - 0.25 * (h_x * v + h * v_x) * b_x * v) + @.. tmp = (0.5 * h * (h * v_x + h_x * v) * v_x_upwind + - + 0.25 * (h_x * v + h * v_x) * b_x * v) if equations.bathymetry_type isa BathymetryVariable - @. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v + @.. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v end - @. p_h = p_h + tmp - @. tmp = tmp * h + @.. p_h = p_h + tmp + @.. tmp = tmp * h mul!(p_x, D1_upwind.plus, tmp, 1.0, 1.0) # Plain: h_t + (h v)_x = 0 # # Split form for energy conservation: # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) + @.. dh = -(h_x * v + h * v_x) # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -(g * h_hpb_x - g * eta * h_x - + - 0.5 * h * v2_x - - - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - - - 0.5 * h * v * v_x - + p_x - + 1.5 * p_h * b_x) + @.. tmp = -(g * h_hpb_x - g * eta * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) if equations.bathymetry_type isa BathymetryVariable - @. tmp = tmp - psi * b_x + @.. tmp = tmp - psi * b_x end end @@ -857,8 +857,8 @@ function energy_total_modified(q_global, # the total water height `eta = h + b` and the velocity `v`. eta, v = q_global.x let D = q_global.x[3] - @. b = equations.eta0 - D - @. h = eta - b + @.. b = equations.eta0 - D + @.. h = eta - b end N = length(v) @@ -877,7 +877,7 @@ function energy_total_modified(q_global, fill!(b_x, zero(eltype(b_x))) else (; b, b_x) = cache - @. b = equations.eta0 - q_global.x[3] + @.. b = equations.eta0 - q_global.x[3] if D1 isa PeriodicUpwindOperators mul!(b_x, D1.central, b) else @@ -885,9 +885,9 @@ function energy_total_modified(q_global, end end - @. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 + @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 if equations.bathymetry_type isa BathymetryVariable - @. e += 1 / 8 * h * (v * b_x)^2 + @.. e += 1 / 8 * h * (v * b_x)^2 end return e diff --git a/test/test_hyperbolic_serre_green_naghdi_1d.jl b/test/test_hyperbolic_serre_green_naghdi_1d.jl index dbcf8ed8..91fb2160 100644 --- a/test/test_hyperbolic_serre_green_naghdi_1d.jl +++ b/test/test_hyperbolic_serre_green_naghdi_1d.jl @@ -186,6 +186,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "hyperbolic_serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=1_000) end + + @trixi_testset "hyperbolic_serre_green_naghdi_conservation.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "hyperbolic_serre_green_naghdi_conservation.jl"), + l2=[ + 1.3601939454962908, + 2.367812627978776, + 3.565537895102537e-14, + 0.8095219703329345, + 1.3613680507825028, + ], + linf=[ + 1.0010230791351136, + 0.7870081638593119, + 9.325873406851315e-15, + 0.24169078965789928, + 1.0010090485188015, + ], + cons_error=[6.230038707144558e-11, + 0.00029618251641450044, + 4.547473508864641e-13, + 0.00759853016721156, + 0.0014395235644997229], + change_entropy=-0.20915006380346313, + change_entropy_modified=-0.09431070096138683, + atol=1e-11, # to make CI pass + atol_ints=4e-9) # to make CI pass + + @test_allocations(semi, sol, allocs=1_000) + end end end # module