Merge wet/dry capability to main (#1501)

* add dummy commit in order to open a dev to main PR * [WIP] Wet/dry capabilities for 2D shallow water equations (#1340) * HR of Chen and Noelle (1D) and edit SWE struct * Overload limiter (SWE 1D) to cut off waterheight * New indicatorHG (SWE 1D) to apply FV on dry cells * Threshold in rhs! before calculation (SWE 1D) * New lake_at_rest_error for SWE 1D * New wet/dry elixirs for testing scheme for SWE 1D * HR of Chen and Noelle (2D) and edit SWE struct * Overload limiter (SWE 2D) to cut off waterheight * New indicatorHG (SWE 2D) to apply FV on dry cells * Threshold in rhs! before calculation (SWE 2D) * New lake_at_rest_error for SWE 2D * New wet/dry elixirs for testing scheme for SWE 2D * Elixir SWE 2D: 3 mounds, problem with boundaries * Fixed MethodError; apply_thresholds! too strict * Fixed MethodError; apply_thresholds! too strict * Move threshold on volume integral in stage_limiter * Indentation, spacing and comments adjustment * Renaming numerical HLL type flux (SWE 1D) * Move threshold on volume integral in stage_limiter * Renaming numerical HLL type flux (SWE 2D) * Indentation, spacing and comments adjustment * Describing docs for Chen and Noelle HR (SWE 1D) * Edit SWE 1D elixirs, error-based solver and docs * Including tests on new SWE 1D elixirs * Describing docs for Chen and Noelle HR (SWE 2D) * Edit SWE 2D elixirs, error-based solver and docs * Including tests on new SWE 2D elixirs * New/reorganize positivity limiter (SWE 2D) * New/reorganize positivity limiter (SWE 1D) * Editing docs SWE 1D * Editing docs SWE 2D * Rearrange cut off at interfaces, edit tests SWE 1D * Edit docs, add Ref * Edit docs and indenting (SWE 2D) * Rearrange cut off at interfaces, edit tests SWE 2D * Remove tree/structured mesh elixir from repo SWE2D * Create unstructured mesh elixir SWE 2D * Add 1D lake-at-rest-error logic to pass 1D tests * Add 2D lake-at-rest-error logic to pass 2D tests * Fixed typo. Confusing name, but correct math * Correction of comments and docstrings * Correction of comments and docstrings * Rename mesh file in elixir for UnstructuredMesh * Update test_unstructured_2d.jl forgot an end statement for the new test * Fixing typos * fix dispatching error on new lake-at-rest error calculation. See if this fixes broken tests * Editing initial condition in parabolic bowl elixir * Delete unnecessary variable in elixir * adjust lake-at-rest error computation strategy. move specialized version of error into the wet-dry elixir as the new functionality was only needed in this speacial case. update corresponding test values as the bottom is now truly discontinuous * update structured mesh version of the wet-dry well-balancedness test * fix typos * update values in parabolic bowl test on StructuredMesh * update parabolic bowl test on TreeMesh * revert the 1D computation of the lake-at-rest error to the standard way. This will change once the 1D wet/dry merges * Reset lake-at-rest error computation strategy. New version of error only in wet-dry elixir (special case) Update test values as the bottom is now truly discontinuous * Fix typo * Shorten test run for parabolic bowl 1D * Choose lower resolution for parabolic bowl and update test values * Further reduce resolution for parabolic bowl and update test values * adjust special initial conditions and well-balancedness error routines to avoid the need of element IDs * Remove MPI from well-balanced test * simplify workaround to set discontinuous initial data * Simplify workaround to set discontinuity * Change structure of Chen&Noelle flux * Fix typos and indenting * Adjust call of solve and use ode_default_options * Edit docstring * Replace boolean with if, remove set_node_vars Shorten test runs on TreeMesh and UnstructuredMesh * Change structure of Chen&Noelle flux * Fix typos and indenting * Adjust call of solve and use ode_default_options * Edit docstring * Replace boolean with if, remove set_node_vars Shorten test runs on TreeMesh and UnstructuredMesh * Update comment regarding H0 for lake-at-rest error * Add the original source to the parabolic bowl test * Update comment regarding H0 for lake-at-rest error * Add the original source to the parabolic bowl test * New sc indicator especially for SWE * Remove threshold parameter from SWE limiter call * update some docstrings * remove type instability in positivty limiter * typo fix * move safety check for dry state in the new positivity limiter into the same element loop * more docstring updates * remove dummy comment added in the dev initial commit * adjust default threshold values to be precision agnostic * update comment on the default threshold value in the new TreeMesh elixirs * update comments for the three new TreeMesh examples * update IC comment for three mound test * update IC comments for new StructuredMesh2D tests * update comment on shallow water constructor * adjust comments in the shallow_water_2d file * adjust comment regarding threshold_limiter in the new elixirs * fix typos found by SpellCheck * Edit docs * Import Printf macros for printing wb error * Remove type instability in Chen & Noelle HR * Change logic for setting SC indicator to one * Change logic for default values of SWE struct * Outsource HG shock capturing indicator for SWE Create different function to compute indicator Edit comments Change wet/dry clipping to if-else logic * Move limiterthreshold into function & edit docs Threshold was a passed variable in elixir before. Now, it is taken right from the SWE struct in the limiter Edit docs * Move new limiter safety check in same element loop * Adjust default threshold values * Remove type instability * Import Printf package for terminal output * Edit docs * Add Printf package to the test/Project.toml Used for printing lake-at-rest error in well-balancedness test * Add Printf package to the test/Project.toml Used for printing lake-at-rest error in well-balancedness test * Typo fix in elixir_shallowwater_well_balanced_wet_dry.jl * Typo fix in elixir_shallowwater_well_balanced_wet_dry.jl * unify new code with required formatting * fix weird formatting and add 'format: noindent' where missing. fix crashing structured mesh run * add unit test for new show routine * apply JuliaFormatter * simplify elixir as we can set discontinuous ICs in 1D. Also update beach test values * dummy commit to check push access * remove dummy comment * typo fix --------- Co-authored-by: Andrew Winters <[email protected]> Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * adjust comments and remove duplicate code * add TODOs for code pieces that should move to TrixiShallowWater package * remove accidentally added file * apply formatter to avoid errors with new comments * move TODO comments to avoid errors in Documentation build * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * remove unnecessary analysis quantities from several new elixirs * rename local threshold variable in new indicator to avoid confusion * update NEWS.md with wetting and drying feature * fix fomartting issue from conflict resolution --------- Co-authored-by: svengoldberg <[email protected]> Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Jul 14, 2023 · 905c8e2 · 905c8e2
1 parent 0816ed0
commit 905c8e2
Show file tree

Hide file tree

Showing 32 changed files with 2,608 additions and 72 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -10,6 +10,7 @@ for human readability.
 
 - Experimental support for 3D parabolic diffusion terms has been added.
 - Capability to set truly discontinuous initial conditions in 1D.
+- Wetting and drying feature and examples for 1D and 2D shallow water equations
 
 #### Changed
 

diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl b/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl
@@ -0,0 +1,113 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+ ###############################################################################
+ # Semidiscretization of the shallow water equations
+#
+# TODO: TrixiShallowWater: wet/dry example elixir
+
+equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=1.4)
+
+"""
+    initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D)
+
+Initial condition for the [`ShallowWaterEquations2D`](@ref) to test the [`hydrostatic_reconstruction_chen_noelle`](@ref)
+and its handling of discontinuous water heights at the start in combination with wetting and
+drying. The bottom topography is given by a conical island in the middle of the domain. Around that
+island, there is a cylindrical water column at t=0 and the rest of the domain is dry. This
+discontinuous water height is smoothed by a logistic function. This simulation uses periodic
+boundary conditions.
+"""
+function initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D)
+  # Set the background values
+
+  v1 = 0.0
+  v2 = 0.0
+
+  x1, x2 = x
+  b = max(0.1, 1.0 - 4.0 * sqrt(x1^2 + x2^2))
+
+  # use a logistic function to transfer water height value smoothly
+  L  = equations.H0    # maximum of function
+  x0 = 0.3   # center point of function
+  k  = -25.0 # sharpness of transfer
+
+  H = max(b, L/(1.0 + exp(-k*(sqrt(x1^2+x2^2) - x0))))
+
+  # It is mandatory to shift the water level at dry areas to make sure the water height h
+  # 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 v1) / 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 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(H, b + equations.threshold_limiter)
+  return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_conical_island
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(4)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max=0.5,
+                                                     alpha_min=0.001,
+                                                     alpha_smooth=true,
+                                                     variable=waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the StructuredMesh and setup a periodic mesh
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0,  1.0)
+
+cells_per_dimension = (16, 16)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,))
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback=callbacks);
+
+summary_callback() # print the timer summary
diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl
@@ -0,0 +1,119 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+#
+# TODO: TrixiShallowWater: wet/dry example elixir
+
+equations = ShallowWaterEquations2D(gravity_constant=9.81)
+
+"""
+    initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D)
+
+Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its
+wet-dry mechanics. This test has an analytical solution. The initial condition is defined by the
+analytical solution at time t=0. The bottom topography defines a bowl and the water level is given
+by an oscillating lake.
+
+The original test and its analytical solution were first presented in
+- William C. Thacker (1981)
+  Some exact solutions to the nonlinear shallow-water wave equations
+  [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882).
+
+The particular setup below is taken from Section 6.2 of
+- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018)
+  An entropy stable discontinuous Galerkin method for the shallow water equations on
+  curvilinear meshes with wet/dry fronts accelerated by GPUs
+  [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038).
+"""
+function initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D)
+  a = 1.0
+  h_0 = 0.1
+  sigma = 0.5
+  ω = sqrt(2 * equations.gravity * h_0) / a
+
+  v1 = -sigma * ω * sin(ω * t)
+  v2 = sigma * ω * cos(ω * t)
+
+  b = h_0 * ((x[1])^2 + (x[2])^2) / a^2
+
+  H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) + 2 * x[2] * sin(ω * t) - sigma) + h_0
+
+  # It is mandatory to shift the water level at dry areas to make sure the water height h
+  # 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 v1) / 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 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(H, b + equations.threshold_limiter)
+  return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_parabolic_bowl
+
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(4)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max=0.6,
+                                                     alpha_min=0.001,
+                                                     alpha_smooth=true,
+                                                     variable=waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+
+###############################################################################
+
+coordinates_min = (-2.0, -2.0)
+coordinates_max = (2.0, 2.0)
+
+cells_per_dimension = (150, 150)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false,
+                                     extra_analysis_integrals=(energy_kinetic,
+                                                               energy_internal))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback=callbacks);
+
+summary_callback() # print the timer summary