diff --git a/NEWS.md b/NEWS.md index ba3b6c7e..d28b6e90 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,9 @@ for human readability. #### Added - Experimental support for well-balanced mortars together with AMR [#45] +- New boundary conditions `BoundaryConditionWaterHeight` and `BoundaryConditionMomentum` now + available for the `ShallowWaterEquationsWetDry` to impose either the water height or the momentum + at the boundary. [#91] #### Changed diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_shock.jl b/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_shock.jl new file mode 100644 index 00000000..25cbb7b0 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_shock.jl @@ -0,0 +1,184 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater +using Roots + +############################################################################### +# semidiscretization of the shallow water equations for a transonic moving water +# steady-state with a standing shock. + +equations = ShallowWaterEquationsWetDry1D(gravity = 9.812, H0 = 3.25) + +""" + inverse_transform(E, hv, sigma, b) + +Inverse transformation from equilibrium variables (E, hv) to conservative variables (h, hv). Besides the +equilibrium variables, which are the total energy `E` and momentum `hv`, the function also depends +on the bottom topography `b` and the flow regime `sigma` (supersonic = 1 , sonic = 0 or subsonic = -1). + +The implementation follows the procedure described in Section 2.1 of the paper: +- Sebastian Noelle, Yulong Xing and Chi-Wang Shu (2007) + High Order Well-balanced Finite Volume WENO Schemes for Shallow Water Equation with Moving Water + [DOI: 10.1016/j.jcp.2007.03.031](https://doi.org/10.1016/j.jcp.2007.03.031). +""" +function inverse_transform(E, hv, sigma, b) + # Extract the gravitational acceleration + g = equations.gravity + + # Compute water height and specific energy at the sonic point + h_0 = 1 / g * (g * abs(hv))^(2 / 3) + phi_0 = 3 / 2 * (g * abs(hv))^(2 / 3) + + # normalized total energy + E_hat = (E - g * b) / phi_0 + + # Check if the state is admissible and compute the water height + # as in equation (2.19) of the reference in the docstring. + if sigma == 0 && E_hat ≈ 1 # sonic state + h_hat = 1 + elseif abs(sigma) == 1 && E_hat > 1 # supersonic / subsonic state + if sigma == 1 # supersonic + h_hat_init = 0.5 # needs to be < 1 + else # subsonic + h_hat_init = 2.0 # needs to be > 1 + end + + # Setup the root finding problem. + f(h_hat) = E_hat - 2 / 3 * ((1 / (2 * h_hat^2)) + h_hat) + D(f) = h_hat -> Trixi.ForwardDiff.derivative(f, float(h_hat)) + + # Solve the root finding problem using Newton's method + h_hat = Roots.newton((f, D(f)), h_hat_init) + else + throw(error("The given state is not admissible: E_hat = $E_hat, sigma = $sigma")) + end + + # Compute and return the water height `h` from the normalized water height `h_hat = h / h_0`. + return h_hat * h_0 +end + +""" + initial_condition_moving_water_transonic(x, t, equations::ShallowWaterEquationsWetDry1D) + +Set the initial condition for a transonic moving water steady-state and a quadratic bottom +topography, to test the well-balancedness of the scheme. + +The test parameters are taken from Section 4.1 of the paper: +- Sebastian Noelle, Yulong Xing and Chi-Wang Shu (2007) + High Order Well-balanced Finite Volume WENO Schemes for Shallow Water Equation with Moving Water + [DOI: 10.1016/j.jcp.2007.03.031](https://doi.org/10.1016/j.jcp.2007.03.031). +""" +function initial_condition_moving_water_shock(x, t, + equations::ShallowWaterEquationsWetDry1D) + # Extract the gravitational acceleration + g = equations.gravity + + hv = 0.18 # momentum + + # Set the total energy before and after the shock + if x[1] < 11.665504281554291 + E = 3 / 2 * (g * 0.18)^(2 / 3) + g * 0.2 + else + E = 0.18^2 / (2 * 0.33^2) + g * 0.33 + end + + # # Set the quadratic bottom topography function + if 8 <= x[1] <= 12 + b = 0.2 - 0.05 * (x[1] - 10.0)^2 + else + b = 0.0 + end + + # Set the sign function to label the flow regime (subsonic = -1, sonic = 0, supersonic = 1). + # A small tolerance is required to avoid numerical issues in the inverse_transform function + # close to the sonic point at x = 10. + tol = 1e-12 + if x[1] <= 10.0 - tol || x[1] >= 11.665504281554291 + sigma = -1 # subsonic + elseif 10 - tol < x[1] < 10.0 + tol + sigma = 0 # sonic + elseif 10 + tol <= x[1] < 11.665504281554291 + sigma = 1 # supersonic + end + + h = inverse_transform(E, hv, sigma, b) + + return SVector(h, hv, b) +end + +initial_condition = initial_condition_moving_water_shock + +boundary_condition_inflow = BoundaryConditionMomentum(0.18, equations) +boundary_condition_outflow = BoundaryConditionDirichlet(initial_condition_moving_water_shock) + +boundary_conditions = (x_neg = boundary_condition_inflow, + x_pos = boundary_condition_outflow) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxPlusDissipation(flux_wintermeyer_etal, DissipationLocalLaxFriedrichs()), + flux_nonconservative_wintermeyer_etal) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = 0.0 +coordinates_max = 32.0 # This needs to be a multiple of 2 to match the corners of the bottom topography +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000, + periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# 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, + save_analysis = true) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.7) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_subsonic.jl b/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_subsonic.jl new file mode 100644 index 00000000..e37cac4d --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_subsonic.jl @@ -0,0 +1,155 @@ +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater +using Roots + +############################################################################### +# semidiscretization of the shallow water equations for a subsonic moving water steady-state + +equations = ShallowWaterEquationsWetDry1D(gravity = 9.812, H0 = 3.25) + +""" + inverse_transform(E, hv, sigma, b) + +Inverse transformation from equilibrium variables (E, hv) to conservative variables (h, hv). Besides the +equilibrium variables, which are the total energy `E` and momentum `hv`, the function also depends +on the bottom topography `b` and the flow regime `sigma` (supersonic = 1 , sonic = 0 or subsonic = -1). + +The implementation follows the procedure described in Section 2.1 of the paper: +- Sebastian Noelle, Yulong Xing and Chi-Wang Shu (2007) + High Order Well-balanced Finite Volume WENO Schemes for Shallow Water Equation with Moving Water + [DOI: 10.1016/j.jcp.2007.03.031](https://doi.org/10.1016/j.jcp.2007.03.031). +""" +function inverse_transform(E, hv, sigma, b) + # Extract the gravitational acceleration + g = equations.gravity + + # Compute water height and specific energy at the sonic point + h_0 = 1 / g * (g * abs(hv))^(2 / 3) + phi_0 = 3 / 2 * (g * abs(hv))^(2 / 3) + + # normalized total energy + E_hat = (E - g * b) / phi_0 + + # Check if the state is admissible and compute the water height + # as in equation (2.19) of the reference in the docstring. + if sigma == 0 && E_hat ≈ 1 # sonic state + h_hat = 1 + elseif abs(sigma) == 1 && E_hat > 1 # supersonic / subsonic state + # Pick an initial guess for the root finding problem based on the flow regime + if sigma == 1 # supersonic + h_hat_init = 0.5 # needs to be < 1 + else # subsonic + h_hat_init = 2 # needs to be > 1 + end + + # Setup the root finding problem. + f(h_hat) = E_hat - 2 / 3 * ((1 / (2 * h_hat^2)) + h_hat) + D(f) = h_hat -> Trixi.ForwardDiff.derivative(f, float(h_hat)) + + # Solve the root finding problem using Newton's method + h_hat = Roots.newton((f, D(f)), h_hat_init) + else + throw(error("The given state is not admissible: E_hat = $E_hat, sigma = $sigma")) + end + + # Compute and return the water height `h` from the normalized water height `h_hat = h / h_0`. + return h_hat * h_0 +end + +""" + initial_condition_moving_water_subsonic(x, t, equations::ShallowWaterEquations1D) + +Set the initial condition for a subsonic moving water steady-state and quadratic bottom topography, +to test the well-balancedness of the scheme. + +The test setup is taken from Section 4.1 of the paper: +- Sebastian Noelle, Yulong Xing and Chi-Wang Shu (2007) + High Order Well-balanced Finite Volume WENO Schemes for Shallow Water Equation with Moving Water + [DOI: 10.1016/j.jcp.2007.03.031](https://doi.org/10.1016/j.jcp.2007.03.031). +""" +function initial_condition_moving_water_subsonic(x, t, + equations::ShallowWaterEquationsWetDry1D) + # Set initial conditions + hv = 4.42 # momentum + E = 22.06605 # total energy + + # Set the quadratic bottom topography function + if 8 <= x[1] <= 12 + b = 0.2 - 0.05 * (x[1] - 10.0)^2 + else + b = 0.0 + end + + sigma = -1 # sign function to label the flow regime (subsonic = -1, sonic = 0, supersonic = 1) + + # Compute the water height using the inverse transformation + h = inverse_transform(E, hv, sigma, b) + + return SVector(h, hv, b) +end + +initial_condition = initial_condition_moving_water_subsonic + +boundary_condition_inflow = BoundaryConditionMomentum(4.42, equations) +boundary_condition_outflow = BoundaryConditionWaterHeight(2.0, equations) + +boundary_conditions = (x_neg = boundary_condition_inflow, + x_pos = boundary_condition_outflow) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) + +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = 0.0 +coordinates_max = 32.0 # This needs to be a multiple of 2 to match the corners of the bottom topography +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000, + periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# 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, + save_analysis = true) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_transonic.jl b/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_transonic.jl new file mode 100644 index 00000000..17553f75 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_moving_water_transonic.jl @@ -0,0 +1,166 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater +using Roots + +############################################################################### +# semidiscretization of the shallow water equations for a transonic moving water steady-state + +equations = ShallowWaterEquationsWetDry1D(gravity = 9.812, H0 = 3.25) + +""" + inverse_transform(E, hv, sigma, b) + +Inverse transformation from equilibrium variables (E, hv) to conservative variables (h, hv). Besides the +equilibrium variables, which are the total energy `E` and momentum `hv`, the function also depends +on the bottom topography `b` and the flow regime `sigma` (supersonic = 1 , sonic = 0 or subsonic = -1). + +The implementation follows the procedure described in Section 2.1 of the paper: +- Sebastian Noelle, Yulong Xing and Chi-Wang Shu (2007) + High Order Well-balanced Finite Volume WENO Schemes for Shallow Water Equation with Moving Water + [DOI: 10.1016/j.jcp.2007.03.031](https://doi.org/10.1016/j.jcp.2007.03.031). +""" +function inverse_transform(E, hv, sigma, b) + # Extract the gravitational acceleration + g = equations.gravity + + # Compute water height and specific energy at the sonic point + h_0 = 1 / g * (g * abs(hv))^(2 / 3) + phi_0 = 3 / 2 * (g * abs(hv))^(2 / 3) + + # normalized total energy + E_hat = (E - g * b) / phi_0 + + # Check if the state is admissible and compute the water height + # as in equation (2.19) of the reference in the docstring. + if sigma == 0 && E_hat ≈ 1 # sonic state + h_hat = 1 + elseif abs(sigma) == 1 && E_hat > 1 # supersonic / subsonic state + # Pick an initial guess for the root finding problem based on the flow regime + if sigma == 1 # supersonic + h_hat_init = 0.5 # needs to be < 1 + else # subsonic + h_hat_init = 2.0 # needs to be > 1 + end + + # Setup the root finding problem. + f(h_hat) = E_hat - 2 / 3 * ((1 / (2 * h_hat^2)) + h_hat) + D(f) = h_hat -> Trixi.ForwardDiff.derivative(f, float(h_hat)) + + # Solve the root finding problem using Newton's method + h_hat = Roots.newton((f, D(f)), h_hat_init) + else + throw(error("The given state is not admissible: E_hat = $E_hat, sigma = $sigma")) + end + + # Compute and return the water height `h` from the normalized water height `h_hat = h / h_0`. + return h_hat * h_0 +end + +""" + initial_condition_moving_water_transonic(x, t, equations::ShallowWaterEquationsWetDry1D) + +Set the initial condition for a smooth transonic moving water steady-state and a quadratic bottom +topography, to test the well-balancedness of the scheme. + +The test parameters are taken from Section 5 of the paper: +- Christian Klingenberg, Alexander Kurganov, Yongle Liu, and Markus Zenk + Moving-Water Equilibria Preserving HLL-Type Schemes for the Shallow Water Equations + [DOI: 10.4208/cmr.2020-0013](https://doi.org/10.4208/cmr.2020-0013) +""" +function initial_condition_moving_water_transonic(x, t, + equations::ShallowWaterEquationsWetDry1D) + # Set initial conditions + hv = 1.53 # momentum + E = 11.090714039778195 # total energy + + # Set the quadratic bottom topography function + if 8 <= x[1] <= 12 + b = 0.2 - 0.05 * (x[1] - 10.0)^2 + else + b = 0.0 + end + + # Set the sign function to label the flow regime (subsonic = -1, sonic = 0, supersonic = 1). + # A small tolerance is required to avoid numerical issues in the inverse_transform function + # close to the sonic point at x = 10. + tol = 1e-10 + if x[1] <= 10.0 - tol + sigma = -1 # subsonic + elseif 10 - tol < x[1] < 10 + tol + sigma = 0 # sonic + else + sigma = 1 # supersonic + end + + # Compute the water height using the inverse transformation + h = inverse_transform(E, hv, sigma, b) + + return SVector(h, hv, b) +end + +initial_condition = initial_condition_moving_water_transonic + +boundary_condition_inflow = BoundaryConditionMomentum(1.53, equations) +boundary_condition_outflow = BoundaryConditionDirichlet(initial_condition_moving_water_transonic) + +boundary_conditions = (x_neg = boundary_condition_inflow, + x_pos = boundary_condition_outflow) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) + +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = 0.0 +coordinates_max = 32.0 # This needs to be a multiple of 2 to match the corners of the bottom topography +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000, + periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# 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, + save_analysis = true) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_inflow_outflow.jl b/examples/tree_2d_dgsem/elixir_shallowwater_inflow_outflow.jl new file mode 100644 index 00000000..c7f104f1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_inflow_outflow.jl @@ -0,0 +1,94 @@ +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations to test inflow/outflow boundary conditions + +equations = ShallowWaterEquationsWetDry2D(gravity = 9.81) + +# Setup initial conditions for a smooth channel flow with constant water height and velocity +function initial_condition_channel_flow(x, t, equations::ShallowWaterEquationsWetDry2D) + H = 1.0 + v1 = -0.1 + v2 = -0.1 + b = 0.0 + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_channel_flow + +# Setup boundary conditions. +# At the inlet, we prescribe the momentum, starting with a negative value to simulate outflow, +# and gradually transitioning to a positive value to simulate inflow. At the outlet, we prescribe +# the water height as a time-dependent cosine wave. This setup is designed to test the behavior +# of the boundary conditions under both inflow and outflow scenarios. +boundary_condition_inflow = BoundaryConditionMomentum(t -> -0.1 + 0.05 * t, + t -> -0.1 + 0.05 * t, + equations) +boundary_condition_outflow = BoundaryConditionWaterHeight(t -> 1.0 + 0.1 * cos(π / 2 * t), + equations) +boundary_conditions = (x_neg = boundary_condition_inflow, + x_pos = boundary_condition_outflow, + y_neg = boundary_condition_inflow, + y_pos = boundary_condition_outflow) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxPlusDissipation(flux_wintermeyer_etal, DissipationLocalLaxFriedrichs()), + flux_nonconservative_wintermeyer_etal) +polydeg = 3 +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Create a non-periodic TreeMesh +coordinates_min = (-10.0, -10.0) +coordinates_max = (10.0, 10.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# 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, + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 1.0, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_inflow_outflow.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_inflow_outflow.jl new file mode 100644 index 00000000..dec0b8d2 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_inflow_outflow.jl @@ -0,0 +1,94 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations to test inflow/outflow boundary conditions + +equations = ShallowWaterEquationsWetDry2D(gravity = 9.81) + +# Setup initial conditions for a smooth channel flow with constant water height and velocity +function initial_condition_channel_flow(x, t, equations::ShallowWaterEquationsWetDry2D) + H = 1.0 + v1 = 0.4 + v2 = 0.4 + b = 0.0 + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_channel_flow + +# Setup boundary conditions. +# At the inlet, we prescribe constant momentum to simulate inflow. +# At the outlet, we prescribe the water height as a time-dependent cosine wave. +boundary_condition_inflow = BoundaryConditionMomentum(0.4, 0.4, equations) +boundary_condition_outflow = BoundaryConditionWaterHeight(t -> 1.0 + 0.1 * cos(π / 2 * t), + equations) + +boundary_conditions = Dict(:Bottom => boundary_condition_inflow, + :Top => boundary_condition_outflow, + :Left => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxPlusDissipation(flux_wintermeyer_etal, DissipationLocalLaxFriedrichs()), + flux_nonconservative_wintermeyer_etal) +polydeg = 3 +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "mesh_wobbly_channel.mesh") +isfile(default_mesh_file) || + Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/431baa423ce86aadba70d38a3194947b/raw/50914cb30e72e9a58d4723e161476435c6dea182/mesh_wobbly_channel.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# 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, + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.4, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index fae82ea3..6ac3ff8f 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -6,11 +6,12 @@ module TrixiShallowWater # https://github.com/trixi-framework/TrixiShallowWater.jl/pull/10#discussion_r1433720559 using Trixi # Import additional symbols that are not exported by Trixi.jl -using Trixi: get_node_vars, set_node_vars!, waterheight +using Trixi: get_node_vars, set_node_vars! using MuladdMacro: @muladd using StaticArrays: SVector, @SMatrix, MVector using Static: True, False using LinearAlgebra: norm +using Roots: Order2, solve, ZeroProblem include("equations/equations.jl") include("equations/numerical_fluxes.jl") @@ -33,6 +34,8 @@ export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, export ManningFriction, MeyerPeterMueller, GrassModel, ShieldsStressModel, dissipation_roe, water_sediment_height, source_term_bottom_friction +export BoundaryConditionWaterHeight, BoundaryConditionMomentum + export nlayers, eachlayer export PositivityPreservingLimiterShallowWater diff --git a/src/equations/equations.jl b/src/equations/equations.jl index b215b9ae..8e4540c6 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -5,6 +5,16 @@ @muladd begin #! format: noindent +# Struct used for multiple dispatch on boundary conditions that set the water height at the boundary. +struct BoundaryConditionWaterHeight{F <: Function} + h_boundary::F +end + +# Struct used for multiple dispatch on boundary conditions that set the momentum at the boundary. +struct BoundaryConditionMomentum{F <: Function} + hv_boundary::F +end + #################################################################################################### # Include files with actual implementations for different systems of equations. diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index 25fde5d4..0b3db28b 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -217,6 +217,159 @@ For details see Section 9.2.5 of the book: return flux, noncons_flux end +""" + BoundaryConditionWaterHeight(h_boundary, equations::ShallowWaterEquationsWetDry1D) + +Create a boundary condition that uses `h_boundary` to specify a fixed water height at the +boundary and extrapolates the velocity from the incoming Riemann invariant. + +The external water height `h_boundary` can be specified as a constant value or as a function of time, e.g. +```julia + BoundaryConditionWaterHeight(h_boundary, equations)) + BoundaryConditionWaterHeight(t -> h_boundary(t), equations)) +``` + +More details can be found in the paper: +- Lixiang Song, Jianzhong Zhou, Jun Guo, Qiang Zou, Yi Liu (2011) + A robust well-balanced finite volume model for shallow water flows + with wetting and drying over irregular terrain + [doi: 10.1016/j.advwatres.2011.04.017](https://doi.org/10.1016/j.advwatres.2011.04.017) + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +function BoundaryConditionWaterHeight(h_boundary::Real, + equations::ShallowWaterEquationsWetDry1D{RealT}) where {RealT} + # Convert function output to the correct type + h_boundary = convert(RealT, h_boundary) + return BoundaryConditionWaterHeight(t -> h_boundary) +end + +function BoundaryConditionWaterHeight(h_boundary::Function, + equations::ShallowWaterEquationsWetDry1D{RealT}) where {RealT} + # Check if the function output is of the correct type + if !(typeof(h_boundary(one(RealT))) == RealT) + throw(ArgumentError("Boundary value functions must return a value of type $(RealT)")) + end + return BoundaryConditionWaterHeight(t -> h_boundary(t)) +end + +# Version for `TreeMesh` +function (boundary_condition::BoundaryConditionWaterHeight)(u_inner, + orientation, + direction, x, t, + surface_flux_functions, + equations::ShallowWaterEquationsWetDry1D) + # Extract the gravitational acceleration + g = equations.gravity + + # Get the water height and velocity from the inner state + h_inner = waterheight(u_inner, equations) + v_inner = velocity(u_inner, equations) + + # Extract the external water height from the boundary condition + h_boundary = boundary_condition.h_boundary(t) + + # Calculate the boundary state based on the direction. + # To extrapolate the external velocity assume that the Riemann invariant remains constant across + # the incoming characteristic. In the case of inflow we assume that the tangential velocity at + # the boundary is zero. + if direction == 1 # x- + v_boundary = v_inner + 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + u_boundary = SVector(h_boundary, h_boundary * v_boundary, u_inner[3]) + elseif direction == 2 # x+ + v_boundary = v_inner - 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + u_boundary = SVector(h_boundary, h_boundary * v_boundary, u_inner[3]) + end + + # Evaluate the conservative flux at the boundary + flux = Trixi.flux(u_boundary, orientation, equations) + + # Return the conservative and nonconservative fluxes. + # The nonconservative part is zero as we assume a constant bottom topography at the boundary. + return (flux, zero(u_inner)) +end + +""" + BoundaryConditionMomentum(hv_boundary, equations::ShallowWaterEquationsWetDry1D) + +Create a boundary condition that sets a fixed momentum `hv_boundary` at the boundary and +extrapolates the water height `h_boundary` from the incoming Riemann invariant. + +The external momentum can be specified as a constant value or as a function of time, e.g. +```julia + BoundaryConditionMomentum(hv_boundary, equations) + BoundaryConditionMomentum(t -> hv_boundary(t), equations) +``` + +More details can be found in the paper: +- Lixiang Song, Jianzhong Zhou, Jun Guo, Qiang Zou, Yi Liu (2011) + A robust well-balanced finite volume model for shallow water flows + with wetting and drying over irregular terrain + [doi: 10.1016/j.advwatres.2011.04.017](https://doi.org/10.1016/j.advwatres.2011.04.017) + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +function BoundaryConditionMomentum(hv_boundary::Real, + equations::ShallowWaterEquationsWetDry1D{RealT}) where {RealT} + # Convert function output to the correct type + hv_boundary = convert(RealT, hv_boundary) + return BoundaryConditionMomentum(t -> hv_boundary) +end + +function BoundaryConditionMomentum(hv_boundary::Function, + equations::ShallowWaterEquationsWetDry1D{RealT}) where {RealT} + # Check if the function output is of the correct type + if !(typeof(hv_boundary(one(RealT))) == RealT) + throw(ArgumentError("Boundary value functions must return a value of type $(RealT)")) + end + return BoundaryConditionMomentum(t -> hv_boundary(t)) +end + +# Version for `TreeMesh` +function (boundary_condition::BoundaryConditionMomentum)(u_inner, + orientation, + direction, x, t, + surface_flux_functions, + equations::ShallowWaterEquationsWetDry1D) + # Extract the gravitational acceleration + g = equations.gravity + + # Get the water height and velocity from the inner state + h_inner = waterheight(u_inner, equations) + v_inner = velocity(u_inner, equations) + + # Extract the external momentum from the boundary condition + hv_boundary = boundary_condition.hv_boundary(t) + + # Calculate the boundary state based on the direction. + # To extrapolate the external water height `h_boundary` assume that the Riemann invariant remains + # constant across the incoming characteristic. + # Requires one to solve for the roots of a nonlinear function, see Eq. (52) in the reference above. + # For convenience we substitute x = h_boundary and solve for x, using the Steffensen method. + if direction == 1 # x- + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) + + (v_inner - 2 * sqrt(g * h_inner)) * x - hv_boundary, + h_inner) + h_boundary = solve(fx, Order2()) + u_boundary = SVector(h_boundary, hv_boundary, u_inner[3]) + elseif direction == 2 # x+ + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) - + (v_inner + 2 * sqrt(g * h_inner)) * x + hv_boundary, + h_inner) + h_boundary = solve(fx, Order2()) + u_boundary = SVector(h_boundary, hv_boundary, u_inner[3]) + end + + # Evaluate the conservative flux at the boundary + flux = Trixi.flux(u_boundary, orientation, equations) + + # Return the conservative and nonconservative fluxes. + # The nonconservative part is zero as we assume a constant bottom topography at the boundary. + return (flux, zero(u_inner)) +end + # Calculate 1D flux for a single point # Note, the bottom topography has no flux @inline function Trixi.flux(u, orientation::Integer, diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index 523baceb..00bdecbf 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -4,7 +4,6 @@ # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin #! format: noindent - @doc raw""" ShallowWaterEquationsWetDry2D(; gravity, H0 = 0, threshold_limiter = nothing, threshold_wet = nothing) @@ -243,6 +242,304 @@ Should be used together with [`Trixi.TreeMesh`](@extref). return flux, noncons_flux end +""" + BoundaryConditionWaterHeight(h_boundary, equations::ShallowWaterEquationsWetDry2D) + +Create a boundary condition that uses `h_boundary` to specify a fixed water height at the +boundary and extrapolates the velocity from the incoming Riemann invariant. + +The external water height `h_boundary` can be specified as a constant value or as a function of time, e.g. +```julia + BoundaryConditionWaterHeight(h_boundary, equations)) + BoundaryConditionWaterHeight(t -> h_boundary(t), equations)) +``` + +More details can be found in the paper: +- Lixiang Song, Jianzhong Zhou, Jun Guo, Qiang Zou, Yi Liu (2011) + A robust well-balanced finite volume model for shallow water flows + with wetting and drying over irregular terrain + [doi: 10.1016/j.advwatres.2011.04.017](https://doi.org/10.1016/j.advwatres.2011.04.017) + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +function BoundaryConditionWaterHeight(h_boundary::Real, + equations::ShallowWaterEquationsWetDry2D{RealT}) where {RealT} + # Convert function output to the correct type + h_boundary = convert(RealT, h_boundary) + return BoundaryConditionWaterHeight(t -> h_boundary) +end + +function BoundaryConditionWaterHeight(h_boundary::Function, + equations::ShallowWaterEquationsWetDry2D{RealT}) where {RealT} + # Check if the function output is of the correct type + if !(typeof(h_boundary(one(RealT))) == RealT) + throw(ArgumentError("Boundary value functions must return a value of type $(RealT)")) + end + return BoundaryConditionWaterHeight(t -> h_boundary(t)) +end + +# Version for `TreeMesh` +function (boundary_condition::BoundaryConditionWaterHeight)(u_inner, + orientation, + direction, x, t, + surface_flux_functions, + equations::ShallowWaterEquationsWetDry2D) + # Extract the gravitational acceleration + g = equations.gravity + + # Get the water height and velocity from the inner state + h_inner = waterheight(u_inner, equations) + v1_inner, v2_inner = velocity(u_inner, equations) + + # Extract the external water height from the boundary condition + h_boundary = boundary_condition.h_boundary(t) + + # Calculate the boundary state based on the direction. + # To extrapolate the external velocity assume that the Riemann invariant remains constant across + # the incoming characteristic. In the case of inflow we assume that the tangential velocity at + # the boundary is zero. + if direction == 1 # x- + v1_boundary = v1_inner + + 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + v1_boundary > 0 ? hv2_boundary = zero(u_inner[3]) : hv2_boundary = u_inner[3] + u_boundary = SVector(h_boundary, h_boundary * v1_boundary, hv2_boundary, + u_inner[4]) + elseif direction == 2 # x+ + v1_boundary = v1_inner - + 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + v1_boundary < 0 ? hv2_boundary = zero(u_inner[3]) : hv2_boundary = u_inner[3] + u_boundary = SVector(h_boundary, h_boundary * v1_boundary, hv2_boundary, + u_inner[4]) + elseif direction == 3 # y- + v2_boundary = v2_inner + + 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + v2_boundary > 0 ? hv1_boundary = zero(u_inner[2]) : hv1_boundary = u_inner[2] + u_boundary = SVector(h_boundary, hv1_boundary, h_boundary * v2_boundary, + u_inner[4]) + elseif direction == 4 # y+ + v2_boundary = v2_inner - + 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + v2_boundary < 0 ? hv1_boundary = zero(u_inner[2]) : hv1_boundary = u_inner[2] + u_boundary = SVector(h_boundary, hv1_boundary, h_boundary * v2_boundary, + u_inner[4]) + end + + # Evaluate the conservative flux at the boundary + flux = Trixi.flux(u_boundary, orientation, equations) + + # Return the conservative and nonconservative fluxes. + # The nonconservative part is zero as we assume a constant bottom topography at the boundary. + return (flux, zero(u_inner)) +end + +# Version for `UnstructuredMesh2D` and `P4estMesh` +function (boundary_condition::BoundaryConditionWaterHeight)(u_inner, + normal_direction, + x, t, + surface_flux_functions, + equations::ShallowWaterEquationsWetDry2D) + # Extract the gravitational acceleration + g = equations.gravity + + # Normalized normal vector + normal = normal_direction / Trixi.norm(normal_direction) + + # Apply the rotation that maps `normal` onto the x-axis to `u_inner`. + u_rotated = Trixi.rotate_to_x(u_inner, normal, equations) + + # Get the water height and velocity from the inner state + h_inner = waterheight(u_rotated, equations) + v_inner_normal, _ = velocity(u_rotated, equations) + + # Extract the external water height from the boundary condition + h_boundary = boundary_condition.h_boundary(t) + + # Calculate the boundary state in the rotated coordinate system. + # To extrapolate the external velocity assume that the Riemann invariant remains constant across + # the incoming characteristic. In the case of inflow we assume that the tangential velocity at + # the boundary is zero. + v_boundary_normal = v_inner_normal - 2 * (sqrt(g * h_boundary) - sqrt(g * h_inner)) + v_boundary_normal < 0 ? hv_boundary_tangential = zero(u_rotated[3]) : + hv_boundary_tangential = u_rotated[3] + + u_boundary = SVector(h_boundary, h_boundary * v_boundary_normal, + hv_boundary_tangential, u_rotated[4]) + + # Compute the boundary flux in the rotated coordinate system. + flux = Trixi.flux(u_boundary, 1, equations) + + # Apply the back-rotation that maps the x-axis onto `normal` to the boundary flux. + flux = Trixi.rotate_from_x(flux, normal, equations) * Trixi.norm(normal_direction) + + # Return the conservative and nonconservative fluxes. The nonconservative part is zero as we assume + # a constant bottom topography at the boundary. + return (flux, zero(u_inner)) +end + +""" + BoundaryConditionMomentum(hv1_boundary, hv2_boundary, equations::ShallowWaterEquationsWetDry2D) + +Create a boundary condition that sets a fixed momentum in x- and y- directions, `hv1_boundary` +and `hv2_boundary`, at the boundary and extrapolates the water height `h_boundary` from the incoming +Riemann invariant. + +The external momentum can be specified as a constant value or as a function of time, e.g. +```julia + BoundaryConditionMomentum(hv1_boundary, hv2_boundary, equations) + BoundaryConditionMomentum(t -> hv1_boundary(t), t -> hv2_boundary(t), equations) +``` + +More details can be found in the paper: +- Lixiang Song, Jianzhong Zhou, Jun Guo, Qiang Zou, Yi Liu (2011) + A robust well-balanced finite volume model for shallow water flows + with wetting and drying over irregular terrain + [doi: 10.1016/j.advwatres.2011.04.017](https://doi.org/10.1016/j.advwatres.2011.04.017) + +!!! warning "Experimental code" + This is an experimental feature and can change any time. +""" +function BoundaryConditionMomentum(hv1_boundary::Real, hv2_boundary::Real, + equations::ShallowWaterEquationsWetDry2D{RealT}) where {RealT} + # Convert function output to the correct type + hv1_boundary = convert(RealT, hv1_boundary) + hv2_boundary = convert(RealT, hv2_boundary) + return BoundaryConditionMomentum((t -> (hv1_boundary, hv2_boundary))) +end + +function BoundaryConditionMomentum(hv1_boundary::Function, hv2_boundary::Function, + equations::ShallowWaterEquationsWetDry2D{RealT}) where {RealT} + # Check if the function output is of the correct type + if !(typeof(hv1_boundary(one(RealT))) == RealT && + typeof(hv2_boundary(one(RealT))) == RealT) + throw(ArgumentError("Boundary value functions must return a value of type $(RealT)")) + end + return BoundaryConditionMomentum(t -> (hv1_boundary(t), hv2_boundary(t))) +end + +# Version for `TreeMesh` +function (boundary_condition::BoundaryConditionMomentum)(u_inner, + orientation, + direction, x, t, + surface_flux_functions, + equations::ShallowWaterEquationsWetDry2D) + # Extract the gravitational acceleration + g = equations.gravity + + # Get the water height and velocity from the inner state + h_inner = waterheight(u_inner, equations) + v1_inner, v2_inner = velocity(u_inner, equations) + + # Extract the external momentum from the boundary condition + hv1_boundary, hv2_boundary = boundary_condition.hv_boundary(t) + + # Calculate the boundary state based on the direction. + # To extrapolate the external water height `h_boundary` assume that the Riemann invariant remains + # constant across the incoming characteristic. + # Requires one to solve for the roots of a nonlinear function, see Eq. (52) in the reference above. + # For convenience we substitute x = h_boundary and solve for x, using the Steffensen method. + if direction == 1 # x- + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) + + (v1_inner - 2 * sqrt(g * h_inner)) * x - hv1_boundary, + h_inner) + h_boundary = solve(fx, Order2()) + if hv1_boundary > 0 + u_boundary = SVector(h_boundary, hv1_boundary, hv2_boundary, u_inner[4]) + else + u_boundary = SVector(h_boundary, hv1_boundary, u_inner[3], u_inner[4]) + end + elseif direction == 2 # x+ + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) - + (v1_inner + 2 * sqrt(g * h_inner)) * x + hv1_boundary, + h_inner) + h_boundary = solve(fx, Order2()) + if hv1_boundary < 0 + u_boundary = SVector(h_boundary, hv1_boundary, hv2_boundary, u_inner[4]) + else + u_boundary = SVector(h_boundary, hv1_boundary, u_inner[3], u_inner[4]) + end + elseif direction == 3 # y- + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) + + (v2_inner - 2 * sqrt(g * h_inner)) * x - hv2_boundary, + h_inner) + h_boundary = solve(fx, Order2()) + if hv2_boundary > 0 + u_boundary = SVector(h_boundary, hv1_boundary, hv2_boundary, u_inner[4]) + else + u_boundary = SVector(h_boundary, u_inner[2], hv2_boundary, u_inner[4]) + end + elseif direction == 4 # y+ + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) - + (v2_inner + 2 * sqrt(g * h_inner)) * x + hv2_boundary, + h_inner) + h_boundary = solve(fx, Order2()) + if hv2_boundary < 0 + u_boundary = SVector(h_boundary, hv1_boundary, hv2_boundary, u_inner[4]) + else + u_boundary = SVector(h_boundary, u_inner[2], hv2_boundary, u_inner[4]) + end + end + + # Evaluate the conservative flux at the boundary + flux = Trixi.flux(u_boundary, orientation, equations) + + # Return the conservative and nonconservative fluxes. + # The nonconservative part is zero as we assume a constant bottom topography at the boundary. + return (flux, zero(u_inner)) +end + +# Version for `UnstructuredMesh2D` and `P4estMesh` +function (boundary_condition::BoundaryConditionMomentum)(u_inner, + normal_direction, + x, t, + surface_flux_functions, + equations::ShallowWaterEquationsWetDry2D) + + # Extract the gravitational acceleration + g = equations.gravity + + # Normalized normal vector + normal = normal_direction / Trixi.norm(normal_direction) + + # Apply the rotation that maps `normal` onto the x-axis to `u_inner` and `hv_boundary``. + u_rotated = Trixi.rotate_to_x(u_inner, normal, equations) + + # Get the water height and velocity from the inner state + h_inner = waterheight(u_rotated, equations) + v_inner_normal, _ = velocity(u_rotated, equations) + + # Extract the external momentum from the boundary condition + hv1_boundary, hv2_boundary = boundary_condition.hv_boundary(t) + + hv_boundary_normal = hv1_boundary * normal[1] + hv2_boundary * normal[2] + hv_boundary_tangential = -hv1_boundary * normal[2] + hv2_boundary * normal[1] + + # Calculate the boundary state in the rotated coordinate system. + # To extrapolate the external water height `h_boundary` assume that the Riemann invariant remains + # constant across the incoming characteristic. + # Requires one to solve for the roots of a nonlinear function, see Eq. (52) in the reference above. + # For convenience we substitute x = h_boundary and solve for x. + fx = ZeroProblem(x -> 2 * sqrt(g) * x^(3 / 2) - + (v_inner_normal + 2 * sqrt(g * h_inner)) * x + + hv_boundary_normal, h_inner) + h_boundary = solve(fx, Order2()) + + hv_boundary_normal < 0 ? nothing : hv_boundary_tangential = u_rotated[3] + + u_boundary = SVector(h_boundary, hv_boundary_normal, hv_boundary_tangential, + u_inner[4]) + + # Compute the boundary flux in the rotated coordinate system. + flux = Trixi.flux(u_boundary, 1, equations) + + # Apply the back-rotation that maps the x-axis onto `normal` to the boundary flux. + flux = Trixi.rotate_from_x(flux, normal, equations) * Trixi.norm(normal_direction) + + # Return the conservative and nonconservative fluxes. The nonconservative part is zero as we assume + # a constant bottom topography at the boundary. + return (flux, zero(u_inner)) +end + # Calculate 1D flux for a single point # Note, the bottom topography has no flux @inline function Trixi.flux(u, orientation::Integer, @@ -818,6 +1115,46 @@ end equations.basic_swe) end +@inline function Trixi.rotate_to_x(u, normal_vector, + equations::ShallowWaterEquationsWetDry2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 n_2 0; + # 0 t_1 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] + s * u[3], + -s * u[2] + c * u[3], + u[4]) +end + +@inline function Trixi.rotate_from_x(u, normal_vector, + equations::ShallowWaterEquationsWetDry2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 t_1 0; + # 0 n_2 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] - s * u[3], + s * u[2] + c * u[3], + u[4]) +end + # Helper function to extract the velocity vector from the conservative variables @inline function Trixi.velocity(u, equations::ShallowWaterEquationsWetDry2D) return Trixi.velocity(u, equations.basic_swe) @@ -858,7 +1195,7 @@ end end @inline function Trixi.waterheight_pressure(u, equations::ShallowWaterEquationsWetDry2D) - return Trixi.waterheight(u, equations) * Trixi.pressure(u, equations) + return waterheight(u, equations) * Trixi.pressure(u, equations) end # Entropy function for the shallow water equations is the total energy diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index c8f81cc2..494267d7 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -361,6 +361,104 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @trixi_testset "elixir_shallowwater_moving_water_subsonic" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_moving_water_subsonic.jl"), + l2=[ + 8.417228364425849e-8, + 1.2192883088824584e-7, + 5.18813898298437e-17 + ], + linf=[ + 1.0690909364452494e-6, + 1.031855541455684e-6, + 5.195496810550537e-16 + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_moving_water_subsonic (pressure inflow / momentum outflow)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_moving_water_subsonic.jl"), + l2=[ + 8.417228365380916e-8, + 1.2192883078461691e-7, + 5.18813898298437e-17 + ], + linf=[ + 1.0690909364452494e-6, + 1.031855541455684e-6, + 5.195496810550537e-16 + ], + boundary_conditions=(x_neg = boundary_condition_outflow, + x_pos = boundary_condition_inflow), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_moving_water_transonic" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_moving_water_transonic.jl"), + l2=[ + 1.4350478061746362e-8, + 2.7182491041365426e-8, + 5.18813898298437e-17 + ], + linf=[ + 2.0158945446269172e-7, + 2.638769658336315e-7, + 5.195496810550537e-16 + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_moving_water_shock" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_moving_water_shock.jl"), + l2=[ + 0.005085389784760944, + 0.0035188809785392672, + 5.18813898298437e-17 + ], + linf=[ + 0.1072385165561712, + 0.045157748898879274, + 5.195496810550537e-16 + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "elixir_shallowwater_shock_capturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_shock_capturing.jl"), diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index 9c59ea1b..468e76cf 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -390,6 +390,72 @@ isdir(outdir) && rm(outdir, recursive = true) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_inflow_outflow.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_inflow_outflow.jl"), + l2=[ + 0.164716617086721, + 0.5257126140039803, + 0.5257126140039803, + 0.0 + ], + linf=[ + 0.5595760580954796, + 1.3874204364467229, + 1.3874204364467246, + 0.0 + ], + # Increase iterations for coverage testing to trigger both inflow and outflow conditions + coverage_override=(maxiters = 65, tspan = (0.0, 2.5))) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_inflow_outflow_reverse.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_inflow_outflow.jl"), + l2=[ + 0.16471661708672095, + 0.52571261400398, + 0.5257126140039801, + 0.0 + ], + linf=[ + 0.5595760580954816, + 1.3874204364467226, + 1.3874204364467244, + 0.0 + ], + v1=0.1, v2=0.1, + boundary_condition_inflow=BoundaryConditionMomentum(t -> 0.1 - + 0.05 * + t, + t -> 0.1 - + 0.05 * + t, + equations), + boundary_conditions=(x_neg = boundary_condition_outflow, + x_pos = boundary_condition_inflow, + y_neg = boundary_condition_outflow, + y_pos = boundary_condition_inflow), + # Increase iterations for coverage testing to trigger both inflow and outflow conditions + coverage_override=(maxiters = 65, tspan = (0.0, 2.5))) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end end # SWE @testset "Two-Layer Shallow Water" begin diff --git a/test/test_unit.jl b/test/test_unit.jl index ed14318d..3aed0fcd 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -399,6 +399,70 @@ end @test_throws ArgumentError TrixiShallowWater.default_threshold_partially_wet(Int64) @test_throws ArgumentError TrixiShallowWater.default_threshold_desingularization(Int64) end + +@timed_testset "Consistency check for boundary condition arguments" begin + let + equations = ShallowWaterEquationsWetDry2D(gravity = 9.81) + u_inner = SVector(1.0, 0.3, 0.3, 0.1) + t = 1.0 + RealT = typeof(equations.gravity) + + # Check that the boundary condition functions have the correct output type + boundary_condition = BoundaryConditionWaterHeight(1.0, equations) + @test typeof(boundary_condition.h_boundary(t)) == RealT + @test BoundaryConditionWaterHeight(t -> 1.0, equations).h_boundary(t) == + boundary_condition.h_boundary(t) + @test BoundaryConditionWaterHeight(1, equations).h_boundary(t) == + boundary_condition.h_boundary(t) + @test BoundaryConditionWaterHeight(1.0f0, equations).h_boundary(t) == + boundary_condition.h_boundary(t) + + @test_throws ArgumentError BoundaryConditionWaterHeight(t -> 1, equations) + + boundary_condition = BoundaryConditionMomentum(0.3, 0.1, equations) + @test typeof(boundary_condition.hv_boundary(t)) == Tuple{RealT, RealT} + @test BoundaryConditionMomentum(t -> 0.3, t -> 0.1, equations).hv_boundary(t) == + boundary_condition.hv_boundary(t) + # Here we only check the type since 0.1f0 != 0.1 + @test typeof(BoundaryConditionMomentum(0.3, 0.1f0, equations).hv_boundary(t)) == + typeof(boundary_condition.hv_boundary(t)) + + @test_throws ArgumentError BoundaryConditionMomentum(t -> 0.3 * t, t -> 1, + equations) + @test_throws ArgumentError BoundaryConditionMomentum(t -> 0.3 * t, t -> 1.0f0, + equations) + end + + let + equations = ShallowWaterEquationsWetDry1D(gravity = 9.81) + u_inner = SVector(1.0, 0.3, 0.1) + t = 1.0 + RealT = typeof(equations.gravity) + + # Check that the boundary condition functions have the correct output type + boundary_condition = BoundaryConditionWaterHeight(1.0, equations) + @test typeof(boundary_condition.h_boundary(t)) == RealT + @test BoundaryConditionWaterHeight(t -> 1.0, equations).h_boundary(t) == + boundary_condition.h_boundary(t) + @test BoundaryConditionWaterHeight(1, equations).h_boundary(t) == + boundary_condition.h_boundary(t) + @test BoundaryConditionWaterHeight(1.0f0, equations).h_boundary(t) == + boundary_condition.h_boundary(t) + + @test_throws ArgumentError BoundaryConditionWaterHeight(t -> 1, equations) + + boundary_condition = BoundaryConditionMomentum(0.3, equations) + @test typeof(boundary_condition.hv_boundary(t)) == RealT + @test BoundaryConditionMomentum(t -> 0.3, equations).hv_boundary(t) == + boundary_condition.hv_boundary(t) + # Here we only check the type since 0.1f0 != 0.1 + @test typeof(BoundaryConditionMomentum(0.3f0, equations).hv_boundary(t)) == + typeof(boundary_condition.hv_boundary(t)) + + @test_throws ArgumentError BoundaryConditionMomentum(t -> 0.3f0, + equations) + end +end end # Unit tests end # module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index cf71add2..dc413f83 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -347,6 +347,32 @@ isdir(outdir) && rm(outdir, recursive = true) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_inflow_outflow.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_inflow_outflow.jl"), + l2=[ + 0.05395240720695749, + 0.17993463107947336, + 0.17890511959110378, + 0.0 + ], + linf=[ + 0.10340479800696867, + 0.3998099631746912, + 0.4041504686758526, + 0.0 + ], + tspan=(0.0, 3.0)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end end # SWE @testset "Two-Layer Shallow Water" begin