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 <[email protected]>

---------

Co-authored-by: Joshua Lampert <[email protected]>

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"

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)

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

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

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,36 +377,29 @@ 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
         #
         # 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