User defined RHS Splitting for IMEX

examples/p4est_2d_dgsem/elixir_euler_imex_warm_bubble.jl

@@ -0,0 +1,147 @@
+using Trixi
+using OrdinaryDiffEqBDF
+using SparseDiffTools ## This is needed to force 'autodiff = AutoFiniteDiff()' in the ODE solver.
+
+function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
+    g = 9.81
+    c_p = 1004.0
+    c_v = 717.0
+
+    # center of perturbation
+    center_x = 10000.0
+    center_z = 2000.0
+    # radius of perturbation
+    radius = 2000.0
+    # distance of current x to center of perturbation
+    r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2)
+
+    # perturbation in potential temperature
+    potential_temperature_ref = 300.0
+    potential_temperature_perturbation = 0.0
+    if r <= radius
+        potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2
+    end
+    potential_temperature = potential_temperature_ref + potential_temperature_perturbation
+
+    # Exner pressure, solves hydrostatic equation for x[2]
+    exner = 1 - g / (c_p * potential_temperature) * x[2]
+
+    # pressure
+    p_0 = 100_000.0  # reference pressure
+    R = c_p - c_v    # gas constant (dry air)
+    p = p_0 * exner^(c_p / R)
+
+    # temperature
+    T = potential_temperature * exner
+
+    # density
+    rho = p / (R * T)
+
+    v1 = 20.0
+    v2 = 0.0
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+@inline function boundary_condition_slip_wall_vel(u_inner, normal_direction::AbstractVector,
+                                                  x, t,
+                                                  surface_flux_function,
+                                                  equations::CompressibleEulerEquations2D)
+    # normalize the outward pointing direction
+    normal = normal_direction / Trixi.norm(normal_direction)
+
+    # compute the normal velocity
+    u_normal = normal[1] * u_inner[2] + normal[2] * u_inner[3]
+
+    # create the "external" boundary solution state
+    u_boundary = SVector(u_inner[1],
+                         u_inner[2] - 2 * u_normal * normal[1],
+                         u_inner[3] - 2 * u_normal * normal[2],
+                         u_inner[4])
+
+    # calculate the boundary flux
+    flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+
+    return flux
+end
+
+@inline function boundary_condition_zero(u_inner, normal_direction::AbstractVector,
+                                         x, t,
+                                         surface_flux_function,
+                                         equations::CompressibleEulerEquations2D)
+    flux = flux_zero(u_inner, u_inner, normal_direction, equations)
+
+    return flux
+end
+
+@inline function flux_zero(u_ll, u_rr, normal_direction,
+                           equations::CompressibleEulerEquations2D)
+    return SVector(0.0, 0.0, 0.0, 0.0)
+end
+
+@inline function source_terms_gravity(u, x, t, equations::CompressibleEulerEquations2D)
+    g = 9.81
+    rho, _, rho_v2, _ = u
+    return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2)
+end
+
+equations = CompressibleEulerEquations2D(1004 / 717)
+
+polydeg = 2
+basis = LobattoLegendreBasis(polydeg)
+
+volume_integral_explicit = VolumeIntegralFluxDifferencing(flux_kennedy_gruber)
+solver_explicit = DGSEM(basis, FluxLMARS(340.0), volume_integral_explicit)
+
+volume_integral_implicit = VolumeIntegralFluxDifferencing(flux_zero)
+solver_implicit = DGSEM(basis, flux_zero, volume_integral_implicit)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (20_000.0, 10_000.0)
+
+trees_per_dimension = (16, 16)
+
+mesh = P4estMesh(trees_per_dimension, polydeg = polydeg,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (true, false), initial_refinement_level = 0)
+
+boundary_conditions_explicit = Dict(:y_neg => boundary_condition_slip_wall_vel,
+                                    :y_pos => boundary_condition_slip_wall_vel)
+boundary_conditions_implicit = Dict(:y_neg => boundary_condition_zero,
+                                    :y_pos => boundary_condition_zero)
+
+initial_condition = initial_condition_warm_bubble
+
+semi = SemidiscretizationHyperbolicSplit(mesh,
+                                         (equations, equations),
+                                         initial_condition,
+                                         solver_implicit,
+                                         solver_explicit;
+                                         boundary_conditions = (boundary_conditions_implicit,
+                                                                boundary_conditions_explicit),
+                                         source_terms = (source_terms_gravity, nothing),)
+T = 10.0
+tspan = (0.0, T)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100, solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode,
+            SBDF2(autodiff = AutoFiniteDiff());
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false,
+            callback = callbacks,);

