Move SWEs from Trixi.jl to TrixiShallowWater.jl

.github/workflows/Documenter.yml

@@ -31,7 +31,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1'
+          version: '1.10'
           show-versioninfo: true
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1

.github/workflows/ci.yml

@@ -51,6 +51,8 @@ jobs:
           - tree_2d
           - structured_2d
           - unstructured_2d
+          - t8code_2d
+          - dgmulti
           - p4est_2d
           - unit
           - upstream
@@ -99,6 +101,8 @@ jobs:
           - tree_2d
           - structured_2d
           - unstructured_2d
+          - t8code_2d
+          - dgmulti
           - p4est_2d
           - unit
     steps:

NEWS.md

@@ -5,6 +5,20 @@ TrixiShallowWater.jl follows the interpretation of
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.2 from v0.1.x
+
+#### Added
+- New equation `ShallowWaterEquationsQuasi1D` and functions `FluxHydrostaticReconstruction`, 
+  `flux_nonconservative_audusse_etal`, and `hydrostatic_reconstruction_audusse_etal` are now available
+  through TrixiShallowWater.jl instead of Trixi.jl. ([#96])
+
+#### Changed
+- `ShallowWaterEquationsWetDry` have been renamed to `ShallowWaterEquations`. The source code
+  for these equations is now implemented directly in TrixiShallowWater.jl ([#96]).
+
+#### Deprecated
+
+#### Removed
 
 ## Changes in the v0.1 lifecycle
 

Project.toml

@@ -19,5 +19,5 @@ Printf = "1"
 Roots = "2.1.6"
 Static = "1.1.1"
 StaticArrays = "1.9"
-Trixi = "0.11.10"
-julia = "1.10"
+Trixi = "0.12"
+julia = "1.10"

docs/src/tutorials/elixir_shallowwater_dam_break_triangular.jl

@@ -13,12 +13,12 @@
 # - Set up a SWE solver for wet/dry transitions
 # - Create custom initial conditions and source terms
 # - Save solution data at gauge points
-# - Visualization with [`Makie.jl`](https://docs.makie.org/dev/)
+# - Visualization with [Makie.jl](https://docs.makie.org/dev/)
 
 # ## Load required packages
 # Before we start, we need to load the required packages. Besides TrixiShallowWater.jl, we require
-# [`Trixi.jl`](@extref Trixi.jl) for the spatial discretization and [`OrdinaryDiffEqSSPRK.jl`](https://docs.sciml.ai/OrdinaryDiffEq/stable/) for time integration.
-# In addition to that [`CairoMakie.jl`](https://docs.makie.org/dev/) is used for visualization and [`CSV.jl`](https://csv.juliadata.org/stable/) and [`DataFrames.jl`](https://dataframes.juliadata.org/stable/) will
+# [Trixi.jl](@extref Trixi.jl) for the spatial discretization and [OrdinaryDiffEqSSPRK.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) for time integration.
+# In addition to that [CairoMakie.jl](https://docs.makie.org/dev/) is used for visualization and [CSV.jl](https://csv.juliadata.org/stable/) and [DataFrames.jl](https://dataframes.juliadata.org/stable/) will
 # be used to load the experimental data.
 
 # Standard packages
@@ -33,16 +33,16 @@ using CSV
 # ## Prepare and run the problem setup
 
 # In the first step we will set up the equation system. In this example we want to solve the
-# one-dimensional shallow water equations, so we will use the [`ShallowWaterEquationsWetDry1D`](@ref ShallowWaterEquationsWetDry1D)
+# one-dimensional shallow water equations, so we will use the [`ShallowWaterEquations1D`](@ref)
 # and specify the gravitational acceleration to `gravity = 9.812`. In contrast to the
-# [`Trixi.ShallowWaterEquations1D`](@extref Trixi.ShallowWaterEquations1D) type, this equation type
+# [`ShallowWaterEquations1D`](@ref) type, this equation type
 # contains additional parameters and methods that are needed to handle wetting and drying.
-equations = ShallowWaterEquationsWetDry1D(gravity = 9.812)
+equations = ShallowWaterEquations1D(gravity = 9.812)
 
 # We then create a function to supply the initial condition for the simulation. Note, in the last
 # step of this function the water height needs to be shifted by a small value to avoid division by zero.
 function initial_condition_dam_break_triangular(x, t,
-                                                equations::ShallowWaterEquationsWetDry1D)
+                                                equations::ShallowWaterEquations1D)
     b = 0.0  # Bottom topography
     h = 0.0  # Water height
     v = 0.0  # Velocity
@@ -75,7 +75,7 @@ initial_condition = initial_condition_dam_break_triangular;
 # As we want to compare the results to experimental data, we also need to account for bottom friction.
 # For this we create a new source term, which adds a Manning friction term to the momentum equation.
 @inline function source_term_manning_friction(u, x, t,
-                                              equations::ShallowWaterEquationsWetDry1D)
+                                              equations::ShallowWaterEquations1D)
     h, hv, _ = u
 
     n = 0.0125  # friction coefficient

docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl

@@ -30,12 +30,12 @@
 
 # ## Load required packages
 # The core solver component is TrixiShallowWater.jl,
-# which requires [`Trixi.jl`](@extref Trixi.jl) for the underlying spatial discretization
-# and `OrdinaryDiffEqSSPRK.jl` for time integration.
-# `HOHQMesh.jl` is needed to generate an unstructured mesh for this problem.
-# `TrixiBottomTopography.jl` is needed to create a bathymetry approximation that is directly
+# which requires [Trixi.jl](@extref Trixi.jl) for the underlying spatial discretization
+# and OrdinaryDiffEqSSPRK.jl for time integration.
+# HOHQMesh.jl is needed to generate an unstructured mesh for this problem.
+# TrixiBottomTopography.jl is needed to create a bathymetry approximation that is directly
 # usable by Trixi.jl.
-# Finally, we include [`CairoMakie.jl`](https://docs.makie.org/stable/) for insitu visualization and `Trixi2Vtk.jl` for postprocessing.
+# Finally, we include [CairoMakie.jl](https://docs.makie.org/stable/) for insitu visualization and Trixi2Vtk.jl for postprocessing.
 using HOHQMesh
 using OrdinaryDiffEqSSPRK
 using Trixi
@@ -152,12 +152,12 @@ generate_mesh(monai);
 # for the tsunami runup problem.
 
 # For this example we solve the two-dimensional shallow water equations,
-# so we use the [`ShallowWaterEquationsWetDry2D`](@ref ShallowWaterEquationsWetDry2D)
+# so we use the [`ShallowWaterEquations2D`](@ref)
 # and specify the gravitational acceleration to `gravity = 9.812`
 # as well as a background water height `H0 = 0.0`.
-# In contrast to the [`Trixi.ShallowWaterEquations2D`](@extref Trixi.ShallowWaterEquations2D) type,
-# this equation type contains additional parameters and methods needed to handle wetting and drying.
-equations = ShallowWaterEquationsWetDry2D(gravity = 9.81, H0 = 0.0)
+# In contrast to the [`ShallowWaterEquations2D`](@ref) type,
+# this equation type allows contains additional parameters and methods needed to handle wetting and drying.
+equations = ShallowWaterEquations2D(gravity = 9.81, H0 = 0.0)
 
 # Next, we construct an approximation to the bathymetry with TrixiBottomTopography.jl using
 # a [`BicubicBSpline`](https://trixi-framework.github.io/TrixiBottomTopography.jl/stable/reference/#TrixiBottomTopography.BicubicBSpline)
@@ -176,7 +176,7 @@ bathymetry(x::Float64, y::Float64) = spline_interpolation(bath_spline_struct, x,
 
 # We then create a function to supply the initial condition for the simulation.
 @inline function initial_condition_monai_tsunami(x, t,
-                                                 equations::ShallowWaterEquationsWetDry2D)
+                                                 equations::ShallowWaterEquations2D)
     ## Initially water is at rest
     v1 = 0.0
     v2 = 0.0
@@ -221,7 +221,7 @@ H_from_wave_maker(t::Float64) = spline_interpolation(h_spline_struct, t);
 # wave maker.
 @inline function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector,
                                                x, t, surface_flux_functions,
-                                               equations::ShallowWaterEquationsWetDry2D)
+                                               equations::ShallowWaterEquations2D)
     ## Extract the numerical flux functions to compute the conservative and nonconservative
     ## pieces of the approximation
     surface_flux_function, nonconservative_flux_function = surface_flux_functions
