Add auxiliary variables, part 1

NEWS.md

@@ -8,6 +8,18 @@ for human readability.
 ## Changes in the v0.16 lifecycle
 
 #### Added
+- Experimental support for auxiliary variables ([#2348]).
+  An additional container `aux_vars` in `cache` is made available in central functions like
+  flux computations. Possible applications are steady background states, variable velocity
+  fields, geometrical information, or any other pointwise, passive (constant in time)
+  quantity that is required in addition to the unknows in the governing equations. The
+  auxiliary variables are set up by supplying a function to the
+  `SemidiscretizationHyperbolic` constructor via the keyword argument `aux_field`.
+  The current `equations` need to set `have_aux_node_vars to `True()` and `n_aux_node_vars`
+  to the number of auxiliary variables per node.
+  So far, a simplifying continuity assumption is made for the auxiliary variables, which
+  e.g. allows to directly compute values at neighboring (mortar) interfaces instead of
+  MPI-communicating their values.
 - Added support for plotting 1D solutions with Makie.jl, matching the existing Plots.jl interface ([#3035]).
 - `VolumeIntegralAdaptive` is now also available with `VolumeIntegralSubcellLimiting` for `TreeMesh` in 2D and 3D using the heuristic a-priori indicator `IndicatorHennemannGassner` ([#2924], [#2986]).
 - A new EOS type `AbstractHelmholtzEOS`, with concrete implementation `HelmholtzIdealGas`. This implementation roughly follows Klein et al.'s approach in

examples/tree_2d_dgsem/elixir_acoustics_convergence_auxvars.jl

@@ -0,0 +1,66 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the acoustic perturbation equations
+
+function auxiliary_variables_mean_values(x, equations)
+    # constant auxiliary variables (mean state)
+    return global_mean_vars(equations)
+end
+
+equations = AcousticPerturbationEquations2DAuxVars(v_mean_global = (0.5, 0.3),
+                                                   c_mean_global = 2.0,
+                                                   rho_mean_global = 0.9)
+
+initial_condition = initial_condition_convergence_test
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (2.0, 2.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000, periodicity = true) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test,
+                                    boundary_conditions = boundary_condition_periodic,
+                                    aux_field = auxiliary_variables_mean_values)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+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)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+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, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

examples/tree_2d_dgsem/elixir_acoustics_gauss_wall_amr_auxvars.jl

@@ -0,0 +1,95 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+using Plots
+
+###############################################################################
+# semidiscretization of the acoustic perturbation equations
+
+equations = AcousticPerturbationEquations2DAuxVars(v_mean_global = (0.5, 0.0),
+                                                   c_mean_global = 1.0,
+                                                   rho_mean_global = 1.0)
+
+# Create DG solver with polynomial degree = 5 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-100.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (100.0, 200.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                n_cells_max = 100_000,
+                periodicity = false)
+
+"""
+    initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2DAuxVars)
+
+A Gaussian pulse, used in the `gauss_wall` example elixir in combination with
+[`boundary_condition_wall`](@ref). Uses the global mean values from `equations`.
+"""
+function initial_condition_gauss_wall(x, t,
+                                      equations::AcousticPerturbationEquations2DAuxVars)
+    RealT = eltype(x)
+    v1_prime = 0
+    v2_prime = 0
+    p_prime = exp(-log(convert(RealT, 2)) * (x[1]^2 + (x[2] - 25)^2) / 25)
+    p_prime_scaled = p_prime / equations.c_mean_global^2
+
+    return SVector(v1_prime, v2_prime, p_prime_scaled)
+end
+initial_condition = initial_condition_gauss_wall
+
+function auxiliary_variables_mean_values(x, equations)
+    # constant auxiliary variables (mean state)
+    return global_mean_vars(equations)
+end
+
+@inline function p_prime(u, equations::AcousticPerturbationEquations2DAuxVars)
+    return u[3]
+end
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition_wall,
+                                    aux_field = auxiliary_variables_mean_values)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 30.0
+tspan = (0.0, 30.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+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 = cons2state)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = p_prime),
+                                      base_level = 2,
+                                      med_level = 4, med_threshold = 0.1,
+                                      max_level = 5, max_threshold = 0.2)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+# 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,
+                        amr_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks)

examples/tree_2d_dgsem/elixir_acoustics_gauss_wall_auxvars.jl

@@ -0,0 +1,81 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the acoustic perturbation equations
+
+equations = AcousticPerturbationEquations2DAuxVars(v_mean_global = (0.5, 0.0),
+                                                   c_mean_global = 1.0,
+                                                   rho_mean_global = 1.0)
+
+# Create DG solver with polynomial degree = 5 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-100.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (100.0, 200.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 100_000,
+                periodicity = false)
+
+"""
+    initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2DAuxVars)
+
+A Gaussian pulse, used in the `gauss_wall` example elixir in combination with
+[`boundary_condition_wall`](@ref). Uses the global mean values from `equations`.
+"""
+function initial_condition_gauss_wall(x, t,
+                                      equations::AcousticPerturbationEquations2DAuxVars)
+    RealT = eltype(x)
+    v1_prime = 0
+    v2_prime = 0
+    p_prime = exp(-log(convert(RealT, 2)) * (x[1]^2 + (x[2] - 25)^2) / 25)
+    p_prime_scaled = p_prime / equations.c_mean_global^2
+
+    return SVector(v1_prime, v2_prime, p_prime_scaled)
+end
+initial_condition = initial_condition_gauss_wall
+
+function auxiliary_variables_mean_values(x, equations)
+    # constant auxiliary variables (mean state)
+    return global_mean_vars(equations)
+end
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition_wall,
+                                    aux_field = auxiliary_variables_mean_values)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 30.0
+tspan = (0.0, 30.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+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 = cons2state)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+# 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, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks)

examples/tree_2d_dgsem/elixir_advection_variable_swirling_flow.jl

@@ -0,0 +1,126 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+using Plots
+
+# Advection test case following
+# Randall J. LeVeque
+# High-Resolution Conservative Algorithms for Advection in Incompressible Flow
+# https://doi.org/10.1137/0733033
+
+############################################################################################
+# initial condition
+
+equations = LinearVariableScalarAdvectionEquation2D()
+
+function initial_condition_advected_objects(x, t,
+                                            equations::LinearVariableScalarAdvectionEquation2D)
+    RealT = eltype(x)
+
+    # smooth hump
+    x_0, y_0, r_0 = 0.25f0, 0.5f0, convert(RealT, 0.15)
+    r = sqrt((x[1] - x_0)^2 + (x[2] - y_0)^2)
+    r = min(r, r_0) / r_0
+    hump = 0.25f0 * (1 + cospi(r))
+
+    # cone
+    x_1, y_1, r_1 = 0.5f0, 0.25f0, convert(RealT, 0.15)
+    r = sqrt((x[1] - x_1)^2 + (x[2] - y_1)^2)
+    cone = 1.0f0 - min(r, r_1) / r_1
+
+    # slotted disc
+    x_2, y_2, r_2 = 0.5f0, 0.75f0, convert(RealT, 0.15)
+    w, l = 0.05f0, convert(RealT, 0.25)
+    r = sqrt((x[1] - x_2)^2 + (x[2] - y_2)^2)
+    disc = 0
+    if r <= r_2 && (x[2] >= y_2 - r_2 + l || abs(x[1] - x_2) >= w)
+        disc = 1.0f0
+    end
+
+    return SVector(hump + cone + disc)
+end
+
+# velocity field at time t = 0 (see reference)
+@inline function velocity_swirling(x, equations)
+    u = sinpi(x[1])^2 * sinpi(2 * x[2])
+    v = -sinpi(x[2])^2 * sinpi(2 * x[1])
+    return SVector(u, v)
+end
+
+############################################################################################
+# semidiscretization
+polydeg = 3
+
+initial_condition = initial_condition_advected_objects
+
+surface_flux = flux_lax_friedrichs
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                periodicity = false)
+
+# boundary conditions
+boundary_condition_dirichlet = BoundaryConditionDirichlet(initial_condition_advected_objects)
+boundary_conditions = (; x_neg = boundary_condition_dirichlet,
+                       x_pos = boundary_condition_dirichlet,
+                       y_neg = boundary_condition_dirichlet,
+                       y_pos = boundary_condition_dirichlet)
+
+# the velocity is passed as auxiliary_field into the cache
+semi = SemidiscretizationHyperbolic(mesh,
+                                    equations,
+                                    initial_condition,
+                                    solver,
+                                    boundary_conditions = boundary_conditions,
+                                    aux_field = velocity_swirling)
+
+##############################################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+solution_variables = cons2prim
+
+analysis_callback = AnalysisCallback(semi,
+                                     interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     output_directory = "out_swirling",
+                                     solution_variables = solution_variables)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 2,
+                                      med_level = 4, med_threshold = 0.3,
+                                      max_level = 6, max_threshold = 0.8)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 20,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+visualization = VisualizationCallback(semi; interval = 20) #, show_mesh = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback,
+                        amr_callback,
+                        #visualization,
+                        save_solution)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, ode_default_options()..., callback = callbacks);