examples/p4est_3d_dgsem/elixir_mhd_imex_shockcapturing_amr.jl

@@ -0,0 +1,150 @@
+using Trixi
+using OrdinaryDiffEqBDF
+using SparseDiffTools ## This is needed to force 'autodiff = AutoFiniteDiff()' in the ODE solver.
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations3D(1.4)
+function flux_zero(u_ll, u_rr, normal_direction, equations::IdealGlmMhdEquations3D)
+    return zero(u_ll)
+end
+"""
+    initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
+
+Weak magnetic blast wave setup taken from Section 6.1 of the paper:
+- A. M. Rueda-Ramírez, S. Hennemann, F. J. Hindenlang, A. R. Winters, G. J. Gassner (2021)
+  An entropy stable nodal discontinuous Galerkin method for the resistive MHD
+  equations. Part II: Subcell finite volume shock capturing
+  [doi: 10.1016/j.jcp.2021.110580](https://doi.org/10.1016/j.jcp.2021.110580)
+"""
+function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
+    # Center of the blast wave is selected for the domain [0, 3]^3
+    inicenter = (1.5, 1.5, 1.5)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    delta_0 = 0.1
+    r_0 = 0.3
+    lambda = exp(5.0 / delta_0 * (r - r_0))
+
+    prim_inner = SVector(1.2, 0.1, 0.0, 0.1, 0.9, 1.0, 1.0, 1.0, 0.0)
+    prim_outer = SVector(1.2, 0.2, -0.4, 0.2, 0.3, 1.0, 1.0, 1.0, 0.0)
+    prim_vars = (prim_inner + lambda * prim_outer) / (1.0 + lambda)
+
+    return prim2cons(prim_vars, equations)
+end
+initial_condition = initial_condition_blast_wave
+
+# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
+# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
+# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
+# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
+# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
+# `StepsizeCallback` (CFL-Condition) and less diffusion.
+surface_flux = (FluxLaxFriedrichs(max_abs_speed_naive), flux_nonconservative_powell)
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver_explicit = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+                        volume_integral = volume_integral)
+solver_implicit = DGSEM(polydeg = polydeg, surface_flux = (flux_zero, flux_zero),
+                        volume_integral = VolumeIntegralFluxDifferencing((flux_zero,
+                                                                          flux_zero)))
+
+# Mapping as described in https://arxiv.org/abs/2012.12040 but with slightly less warping.
+# The mapping will be interpolated at tree level, and then refined without changing
+# the geometry interpolant.
+function mapping(xi_, eta_, zeta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    xi = 1.5 * xi_ + 1.5
+    eta = 1.5 * eta_ + 1.5
+    zeta = 1.5 * zeta_ + 1.5
+
+    y = eta +
+        3 / 11 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+         cos(0.5 * pi * (2 * eta - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    x = xi +
+        3 / 11 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+         cos(2 * pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    z = zeta +
+        3 / 11 * (cos(0.5 * pi * (2 * x - 3) / 3) *
+         cos(pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    return SVector(x, y, z)
+end
+
+trees_per_dimension = (1, 1, 1)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 2,
+                 mapping = mapping,
+                 initial_refinement_level = 2,
+                 periodicity = true)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolicSplit(mesh,
+                                         (equations, equations),
+                                         initial_condition,
+                                         solver_implicit,
+                                         solver_explicit;
+                                         source_terms = (nothing, nothing),)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi,
+                                variable = density_pressure)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 2,
+                                      max_level = 4, max_threshold = 0.15)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+cfl = 1.4
+stepsize_callback = StepsizeCallback(cfl = cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode,
+            SBDF(order = 1, autodiff = AutoFiniteDiff());
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -147,6 +147,7 @@ include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
 include("semidiscretization/semidiscretization_coupled.jl")
+include("semidiscretization/semidiscretization_split.jl")
 include("time_integration/time_integration.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
@@ -276,6 +277,8 @@ export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integ
 
 export SemidiscretizationHyperbolicParabolic
 
+export SemidiscretizationHyperbolicSplit
+
 export SemidiscretizationEulerAcoustics
 
 export SemidiscretizationEulerGravity, ParametersEulerGravity,