@@ -263,7 +263,7 @@ boundary_condition = Dict(:Bottom => boundary_condition_slip_wall,
 # For this application, we also need to model the bottom friction.
 # Thus, we create a new source term, which adds a Manning friction term to the momentum equations.
 @inline function source_terms_manning_friction(u, x, t,
-                                               equations::ShallowWaterEquationsWetDry2D)
+                                               equations::ShallowWaterEquations2D)
     h, hv_1, hv_2, _ = u
 
     n = 0.001 # friction coefficient

examples/dgmulti_1d/elixir_shallow_water_quasi_1d.jl

@@ -0,0 +1,49 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the quasi 1d shallow water equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = ShallowWaterEquationsQuasi1D(gravity = 9.81)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = (FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs()),
+                flux_nonconservative_chan_etal)
+
+dg = DGMulti(polydeg = 4, element_type = Line(), approximation_type = SBP(),
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+cells_per_dimension = (8,)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (0.0,), coordinates_max = (sqrt(2),),
+                   periodicity = true)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8,
+            ode_default_options()..., callback = callbacks)

examples/dgmulti_2d/elixir_shallowwater_source_terms.jl

@@ -0,0 +1,46 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations2D(gravity = 9.81)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal)
+dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+cells_per_dimension = (8, 8)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (0.0, 0.0), coordinates_max = (sqrt(2), sqrt(2)),
+                   periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-7, reltol = 1.0e-7,
+            ode_default_options()..., callback = callbacks);

