diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl b/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl index b0eb2a2763..ec3b767caf 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl @@ -4,6 +4,8 @@ using Trixi ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=1.4) diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl index cabf3e18ba..a1734705d2 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -4,6 +4,8 @@ using Trixi ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations2D(gravity_constant=9.81) diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index c75b2e9d12..b18b02e0b4 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -5,6 +5,8 @@ using Printf: @printf, @sprintf ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations2D(gravity_constant=9.812) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index c1cfa350ef..a4c65e21b1 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -4,6 +4,8 @@ using Trixi ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations1D(gravity_constant=9.812) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 7836b02ef6..4b641df41c 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -4,6 +4,8 @@ using Trixi ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations1D(gravity_constant=9.81) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 550c7db390..8de46c6179 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -5,6 +5,8 @@ using Printf: @printf, @sprintf ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations1D(gravity_constant=9.812) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl b/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl index f7fe4edcd9..41d0419d39 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl @@ -4,6 +4,8 @@ using Trixi ############################################################################### # semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=1.4) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 4302329623..ac5c825263 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -4,6 +4,8 @@ using Trixi ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations2D(gravity_constant=9.81) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 8b2070300f..6fede2fa4e 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -5,6 +5,8 @@ using Printf: @printf, @sprintf ############################################################################### # Semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir equations = ShallowWaterEquations2D(gravity_constant=9.812) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl index b9b932ddb2..277f90b025 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl @@ -5,6 +5,9 @@ using Trixi ############################################################################### # semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir + equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=1.875, threshold_limiter=1e-12, threshold_wet=1e-14) diff --git a/src/Trixi.jl b/src/Trixi.jl index e498034e7c..bc24e3c5ae 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -162,6 +162,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, + # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, flux_hll_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, @@ -217,6 +218,7 @@ export DG, VolumeIntegralFluxDifferencing, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, + # TODO: TrixiShallowWater: move new indicator IndicatorHennemannGassnerShallowWater, VolumeIntegralUpwind, SurfaceIntegralWeakForm, SurfaceIntegralStrongForm, diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index 01dd92b2fb..ab0f34efb7 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -6,5 +6,6 @@ #! format: noindent include("positivity_zhang_shu.jl") +# TODO: TrixiShallowWater: move specific limiter file include("positivity_shallow_water.jl") end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index febc027867..36276026fe 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# TODO: TrixiShallowWater: generic wet/dry limiter + """ PositivityPreservingLimiterShallowWater(; variables) diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index ff8be26bef..13c6866e89 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# TODO: TrixiShallowWater: 1D wet/dry limiter should move + function limiter_shallow_water!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations::ShallowWaterEquations1D, diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl index acebf24840..da3a25fdcf 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# TODO: TrixiShallowWater: 2D wet/dry limiter should move + function limiter_shallow_water!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations::ShallowWaterEquations2D, dg::DGSEM, cache) diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index f17299cbf9..8e452c7dc7 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -263,6 +263,8 @@ See [`FluxHLL`](@ref). """ const flux_hll = FluxHLL() +# TODO: TrixiShallowWater: move the chen_noelle flux structure to the new package + # An empty version of the `min_max_speed_chen_noelle` function is declared here # in order to create a dimension agnostic version of `flux_hll_chen_noelle`. # The full description of this wave speed estimate can be found in the docstrings diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index e06f9b2560..8112d18098 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -50,6 +50,8 @@ References for the SWE are many but a good introduction is available in Chapter Finite Volume Methods for Hyperbolic Problems [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) """ +# TODO: TrixiShallowWater: where should the `threshold_limiter` and `threshold_wet` live? +# how to "properly" export these constants across the two packages? struct ShallowWaterEquations1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 3} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height @@ -347,6 +349,7 @@ Further details on the hydrostatic reconstruction and its motivation can be foun A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ +# TODO: TrixiShallowWater: move wet/dry specific routine @inline function flux_nonconservative_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) @@ -495,6 +498,7 @@ Further details on this hydrostatic reconstruction and its motivation can be fou A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ +# TODO: TrixiShallowWater: move wet/dry specific routine @inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations::ShallowWaterEquations1D) # Unpack left and right water heights and bottom topographies @@ -617,6 +621,7 @@ Further details on this hydrostatic reconstruction and its motivation can be fou A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ +# TODO: TrixiShallowWater: move wet/dry specific routine @inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) # Get the velocity quantities @@ -736,6 +741,9 @@ end # Calculate the error for the "lake-at-rest" test case where H = h+b should # be a constant value over time. Note, assumes there is a single reference # water height `H0` with which to compare. +# +# TODO: TrixiShallowWater: where should `threshold_limiter` live? May need +# to modify or have different versions of the `lake_at_rest_error` function @inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D) h, _, b = u diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index eaeb48b999..e326cc54cb 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -53,6 +53,8 @@ References for the SWE are many but a good introduction is available in Chapter Finite Volume Methods for Hyperbolic Problems [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) """ +# TODO: TrixiShallowWater: where should the `threshold_limiter` and `threshold_wet` live? +# how to "properly" export these constants across the two packages? struct ShallowWaterEquations2D{RealT <: Real} <: AbstractShallowWaterEquations{2, 4} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height @@ -469,6 +471,7 @@ Further details on this hydrostatic reconstruction and its motivation can be fou A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ +# TODO: TrixiShallowWater: move wet/dry specific routine @inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations::ShallowWaterEquations2D) # Unpack left and right water heights and bottom topographies @@ -622,6 +625,7 @@ Further details on the hydrostatic reconstruction and its motivation can be foun A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ +# TODO: TrixiShallowWater: move wet/dry specific routine @inline function flux_nonconservative_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) # Pull the water height and bottom topography on the left @@ -962,6 +966,7 @@ the reference below. The definition of the wave speeds are given in Equation (2. A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ +# TODO: TrixiShallowWater: move wet/dry specific routine @inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) h_ll = waterheight(u_ll, equations) @@ -1106,6 +1111,9 @@ end # Calculate the error for the "lake-at-rest" test case where H = h+b should # be a constant value over time. Note, assumes there is a single reference # water height `H0` with which to compare. +# +# TODO: TrixiShallowWater: where should `threshold_limiter` live? May need +# to modify or have different versions of the `lake_at_rest_error` function @inline function lake_at_rest_error(u, equations::ShallowWaterEquations2D) h, _, _, b = u diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index e126eec7c2..4b64481cca 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# TODO: TrixiShallowWater: 1D two layer equations should move to new package + @doc raw""" ShallowWaterTwoLayerEquations1D(gravity, H0, rho_upper, rho_lower) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index a54831c711..87249e9194 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -5,48 +5,50 @@ @muladd begin #! format: noindent +# TODO: TrixiShallowWater: 2D two layer equations should move to new package + @doc raw""" ShallowWaterTwoLayerEquations2D(gravity, H0, rho_upper, rho_lower) Two-Layer Shallow water equations (2LSWE) in two space dimension. The equations are given by ```math \begin{alignat*}{8} -&\frac{\partial}{\partial t}h_{upper} +&\frac{\partial}{\partial t}h_{upper} &&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) -&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}\right) \quad +&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}\right) \quad &&= \quad 0 \\ -&\frac{\partial}{\partial t}\left(h_{upper} v_{1,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}^2 + \frac{gh_{upper}^2}{2}\right) -&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{1,upper} v_{2,upper}\right) \quad +&\frac{\partial}{\partial t}\left(h_{upper} v_{1,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}^2 + \frac{gh_{upper}^2}{2}\right) +&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{1,upper} v_{2,upper}\right) \quad &&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right) \\ -&\frac{\partial}{\partial t}\left(h_{upper} v_{2,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper} v_{2,upper}\right) -&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}^2 + \frac{gh_{upper}^2}{2}\right) +&\frac{\partial}{\partial t}\left(h_{upper} v_{2,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper} v_{2,upper}\right) +&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}^2 + \frac{gh_{upper}^2}{2}\right) &&= -gh_{upper}\frac{\partial}{\partial y}\left(b+h_{lower}\right)\\ -&\frac{\partial}{\partial t}h_{lower} -&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}\right) -&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}\right) +&\frac{\partial}{\partial t}h_{lower} +&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}\right) +&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}\right) &&= \quad 0 \\ -&\frac{\partial}{\partial t}\left(h_{lower} v_{1,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}^2 + \frac{gh_{lower}^2}{2}\right) -&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{1,lower} v_{2,lower}\right) +&\frac{\partial}{\partial t}\left(h_{lower} v_{1,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}^2 + \frac{gh_{lower}^2}{2}\right) +&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{1,lower} v_{2,lower}\right) &&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\frac{\rho_{upper}}{\rho_{lower}} h_{upper}\right)\\ -&\frac{\partial}{\partial t}\left(h_{lower} v_{2,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower} v_{2,lower}\right) -&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}^2 + \frac{gh_{lower}^2}{2}\right) +&\frac{\partial}{\partial t}\left(h_{lower} v_{2,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower} v_{2,lower}\right) +&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}^2 + \frac{gh_{lower}^2}{2}\right) &&= -gh_{lower}\frac{\partial}{\partial y}\left(b+\frac{\rho_{upper}}{\rho_{lower}} h_{upper}\right) \end{alignat*} ``` -The unknown quantities of the 2LSWE are the water heights of the lower layer ``h_{lower}`` and the -upper +The unknown quantities of the 2LSWE are the water heights of the lower layer ``h_{lower}`` and the +upper layer ``h_{upper}`` and the respective velocities in x-direction ``v_{1,lower}`` and ``v_{1,upper}`` and in y-direction -``v_{2,lower}`` and ``v_{2,upper}``. The gravitational constant is denoted by `g`, the layer densitites by -``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable bottom topography function by ``b(x)``. -Conservative variable water height ``h_{lower}`` is measured from the bottom topography ``b`` and ``h_{upper}`` -relative to ``h_{lower}``, therefore one also defines the total water heights as ``H_{lower} = h_{lower} + b`` and +``v_{2,lower}`` and ``v_{2,upper}``. The gravitational constant is denoted by `g`, the layer densitites by +``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable bottom topography function by ``b(x)``. +Conservative variable water height ``h_{lower}`` is measured from the bottom topography ``b`` and ``h_{upper}`` +relative to ``h_{lower}``, therefore one also defines the total water heights as ``H_{lower} = h_{lower} + b`` and ``H_{upper} = h_{upper} + h_{lower} + b``. -The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid +The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid ``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the upper layer. The additional quantity ``H_0`` is also available to store a reference value for the total water @@ -55,13 +57,13 @@ height that is useful to set initial conditions or test the "lake-at-rest" well- The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. -In addition to the unknowns, Trixi currently stores the bottom topography values at the -approximation points despite being fixed in time. This is done for convenience of computing the -bottom topography gradients on the fly during the approximation as well as computing auxiliary +In addition to the unknowns, Trixi currently stores the bottom topography values at the +approximation points despite being fixed in time. This is done for convenience of computing the +bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H`` or the entropy variables. This affects the implementation and use of these equations in various ways: * The flux values corresponding to the bottom topography must be zero. -* The bottom topography values must be included when defining initial conditions, boundary +* The bottom topography values must be included when defining initial conditions, boundary conditions or source terms. * [`AnalysisCallback`](@ref) analyzes this variable. * Trixi's visualization tools will visualize the bottom topography by default. @@ -113,7 +115,7 @@ end initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations2D) A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref). Constants must be set to ``rho_{upper} = 0.9``, +[`source_terms_convergence_test`](@ref). Constants must be set to ``rho_{upper} = 0.9``, ``rho_{lower} = 1.0``, ``g = 10.0``. """ function initial_condition_convergence_test(x, t, @@ -141,7 +143,7 @@ Source terms used for convergence tests in combination with """ @inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations2D) - # Same settings as in `initial_condition_convergence_test`. + # Same settings as in `initial_condition_convergence_test`. # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2] ω = 2.0 * pi * sqrt(2.0) @@ -325,7 +327,7 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version +additional term that couples the momentum of both layers. This is a slightly modified version to account for the additional source term compared to the standard SWE described in the paper. Further details are available in the paper: @@ -345,7 +347,7 @@ Further details are available in the paper: z = zero(eltype(u_ll)) # Bottom gradient nonconservative term: (0, g*h_upper*(b + h_lower)_x, g*h_upper*(b + h_lower)_y , - # 0, g*h_lower*(b + r*h_upper)_x, + # 0, g*h_lower*(b + r*h_upper)_x, # g*h_lower*(b + r*h_upper)_y, 0) if orientation == 1 f = SVector(z, @@ -397,8 +399,8 @@ end !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both +Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains +the gradients of the bottom topography and an additional term that couples the momentum of both layers [`ShallowWaterTwoLayerEquations2D`](@ref). Further details are available in the paper: @@ -506,13 +508,13 @@ end flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations2D) -Total energy conservative (mathematical entropy for two-layer shallow water equations). When the -bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will +Total energy conservative (mathematical entropy for two-layer shallow water equations). When the +bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). Details are available in Eq. (4.1) in the paper: - Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) - Well-balanced and energy stable schemes for the shallow water equations with discontinuous + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) and the application to two layers is shown in the paper: - Ulrik Skre Fjordholm (2012) @@ -606,11 +608,11 @@ end """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations2D) - + Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -696,9 +698,9 @@ end flux_es_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterTwoLayerEquations1D) -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative [`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy -variables. +variables. Further details are available in the paper: - Ulrik Skre Fjordholm (2012) @@ -723,7 +725,7 @@ formulation. q_rr = cons2entropy(u_rr, equations) q_ll = cons2entropy(u_ll, equations) - # Average values from left and right + # Average values from left and right u_avg = (u_ll + u_rr) / 2 # Introduce variables for better readability @@ -791,10 +793,10 @@ formulation. end # Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the -# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate -# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them -# analytically. -# +# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate +# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them +# analytically. +# # A good overview of the derivation is given in: # - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) # Open boundary conditions for nonlinear channel Flows @@ -914,7 +916,7 @@ end # Convert conservative variables to entropy variables # Note, only the first four are the entropy variables, the fifth entry still just carries the bottom -# topography values for convenience. +# topography values for convenience. # In contrast to general usage the entropy variables are denoted with q instead of w, because w is # already used for velocity in y-Direction @inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations2D) diff --git a/src/equations/untitled_upwind.jl b/src/equations/untitled_upwind.jl new file mode 100644 index 0000000000..ecfcf62a74 --- /dev/null +++ b/src/equations/untitled_upwind.jl @@ -0,0 +1,115 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +# Free-stream initial condition +initial_condition = initial_condition_constant + +# Boundary conditions for free-stream testing +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) + +# boundary_conditions = Dict( :Body => boundary_condition_free_stream, +# :Button1 => boundary_condition_free_stream, +# :Button2 => boundary_condition_free_stream, +# :Eye1 => boundary_condition_free_stream, +# :Eye2 => boundary_condition_free_stream, +# :Smile => boundary_condition_free_stream, +# :Bowtie => boundary_condition_free_stream ) + +# for a warped style mesh +# boundary_conditions = Dict( :Top => boundary_condition_free_stream, +# :Bottom => boundary_condition_free_stream, +# :Right => boundary_condition_free_stream, +# :Left => boundary_condition_free_stream ) + +# for a warped style mesh +boundary_conditions = Dict(:Top => boundary_condition_free_stream, + :Bottom => boundary_condition_free_stream, + :Right => boundary_condition_free_stream, + :Slant => boundary_condition_free_stream, + :Bezier => boundary_condition_free_stream) + +# for the outer circle mesh +# boundary_conditions = Dict(:OuterCircle => boundary_condition_free_stream) + +# for 3 circles mesh +# boundary_conditions = Dict(:outer => boundary_condition_free_stream, +# :Circle1 => boundary_condition_free_stream, +# :Circle2 => boundary_condition_free_stream, +# :Circle3 => boundary_condition_free_stream) + +############################################################################### +# Get the Upwind FDSBP approximation space + +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order = 1, + accuracy_order = 4, + xmin = -1.0, xmax = 1.0, + N = 18) + +flux_splitting = splitting_lax_friedrichs +solver = FDSBP(D_upw, + surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)), + # volume_integral=VolumeIntegralStrongForm()) + volume_integral = VolumeIntegralUpwind(flux_splitting)) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "mesh_multiple_flips.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh", + default_mesh_file) + +# MADE A NEW FLIPPED ORIENTAION MESH +#default_mesh_file = joinpath(@__DIR__, "mesh_new_flip.mesh") +# default_mesh_file = joinpath(@__DIR__, "mesh_new_flip_corners.mesh") +# default_mesh_file = joinpath(@__DIR__, "WarpedMesh_00.mesh") + +# Circular mesh with flips +default_mesh_file = joinpath(@__DIR__, "mesh_docs.mesh") +# default_mesh_file = joinpath(@__DIR__, "mesh_outer_circle.mesh") +# default_mesh_file = joinpath(@__DIR__, "Circles3.mesh") + +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, + alive_callback, save_solution) + +############################################################################### +# run the simulation + +# set small tolerances for the free-stream preservation test +sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12, + save_everystep = false, callback = callbacks) +summary_callback() # print the timer summary diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index f62babe38d..d076181b5a 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -123,6 +123,8 @@ See also [`VolumeIntegralShockCapturingHG`](@ref). "A provably entropy stable subcell shock capturing approach for high order split form DG" [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) """ +# +# TODO: TrixiShallowWater: move the new indicator and all associated routines to the new package struct IndicatorHennemannGassnerShallowWater{RealT <: Real, Variable, Cache} <: AbstractIndicator alpha_max::RealT diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 04252b68bd..fd87a87338 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -27,6 +27,8 @@ end # Modified indicator for ShallowWaterEquations1D to apply full FV method on cells # containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a # full FV element. +# +# TODO: TrixiShallowWater: move new indicator type function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, 3 }, mesh, diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index d97a5eb516..fec42c2621 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -31,6 +31,8 @@ end # Modified indicator for ShallowWaterEquations2D to apply full FV method on cells # containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a # full FV element. +# +# TODO: TrixiShallowWater: move new indicator type function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, 4 }, mesh, diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index a8a0b41d12..75937ba82a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -1,5 +1,7 @@ module TestExamplesStructuredMesh2D +# TODO: TrixiShallowWater: move any wet/dry tests to new package + using Test using Trixi diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 573d02f363..cafa17edd4 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -1,5 +1,7 @@ module TestExamples1DShallowWater +# TODO: TrixiShallowWater: move any wet/dry tests to new package + using Test using Trixi diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index 0d8a83806f..8372d0d467 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -1,5 +1,7 @@ module TestExamples1DShallowWaterTwoLayer +# TODO: TrixiShallowWater: move two layer tests to new package + using Test using Trixi diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 4c7ada723a..7670d28f43 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -1,5 +1,7 @@ module TestExamples2DShallowWater +# TODO: TrixiShallowWater: move any wet/dry tests to new package + using Test using Trixi diff --git a/test/test_tree_2d_shallowwater_twolayer.jl b/test/test_tree_2d_shallowwater_twolayer.jl index 4bb4506471..7ad5b0f731 100644 --- a/test/test_tree_2d_shallowwater_twolayer.jl +++ b/test/test_tree_2d_shallowwater_twolayer.jl @@ -1,5 +1,7 @@ module TestExamples2DShallowWaterTwoLayer +# TODO: TrixiShallowWater: move two layer tests to new package + using Test using Trixi @@ -19,10 +21,10 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @trixi_testset "elixir_shallowwater_twolayer_convergence.jl with flux_es_fjordholm_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_convergence.jl"), - l2 = [0.00024709443131137236, 0.0019215286339769443, 0.0023833298173254447, + l2 = [0.00024709443131137236, 0.0019215286339769443, 0.0023833298173254447, 0.00021258247976270914, 0.0011299428031136195, 0.0009191313765262401, - 8.873630921431545e-6], - linf = [0.0016099763244645793, 0.007659242165565017, 0.009123320235427057, + 8.873630921431545e-6], + linf = [0.0016099763244645793, 0.007659242165565017, 0.009123320235427057, 0.0013496983982568267, 0.0035573687287770994, 0.00296823235874899, 3.361991620143279e-5], surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal), @@ -31,19 +33,19 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_well_balanced.jl"), - l2 = [3.2935164267930016e-16, 4.6800825611195103e-17, 4.843057532147818e-17, - 0.0030769233188015013, 1.4809161150389857e-16, 1.509071695038043e-16, + l2 = [3.2935164267930016e-16, 4.6800825611195103e-17, 4.843057532147818e-17, + 0.0030769233188015013, 1.4809161150389857e-16, 1.509071695038043e-16, 0.0030769233188014935], - linf = [2.248201624865942e-15, 2.346382070278936e-16, 2.208565017494899e-16, - 0.026474051138910493, 9.237568031609006e-16, 7.520758026187046e-16, + linf = [2.248201624865942e-15, 2.346382070278936e-16, 2.208565017494899e-16, + 0.026474051138910493, 9.237568031609006e-16, 7.520758026187046e-16, 0.026474051138910267], tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_twolayer_well_balanced with flux_lax_friedrichs.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_well_balanced.jl"), - l2 = [2.0525741072929735e-16, 6.000589392730905e-17, 6.102759428478984e-17, - 0.0030769233188014905, 1.8421386173122792e-16, 1.8473184927121752e-16, + l2 = [2.0525741072929735e-16, 6.000589392730905e-17, 6.102759428478984e-17, + 0.0030769233188014905, 1.8421386173122792e-16, 1.8473184927121752e-16, 0.0030769233188014935], linf = [7.355227538141662e-16, 2.960836949170518e-16, 4.2726562436938764e-16, 0.02647405113891016, 1.038795478061861e-15, 1.0401789378532516e-15, diff --git a/test/test_unit.jl b/test/test_unit.jl index 6392da94fd..2012f5d289 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -402,6 +402,7 @@ isdir(outdir) && rm(outdir, recursive=true) indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) + # TODO: TrixiShallowWater: move unit test indicator_hg_swe = IndicatorHennemannGassnerShallowWater(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg_swe) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 9c48d8ba7a..fbe88a2a0a 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -1,5 +1,7 @@ module TestExamplesUnstructuredMesh2D +# TODO: TrixiShallowWater: move any wet/dry and two layer tests + using Test using Trixi