examples/p4est_2d_dgsem/elixir_shallowwater_perturbation_amr.jl

@@ -7,9 +7,9 @@ using TrixiShallowWater
 # semidiscretization of the shallow water equations with a continuous
 # bottom topography function and a perturbation in the water height
 
-equations = ShallowWaterEquationsWetDry2D(gravity = 9.812, H0 = 2.1)
+equations = ShallowWaterEquations2D(gravity = 9.812, H0 = 2.1)
 
-function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWetDry2D)
+function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D)
     # Calculate primitive variables
     H = equations.H0
     v1 = 0.0
@@ -111,7 +111,7 @@ save_solution = SaveSolutionCallback(dt = 0.04,
                                      save_final_solution = true)
 
 # Define the perturbation of water height as a  variable to use in the AMR indicator
-@inline function waterheight_total(u, equations::ShallowWaterEquationsWetDry2D)
+@inline function waterheight_total(u, equations::ShallowWaterEquations2D)
     return u[1] + u[4]
 end
 

examples/p4est_2d_dgsem/elixir_shallowwater_perturbation_wet_dry_amr.jl

@@ -17,11 +17,11 @@ using TrixiShallowWater
 # semidiscretization of the shallow water equations with a discontinuous
 # bottom topography function for a perturbed water height on a nonconforming mesh with AMR
 
-equations = ShallowWaterEquationsWetDry2D(gravity = 9.812, H0 = 1.235,
-                                          threshold_limiter = 1e-12,
-                                          threshold_desingularization = 1e-4)
+equations = ShallowWaterEquations2D(gravity = 9.812, H0 = 1.235,
+                                    threshold_limiter = 1e-12,
+                                    threshold_desingularization = 1e-4)
 
-function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWetDry2D)
+function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D)
     # Calculate primitive variables
     H = equations.H0
     v1 = 0.0
@@ -38,7 +38,7 @@ function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWe
     # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
     # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
     # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above
-    # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0.
+    # for the ShallowWaterEquations and added to the initial condition if h = 0.
     # This default value can be changed within the constructor call depending on the simulation setup.
     H = max(equations.H0, b + equations.threshold_limiter)
 
@@ -124,7 +124,7 @@ tspan = (0.0, 4.0)
 ode = semidiscretize(semi, tspan)
 
 function initial_condition_discontinuous_perturbation(x, t, element_id,
-                                                      equations::ShallowWaterEquationsWetDry2D)
+                                                      equations::ShallowWaterEquations2D)
     # Set the background values for velocity
     H = equations.H0
     v1 = zero(H)
@@ -158,7 +158,7 @@ function initial_condition_discontinuous_perturbation(x, t, element_id,
 
     # Avoid division by zero by adjusting the initial condition with a small dry state threshold
     # that defaults to 500*eps() ≈ 1e-13 in double precision and is set in the constructor above
-    # for the ShallowWaterEquationsWetDry struct.
+    # for the ShallowWaterEquations struct.
     H = max(H, b + equations.threshold_limiter)
     return prim2cons(SVector(H, v1, v2, b), equations)
 end
@@ -193,7 +193,7 @@ save_solution = SaveSolutionCallback(dt = 0.2,
 
 # Another possible AMR indicator function could be the velocity, such that it only fires
 # in regions where the water is moving
-# @inline function velocity_norm(u, equations::ShallowWaterEquationsWetDry2D)
+# @inline function velocity_norm(u, equations::ShallowWaterEquations2D)
 #    v1, v2 = velocity(u, equations)
 #    return sqrt(v1^2 + v2^2)
 # end

examples/p4est_2d_dgsem/elixir_shallowwater_source_terms.jl

@@ -0,0 +1,61 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations2D(gravity = 9.81)
+
+initial_condition = initial_condition_convergence_test # MMS EOC test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the P4estMesh and setup a periodic mesh
+
+coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (sqrt(2.0), sqrt(2.0))  # maximum coordinates (max(x), max(y))
+
+# Create P4estMesh with 8 x 8 trees and 16 x 16 elements
+trees_per_dimension = (8, 8)
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 initial_refinement_level = 1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 200,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8,
+            ode_default_options()..., callback = callbacks);