diff --git a/NEWS.md b/NEWS.md index 05c246236b2..919cebfdd63 100644 --- a/NEWS.md +++ b/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 diff --git a/examples/tree_2d_dgsem/elixir_acoustics_convergence_auxvars.jl b/examples/tree_2d_dgsem/elixir_acoustics_convergence_auxvars.jl new file mode 100644 index 00000000000..4e5a17f79b9 --- /dev/null +++ b/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); diff --git a/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall_amr_auxvars.jl b/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall_amr_auxvars.jl new file mode 100644 index 00000000000..6377e452910 --- /dev/null +++ b/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) diff --git a/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall_auxvars.jl b/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall_auxvars.jl new file mode 100644 index 00000000000..c22a637fcaa --- /dev/null +++ b/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) diff --git a/examples/tree_2d_dgsem/elixir_advection_variable_swirling_flow.jl b/examples/tree_2d_dgsem/elixir_advection_variable_swirling_flow.jl new file mode 100644 index 00000000000..f1ce704de24 --- /dev/null +++ b/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); diff --git a/src/Trixi.jl b/src/Trixi.jl index 0d758c31ce0..f8b4a4574d5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -201,7 +201,7 @@ include("visualization/visualization.jl") # export types/functions that define the public API of Trixi.jl -export AcousticPerturbationEquations2D, +export AcousticPerturbationEquations2D, AcousticPerturbationEquations2DAuxVars, CompressibleEulerEquations1D, CompressibleEulerEquations2D, CompressibleEulerEquations3D, CompressibleEulerMulticomponentEquations1D, @@ -214,6 +214,7 @@ export AcousticPerturbationEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D, + LinearVariableScalarAdvectionEquation2D, InviscidBurgersEquation1D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, LinearizedEulerEquations1D, LinearizedEulerEquations2D, LinearizedEulerEquations3D, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 976954d8359..cf948b7f618 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -460,7 +460,19 @@ function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::Union{DGSEM, FDSBP}, cache) - # Calculate + # The entropy_timederivative may depend on auxiliary variables. + return analyze(entropy_timederivative, du, u, t, + mesh, + have_aux_node_vars(equations), equations, + dg, cache) +end + +@inline function analyze(::typeof(entropy_timederivative), du, u, t, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + have_aux_node_vars::False, equations, + dg::Union{DGSEM, FDSBP}, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, du) do u, i, j, element, equations, dg, du @@ -470,6 +482,23 @@ function analyze(::typeof(entropy_timederivative), du, u, t, end end +@inline function analyze(::typeof(entropy_timederivative), du, u, t, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + have_aux_node_vars::True, equations, + dg::Union{DGSEM, FDSBP}, cache) + @unpack aux_node_vars = cache.aux_vars + # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ + integrate_via_indices(u, mesh, equations, dg, cache, + du) do u, i, j, element, equations, dg, du + u_node = get_node_vars(u, equations, dg, i, j, element) + aux_node = get_aux_node_vars(aux_node_vars, equations, dg, i, j, element) + du_node = get_node_vars(du, equations, dg, i, j, element) + return dot(cons2entropy(u_node, aux_node, equations), du_node) + end +end + function analyze(::Val{:l2_divb}, du, u, t, mesh::TreeMesh{2}, equations, dg::DGSEM, cache) diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 0e91e75803f..7e8469f6df7 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -5,6 +5,34 @@ @muladd begin #! format: noindent +@inline function convert_to_solution_variables(u, solution_variables, cache, + have_aux_node_vars::False, + equations) + # Reinterpret the solution array as an array of conservative variables, + # compute the solution variables via broadcasting, and reinterpret the + # result as a plain array of floating point numbers + return Array(reinterpret(eltype(u), + solution_variables.(reinterpret(SVector{nvariables(equations), + eltype(u)}, u), + Ref(equations)))) +end + +@inline function convert_to_solution_variables(u, solution_variables, cache, + have_aux_node_vars::True, + equations) + @unpack aux_node_vars = cache.aux_vars + # Reinterpret the solution array as an array of conservative variables, + # compute the solution variables via broadcasting, and reinterpret the + # result as a plain array of floating point numbers + return Array(reinterpret(eltype(u), + solution_variables.(reinterpret(SVector{nvariables(equations), + eltype(u)}, u), + reinterpret(SVector{n_aux_node_vars(equations), + eltype(aux_node_vars)}, + aux_node_vars), + Ref(equations)))) +end + function save_solution_file(u, time, dt, timestep, mesh::Union{TreeMeshSerial, StructuredMesh, StructuredMeshView, @@ -31,14 +59,9 @@ function save_solution_file(u, time, dt, timestep, data = u n_vars = nvariables(equations) else - # Reinterpret the solution array as an array of conservative variables, - # compute the solution variables via broadcasting, and reinterpret the - # result as a plain array of floating point numbers - data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{nvariables(equations), - eltype(u)}, u), - Ref(equations)))) - + data = convert_to_solution_variables(u, solution_variables, cache, + have_aux_node_vars(equations), + equations) # Find out variable count by looking at output from `solution_variables` function n_vars = size(data, 1) end @@ -206,14 +229,9 @@ function save_solution_file(u, time, dt, timestep, data = u n_vars = nvariables(equations) else - # Reinterpret the solution array as an array of conservative variables, - # compute the solution variables via broadcasting, and reinterpret the - # result as a plain array of floating point numbers - data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{nvariables(equations), - eltype(u)}, u), - Ref(equations)))) - + data = convert_to_solution_variables(u, solution_variables, cache, + have_aux_node_vars(equations), + equations) # Find out variable count by looking at output from `solution_variables` function n_vars = size(data, 1) end diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index e32221131c4..9b3e3d74adc 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -149,7 +149,8 @@ function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, u = wrap_array(u_ode, mesh, equations, solver, cache) return cfl_hyperbolic(t) * max_dt(u, t, mesh, - have_constant_speed(equations), equations, + have_constant_speed(equations), + have_aux_node_vars(equations), equations, solver, cache) end @@ -171,7 +172,8 @@ function calculate_dt(u_ode, t, cfl_hyperbolic::Real, cfl_parabolic::Real, u = wrap_array(u_ode, mesh, equations, solver, cache) return cfl_hyperbolic * max_dt(u, t, mesh, - have_constant_speed(equations), equations, + have_constant_speed(equations), + have_aux_node_vars(equations), equations, solver, cache) end @@ -184,7 +186,8 @@ function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, u = wrap_array(u_ode, mesh, equations, solver, cache) dt_hyperbolic = cfl_hyperbolic(t) * max_dt(u, t, mesh, - have_constant_speed(equations), equations, + have_constant_speed(equations), + have_aux_node_vars(equations), equations, solver, cache) cfl_para = cfl_parabolic(t) @@ -199,14 +202,14 @@ function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, end end -function calc_max_scaled_speed(backend::Nothing, u, mesh, constant_speed, equations, dg, - cache) +function calc_max_scaled_speed(backend::Nothing, u, mesh, have_constant_speed, + have_aux_node_vars, equations, dg, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements max_scaled_speed = zero(eltype(u)) @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache) - max_lambda = max_scaled_speed_per_element(u, typeof(mesh), constant_speed, - equations, dg, + max_lambda = max_scaled_speed_per_element(u, typeof(mesh), have_constant_speed, + have_aux_node_vars, equations, dg, contravariant_vectors, inverse_jacobian, element) @@ -217,9 +220,8 @@ function calc_max_scaled_speed(backend::Nothing, u, mesh, constant_speed, equati return max_scaled_speed end -function calc_max_scaled_speed(backend::Backend, u, ::MeshT, constant_speed, equations, - dg, - cache) where {MeshT} +function calc_max_scaled_speed(backend::Backend, u, ::MeshT, have_constant_speed, + have_aux_node_vars, equations, dg, cache) where {MeshT} @unpack contravariant_vectors, inverse_jacobian = cache.elements num_elements = nelements(dg, cache) @@ -228,7 +230,7 @@ function calc_max_scaled_speed(backend::Backend, u, ::MeshT, constant_speed, equ # Provide a custom neutral and init element since we "reduce" over 1:num_elements max_scaled_speed = AcceleratedKernels.mapreduce(Base.max, 1:num_elements, backend; init, neutral) do element - max_scaled_speed_per_element(u, MeshT, constant_speed, + max_scaled_speed_per_element(u, MeshT, have_constant_speed, have_aux_node_vars, equations, dg, contravariant_vectors, inverse_jacobian, diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl index 8a029543575..1a58fcda407 100644 --- a/src/callbacks_step/stepsize_dg1d.jl +++ b/src/callbacks_step/stepsize_dg1d.jl @@ -6,8 +6,8 @@ #! format: noindent function max_dt(u, t, mesh::TreeMesh{1}, - constant_speed::False, equations, - dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere max_scaled_speed = nextfloat(zero(t)) @@ -54,7 +54,8 @@ function max_dt(u, t, mesh::TreeMesh{1}, end function max_dt(u, t, mesh::TreeMesh{1}, - constant_speed::True, equations, + have_constant_speed::True, + have_aux_node_vars::False, equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -96,8 +97,8 @@ function max_dt(u, t, mesh::TreeMesh, # for all dimensions end function max_dt(u, t, mesh::StructuredMesh{1}, - constant_speed::False, equations, - dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere max_scaled_speed = nextfloat(zero(t)) @@ -123,8 +124,8 @@ function max_dt(u, t, mesh::StructuredMesh{1}, end function max_dt(u, t, mesh::StructuredMesh{1}, - constant_speed::True, equations, - dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 201b112b65e..3d03c27664b 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -6,7 +6,8 @@ #! format: noindent function max_dt(u, t, mesh::TreeMesh{2}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -29,6 +30,32 @@ function max_dt(u, t, mesh::TreeMesh{2}, return 2 / (nnodes(dg) * max_scaled_speed) end +function max_dt(u, t, mesh::TreeMesh{2}, + have_constant_speed::False, have_aux_node_vars::True, + equations, dg::DG, cache) + @unpack aux_node_vars = cache.aux_vars + # to avoid a division by zero if the speed vanishes everywhere, + # e.g. for steady-state linear advection + max_scaled_speed = nextfloat(zero(t)) + + @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache) + max_lambda1 = max_lambda2 = zero(max_scaled_speed) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + aux_node = get_aux_node_vars(aux_node_vars, + equations, dg, i, j, element) + lambda1, lambda2 = max_abs_speeds(u_node, aux_node, equations) + max_lambda1 = max(max_lambda1, lambda1) + max_lambda2 = max(max_lambda2, lambda2) + end + inv_jacobian = cache.elements.inverse_jacobian[element] + max_scaled_speed = max(max_scaled_speed, + inv_jacobian * (max_lambda1 + max_lambda2)) + end + + return 2 / (nnodes(dg) * max_scaled_speed) +end + function max_dt(u, t, mesh::TreeMesh{2}, constant_diffusivity::False, equations, equations_parabolic::AbstractEquationsParabolic, @@ -54,9 +81,9 @@ function max_dt(u, t, mesh::TreeMesh{2}, end function max_dt(u, t, mesh::TreeMesh{2}, - constant_speed::True, equations, dg::DG, cache) - # Avoid division by zero if the speed vanishes everywhere, - # e.g. for steady-state linear advection + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) + # Avoid division by zero if the speed vanishes everywhere max_scaled_speed = nextfloat(zero(t)) max_lambda1, max_lambda2 = max_abs_speeds(equations) @@ -73,16 +100,37 @@ function max_dt(u, t, mesh::TreeMesh{2}, end function max_dt(u, t, mesh::TreeMeshParallel{2}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) + # call the method accepting a general `mesh::TreeMesh{2}` + # TODO: MPI, we should improve this; maybe we should dispatch on `u` + # and create some MPI array type, overloading broadcasting and mapreduce etc. + # Then, this specific array type should also work well with DiffEq etc. + dt = invoke(max_dt, + Tuple{typeof(u), typeof(t), TreeMesh{2}, + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) + # Base.min instead of min needed, see comment in src/auxiliary/math.jl + dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] + + return dt +end + +function max_dt(u, t, mesh::TreeMeshParallel{2}, + have_constant_speed::False, have_aux_node_vars::True, + equations, dg::DG, cache) # call the method accepting a general `mesh::TreeMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), TreeMesh{2}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -90,16 +138,18 @@ function max_dt(u, t, mesh::TreeMeshParallel{2}, end function max_dt(u, t, mesh::TreeMeshParallel{2}, - constant_speed::True, equations, dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::TreeMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), TreeMesh{2}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -113,11 +163,12 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}, StructuredMeshView{2}}, - constant_speed, equations, dg::DG, cache) + have_constant_speed, have_aux_node_vars, equations, dg::DG, cache) backend = trixi_backend(u) - max_lambda = calc_max_scaled_speed(backend, u, mesh, constant_speed, equations, dg, - cache) + max_lambda = calc_max_scaled_speed(backend, u, mesh, have_constant_speed, + have_aux_node_vars, + equations, dg, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -132,7 +183,9 @@ end P4estMesh{2}, T8codeMesh{2}, StructuredMeshView{2}}}, - constant_speed::False, equations, dg::DG, + have_constant_speed::False, + have_aux_node_vars::False, equations, + dg::DG, contravariant_vectors, inverse_jacobian, element) max_lambda1 = max_lambda2 = zero(eltype(u)) @@ -209,7 +262,9 @@ end P4estMeshView{2}, T8codeMesh{2}, StructuredMeshView{2}}}, - constant_speed::True, equations, dg::DG, + have_constant_speed::True, + have_aux_node_vars::False, equations, + dg::DG, contravariant_vectors, inverse_jacobian, element) max_scaled_speed = zero(eltype(u)) @@ -279,16 +334,18 @@ function max_dt(u, t, end function max_dt(u, t, mesh::P4estMeshParallel{2}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::P4estMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), P4estMesh{2}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -296,16 +353,18 @@ function max_dt(u, t, mesh::P4estMeshParallel{2}, end function max_dt(u, t, mesh::P4estMeshParallel{2}, - constant_speed::True, equations, dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::P4estMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), P4estMesh{2}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -313,16 +372,18 @@ function max_dt(u, t, mesh::P4estMeshParallel{2}, end function max_dt(u, t, mesh::T8codeMeshParallel{2}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::T8codeMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), T8codeMesh{2}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -330,16 +391,18 @@ function max_dt(u, t, mesh::T8codeMeshParallel{2}, end function max_dt(u, t, mesh::T8codeMeshParallel{2}, - constant_speed::True, equations, dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::T8codeMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), T8codeMesh{2}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index ca434918f53..5fa8bace5ee 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -6,7 +6,8 @@ #! format: noindent function max_dt(u, t, mesh::TreeMesh{3}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -57,7 +58,8 @@ function max_dt(u, t, mesh::TreeMesh{3}, end function max_dt(u, t, mesh::TreeMesh{3}, - constant_speed::True, equations, dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -78,11 +80,13 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, - constant_speed, equations, dg::DG, cache) + have_constant_speed, have_aux_node_vars::False, equations, dg::DG, + cache) backend = trixi_backend(u) - max_lambda = calc_max_scaled_speed(backend, u, mesh, constant_speed, equations, dg, - cache) + max_lambda = calc_max_scaled_speed(backend, u, mesh, have_constant_speed, + have_aux_node_vars, + equations, dg, cache) # Avoid division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -95,7 +99,8 @@ end ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - constant_speed::False, equations, dg, + have_constant_speed::False, + have_aux_node_vars::False, equations, dg, contravariant_vectors, inverse_jacobian, element) max_lambda1 = max_lambda2 = max_lambda3 = zero(eltype(u)) @@ -178,7 +183,9 @@ end ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - constant_speed::True, equations, dg::DG, + have_constant_speed::True, + have_aux_node_vars::False, equations, + dg::DG, contravariant_vectors, inverse_jacobian, element) max_scaled_speed = zero(eltype(u)) @@ -262,16 +269,18 @@ function max_dt(u, t, end function max_dt(u, t, mesh::P4estMeshParallel{3}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::P4estMesh{3}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), P4estMesh{3}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -279,16 +288,18 @@ function max_dt(u, t, mesh::P4estMeshParallel{3}, end function max_dt(u, t, mesh::P4estMeshParallel{3}, - constant_speed::True, equations, dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::P4estMesh{3}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), P4estMesh{3}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -296,16 +307,18 @@ function max_dt(u, t, mesh::P4estMeshParallel{3}, end function max_dt(u, t, mesh::T8codeMeshParallel{3}, - constant_speed::False, equations, dg::DG, cache) + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::T8codeMesh{3}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), T8codeMesh{3}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] @@ -313,16 +326,18 @@ function max_dt(u, t, mesh::T8codeMeshParallel{3}, end function max_dt(u, t, mesh::T8codeMeshParallel{3}, - constant_speed::True, equations, dg::DG, cache) + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DG, cache) # call the method accepting a general `mesh::T8codeMesh{3}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. dt = invoke(max_dt, Tuple{typeof(u), typeof(t), T8codeMesh{3}, - typeof(constant_speed), typeof(equations), typeof(dg), - typeof(cache)}, - u, t, mesh, constant_speed, equations, dg, cache) + typeof(have_constant_speed), typeof(have_aux_node_vars), + typeof(equations), typeof(dg), typeof(cache)}, + u, t, mesh, have_constant_speed, have_aux_node_vars, equations, dg, + cache) # Base.min instead of min needed, see comment in src/auxiliary/math.jl dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[] diff --git a/src/equations/acoustic_perturbation_2d_aux_vars.jl b/src/equations/acoustic_perturbation_2d_aux_vars.jl new file mode 100644 index 00000000000..7de4780619b --- /dev/null +++ b/src/equations/acoustic_perturbation_2d_aux_vars.jl @@ -0,0 +1,215 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + AcousticPerturbationEquations2DAuxVars(v_mean_global, c_mean_global, rho_mean_global) + +Alternative implementation of [`AcousticPerturbationEquations2D`](@ref) using auxiliary +variables. +""" +struct AcousticPerturbationEquations2DAuxVars{RealT <: Real} <: + AbstractAcousticPerturbationEquations{2, 3} + v_mean_global::SVector{2, RealT} + c_mean_global::RealT + rho_mean_global::RealT +end + +function AcousticPerturbationEquations2DAuxVars(; v_mean_global::NTuple{2, <:Real}, + c_mean_global::Real, + rho_mean_global::Real) + return AcousticPerturbationEquations2DAuxVars(SVector(v_mean_global), c_mean_global, + rho_mean_global) +end + +have_aux_node_vars(::AcousticPerturbationEquations2DAuxVars) = True() +n_aux_node_vars(::AcousticPerturbationEquations2DAuxVars) = 4 + +""" + global_mean_vars(equations::AcousticPerturbationEquations2DAuxVars) + +Returns the global mean variables stored in `equations`. This makes it easier +to define flexible initial conditions for problems with constant mean flow. +""" +function global_mean_vars(equations::AcousticPerturbationEquations2DAuxVars) + return equations.v_mean_global[1], equations.v_mean_global[2], + equations.c_mean_global, + equations.rho_mean_global +end + +""" + initial_condition_convergence_test(x, t, equations::AcousticPerturbationEquations2DAuxVars) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref). Uses the global mean values from `equations`. +""" +function initial_condition_convergence_test(x, t, + equations::AcousticPerturbationEquations2DAuxVars) + RealT = eltype(x) + a = 1 + c = 2 + L = 2 + f = 2.0f0 / L + A = convert(RealT, 0.2) + omega = 2 * convert(RealT, pi) * f + init = c + A * sin(omega * (x[1] + x[2] - a * t)) + + v1_prime = init + v2_prime = init + p_prime = init^2 + p_prime_scaled = p_prime / equations.c_mean_global^2 + + return SVector(v1_prime, v2_prime, p_prime_scaled) +end + +""" + source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2DAuxVars) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref). +""" +function source_terms_convergence_test(u, aux, x, t, + equations::AcousticPerturbationEquations2DAuxVars) + v1_mean, v2_mean, c_mean, rho_mean = aux + + RealT = eltype(u) + a = 1 + c = 2 + L = 2 + f = 2.0f0 / L + A = convert(RealT, 0.2) + omega = 2 * convert(RealT, pi) * f + + si, co = sincos(omega * (x[1] + x[2] - a * t)) + tmp = v1_mean + v2_mean - a + + du1 = du2 = A * omega * co * (2 * c / rho_mean + tmp + 2 / rho_mean * A * si) + du3 = A * omega * co * (2 * c_mean^2 * rho_mean + 2 * c * tmp + 2 * A * tmp * si) / + c_mean^2 + + return SVector(du1, du2, du3) +end + +""" + boundary_condition_wall(u_inner, aux_inner, orientation, direction, x, t, + surface_flux_function, + equations::AcousticPerturbationEquations2DAuxVars) + +Boundary conditions for a solid wall. +""" +function boundary_condition_wall(u_inner, aux_inner, orientation, direction, x, t, + surface_flux_function, + equations::AcousticPerturbationEquations2DAuxVars) + # Boundary state is equal to the inner state except for the perturbed velocity. For boundaries + # in the -x/+x direction, we multiply the perturbed velocity in the x direction by -1. + # Similarly, for boundaries in the -y/+y direction, we multiply the perturbed velocity in the + # y direction by -1 + if direction in (1, 2) # x direction + u_boundary = SVector(-u_inner[1], u_inner[2], u_inner[3]) + else # y direction + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3]) + end + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, aux_inner, aux_inner, + orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, aux_inner, aux_inner, + orientation, equations) + end + + return flux +end + +# Calculate 1D flux for a single point +@inline function flux(u, aux, orientation::Integer, + equations::AcousticPerturbationEquations2DAuxVars) + v1_prime, v2_prime, p_prime_scaled = u + v1_mean, v2_mean, c_mean, rho_mean = aux + + # Calculate flux for conservative state variables + RealT = eltype(u) + if orientation == 1 + f1 = v1_mean * v1_prime + v2_mean * v2_prime + + c_mean^2 * p_prime_scaled / rho_mean + f2 = zero(RealT) + f3 = rho_mean * v1_prime + v1_mean * p_prime_scaled + else + f1 = zero(RealT) + f2 = v1_mean * v1_prime + v2_mean * v2_prime + + c_mean^2 * p_prime_scaled / rho_mean + f3 = rho_mean * v2_prime + v2_mean * p_prime_scaled + end + + return SVector(f1, f2, f3) +end + +# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive` +@inline function max_abs_speed(u_ll, u_rr, aux_ll, aux_rr, orientation::Integer, + equations::AcousticPerturbationEquations2DAuxVars) + # Calculate v = v_prime + v_mean + v_prime_ll = u_ll[orientation] + v_prime_rr = u_rr[orientation] + v_mean_ll = aux_ll[orientation] + v_mean_rr = aux_rr[orientation] + + v_ll = v_prime_ll + v_mean_ll + v_rr = v_prime_rr + v_mean_rr + + c_mean_ll = aux_ll[3] + c_mean_rr = aux_rr[3] + + return max(abs(v_ll) + c_mean_ll, abs(v_rr) + c_mean_rr) +end + +@inline have_constant_speed(::AcousticPerturbationEquations2DAuxVars) = False() + +@inline function max_abs_speeds(u, aux, + equations::AcousticPerturbationEquations2DAuxVars) + v1_mean, v2_mean, c_mean, _ = aux + + return abs(v1_mean) + c_mean, abs(v2_mean) + c_mean +end + +function varnames(::typeof(cons2cons), ::AcousticPerturbationEquations2DAuxVars) + ("v1_prime", "v2_prime", "p_prime_scaled") +end + +# Convenience functions for retrieving state variables and mean variables +function cons2state(u, aux, ::AcousticPerturbationEquations2DAuxVars) + return SVector(u[1], u[2], u[3]) +end + +function varnames(::typeof(cons2state), ::AcousticPerturbationEquations2DAuxVars) + ("v1_prime", "v2_prime", "p_prime_scaled") +end + +function cons2aux(u, aux, ::AcousticPerturbationEquations2DAuxVars) + return SVector(aux[1], aux[2], aux[3], aux[4]) +end + +function varnames(::typeof(cons2aux), ::AcousticPerturbationEquations2DAuxVars) + ("v1_mean", "v2_mean", "c_mean", "rho_mean") +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, aux, equations::AcousticPerturbationEquations2DAuxVars) + p_prime_scaled = u[3] + c_mean = aux[3] + p_prime = p_prime_scaled * c_mean^2 + + return SVector(u[1], u[2], p_prime) +end + +function varnames(::typeof(cons2prim), ::AbstractAcousticPerturbationEquations{2}) + ("v1_prime", "v2_prime", "p_prime", + "v1_mean", "v2_mean", "c_mean", "rho_mean") +end + +# Convert conservative variables to entropy variables +@inline cons2entropy(u, aux, equations::AcousticPerturbationEquations2DAuxVars) = u +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index d1ffe6730f9..95ee72e1259 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -81,6 +81,14 @@ function Base.show(io::IO, ::MIME"text/plain", equations::AbstractEquations) "variable " * string(variable), varnames(cons2cons, equations)[variable]) end + if have_aux_node_vars(equations) == True() + summary_line(io, "#auxiliary variables", n_aux_node_vars(equations)) + for variable in eachauxvariable(equations) + summary_line(increment_indent(io), + "variable " * string(variable), + varnames(cons2aux, equations)[variable]) + end + end summary_footer(io) end end @@ -181,6 +189,27 @@ end return flux end +@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, aux_inner, + orientation_or_normal, + direction, + x, t, + surface_flux_function, + equations) + u_boundary = boundary_condition.boundary_value_function(x, t, equations) + # So far, there are no separate auxiliary variables on boundaries + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, aux_inner, aux_inner, + orientation_or_normal, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, aux_inner, aux_inner, + orientation_or_normal, equations) + end + + return flux +end + # Dirichlet-type boundary condition for use with TreeMesh or StructuredMesh # passing a tuple of surface flux functions for nonconservative terms @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, @@ -375,6 +404,24 @@ This is the default fallback for nonlinear equations. """ have_constant_speed(::AbstractEquations) = False() +""" + have_aux_node_vars(equations) +Trait function determining whether auxiliary variables should be stored at each node +alongside the solution variables and used to compute the flux for `equations`. When +`True()`, the signature of [`flux`](@ref) becomes the following: + +The return value will be `True()` or `False()` to allow dispatching on the return type. +""" +have_aux_node_vars(::AbstractEquations) = False() +""" + n_aux_node_vars(equations) + +Number of auxiliary variables used by `equations`. This function needs to be specialized +only if `equations` has auxiliary variables. +""" +function n_aux_node_vars end +@inline eachauxvariable(equations::AbstractEquations) = Base.OneTo(n_aux_node_vars(equations)) + """ default_analysis_errors(equations) @@ -630,6 +677,8 @@ include("linear_scalar_advection_1d.jl") include("linear_scalar_advection_2d.jl") include("linear_scalar_advection_3d.jl") +include("linear_variable_scalar_advection_2d.jl") + # Inviscid Burgers abstract type AbstractInviscidBurgersEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end @@ -768,6 +817,7 @@ include("lattice_boltzmann_3d.jl") abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("acoustic_perturbation_2d.jl") +include("acoustic_perturbation_2d_aux_vars.jl") # Linearized Euler equations abstract type AbstractLinearizedEulerEquations{NDIMS, NVARS} <: diff --git a/src/equations/linear_variable_scalar_advection_2d.jl b/src/equations/linear_variable_scalar_advection_2d.jl new file mode 100644 index 00000000000..10e6001645b --- /dev/null +++ b/src/equations/linear_variable_scalar_advection_2d.jl @@ -0,0 +1,75 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +struct LinearVariableScalarAdvectionEquation2D{} <: + AbstractLinearScalarAdvectionEquation{2} end + +have_aux_node_vars(::LinearVariableScalarAdvectionEquation2D) = True() +n_aux_node_vars(::LinearVariableScalarAdvectionEquation2D) = 2 + +@inline function flux(u, aux_vars, orientation::Integer, + equations::LinearVariableScalarAdvectionEquation2D) + a = aux_vars[orientation] + return a * u +end + +@inline function flux(u, aux_vars, normal_direction::AbstractVector, + equation::LinearVariableScalarAdvectionEquation2D) + a = dot(aux_vars, normal_direction) # velocity in normal direction + return a * u +end + +function flux_godunov(u_ll, u_rr, aux_ll, aux_rr, normal_direction::AbstractVector, + equation::LinearVariableScalarAdvectionEquation2D) + # velocity in normal direction + v_ll = dot(aux_ll, normal_direction) + v_rr = dot(aux_rr, normal_direction) + + a_normal = 0.5f0 * (v_ll + v_rr) + if a_normal >= 0 + return v_ll * u_ll + else + return v_rr * u_rr + end +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation +@inline function max_abs_speed(u_ll, u_rr, aux_ll, aux_rr, + orientation::Integer, + equation::LinearVariableScalarAdvectionEquation2D) + v_ll = aux_ll[orientation] + v_rr = aux_rr[orientation] + return max(abs(v_ll), abs(v_rr)) +end + +@inline function max_abs_speed(u_ll, u_rr, aux_ll, aux_rr, + normal_direction::AbstractVector, + equation::LinearVariableScalarAdvectionEquation2D) + # velocity in normal direction + v_ll = dot(aux_ll, normal_direction) + v_rr = dot(aux_rr, normal_direction) + return max(abs(v_ll), abs(v_rr)) +end + +# Maximum wave speeds in each direction for CFL calculation +@inline function Trixi.max_abs_speeds(u, aux_vars, + equations::LinearVariableScalarAdvectionEquation2D) + return abs.(aux_vars) +end + +@inline cons2entropy(u, aux, ::LinearVariableScalarAdvectionEquation2D) = u +@inline cons2prim(u, aux, ::LinearVariableScalarAdvectionEquation2D) = SVector(u[1], + aux[1], + aux[2]) +@inline cons2aux(u, aux, ::LinearVariableScalarAdvectionEquation2D) = SVector(aux[1], + aux[2]) + +varnames(::typeof(cons2cons), ::LinearVariableScalarAdvectionEquation2D) = ("scalar",) +varnames(::typeof(cons2prim), ::LinearVariableScalarAdvectionEquation2D) = ("scalar", + "v1", "v2") +varnames(::typeof(cons2aux), ::LinearVariableScalarAdvectionEquation2D) = ("v1", "v2") +end diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 545155b4889..56aa80a86c2 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -24,6 +24,18 @@ DG method (except floating point errors). return 0.5f0 * (f_ll + f_rr) end +# central flux for equations with auxiliary variables +@inline function flux_central(u_ll, u_rr, aux_ll, aux_rr, + orientation_or_normal_direction, + equations::AbstractEquations) + # Calculate regular 1D fluxes + f_ll = flux(u_ll, aux_ll, orientation_or_normal_direction, equations) + f_rr = flux(u_rr, aux_rr, orientation_or_normal_direction, equations) + + # Average regular fluxes + return 0.5f0 * (f_ll + f_rr) +end + """ FluxPlusDissipation(numerical_flux, dissipation) @@ -44,6 +56,18 @@ end dissipation(u_ll, u_rr, orientation_or_normal_direction, equations)) end +@inline function (numflux::FluxPlusDissipation)(u_ll, u_rr, aux_ll, aux_rr, + orientation_or_normal_direction, + equations) + @unpack numerical_flux, dissipation = numflux + + return (numerical_flux(u_ll, u_rr, aux_ll, aux_rr, + orientation_or_normal_direction, equations) + + + dissipation(u_ll, u_rr, aux_ll, aux_rr, + orientation_or_normal_direction, equations)) +end + function Base.show(io::IO, f::FluxPlusDissipation) print(io, "FluxPlusDissipation(", f.numerical_flux, ", ", f.dissipation, ")") return nothing @@ -177,6 +201,16 @@ DissipationLocalLaxFriedrichs() = DissipationLocalLaxFriedrichs(max_abs_speed) return -0.5f0 * λ * (u_rr - u_ll) end +# same as above for equations with auxiliary variables +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + aux_ll, aux_rr, + orientation_or_normal_direction, + equations) + λ = dissipation.max_abs_speed(u_ll, u_rr, aux_ll, aux_rr, + orientation_or_normal_direction, equations) + return -0.5f0 * λ * (u_rr - u_ll) +end + function Base.show(io::IO, d::DissipationLocalLaxFriedrichs) print(io, "DissipationLocalLaxFriedrichs(", d.max_abs_speed, ")") return nothing diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index 661fcdfd5ad..b0a278479b1 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -313,8 +313,8 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode) while !finalstep dtau = @trixi_timeit timer() "calculate dtau" begin cfl * max_dt(u_gravity, tau, semi_gravity.mesh, - have_constant_speed(equations), equations, - semi_gravity.solver, semi_gravity.cache) + have_constant_speed(equations), have_aux_node_vars(equations), + equations, semi_gravity.solver, semi_gravity.cache) end # evolve solution by one pseudo-time step diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 8d5132bbdd8..a455ba8ac9c 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -40,22 +40,40 @@ end source_terms=nothing, boundary_conditions, RealT=real(solver), - uEltype=RealT) + uEltype=RealT, + aux_field = nothing) Construct a semidiscretization of a hyperbolic PDE. Boundary conditions must be provided explicitly either as a `NamedTuple` or as a single boundary condition that gets applied to all boundaries. + +`aux_field` is an optional function taking a coordinate vector `x` and the current +`equations` as arguments. It is used to fill an additional field `aux_vars` in the `cache`, +which will be available, e.g., in flux computations. 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. Upon refinement, `aux_field` will be called again to recompute the auxiliary variables. +!!! warning "Experimental implementation" + The use of auxiliary variables is an experimental feature and may change in future + releases. Currently, only 2D `TreeMesh` is supported. """ function SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; source_terms = nothing, boundary_conditions, # `RealT` is used as real type for node locations etc. # while `uEltype` is used as element type of solutions etc. - RealT = real(solver), uEltype = RealT) + RealT = real(solver), uEltype = RealT, + aux_field = nothing) @assert ndims(mesh) == ndims(equations) cache = create_cache(mesh, equations, solver, RealT, uEltype) + + # Add auxiliary node variables cache + if have_aux_node_vars(equations) == True() + cache = (; cache..., + create_cache_aux(mesh, equations, solver, cache, aux_field)...) + end + _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) @@ -97,6 +115,12 @@ function remake(semi::SemidiscretizationHyperbolic; uEltype = real(semi.solver), source_terms, boundary_conditions, uEltype) end +# auxiliary variables cache +function create_cache_aux(mesh, equations, solver, cache, aux_field) + aux_vars = init_aux_vars(mesh, equations, solver, cache, aux_field) + return (; aux_vars) +end + # general fallback function digest_boundary_conditions(boundary_conditions, mesh, solver, cache) return boundary_conditions @@ -481,6 +505,10 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli print_boundary_conditions(io, semi) summary_line(io, "source terms", semi.source_terms) + if have_aux_node_vars(semi.equations) == True() + summary_line(io, "auxiliary variables", + semi.cache.aux_vars.aux_field) + end summary_line(io, "solver", semi.solver |> typeof |> nameof) summary_line(io, "total #DOFs per field", ndofsglobal(semi)) summary_footer(io) diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 83fc01c232b..78448b77725 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -1071,6 +1071,13 @@ https://docs.julialang.org/en/v1/manual/functions/#Varargs-Functions return SVector(ntuple(@inline(v->u[v, indices...]), Val(nvariables(equations)))) end +# Return the auxiliary variables at a given volume node index +@inline function get_aux_node_vars(aux_node_vars, equations, ::DG, + indices...) + return SVector(ntuple(@inline(v->aux_node_vars[v, indices...]), + Val(n_aux_node_vars(equations)))) +end + @inline function get_surface_node_vars(u, equations, solver::DG, indices...) # There is a cut-off at `n == 10` inside of the method # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17 @@ -1082,6 +1089,16 @@ end return u_ll, u_rr end +# Return the auxiliary variables at a given surface node index +@inline function get_aux_surface_node_vars(aux_surface_node_vars, equations, ::DG, + indices...) + aux_vars_ll = SVector(ntuple(@inline(v->aux_surface_node_vars[1, v, indices...]), + Val(n_aux_node_vars(equations)))) + aux_vars_rr = SVector(ntuple(@inline(v->aux_surface_node_vars[2, v, indices...]), + Val(n_aux_node_vars(equations)))) + return aux_vars_ll, aux_vars_rr +end + # As above but dispatches on an type argument @inline function get_surface_node_vars(u, equations, ::Type{<:DG}, indices...) u_ll = SVector(ntuple(@inline(v->u[1, v, indices...]), Val(nvariables(equations)))) @@ -1096,6 +1113,14 @@ end return nothing end +@inline function set_aux_node_vars!(aux_node_vars, aux_node, equations, + solver::DG, indices...) + for v in 1:n_aux_node_vars(equations) + aux_node_vars[v, indices...] = aux_node[v] + end + return nothing +end + @inline function add_to_node_vars!(u, u_node, equations, solver::DG, indices...) for v in eachvariable(equations) u[v, indices...] += u_node[v] diff --git a/src/solvers/dg_gpu.jl b/src/solvers/dg_gpu.jl index 5a16ade0a7b..891f9077985 100644 --- a/src/solvers/dg_gpu.jl +++ b/src/solvers/dg_gpu.jl @@ -13,21 +13,23 @@ using Atomix: @atomic end function calc_volume_integral!(backend::Backend, du, u, mesh, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral, dg::DGSEM, cache) nelements(dg, cache) == 0 && return nothing kernel! = volume_integral_KAkernel!(backend) kernel_cache = kernel_filter_cache(cache) - kernel!(du, u, typeof(mesh), have_nonconservative_terms, equations, + kernel!(du, u, typeof(mesh), have_nonconservative_terms, have_aux_node_vars, equations, volume_integral, dg, kernel_cache, ndrange = nelements(dg, cache)) return nothing end @kernel function volume_integral_KAkernel!(du, u, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral, dg::DGSEM, cache) element = @index(Global) volume_integral_kernel!(du, u, element, MeshT, have_nonconservative_terms, - equations, volume_integral, dg, cache) + have_aux_node_vars, equations, volume_integral, dg, cache) end diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 5c55c7e5006..e4986e09dd3 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -351,7 +351,8 @@ end # for the stepsize callback function max_dt(u, t, mesh::DGMultiMesh, - constant_speed::False, equations, dg::DGMulti{NDIMS}, + have_constant_speed::False, have_aux_node_vars::False, + equations, dg::DGMulti{NDIMS}, cache) where {NDIMS} @unpack md = mesh rd = dg.basis @@ -374,7 +375,8 @@ function max_dt(u, t, mesh::DGMultiMesh, end function max_dt(u, t, mesh::DGMultiMesh, - constant_speed::True, equations, dg::DGMulti{NDIMS}, + have_constant_speed::True, have_aux_node_vars::False, + equations, dg::DGMulti{NDIMS}, cache) where {NDIMS} @unpack md = mesh rd = dg.basis diff --git a/src/solvers/dgsem/calc_volume_integral.jl b/src/solvers/dgsem/calc_volume_integral.jl index 32f601f00f7..88634c29345 100644 --- a/src/solvers/dgsem/calc_volume_integral.jl +++ b/src/solvers/dgsem/calc_volume_integral.jl @@ -9,44 +9,48 @@ # dimension and meshtype agnostic, i.e., valid for all 1D, 2D, and 3D meshes. @inline function volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralWeakForm, dg, cache, alpha = true) weak_form_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, dg, cache, alpha) return nothing end @inline function volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralFluxDifferencing, dg, cache, alpha = true) @unpack volume_flux = volume_integral # Volume integral specific data flux_differencing_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux, dg, cache, alpha) return nothing end @inline function volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DGSEM, cache, alpha = true) @unpack volume_flux_fv = volume_integral # Volume integral specific data fv_kernel!(du, u, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg, cache, element, alpha) return nothing end @inline function volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralPureLGLFiniteVolumeO2, dg::DGSEM, cache, alpha = true) # Unpack volume integral specific data @@ -54,7 +58,7 @@ end cons2recon, recon2cons) = volume_integral fvO2_kernel!(du, u, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg, cache, element, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons, @@ -64,14 +68,15 @@ end end @inline function volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralAdaptive{<:IndicatorEntropyChange}, dg::DGSEM, cache) @unpack volume_integral_default, volume_integral_stabilized, indicator = volume_integral @unpack maximum_entropy_increase = indicator volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_integral_default, dg, cache) # Compute entropy production of the default volume integral. @@ -93,7 +98,8 @@ end du[.., element] .= zero(eltype(du)) volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_stabilized, dg, cache) end @@ -101,7 +107,8 @@ end end @inline function volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralEntropyCorrection, dg::DGSEM, cache) @unpack volume_integral_default, volume_integral_stabilized, indicator = volume_integral @@ -111,7 +118,7 @@ end # run default volume integral volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_integral_default, dg, cache) # Check entropy production of "high order" volume integral. @@ -148,7 +155,8 @@ end # Calculate entropy stable volume integral contribution volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_stabilized, dg, cache) dS_volume_integral_stabilized = -entropy_change_reference_element(du, u, @@ -178,12 +186,14 @@ end end function calc_volume_integral!(backend::Nothing, du, u, mesh, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral, dg::DGSEM, cache) MeshT = typeof(mesh) @threaded for element in eachelement(dg, cache) volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral, dg, cache) end @@ -191,7 +201,8 @@ function calc_volume_integral!(backend::Nothing, du, u, mesh, end @inline function calc_volume_integral!(backend::Nothing, du, u, mesh, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralAdaptive{<:IndicatorHennemannGassner}, dg::DGSEM, cache) @unpack volume_integral_default, volume_integral_stabilized, indicator = volume_integral @@ -212,11 +223,13 @@ end if default_volume_integral volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_default, dg, cache) else volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_stabilized, dg, cache) end end @@ -225,7 +238,8 @@ end end function calc_volume_integral!(backend::Nothing, du, u, mesh, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral::VolumeIntegralShockCapturingHGType, dg::DGSEM, cache) @unpack (indicator, volume_integral_default, @@ -247,18 +261,21 @@ function calc_volume_integral!(backend::Nothing, du, u, mesh, if dg_only volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_default, dg, cache) else # Calculate DG volume integral contribution volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_blend_high_order, dg, cache, 1 - alpha_element) # Calculate FV volume integral contribution volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_blend_low_order, dg, cache, alpha_element) end @@ -268,7 +285,8 @@ function calc_volume_integral!(backend::Nothing, du, u, mesh, end function calc_volume_integral!(backend::Nothing, du, u, mesh, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral::VolumeIntegralEntropyCorrectionShockCapturingCombined, dg::DGSEM, cache) (; volume_integral_default, volume_integral_stabilized, indicator) = volume_integral @@ -288,7 +306,8 @@ function calc_volume_integral!(backend::Nothing, du, u, mesh, @threaded for element in eachelement(dg, cache) # run default volume integral volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_default, dg, cache) # Check entropy production of "high order" volume integral. @@ -325,7 +344,8 @@ function calc_volume_integral!(backend::Nothing, du, u, mesh, # Calculate entropy stable volume integral contribution volume_integral_kernel!(du, u, element, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, volume_integral_stabilized, dg, cache) dS_volume_integral_stabilized = -entropy_change_reference_element(du, u, diff --git a/src/solvers/dgsem/indicators.jl b/src/solvers/dgsem/indicators.jl index 5399e077ffb..2c756fbfea4 100644 --- a/src/solvers/dgsem/indicators.jl +++ b/src/solvers/dgsem/indicators.jl @@ -137,7 +137,8 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh, equations, dg::DGSEM # Otherwise, `@threaded` does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. calc_indicator_hennemann_gassner!(indicator_hg, threshold, parameter_s, u, - element, mesh, equations, dg, cache) + element, mesh, have_aux_node_vars(equations), + equations, dg, cache) end if alpha_smooth diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 36282f6f645..c3b7409d57c 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -263,7 +263,7 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - have_nonconservative_terms, + have_nonconservative_terms, have_aux_node_vars, equations, surface_integral, dg::DGSEM{<:LobattoLegendreBasis}, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -357,7 +357,7 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, - have_nonconservative_terms, + have_nonconservative_terms, have_aux_node_vars, equations, surface_integral, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -869,7 +869,7 @@ end function calc_mortar_flux!(surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.mortars @@ -1202,7 +1202,8 @@ function rhs!(du, u, t, u_parent, semis, # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.volume_integral, dg, cache) end @@ -1214,7 +1215,8 @@ function rhs!(du, u, t, u_parent, semis, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.surface_integral, dg, cache) end @@ -1239,7 +1241,8 @@ function rhs!(du, u, t, u_parent, semis, # Calculate mortar fluxes @trixi_timeit timer() "mortar flux" begin calc_mortar_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end @@ -1255,7 +1258,8 @@ function rhs!(du, u, t, u_parent, semis, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, + have_aux_node_vars(equations), equations, dg, cache) end return nothing diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl index a3c7e5ac094..ba1ad9f6600 100644 --- a/src/solvers/dgsem_p4est/dg_2d_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl @@ -64,7 +64,7 @@ end function calc_interface_flux!(backend::Backend, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - have_nonconservative_terms, + have_nonconservative_terms, have_aux_node_vars, equations, surface_integral, dg::DGSEM{<:LobattoLegendreBasis}, cache) ninterfaces(cache.interfaces) == 0 && return nothing diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 2cefabd9539..97b0f3fe425 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -46,7 +46,7 @@ end function calc_mpi_interface_flux!(surface_flux_values, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - have_nonconservative_terms, + have_nonconservative_terms, have_aux_node_vars, equations, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements @@ -89,7 +89,7 @@ function calc_mpi_interface_flux!(surface_flux_values, calc_mpi_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, - equations, + have_aux_node_vars, equations, surface_integral, dg, cache, interface, normal_direction, node, local_side, @@ -111,7 +111,8 @@ end @inline function calc_mpi_interface_flux!(surface_flux_values, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, interface_node_index, local_side, @@ -140,7 +141,8 @@ end @inline function calc_mpi_interface_flux!(surface_flux_values, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, interface_node_index, local_side, @@ -240,7 +242,8 @@ end function calc_mpi_mortar_flux!(surface_flux_values, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, + equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = cache.mpi_mortars diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index b3bee323edf..b1d9ce8e0f6 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -184,7 +184,7 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - have_nonconservative_terms, + have_nonconservative_terms, have_aux_node_vars, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @unpack contravariant_vectors = cache.elements @@ -754,7 +754,8 @@ end function calc_mortar_flux!(surface_flux_values, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.mortars diff --git a/src/solvers/dgsem_p4est/dg_3d_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_gpu.jl index 5f67bdbb908..1a60d6837bf 100644 --- a/src/solvers/dgsem_p4est/dg_3d_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_3d_gpu.jl @@ -182,7 +182,7 @@ end function calc_interface_flux!(backend::Backend, surface_flux_values, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - have_nonconservative_terms, + have_nonconservative_terms, have_aux_node_vars, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @unpack contravariant_vectors = cache.elements diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl index 3370c157376..f947c1f51d4 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl @@ -36,7 +36,8 @@ function rhs!(du, u, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.volume_integral, dg, cache) end @@ -48,7 +49,8 @@ function rhs!(du, u, t, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.surface_integral, dg, cache) end @@ -72,7 +74,8 @@ function rhs!(du, u, t, # Calculate mortar fluxes @trixi_timeit timer() "mortar flux" begin calc_mortar_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end @@ -107,7 +110,8 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, + have_aux_node_vars(equations), equations, dg, cache) end # Finish to send MPI data diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 5ac27049151..25df809a279 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -50,7 +50,8 @@ function rhs!(du, u, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.volume_integral, dg, cache) end @@ -88,7 +89,8 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, + have_aux_node_vars(equations), equations, dg, cache) end return nothing diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 171c67b626e..2f7035591bc 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -37,7 +37,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -80,7 +81,8 @@ end UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) @unpack derivative_split = dg.basis @unpack contravariant_vectors = cache.elements @@ -144,9 +146,11 @@ end UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) - flux_differencing_kernel!(du, u, element, MeshT, have_nonconservative_terms, + flux_differencing_kernel!(du, u, element, MeshT, + have_nonconservative_terms, have_aux_node_vars, combine_conservative_and_nonconservative_fluxes(volume_flux, equations), equations, @@ -162,6 +166,7 @@ end P4estMesh{2}, T8codeMesh{2}}}, have_nonconservative_terms::True, + have_aux_node_vars::False, combine_conservative_and_nonconservative_fluxes::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) @@ -170,8 +175,9 @@ end symmetric_flux, nonconservative_flux = volume_flux # Apply the symmetric flux as usual - flux_differencing_kernel!(du, u, element, MeshT, False(), equations, symmetric_flux, - dg, cache, alpha) + flux_differencing_kernel!(du, u, element, MeshT, False(), have_aux_node_vars, + equations, + symmetric_flux, dg, cache, alpha) # Calculate the remaining volume terms using the nonsymmetric generalized flux for j in eachnode(dg), i in eachnode(dg) @@ -236,6 +242,7 @@ end P4estMesh{2}, T8codeMesh{2}}}, have_nonconservative_terms::True, + have_aux_node_vars::False, combine_conservative_and_nonconservative_fluxes::True, equations, volume_flux, dg::DGSEM, cache, alpha = true) @@ -305,7 +312,8 @@ end ::Type{<:Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack normal_vectors_1, normal_vectors_2 = cache.normal_vectors @@ -428,7 +436,8 @@ end ::Type{<:Union{StructuredMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack normal_vectors_1, normal_vectors_2 = cache.normal_vectors diff --git a/src/solvers/dgsem_structured/dg_2d_compressible_euler.jl b/src/solvers/dgsem_structured/dg_2d_compressible_euler.jl index 77b5bc51a09..f42374c0092 100644 --- a/src/solvers/dgsem_structured/dg_2d_compressible_euler.jl +++ b/src/solvers/dgsem_structured/dg_2d_compressible_euler.jl @@ -23,6 +23,7 @@ UnstructuredMesh2D, P4estMesh{2}}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations2D, volume_flux::typeof(flux_shima_etal_turbo), dg::DGSEM, cache, alpha) @@ -231,6 +232,7 @@ end UnstructuredMesh2D, P4estMesh{2}}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations2D, volume_flux::typeof(flux_ranocha_turbo), dg::DGSEM, cache, alpha) diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index b6e6129a31a..3c13b716a3f 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -37,7 +37,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 element, ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -95,7 +96,8 @@ end ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -178,9 +180,11 @@ end MeshT::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) flux_differencing_kernel!(du, u, element, MeshT, have_nonconservative_terms, + have_aux_node_vars, combine_conservative_and_nonconservative_fluxes(volume_flux, equations), equations, volume_flux, dg, cache, alpha) @@ -193,6 +197,7 @@ end P4estMesh{3}, T8codeMesh{3}}}, have_nonconservative_terms::True, + have_aux_node_vars::False, combine_conservative_and_nonconservative_fluxes::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) @@ -201,8 +206,9 @@ end symmetric_flux, nonconservative_flux = volume_flux # Apply the symmetric flux as usual - flux_differencing_kernel!(du, u, element, MeshT, False(), equations, symmetric_flux, - dg, cache, alpha) + flux_differencing_kernel!(du, u, element, MeshT, False(), have_aux_node_vars, + equations, + symmetric_flux, dg, cache, alpha) # Calculate the remaining volume terms using the nonsymmetric generalized flux for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) @@ -284,6 +290,7 @@ end P4estMesh{3}, T8codeMesh{3}}}, have_nonconservative_terms::True, + have_aux_node_vars::False, combine_conservative_and_nonconservative_fluxes::True, equations, volume_flux, dg::DGSEM, cache, alpha = true) @@ -379,6 +386,7 @@ end ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @unpack weights, derivative_matrix = dg.basis @@ -440,6 +448,7 @@ end ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @unpack weights, derivative_matrix = dg.basis @@ -535,7 +544,8 @@ end fstar3_L, fstar3_R, u, ::Type{<:Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons) diff --git a/src/solvers/dgsem_structured/dg_3d_compressible_euler.jl b/src/solvers/dgsem_structured/dg_3d_compressible_euler.jl index 2022eb9f3e6..f0da7dc7866 100644 --- a/src/solvers/dgsem_structured/dg_3d_compressible_euler.jl +++ b/src/solvers/dgsem_structured/dg_3d_compressible_euler.jl @@ -22,6 +22,7 @@ MeshT::Type{<:Union{StructuredMesh{3}, P4estMesh{3}}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations3D, volume_flux::typeof(flux_shima_etal_turbo), dg::DGSEM, cache, alpha) @@ -355,6 +356,7 @@ end MeshT::Type{<:Union{StructuredMesh{3}, P4estMesh{3}}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations3D, volume_flux::typeof(flux_ranocha_turbo), dg::DGSEM, cache, alpha) diff --git a/src/solvers/dgsem_tree/containers.jl b/src/solvers/dgsem_tree/containers.jl index 6e57612da5f..caaf845efbe 100644 --- a/src/solvers/dgsem_tree/containers.jl +++ b/src/solvers/dgsem_tree/containers.jl @@ -51,6 +51,226 @@ end @inline nvariables(interfaces::AbstractTreeMPIInterfaceContainer) = size(interfaces.u, 2) +# Container for storing values of auxiliary variables at volume/surface quadrature nodes +mutable struct AuxNodeVarsContainer{NDIMS, uEltype <: Real, NDIMSP2, NDIMSP3, AuxField} + aux_node_vars::Array{uEltype, NDIMSP2} # [var, i, j, element] + aux_surface_node_vars::Array{uEltype, NDIMSP2} # [leftright, var, i, interface] + aux_boundary_node_vars::Array{uEltype, NDIMSP2} # [leftright, var, i, boundary] + aux_mortar_node_vars::Array{uEltype, NDIMSP3} # [leftright, var, pos, i, mortar] + aux_mpiinterface_node_vars::Array{uEltype, NDIMSP2} # [leftright, var, i, interface] + aux_mpimortar_node_vars::Array{uEltype, NDIMSP3} # [leftright, var, pos, i, mortar] + + # internal `resize!`able storage + _aux_node_vars::Vector{uEltype} + _aux_surface_node_vars::Vector{uEltype} + _aux_boundary_node_vars::Vector{uEltype} + _aux_mortar_node_vars::Vector{uEltype} + _aux_mpiinterface_node_vars::Vector{uEltype} + _aux_mpimortar_node_vars::Vector{uEltype} + + # save initialization function + aux_field::AuxField +end + +nvariables(aux_vars::AuxNodeVarsContainer) = size(aux_vars.aux_node_vars, 1) +nnodes(aux_vars::AuxNodeVarsContainer) = size(aux_vars.aux_node_vars, 2) + +# Create auxiliary node variable container +function init_aux_vars(mesh, equations, solver, cache, aux_field) + @unpack elements, interfaces, boundaries, mortars = cache + + n_elements = nelements(elements) + n_interfaces = ninterfaces(interfaces) + n_boundaries = nboundaries(boundaries) + n_mortars = nmortars(mortars) + if mpi_isparallel() + n_mpiinterfaces = nmpiinterfaces(cache.mpi_interfaces) + n_mpimortars = nmpimortars(cache.mpi_mortars) + else + n_mpiinterfaces = 0 + n_mpimortars = 0 + end + NDIMS = ndims(mesh) + uEltype = eltype(elements) + nan_uEltype = convert(uEltype, NaN) + + _aux_node_vars = fill(nan_uEltype, + n_aux_node_vars(equations) * + nnodes(solver)^NDIMS * n_elements) + aux_node_vars = unsafe_wrap(Array, pointer(_aux_node_vars), + (n_aux_node_vars(equations), + ntuple(_ -> nnodes(solver), NDIMS)..., + n_elements)) + _aux_surface_node_vars = fill(nan_uEltype, + 2 * n_aux_node_vars(equations) * + nnodes(solver)^(NDIMS - 1) * + n_interfaces) + aux_surface_node_vars = unsafe_wrap(Array, + pointer(_aux_surface_node_vars), + (2, n_aux_node_vars(equations), + ntuple(_ -> nnodes(solver), + NDIMS - 1)..., + n_interfaces)) + _aux_boundary_node_vars = fill(nan_uEltype, + 2 * n_aux_node_vars(equations) * + nnodes(solver)^(NDIMS - 1) * + n_boundaries) + aux_boundary_node_vars = unsafe_wrap(Array, + pointer(_aux_boundary_node_vars), + (2, n_aux_node_vars(equations), + ntuple(_ -> nnodes(solver), + NDIMS - 1)..., + n_boundaries)) + _aux_mortar_node_vars = fill(nan_uEltype, + 2 * n_aux_node_vars(equations) * + 2^(NDIMS - 1) * nnodes(solver)^(NDIMS - 1) * + n_mortars) + aux_mortar_node_vars = unsafe_wrap(Array, + pointer(_aux_mortar_node_vars), + (2, n_aux_node_vars(equations), + 2^(NDIMS - 1), + ntuple(_ -> nnodes(solver), + NDIMS - 1)..., + n_mortars)) + _aux_mpiinterface_node_vars = fill(nan_uEltype, + 2 * n_aux_node_vars(equations) * + nnodes(solver)^(NDIMS - 1) * + n_mpiinterfaces) + aux_mpiinterface_node_vars = unsafe_wrap(Array, + pointer(_aux_mpiinterface_node_vars), + (2, n_aux_node_vars(equations), + ntuple(_ -> nnodes(solver), + NDIMS - 1)..., + n_mpiinterfaces)) + _aux_mpimortar_node_vars = fill(nan_uEltype, + 2 * n_aux_node_vars(equations) * + 2^(NDIMS - 1) * nnodes(solver)^(NDIMS - 1) * + n_mpimortars) + aux_mpimortar_node_vars = unsafe_wrap(Array, + pointer(_aux_mpimortar_node_vars), + (2, n_aux_node_vars(equations), + 2^(NDIMS - 1), + ntuple(_ -> nnodes(solver), + NDIMS - 1)..., + n_mpimortars)) + + aux_vars = AuxNodeVarsContainer{NDIMS, uEltype, NDIMS + 2, NDIMS + 3, + typeof(aux_field)}(aux_node_vars, + aux_surface_node_vars, + aux_boundary_node_vars, + aux_mortar_node_vars, + aux_mpiinterface_node_vars, + aux_mpimortar_node_vars, + _aux_node_vars, + _aux_surface_node_vars, + _aux_boundary_node_vars, + _aux_mortar_node_vars, + _aux_mpiinterface_node_vars, + _aux_mpimortar_node_vars, + aux_field) + init_aux_vars!(aux_vars, mesh, equations, solver, cache) + return aux_vars +end + +function init_aux_vars!(aux_vars, mesh, equations, solver, cache) + init_aux_node_vars!(aux_vars, mesh, equations, solver, cache) + init_aux_surface_node_vars!(aux_vars, mesh, equations, solver, cache) + init_aux_boundary_node_vars!(aux_vars, mesh, equations, solver, cache) + init_aux_mortar_node_vars!(aux_vars, mesh, equations, solver, cache) + if mpi_isparallel() + init_aux_mpiinterface_node_vars!(aux_vars, mesh, equations, solver, cache) + init_aux_mpimortar_node_vars!(aux_vars, mesh, equations, solver, cache) + end +end + +# Initialize auxiliary node variables (generic implementation) +function init_aux_node_vars!(aux_vars, mesh, equations, solver, cache) + @unpack aux_node_vars, aux_field = aux_vars + @unpack node_coordinates = cache.elements + + # all permutations of nodes indices for arbitrary dimension + node_cis = CartesianIndices(ntuple(i -> nnodes(solver), ndims(mesh))) + + @threaded for element in eachelement(solver, cache) + for node_ci in node_cis + x_local = get_node_coords(node_coordinates, equations, solver, + node_ci, element) + set_aux_node_vars!(aux_node_vars, + aux_field(x_local, equations), + equations, solver, node_ci, element) + end + end + return nothing +end + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(aux_vars::AuxNodeVarsContainer{NDIMS}, + capacity_node_vars, capacity_surface_node_vars, + capacity_boundary_node_vars, + capacity_mortar_node_vars, + capacity_mpiinterface_node_vars, + capacity_mpimortar_node_vars) where {NDIMS} + @unpack _aux_node_vars, _aux_surface_node_vars, _aux_boundary_node_vars, _aux_mortar_node_vars, _aux_mpiinterface_node_vars, _aux_mpimortar_node_vars = aux_vars + n_nodes = nnodes(aux_vars) + n_variables = nvariables(aux_vars) + nan_uEltype = convert(eltype(_aux_node_vars), NaN) + + resize!(_aux_node_vars, n_variables * n_nodes^NDIMS * capacity_node_vars) + aux_vars.aux_node_vars = unsafe_wrap(Array, + pointer(_aux_node_vars), + (n_variables, + ntuple(_ -> n_nodes, NDIMS)..., + capacity_node_vars)) + + resize!(_aux_surface_node_vars, + 2 * n_variables * n_nodes^(NDIMS - 1) * + capacity_surface_node_vars) + aux_vars.aux_surface_node_vars = unsafe_wrap(Array, + pointer(_aux_surface_node_vars), + (2, n_variables, + ntuple(_ -> n_nodes, NDIMS - 1)..., + capacity_surface_node_vars)) + resize!(_aux_boundary_node_vars, + 2 * n_variables * n_nodes^(NDIMS - 1) * + capacity_boundary_node_vars) + aux_vars.aux_boundary_node_vars = unsafe_wrap(Array, + pointer(_aux_boundary_node_vars), + (2, n_variables, + ntuple(_ -> n_nodes, NDIMS - 1)..., + capacity_boundary_node_vars)) + resize!(_aux_mortar_node_vars, + 2 * n_variables * 2^(NDIMS - 1) * n_nodes^(NDIMS - 1) * + capacity_mortar_node_vars) + aux_vars.aux_mortar_node_vars = unsafe_wrap(Array, + pointer(_aux_mortar_node_vars), + (2, n_variables, 2^(NDIMS - 1), + ntuple(_ -> n_nodes, NDIMS - 1)..., + capacity_mortar_node_vars)) + resize!(_aux_mpiinterface_node_vars, + 2 * n_variables * n_nodes^(NDIMS - 1) * + capacity_mpiinterface_node_vars) + aux_vars.aux_mpiinterface_node_vars = unsafe_wrap(Array, + pointer(_aux_mpiinterface_node_vars), + (2, n_variables, + ntuple(_ -> n_nodes, + NDIMS - 1)..., + capacity_mpiinterface_node_vars)) + resize!(_aux_mpimortar_node_vars, + 2 * n_variables * 2^(NDIMS - 1) * n_nodes^(NDIMS - 1) * + capacity_mpimortar_node_vars) + _aux_mpimortar_node_vars .= nan_uEltype + aux_vars.aux_mpimortar_node_vars = unsafe_wrap(Array, + pointer(_aux_mpimortar_node_vars), + (2, n_variables, 2^(NDIMS - 1), + ntuple(_ -> n_nodes, NDIMS - 1)..., + capacity_mpimortar_node_vars)) + return nothing +end + abstract type AbstractTreeBoundaryContainer <: AbstractBoundaryContainer end # Return number of boundaries diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index be76b91c102..5de70bca1fb 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -1452,6 +1452,354 @@ function reinitialize_containers!(mesh::Union{TreeMesh{2}, TreeMesh{3}}, equatio nvariables(equations), nnodes(dg), eltype(elements)) end + # re-initialize auxiliary variables container + if hasproperty(cache, :aux_vars) + @unpack aux_vars = cache + resize!(aux_vars, length(leaf_cell_ids), + count_required_interfaces(mesh, leaf_cell_ids), + count_required_boundaries(mesh, leaf_cell_ids), + count_required_mortars(mesh, leaf_cell_ids), + count_required_mpi_interfaces(mesh, leaf_cell_ids), + count_required_mpi_mortars(mesh, leaf_cell_ids)) + init_aux_vars!(aux_vars, mesh, equations, dg, cache) + end + + return nothing +end + +# Initialize auxiliary surface node variables +# 2D TreeMesh implementation, similar to prolong2interfaces +function init_aux_surface_node_vars!(aux_vars, mesh::TreeMesh2D, equations, solver, + cache) + @unpack aux_node_vars, aux_surface_node_vars = aux_vars + @unpack orientations, neighbor_ids = cache.interfaces + + @threaded for interface in eachinterface(solver, cache) + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] + + if orientations[interface] == 1 + # interface in x-direction + for j in eachnode(solver) + for v in axes(aux_surface_node_vars, 2) + aux_surface_node_vars[1, v, j, interface] = aux_node_vars[v, + nnodes(solver), + j, + left_element] + aux_surface_node_vars[2, v, j, interface] = aux_node_vars[v, + 1, + j, + right_element] + end + end + else # if orientations[interface] == 2 + # interface in y-direction + for i in eachnode(solver) + for v in axes(aux_surface_node_vars, 2) + aux_surface_node_vars[1, v, i, interface] = aux_node_vars[v, + i, + nnodes(solver), + left_element] + aux_surface_node_vars[2, v, i, interface] = aux_node_vars[v, + i, + 1, + right_element] + end + end + end + end return nothing end + +# Initialize auxiliary boundary node variables +# 2D TreeMesh implementation, similar to prolong2boundaries +function init_aux_boundary_node_vars!(aux_vars, mesh::TreeMesh2D, equations, solver, + cache) + @unpack aux_node_vars, aux_boundary_node_vars = aux_vars + @unpack orientations, neighbor_ids, neighbor_sides = cache.boundaries + + @threaded for boundary in eachboundary(solver, cache) + element = neighbor_ids[boundary] + if orientations[boundary] == 1 + # boundary in x-direction + if neighbor_sides[boundary] == 1 + # element in -x direction of boundary + for l in eachnode(solver), v in axes(aux_boundary_node_vars, 2) + aux_boundary_node_vars[1, v, l, boundary] = aux_node_vars[v, + nnodes(solver), + l, + element] + end + else # Element in +x direction of boundary + for l in eachnode(solver), v in axes(aux_boundary_node_vars, 2) + aux_boundary_node_vars[2, v, l, boundary] = aux_node_vars[v, 1, l, + element] + end + end + else # if orientations[boundary] == 2 + # boundary in y-direction + if neighbor_sides[boundary] == 1 + # element in -y direction of boundary + for l in eachnode(solver), v in axes(aux_boundary_node_vars, 2) + aux_boundary_node_vars[1, v, l, boundary] = aux_node_vars[v, l, + nnodes(solver), + element] + end + else + # element in +y direction of boundary + for l in eachnode(solver), v in axes(aux_boundary_node_vars, 2) + aux_boundary_node_vars[2, v, l, boundary] = aux_node_vars[v, l, 1, + element] + end + end + end + end + return nothing +end + +# Initialize auxiliary mortar node variables +# 2D TreeMesh implementation, similar to prolong2mortars +# Each mortar has two sides (identified by first variable of u_upper / u_lower) +# On the side with two small elements, values can be copied from the aux vars field +# On the side with one large element, values are usually interpolated to small elements +# We do this differently here and use the same small element values on both side. This +# assumes that the `aux_field` function computes a field with no jumps. +function init_aux_mortar_node_vars!(aux_vars, mesh::TreeMesh2D, equations, solver, + cache) + @unpack aux_node_vars, aux_mortar_node_vars = aux_vars + + @threaded for mortar in eachmortar(solver, cache) + upper_element = cache.mortars.neighbor_ids[2, mortar] + lower_element = cache.mortars.neighbor_ids[1, mortar] + + # Copy solution small to small + if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + for l in eachnode(solver) + for v in axes(aux_mortar_node_vars, 2) + aux_mortar_node_vars[1, v, 2, l, mortar] = aux_node_vars[v, 1, + l, + upper_element] + aux_mortar_node_vars[1, v, 1, l, mortar] = aux_node_vars[v, 1, + l, + lower_element] + end + end + else + # L2 mortars in y-direction + for l in eachnode(solver) + for v in axes(aux_mortar_node_vars, 2) + aux_mortar_node_vars[1, v, 2, l, mortar] = aux_node_vars[v, l, + 1, + upper_element] + aux_mortar_node_vars[1, v, 1, l, mortar] = aux_node_vars[v, l, + 1, + lower_element] + end + end + end + else # large_sides[mortar] == 2 -> small elements on left side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + for l in eachnode(solver) + for v in axes(aux_mortar_node_vars, 2) + aux_mortar_node_vars[1, v, 2, l, mortar] = aux_node_vars[v, + nnodes(solver), + l, + upper_element] + aux_mortar_node_vars[1, v, 1, l, mortar] = aux_node_vars[v, + nnodes(solver), + l, + lower_element] + end + end + else + # L2 mortars in y-direction + for l in eachnode(solver) + for v in axes(aux_mortar_node_vars, 2) + aux_mortar_node_vars[1, v, 2, l, mortar] = aux_node_vars[v, l, + nnodes(solver), + upper_element] + aux_mortar_node_vars[1, v, 1, l, mortar] = aux_node_vars[v, l, + nnodes(solver), + lower_element] + end + end + end + end + end + return nothing +end + +# Initialize auxiliary MPI interface node variables +# 2D TreeMesh implementation, similar to prolong2mpiinterfaces +# However we directly assign to both sides, assuming the aux field had no jumps. Therefore +# we do not need any exchange. +function init_aux_mpiinterface_node_vars!(aux_vars, mesh::TreeMeshParallel{2}, + equations, + solver, cache) + @unpack aux_node_vars, aux_mpiinterface_node_vars = aux_vars + @unpack mpi_interfaces = cache + + @threaded for interface in eachmpiinterface(solver, cache) + local_element = mpi_interfaces.local_neighbor_ids[interface] + + if mpi_interfaces.orientations[interface] == 1 # interface in x-direction + if mpi_interfaces.remote_sides[interface] == 1 # local element in positive direction + for j in eachnode(solver), v in axes(aux_mpiinterface_node_vars, 2) + aux_mpiinterface_node_vars[:, v, j, interface] .= aux_node_vars[v, + 1, + j, + local_element] + end + else # local element in negative direction + for j in eachnode(solver), v in axes(aux_mpiinterface_node_vars, 2) + aux_mpiinterface_node_vars[:, v, j, interface] .= aux_node_vars[v, + nnodes(solver), + j, + local_element] + end + end + else # interface in y-direction + if mpi_interfaces.remote_sides[interface] == 1 # local element in positive direction + for i in eachnode(solver), v in axes(aux_mpiinterface_node_vars, 2) + aux_mpiinterface_node_vars[:, v, i, interface] .= aux_node_vars[v, + i, + 1, + local_element] + end + else # local element in negative direction + for i in eachnode(solver), v in axes(aux_mpiinterface_node_vars, 2) + aux_mpiinterface_node_vars[:, v, i, interface] .= aux_node_vars[v, + i, + nnodes(solver), + local_element] + end + end + end + end + return nothing +end + +# Initialize auxiliary MPI mortar node variables +# 2D TreeMesh implementation, similar to prolong2mpimortars +# However: - We only assign the small element values (only leftright = 1 is used) +# - These have to be communicated +function init_aux_mpimortar_node_vars!(aux_vars, mesh::TreeMeshParallel{2}, equations, + solver, cache) + @unpack aux_node_vars, aux_mpimortar_node_vars = aux_vars + @unpack mpi_mortars = cache + + @threaded for mortar in eachmpimortar(solver, cache) + local_neighbor_ids = mpi_mortars.local_neighbor_ids[mortar] + local_neighbor_positions = mpi_mortars.local_neighbor_positions[mortar] + for (element, position) in zip(local_neighbor_ids, local_neighbor_positions) + if position in (1, 2) # small element + if mpi_mortars.large_sides[mortar] == 1 # -> small elements on right side + if mpi_mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + for l in eachnode(solver) + for v in axes(aux_mpimortar_node_vars, 2) + aux_mpimortar_node_vars[1, v, position, l, mortar] = aux_node_vars[v, + 1, + l, + element] + end + end + else + # L2 mortars in y-direction + for l in eachnode(solver) + for v in axes(aux_mpimortar_node_vars, 2) + aux_mpimortar_node_vars[1, v, position, l, mortar] = aux_node_vars[v, + l, + 1, + element] + end + end + end + else # large_sides[mortar] == 2 -> small elements on left side + if mpi_mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + for l in eachnode(solver) + for v in axes(aux_mpimortar_node_vars, 2) + aux_mpimortar_node_vars[1, v, position, l, mortar] = aux_node_vars[v, + nnodes(solver), + l, + element] + end + end + else + # L2 mortars in y-direction + for l in eachnode(solver) + for v in axes(aux_mpimortar_node_vars, 2) + aux_mpimortar_node_vars[1, v, position, l, mortar] = aux_node_vars[v, + l, + nnodes(solver), + element] + end + end + end + end + end + end + end + + data_size = nnodes(solver) * n_aux_node_vars(equations) + exchange_aux_mpimortars!(aux_mpimortar_node_vars, cache, data_size) + return nothing +end + +function exchange_aux_mpimortars!(aux_mpimortar_node_vars, cache, data_size) + @unpack mpi_cache = cache + + for rank in 1:length(mpi_cache.mpi_neighbor_ranks) + send_buffer = mpi_cache.mpi_send_buffers[rank] + mortars_data_size = length(mpi_cache.mpi_neighbor_mortars[rank]) * data_size * 2 + send_buffer[1:mortars_data_size] .= NaN + for (index, mortar) in enumerate(mpi_cache.mpi_neighbor_mortars[rank]) + base_index = (index - 1) * 2 * data_size + 1 + for position in cache.mpi_mortars.local_neighbor_positions[mortar] + # Determine whether the data belongs to the left or right side + if position in (1, 2) # small element + first = base_index + (position - 1) * data_size + last = first - 1 + data_size + @views send_buffer[first:last] .= vec(aux_mpimortar_node_vars[1, :, + position, + :, + mortar]) + end + end + end + end + + # Start data exchange + for (index, rank) in enumerate(mpi_cache.mpi_neighbor_ranks) + mpi_cache.mpi_send_requests[index] = MPI.Isend(mpi_cache.mpi_send_buffers[index], + rank, mpi_rank(), mpi_comm()) + mpi_cache.mpi_recv_requests[index] = MPI.Irecv!(mpi_cache.mpi_recv_buffers[index], + rank, rank, mpi_comm()) + end + + data = MPI.Waitany(mpi_cache.mpi_recv_requests) + while data !== nothing + recv_buffer = mpi_cache.mpi_recv_buffers[data] + for (index, mortar) in enumerate(mpi_cache.mpi_neighbor_mortars[data]) + base_index = (index - 1) * 2 * data_size + 1 + for position in 1:2 + first = base_index + (position - 1) * data_size + last = first - 1 + data_size + # Skip if received data for `position` is NaN as no real data has been sent for the + # corresponding element + if !isnan(recv_buffer[first]) + @views vec(aux_mpimortar_node_vars[1, :, position, :, mortar]) .= recv_buffer[first:last] + end + end + end + data = MPI.Waitany(mpi_cache.mpi_recv_requests) + end + + # Finish sending + MPI.Waitall(mpi_cache.mpi_send_requests, MPI.Status) +end end # @muladd diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 4c194c3595a..503b95824a4 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -72,7 +72,8 @@ function rhs!(du, u, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.volume_integral, dg, cache) end @@ -111,7 +112,8 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, + have_aux_node_vars(equations), equations, dg, cache) end return nothing @@ -127,7 +129,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 @inline function weak_form_kernel!(du, u, element, ::Type{<:Union{TreeMesh{1}, StructuredMesh{1}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -149,7 +152,8 @@ end @inline function flux_differencing_kernel!(du, u, element, ::Type{<:Union{TreeMesh{1}, StructuredMesh{1}}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -179,7 +183,8 @@ end @inline function flux_differencing_kernel!(du, u, element, MeshT::Type{<:Union{TreeMesh{1}, StructuredMesh{1}}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -187,8 +192,9 @@ end symmetric_flux, nonconservative_flux = volume_flux # Apply the symmetric flux as usual - flux_differencing_kernel!(du, u, element, MeshT, False(), equations, symmetric_flux, - dg, cache, alpha) + flux_differencing_kernel!(du, u, element, MeshT, False(), have_aux_node_vars, + equations, + symmetric_flux, dg, cache, alpha) # Calculate the remaining volume terms using the nonsymmetric generalized flux for i in eachnode(dg) @@ -214,7 +220,7 @@ end @inline function fv_kernel!(du, u, MeshT::Type{<:Union{TreeMesh{1}, StructuredMesh{1}}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded = cache @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes @@ -223,7 +229,7 @@ end fstar1_L = fstar1_L_threaded[Threads.threadid()] fstar1_R = fstar1_R_threaded[Threads.threadid()] calcflux_fv!(fstar1_L, fstar1_R, u, MeshT, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg, element, cache) # Calculate FV volume integral contribution @@ -240,7 +246,7 @@ end @inline function fvO2_kernel!(du, u, MeshT::Type{<:Union{TreeMesh{1}, StructuredMesh{1}}}, - nonconservative_terms, equations, + nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg::DGSEM, cache, element, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons, @@ -275,6 +281,7 @@ end @inline function calcflux_fv!(fstar1_L, fstar1_R, u, ::Type{<:Union{TreeMesh{1}, StructuredMesh{1}}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) for i in 2:nnodes(dg) u_ll = get_node_vars(u, equations, dg, i - 1, element) @@ -290,6 +297,7 @@ end @inline function calcflux_fv!(fstar1_L, fstar1_R, u, ::Type{<:TreeMesh{1}}, have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) volume_flux, nonconservative_flux = volume_flux_fv for i in 2:nnodes(dg) @@ -447,7 +455,8 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + equations, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -478,7 +487,8 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + equations, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces @@ -767,12 +777,14 @@ end # Need dimension specific version to avoid error at dispatching function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, - equations::AbstractEquations{1}, dg::DG, cache) + have_aux_node_vars::False, equations::AbstractEquations{1}, + dg::DG, cache) return nothing end function calc_sources!(backend::Nothing, du, u, t, source_terms, - equations::AbstractEquations{1}, dg::DG, cache) + have_aux_node_vars::False, equations::AbstractEquations{1}, + dg::DG, cache) @unpack node_coordinates = cache.elements @threaded for element in eachelement(dg, cache) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 7a7c5c7f5a7..acb8f12a00b 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -125,7 +125,8 @@ function rhs!(du, u, t, # Calculate volume integral @trixi_timeit_ext backend timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.volume_integral, dg, cache) end @@ -137,7 +138,8 @@ function rhs!(du, u, t, # Calculate interface fluxes @trixi_timeit_ext backend timer() "interface flux" begin calc_interface_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.surface_integral, dg, cache) end @@ -161,7 +163,8 @@ function rhs!(du, u, t, # Calculate mortar fluxes @trixi_timeit_ext backend timer() "mortar flux" begin calc_mortar_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end @@ -178,7 +181,8 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit_ext backend timer() "source terms" begin - calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, + have_aux_node_vars(equations), equations, dg, cache) end return nothing @@ -193,7 +197,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 =# @inline function weak_form_kernel!(du, u, element, ::Type{<:TreeMesh{2}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -219,8 +224,41 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 return nothing end +@inline function weak_form_kernel!(du, u, + element, ::Type{<:TreeMesh{2}}, + have_nonconservative_terms::False, + have_aux_node_vars::True, equations, + dg::DGSEM, cache, alpha = true) + # true * [some floating point value] == [exactly the same floating point value] + # This can (hopefully) be optimized away due to constant propagation. + @unpack derivative_hat = dg.basis + @unpack aux_node_vars = cache.aux_vars + + # Calculate volume terms in one element + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + aux_node = get_aux_node_vars(aux_node_vars, + equations, dg, i, j, element) + + flux1 = flux(u_node, aux_node, 1, equations) + for ii in eachnode(dg) + multiply_add_to_node_vars!(du, alpha * derivative_hat[ii, i], flux1, + equations, dg, ii, j, element) + end + + flux2 = flux(u_node, aux_node, 2, equations) + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, alpha * derivative_hat[jj, j], flux2, + equations, dg, i, jj, element) + end + end + + return nothing +end + @inline function flux_differencing_kernel!(du, u, element, ::Type{<:TreeMesh{2}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -258,7 +296,8 @@ end end @inline function flux_differencing_kernel!(du, u, element, MeshT::Type{<:TreeMesh{2}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -266,8 +305,9 @@ end symmetric_flux, nonconservative_flux = volume_flux # Apply the symmetric flux as usual - flux_differencing_kernel!(du, u, element, MeshT, False(), equations, symmetric_flux, - dg, cache, alpha) + flux_differencing_kernel!(du, u, element, MeshT, False(), have_aux_node_vars, + equations, + symmetric_flux, dg, cache, alpha) # Calculate the remaining volume terms using the nonsymmetric generalized flux for j in eachnode(dg), i in eachnode(dg) @@ -303,7 +343,7 @@ end MeshT::Type{<:Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg::DGSEM, cache, element, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons, @@ -407,7 +447,7 @@ end MeshT::Type{<:Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes @@ -418,8 +458,8 @@ end fstar1_R = fstar1_R_threaded[Threads.threadid()] fstar2_R = fstar2_R_threaded[Threads.threadid()] calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, MeshT, - have_nonconservative_terms, equations, volume_flux_fv, dg, element, - cache) + have_nonconservative_terms, have_aux_node_vars, equations, + volume_flux_fv, dg, element, cache) # Calculate FV volume integral contribution for j in eachnode(dg), i in eachnode(dg) @@ -441,7 +481,8 @@ end # [arXiv: 2008.12044v2](https://arxiv.org/pdf/2008.12044) @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, ::Type{<:TreeMesh{2}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) for j in eachnode(dg), i in 2:nnodes(dg) u_ll = get_node_vars(u, equations, dg, i - 1, j, element) @@ -464,7 +505,8 @@ end @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, ::Type{<:TreeMesh{2}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) volume_flux, nonconservative_flux = volume_flux_fv @@ -598,7 +640,8 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMesh{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -632,7 +675,47 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMesh{2}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::False, + have_aux_node_vars::True, equations, + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral + @unpack u, neighbor_ids, orientations = cache.interfaces + @unpack aux_surface_node_vars = cache.aux_vars + + @threaded for interface in eachinterface(dg, cache) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface) + aux_ll, aux_rr = get_aux_surface_node_vars(aux_surface_node_vars, + equations, dg, i, + interface) + flux = surface_flux(u_ll, u_rr, aux_ll, aux_rr, + orientations[interface], equations) + # Copy flux to left and right element storage + for v in eachvariable(equations) + surface_flux_values[v, i, left_direction, left_id] = flux[v] + surface_flux_values[v, i, right_direction, right_id] = flux[v] + end + end + end + + return nothing +end + +function calc_interface_flux!(backend::Nothing, surface_flux_values, + mesh::TreeMesh{2}, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces @@ -804,18 +887,22 @@ function calc_boundary_flux!(backend::Nothing, cache, t, # Calc boundary fluxes in each direction calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1], have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, surface_integral, dg, cache, 1, firsts[1], lasts[1]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2], have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, surface_integral, dg, cache, 2, firsts[2], lasts[2]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3], have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, surface_integral, dg, cache, 3, firsts[3], lasts[3]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4], have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, surface_integral, dg, cache, 4, firsts[4], lasts[4]) @@ -825,7 +912,8 @@ end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any, 4}, t, boundary_condition, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache, direction, first_boundary, last_boundary) @unpack surface_flux = surface_integral @@ -861,7 +949,50 @@ end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any, 4}, t, boundary_condition, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::False, + have_aux_node_vars::True, equations, + surface_integral, dg::DG, cache, + direction, first_boundary, last_boundary) + @unpack surface_flux = surface_integral + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + @unpack aux_boundary_node_vars = cache.aux_vars + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, boundary) + aux_ll, aux_rr = get_aux_surface_node_vars(aux_boundary_node_vars, + equations, dg, i, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + u_inner = u_ll + aux_inner = aux_ll + else # Element is on the right, boundary on the left + u_inner = u_rr + aux_inner = aux_rr + end + + x = get_node_coords(node_coordinates, equations, dg, i, boundary) + flux = boundary_condition(u_inner, aux_inner, orientations[boundary], + direction, x, t, surface_flux, equations) + + # Copy flux to left and right element storage + for v in eachvariable(equations) + surface_flux_values[v, i, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + +function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any, 4}, + t, + boundary_condition, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache, direction, first_boundary, last_boundary) @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @@ -997,7 +1128,8 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @@ -1014,13 +1146,55 @@ function calc_mortar_flux!(surface_flux_values, # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, mortar, + calc_fstar!(fstar_primary_upper, equations, + surface_flux, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_primary_lower, equations, + surface_flux, dg, u_lower, mortar, orientation) + calc_fstar!(fstar_secondary_upper, equations, + surface_flux, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_secondary_lower, equations, + surface_flux, dg, u_lower, mortar, orientation) + + mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar_primary_upper, fstar_primary_lower, + fstar_secondary_upper, fstar_secondary_lower) + end + return nothing +end + +function calc_mortar_flux!(surface_flux_values, + mesh::TreeMesh{2}, + have_nonconservative_terms::False, + have_aux_node_vars::True, equations, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral + @unpack u_lower, u_upper, orientations = cache.mortars + @unpack (fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded) = cache + @unpack aux_mortar_node_vars = cache.aux_vars + + @threaded for mortar in eachmortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar_primary_upper = fstar_primary_upper_threaded[Threads.threadid()] + fstar_primary_lower = fstar_primary_lower_threaded[Threads.threadid()] + fstar_secondary_upper = fstar_secondary_upper_threaded[Threads.threadid()] + fstar_secondary_lower = fstar_secondary_lower_threaded[Threads.threadid()] + + # Calculate fluxes + orientation = orientations[mortar] + calc_fstar!(fstar_primary_upper, equations, + surface_flux, dg, u_upper, aux_mortar_node_vars, 2, mortar, orientation) - calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, mortar, + calc_fstar!(fstar_primary_lower, equations, + surface_flux, dg, u_lower, aux_mortar_node_vars, 1, mortar, orientation) - calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, mortar, + calc_fstar!(fstar_secondary_upper, equations, + surface_flux, dg, u_upper, aux_mortar_node_vars, 2, mortar, orientation) - calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, mortar, + calc_fstar!(fstar_secondary_lower, equations, + surface_flux, dg, u_lower, aux_mortar_node_vars, 1, mortar, orientation) mortar_fluxes_to_elements!(surface_flux_values, @@ -1028,13 +1202,13 @@ function calc_mortar_flux!(surface_flux_values, mortar, fstar_primary_upper, fstar_primary_lower, fstar_secondary_upper, fstar_secondary_lower) end - return nothing end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @@ -1051,14 +1225,14 @@ function calc_mortar_flux!(surface_flux_values, # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, mortar, - orientation) - calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, mortar, - orientation) - calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, mortar, - orientation) - calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, mortar, - orientation) + calc_fstar!(fstar_primary_upper, equations, + surface_flux, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_primary_lower, equations, + surface_flux, dg, u_lower, mortar, orientation) + calc_fstar!(fstar_secondary_upper, equations, + surface_flux, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_secondary_lower, equations, + surface_flux, dg, u_lower, mortar, orientation) # Add nonconservative fluxes. # These need to be adapted on the geometry (left/right) since the order of @@ -1134,7 +1308,6 @@ function calc_mortar_flux!(surface_flux_values, mortar, fstar_primary_upper, fstar_primary_lower, fstar_secondary_upper, fstar_secondary_lower) end - return nothing end @@ -1163,6 +1336,24 @@ end return nothing end +@inline function calc_fstar!(destination::AbstractArray{<:Any, 2}, + equations, surface_flux, dg::DGSEM, + u, aux, position, interface, orientation) + for i in eachnode(dg) + # Call pointwise two-point numerical flux function + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface) + aux_ll, _ = get_aux_surface_node_vars(aux, equations, dg, + position, i, interface) + # only leftright = 1 is used + flux = surface_flux(u_ll, u_rr, aux_ll, aux_ll, orientation, equations) + + # Copy flux to left and right element storage + set_node_vars!(destination, flux, equations, dg, i) + end + + return nothing +end + @inline function mortar_fluxes_to_elements!(surface_flux_values, mesh::TreeMesh{2}, equations, mortar_l2::LobattoLegendreMortarL2, @@ -1361,12 +1552,20 @@ end # Need dimension specific version to avoid error at dispatching function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, - equations::AbstractEquations{2}, dg::DG, cache) + have_aux_node_vars::False, equations::AbstractEquations{2}, + dg::DG, cache) + return nothing +end + +function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, + have_aux_node_vars::True, equations::AbstractEquations{2}, + dg::DG, cache) return nothing end function calc_sources!(backend::Nothing, du, u, t, source_terms, - equations::AbstractEquations{2}, dg::DG, cache) + have_aux_node_vars::False, equations::AbstractEquations{2}, + dg::DG, cache) @unpack node_coordinates = cache.elements @threaded for element in eachelement(dg, cache) @@ -1381,4 +1580,25 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms, return nothing end + +function calc_sources!(backend::Nothing, du, u, t, source_terms, + have_aux_node_vars::True, equations::AbstractEquations{2}, + dg::DG, cache) + @unpack node_coordinates = cache.elements + @unpack aux_node_vars = cache.aux_vars + + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + aux_local = get_aux_node_vars(aux_node_vars, equations, dg, + i, j, element) + x_local = get_node_coords(node_coordinates, equations, dg, + i, j, element) + du_local = source_terms(u_local, aux_local, x_local, t, equations) + add_to_node_vars!(du, du_local, equations, dg, i, j, element) + end + end + + return nothing +end end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d_compressible_euler.jl b/src/solvers/dgsem_tree/dg_2d_compressible_euler.jl index 507b48b20ea..40d79803eb1 100644 --- a/src/solvers/dgsem_tree/dg_2d_compressible_euler.jl +++ b/src/solvers/dgsem_tree/dg_2d_compressible_euler.jl @@ -67,6 +67,7 @@ end # muladd @inline function flux_differencing_kernel!(_du::PtrArray, u_cons::PtrArray, element, MeshT::Type{<:TreeMesh{2}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations2D, volume_flux::typeof(flux_shima_etal_turbo), dg::DGSEM, cache, alpha) @@ -229,6 +230,7 @@ end @inline function flux_differencing_kernel!(_du::PtrArray, u_cons::PtrArray, element, MeshT::Type{<:TreeMesh{2}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations2D, volume_flux::typeof(flux_ranocha_turbo), dg::DGSEM, cache, alpha) diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 57ed0033804..ebfc0467e4f 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -482,7 +482,8 @@ function rhs!(du, u, t, # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.volume_integral, dg, cache) end @@ -495,7 +496,8 @@ function rhs!(du, u, t, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.surface_integral, dg, cache) end @@ -519,7 +521,8 @@ function rhs!(du, u, t, # Calculate mortar fluxes @trixi_timeit timer() "mortar flux" begin calc_mortar_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end @@ -531,14 +534,16 @@ function rhs!(du, u, t, # Calculate MPI interface fluxes @trixi_timeit timer() "MPI interface flux" begin calc_mpi_interface_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.surface_integral, dg, cache) end # Calculate MPI mortar fluxes @trixi_timeit timer() "MPI mortar flux" begin calc_mpi_mortar_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, + have_nonconservative_terms(equations), + have_aux_node_vars(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end @@ -554,7 +559,8 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, + have_aux_node_vars(equations), equations, dg, cache) end # Finish to send MPI data @@ -727,7 +733,8 @@ end function calc_mpi_interface_flux!(surface_flux_values, mesh::TreeMeshParallel{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u, local_neighbor_ids, orientations, remote_sides = cache.mpi_interfaces @@ -766,14 +773,101 @@ function calc_mpi_interface_flux!(surface_flux_values, return nothing end +function calc_mpi_interface_flux!(surface_flux_values, + mesh::TreeMeshParallel{2}, + have_nonconservative_terms::False, + have_aux_node_vars::True, equations, + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral + @unpack u, local_neighbor_ids, orientations, remote_sides = cache.mpi_interfaces + @unpack aux_mpiinterface_node_vars = cache.aux_vars + + @threaded for interface in eachmpiinterface(dg, cache) + # Get local neighboring element + element = local_neighbor_ids[interface] + + # Determine interface direction with respect to element: + if orientations[interface] == 1 # interface in x-direction + if remote_sides[interface] == 1 # local element in positive direction + direction = 1 + else # local element in negative direction + direction = 2 + end + else # interface in y-direction + if remote_sides[interface] == 1 # local element in positive direction + direction = 3 + else # local element in negative direction + direction = 4 + end + end + + for i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface) + aux_ll, aux_rr = get_aux_surface_node_vars(aux_mpiinterface_node_vars, + equations, dg, i, + interface) + flux = surface_flux(u_ll, u_rr, aux_ll, aux_rr, + orientations[interface], equations) + + # Copy flux to local element storage + for v in eachvariable(equations) + surface_flux_values[v, i, direction, element] = flux[v] + end + end + end + + return nothing +end + +function calc_mpi_mortar_flux!(surface_flux_values, + mesh::TreeMeshParallel{2}, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral + @unpack u_lower, u_upper, orientations = cache.mpi_mortars + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded, fstar_secondary_upper_threaded, fstar_secondary_lower_threaded = cache + + @threaded for mortar in eachmpimortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar_primary_upper = fstar_primary_upper_threaded[Threads.threadid()] + fstar_primary_lower = fstar_primary_lower_threaded[Threads.threadid()] + fstar_secondary_upper = fstar_secondary_upper_threaded[Threads.threadid()] + fstar_secondary_lower = fstar_secondary_lower_threaded[Threads.threadid()] + + # Because `have_nonconservative_terms` is `False` the primary and secondary fluxes + # are identical. So, we could possibly save on computation and just pass two copies later. + orientation = orientations[mortar] + calc_fstar!(fstar_primary_upper, equations, + surface_flux, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_primary_lower, equations, + surface_flux, dg, u_lower, mortar, orientation) + calc_fstar!(fstar_secondary_upper, equations, + surface_flux, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_secondary_lower, equations, + surface_flux, dg, u_lower, mortar, orientation) + + mpi_mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar_primary_upper, fstar_primary_lower, + fstar_secondary_upper, fstar_secondary_lower) + end + + return nothing +end + function calc_mpi_mortar_flux!(surface_flux_values, mesh::TreeMeshParallel{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::True, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u_lower, u_upper, orientations = cache.mpi_mortars @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded, fstar_secondary_upper_threaded, fstar_secondary_lower_threaded = cache + @unpack aux_mpimortar_node_vars = cache.aux_vars @threaded for mortar in eachmpimortar(dg, cache) # Choose thread-specific pre-allocated container @@ -785,14 +879,14 @@ function calc_mpi_mortar_flux!(surface_flux_values, # Because `have_nonconservative_terms` is `False` the primary and secondary fluxes # are identical. So, we could possibly save on computation and just pass two copies later. orientation = orientations[mortar] - calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, mortar, - orientation) - calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, mortar, - orientation) - calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, mortar, - orientation) - calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, mortar, - orientation) + calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, + aux_mpimortar_node_vars, 2, mortar, orientation) + calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, + aux_mpimortar_node_vars, 1, mortar, orientation) + calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, + aux_mpimortar_node_vars, 2, mortar, orientation) + calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, + aux_mpimortar_node_vars, 1, mortar, orientation) mpi_mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, @@ -859,7 +953,6 @@ end direction = 3 end end - multiply_dimensionwise!(view(surface_flux_values, :, :, direction, element), mortar_l2.reverse_upper, fstar_secondary_upper, mortar_l2.reverse_lower, fstar_secondary_lower) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index c497a9027f7..5531f5ee320 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -65,7 +65,8 @@ end MeshT::Type{<:Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes @@ -90,8 +91,8 @@ end fstar1_R = fstar1_R_threaded[Threads.threadid()] fstar2_R = fstar2_R_threaded[Threads.threadid()] calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, MeshT, - have_nonconservative_terms, equations, volume_flux_fv, dg, element, - cache) + have_nonconservative_terms, have_aux_node_vars, equations, + volume_flux_fv, dg, element, cache) # antidiffusive flux calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 75970fbc3dd..b1744eb44d9 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -132,7 +132,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 =# @inline function weak_form_kernel!(du, u, element, ::Type{<:TreeMesh{3}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -164,7 +165,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 end @inline function flux_differencing_kernel!(du, u, element, ::Type{<:TreeMesh{3}}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -214,7 +216,8 @@ end end @inline function flux_differencing_kernel!(du, u, element, MeshT::Type{<:TreeMesh{3}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @@ -222,8 +225,9 @@ end symmetric_flux, nonconservative_flux = volume_flux # Apply the symmetric flux as usual - flux_differencing_kernel!(du, u, element, MeshT, False(), equations, symmetric_flux, - dg, cache, alpha) + flux_differencing_kernel!(du, u, element, MeshT, False(), have_aux_node_vars, + equations, + symmetric_flux, dg, cache, alpha) # Calculate the remaining volume terms using the nonsymmetric generalized flux for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) @@ -269,7 +273,7 @@ end MeshT::Type{<:Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded = cache @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes @@ -283,7 +287,7 @@ end fstar3_R = fstar3_R_threaded[Threads.threadid()] calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, - MeshT, have_nonconservative_terms, equations, + MeshT, have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg, element, cache) # Calculate FV volume integral contribution @@ -309,7 +313,7 @@ end MeshT::Type{<:Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}}, - have_nonconservative_terms, equations, + have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg::DGSEM, cache, element, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons, @@ -328,7 +332,7 @@ end fstar3_R = fstar3_R_threaded[Threads.threadid()] calcflux_fvO2!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, - MeshT, have_nonconservative_terms, equations, + MeshT, have_nonconservative_terms, have_aux_node_vars, equations, volume_flux_fv, dg, element, cache, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons) @@ -359,7 +363,7 @@ end @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, ::Type{<:TreeMesh{3}}, have_nonconservative_terms::False, - equations, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg) u_ll = get_node_vars(u, equations, dg, i - 1, j, k, element) @@ -391,7 +395,8 @@ end @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, ::Type{<:TreeMesh{3}}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache) volume_flux, nonconservative_flux = volume_flux_fv @@ -456,7 +461,7 @@ end fstar3_L, fstar3_R, u, ::Type{<:TreeMesh{3}}, have_nonconservative_terms::False, - equations, + have_aux_node_vars::False, equations, volume_flux_fv, dg::DGSEM, element, cache, sc_interface_coords, reconstruction_mode, slope_limiter, cons2recon, recon2cons) @@ -568,7 +573,8 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMesh{3}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -603,7 +609,8 @@ end function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMesh{3}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces @@ -959,7 +966,8 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{3}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @@ -1017,7 +1025,8 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{3}, - have_nonconservative_terms::True, equations, + have_nonconservative_terms::True, + have_aux_node_vars::False, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @@ -1424,12 +1433,14 @@ end # Need dimension specific version to avoid error at dispatching function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, - equations::AbstractEquations{3}, dg::DG, cache) + have_aux_node_vars::False, equations::AbstractEquations{3}, + dg::DG, cache) return nothing end function calc_sources!(backend::Nothing, du, u, t, source_terms, - equations::AbstractEquations{3}, dg::DG, cache) + have_aux_node_vars::False, equations::AbstractEquations{3}, + dg::DG, cache) @unpack node_coordinates = cache.elements @threaded for element in eachelement(dg, cache) diff --git a/src/solvers/dgsem_tree/dg_3d_compressible_euler.jl b/src/solvers/dgsem_tree/dg_3d_compressible_euler.jl index 1cdf1ca07e6..045053f7856 100644 --- a/src/solvers/dgsem_tree/dg_3d_compressible_euler.jl +++ b/src/solvers/dgsem_tree/dg_3d_compressible_euler.jl @@ -19,6 +19,7 @@ @inline function flux_differencing_kernel!(_du::PtrArray, u_cons::PtrArray, element, MeshT::Type{<:TreeMesh{3}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations3D, volume_flux::typeof(flux_shima_etal_turbo), dg::DGSEM, cache, alpha) @@ -265,6 +266,7 @@ end @inline function flux_differencing_kernel!(_du::PtrArray, u_cons::PtrArray, element, MeshT::Type{<:TreeMesh{3}}, have_nonconservative_terms::False, + have_aux_node_vars::False, equations::CompressibleEulerEquations3D, volume_flux::typeof(flux_ranocha_turbo), dg::DGSEM, cache, alpha) diff --git a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl index 27e5d13049e..95f739119e1 100644 --- a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl @@ -64,7 +64,8 @@ end # Subcell limiting currently only implemented for certain mesh types @inline function volume_integral_kernel!(du, u, element, MeshT::Type{<:Union{TreeMesh{3}, P4estMesh{3}}}, - nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @unpack inverse_weights = dg.basis # Plays role of DG subcell sizes @@ -80,7 +81,7 @@ end fhat3_L = fhat3_L_threaded[Threads.threadid()] fhat3_R = fhat3_R_threaded[Threads.threadid()] calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, - u, MeshT, nonconservative_terms, equations, volume_flux_dg, + u, MeshT, have_nonconservative_terms, equations, volume_flux_dg, dg, element, cache) # low-order FV fluxes @@ -93,13 +94,13 @@ end fstar3_L = fstar3_L_threaded[Threads.threadid()] fstar3_R = fstar3_R_threaded[Threads.threadid()] calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, - u, MeshT, nonconservative_terms, equations, volume_flux_fv, - dg, element, cache) + u, MeshT, have_nonconservative_terms, have_aux_node_vars, equations, + volume_flux_fv, dg, element, cache) # antidiffusive flux calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, - u, MeshT, nonconservative_terms, equations, limiter, + u, MeshT, have_nonconservative_terms, equations, limiter, dg, element, cache) # Calculate volume integral contribution of low-order FV flux diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index ecf9d6b7a48..5da29069b74 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -26,8 +26,10 @@ end # Otherwise, @threaded does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. @inline function calc_indicator_hennemann_gassner!(indicator_hg, threshold, parameter_s, - u, element, mesh::AbstractMesh{1}, - equations, dg, cache) + u, + element, mesh::AbstractMesh{1}, + have_aux_node_vars, equations, dg, + cache) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index c3ac38870e8..37a8d97c288 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -28,7 +28,9 @@ end # Otherwise, @threaded does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. @inline function calc_indicator_hennemann_gassner!(indicator_hg, threshold, parameter_s, - u, element, mesh::AbstractMesh{2}, + u, + element, mesh::AbstractMesh{2}, + have_aux_node_vars, equations, dg, cache) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, @@ -39,10 +41,8 @@ end modal_tmp1 = modal_tmp1_threaded[Threads.threadid()] # Calculate indicator variables at Gauss-Lobatto nodes - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - indicator[i, j] = indicator_hg.variable(u_local, equations) - end + calc_indicator_inner!(indicator, u, element, mesh, indicator_hg.variable, + have_aux_node_vars, equations, dg, cache) # Convert to modal representation multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, @@ -100,6 +100,16 @@ end return nothing end +@inline function calc_indicator_inner!(indicator, u, element, mesh::AbstractMesh{2}, + indicator_variable, + have_aux_node_vars::False, equations, + solver, cache) + for j in eachnode(solver), i in eachnode(solver) + u_local = get_node_vars(u, equations, solver, i, j, element) + indicator[i, j] = indicator_variable(u_local, equations) + end +end + # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, alpha, alpha_tmp, dg, diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index ff670dfc469..e7fe7601ed3 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -37,8 +37,10 @@ end # Otherwise, @threaded does not work here with Julia ARM on macOS. # See https://github.com/JuliaSIMD/Polyester.jl/issues/88. @inline function calc_indicator_hennemann_gassner!(indicator_hg, threshold, parameter_s, - u, element, mesh::AbstractMesh{3}, - equations, dg, cache) + u, + element, mesh::AbstractMesh{3}, + have_aux_node_vars, equations, dg, + cache) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded, modal_tmp2_threaded = indicator_hg.cache diff --git a/src/solvers/fdsbp_tree/fdsbp.jl b/src/solvers/fdsbp_tree/fdsbp.jl index acf38425b04..4c78f090ee4 100644 --- a/src/solvers/fdsbp_tree/fdsbp.jl +++ b/src/solvers/fdsbp_tree/fdsbp.jl @@ -56,7 +56,8 @@ function prolong2mortars!(cache, u, mesh, equations, mortar::Nothing, end function calc_mortar_flux!(surface_flux_values, mesh, - have_nonconservative_terms, equations, + have_nonconservative_terms, + have_aux_node_vars, equations, mortar::Nothing, surface_integral, dg::DG, cache) @assert isempty(eachmortar(dg, cache)) diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 002093121a7..6b7fb37ac63 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,7 +44,8 @@ end # 2D volume integral contributions for `VolumeIntegralStrongForm` function calc_volume_integral!(backend::Nothing, du, u, mesh::TreeMesh{1}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralStrongForm, dg::FDSBP, cache) D = dg.basis # SBP derivative operator @@ -91,7 +92,8 @@ end # of the flux splitting f^-. function calc_volume_integral!(backend::Nothing, du, u, mesh::TreeMesh{1}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralUpwind, dg::FDSBP, cache) # Assume that diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl index 5d30157a28c..b9ff9fb776d 100644 --- a/src/solvers/fdsbp_tree/fdsbp_2d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_2d.jl @@ -44,7 +44,8 @@ end # 2D volume integral contributions for `VolumeIntegralStrongForm` function calc_volume_integral!(backend::Nothing, du, u, mesh::TreeMesh{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralStrongForm, dg::FDSBP, cache) D = dg.basis # SBP derivative operator @@ -100,7 +101,8 @@ end # of the flux splitting f^-. function calc_volume_integral!(backend::Nothing, du, u, mesh::TreeMesh{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralUpwind, dg::FDSBP, cache) # Assume that @@ -218,7 +220,8 @@ end # flux information at each side of an interface. function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMesh{2}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral::SurfaceIntegralUpwind, dg::FDSBP, cache) @unpack splitting = surface_integral diff --git a/src/solvers/fdsbp_tree/fdsbp_3d.jl b/src/solvers/fdsbp_tree/fdsbp_3d.jl index e737a8f2137..ec6444e3a89 100644 --- a/src/solvers/fdsbp_tree/fdsbp_3d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_3d.jl @@ -44,7 +44,8 @@ end # 3D volume integral contributions for `VolumeIntegralStrongForm` function calc_volume_integral!(backend::Nothing, du, u, mesh::TreeMesh{3}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralStrongForm, dg::FDSBP, cache) D = dg.basis # SBP derivative operator @@ -107,7 +108,8 @@ end # of the flux splitting f^-. function calc_volume_integral!(backend::Nothing, du, u, mesh::TreeMesh{3}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralUpwind, dg::FDSBP, cache) # Assume that @@ -254,7 +256,8 @@ end # flux information at each side of an interface. function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMesh{3}, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, surface_integral::SurfaceIntegralUpwind, dg::FDSBP, cache) @unpack splitting = surface_integral diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index 2d7058b9957..7fe80ef373d 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -31,7 +31,8 @@ end # So it is not provably stable for variable coefficients due to the the metric terms. function calc_volume_integral!(backend::Nothing, du, u, mesh::UnstructuredMesh2D, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralStrongForm, dg::FDSBP, cache) D = dg.basis # SBP derivative operator @@ -94,7 +95,8 @@ end # of the flux splitting f^-. function calc_volume_integral!(backend::Nothing, du, u, mesh::UnstructuredMesh2D, - have_nonconservative_terms::False, equations, + have_nonconservative_terms::False, + have_aux_node_vars::False, equations, volume_integral::VolumeIntegralUpwind, dg::FDSBP, cache) # Assume that diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl index 4dfc44b2089..c66ff238d4c 100644 --- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl +++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl @@ -60,8 +60,9 @@ function calculate_cfl(ode_algorithm::AbstractPairedExplicitRK, ode) u = wrap_array(u_ode, mesh, equations, solver, cache) cfl_number = dt_opt / max_dt(u, t0, mesh, - have_constant_speed(equations), equations, - solver, cache) + have_constant_speed(equations), + have_aux_node_vars(equations), + equations, solver, cache) return cfl_number end diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl index 502afdca901..b1c8c9ac43f 100644 --- a/src/visualization/utilities.jl +++ b/src/visualization/utilities.jl @@ -483,6 +483,23 @@ function get_data_1d(original_nodes, unstructured_data, nvisnodes, reinterpolate vcat(original_nodes[1, 1, :], original_nodes[1, end, end]) end +# Apply the `solution_variables` function to all node values stored in `u`. +@inline function apply_solution_variables(u, solution_variables, + have_aux_node_vars::False, equations, + solver, cache) + n_vars_in = nvariables(equations) + n_vars = length(solution_variables(get_node_vars(u, equations, solver), + equations)) + raw_data = Array{eltype(u)}(undef, n_vars, Base.tail(size(u))...) + reshaped_u = reshape(u, n_vars_in, :) + reshaped_r = reshape(raw_data, n_vars, :) + for idx in axes(reshaped_u, 2) + u_node = get_node_vars(reshaped_u, equations, solver, idx) + reshaped_r[:, idx] = solution_variables(u_node, equations) + end + return raw_data +end + # Change order of dimensions (variables are now last) and convert data to `solution_variables` # # Note: This is a low-level function that is not considered as part of Trixi.jl's interface and may @@ -492,22 +509,11 @@ function get_unstructured_data(u, solution_variables, mesh, equations, solver, c raw_data = u n_vars = size(raw_data, 1) else - # Similar to `save_solution_file` in `callbacks_step/save_solution_dg.jl`. - # However, we cannot use `reinterpret` here as `u` might have a non-bits type. - # See https://github.com/trixi-framework/Trixi.jl/pull/2388 - n_vars_in = nvariables(equations) - n_vars = length(solution_variables(get_node_vars(u, equations, solver), - equations)) - raw_data = Array{eltype(u)}(undef, n_vars, Base.tail(size(u))...) - reshaped_u = reshape(u, n_vars_in, :) - reshaped_r = reshape(raw_data, n_vars, :) - for idx in axes(reshaped_u, 2) - reshaped_r[:, idx] = solution_variables(get_node_vars(reshaped_u, equations, - solver, idx), - equations) - end + raw_data = apply_solution_variables(u, solution_variables, + have_aux_node_vars(equations), + equations, solver, cache) + n_vars = size(raw_data, 1) end - unstructured_data = Array{eltype(raw_data)}(undef, ntuple((d) -> nnodes(solver), ndims(equations))..., diff --git a/test/test_mpi_tree.jl b/test/test_mpi_tree.jl index c6df01b2bfb..bad367a616d 100644 --- a/test/test_mpi_tree.jl +++ b/test/test_mpi_tree.jl @@ -106,6 +106,20 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() linf=[0.0253454486893413]) end + @trixi_testset "elixir_advection_variable_swirling_flow.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_variable_swirling_flow.jl"), + l2=[2.90963554e-01], + linf=[1.31858729e+00]) + end + + @trixi_testset "elixir_acoustics_gauss_wall_amr_auxvars.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_acoustics_gauss_wall_amr_auxvars.jl"), + l2=[0.0194350488, 0.0195225440, 0.04819822538525215], + linf=[0.183057645, 0.190845961, 1.03618809]) + end + # Hyperbolic diffusion if !CI_ON_WINDOWS # see comment on `CI_ON_WINDOWS` in `test/test_mpi.jl` @trixi_testset "elixir_hypdiff_lax_friedrichs.jl" begin diff --git a/test/test_performance_specializations_2d.jl b/test/test_performance_specializations_2d.jl index 7dceea2b6a7..3bf00c3e104 100644 --- a/test/test_performance_specializations_2d.jl +++ b/test/test_performance_specializations_2d.jl @@ -30,11 +30,13 @@ isdir(outdir) && rm(outdir, recursive = true) u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, 1] @@ -44,10 +46,11 @@ isdir(outdir) && rm(outdir, recursive = true) du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, 1] @@ -69,11 +72,13 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, 1] @@ -83,10 +88,11 @@ end du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, 1] @@ -109,11 +115,13 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, 1] @@ -123,10 +131,11 @@ end du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, 1] @@ -148,11 +157,13 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, 1] @@ -162,10 +173,11 @@ end du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, 1] diff --git a/test/test_performance_specializations_3d.jl b/test/test_performance_specializations_3d.jl index 967b0f9cf3e..fefa703d13c 100644 --- a/test/test_performance_specializations_3d.jl +++ b/test/test_performance_specializations_3d.jl @@ -30,11 +30,13 @@ isdir(outdir) && rm(outdir, recursive = true) u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, :, 1] @@ -44,10 +46,11 @@ isdir(outdir) && rm(outdir, recursive = true) du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), - Function, typeof(semi.solver), typeof(semi.cache), Bool}, + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), + typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, :, 1] @@ -69,11 +72,13 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, :, 1] @@ -83,10 +88,11 @@ end du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), - Function, typeof(semi.solver), typeof(semi.cache), Bool}, + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), + typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, :, 1] @@ -109,11 +115,13 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, :, 1] @@ -123,10 +131,11 @@ end du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), - Function, typeof(semi.solver), typeof(semi.cache), Bool}, + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), + typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, :, 1] @@ -148,11 +157,13 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations) + have_aux_node_vars = Trixi.have_aux_node_vars(semi.equations) # Call the optimized default version du .= 0 Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, + semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_specialized = du[:, :, :, :, 1] @@ -162,10 +173,11 @@ end du .= 0 invoke(Trixi.flux_differencing_kernel!, Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)}, - typeof(have_nonconservative_terms), typeof(semi.equations), - Function, typeof(semi.solver), typeof(semi.cache), Bool}, + typeof(have_nonconservative_terms), typeof(have_aux_node_vars), + typeof(semi.equations), Function, typeof(semi.solver), + typeof(semi.cache), Bool}, du, u, 1, typeof(semi.mesh), - have_nonconservative_terms, semi.equations, + have_nonconservative_terms, have_aux_node_vars, semi.equations, semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true) du_baseline = du[:, :, :, :, 1] diff --git a/test/test_tree_2d_acoustics.jl b/test/test_tree_2d_acoustics.jl index 49e183fdc7a..40f46f7caf1 100644 --- a/test/test_tree_2d_acoustics.jl +++ b/test/test_tree_2d_acoustics.jl @@ -35,6 +35,24 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_acoustics_convergence_auxvars.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_acoustics_convergence_auxvars.jl"), + l2=[ + 0.0019921138796370834, + 0.002090394698052287, + 0.0006091925854593805 + ], + linf=[ + 0.00769282588065634, + 0.008276649669227254, + 0.004196479023954813 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_acoustics_gauss.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss.jl"), l2=[ @@ -108,6 +126,28 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_acoustics_gauss_wall_auxvars.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_acoustics_gauss_wall_auxvars.jl"), + l2=[0.019419398248465843, 0.019510701017551826, + 0.04818246051887614], + linf=[0.18193631937316496, 0.1877464607867628, + 1.0355388011792845]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + +@trixi_testset "elixir_acoustics_gauss_wall_amr_auxvars.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_acoustics_gauss_wall_amr_auxvars.jl"), + l2=[0.0194350488, 0.0195225440, 0.04819822538525215], + linf=[0.183057645, 0.190845961, 1.03618809]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_acoustics_gauss_wall.jl (Gauss Legendre)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss_wall.jl"), solver=DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs, diff --git a/test/test_tree_2d_advection.jl b/test/test_tree_2d_advection.jl index 53f3aefe71b..9f02df5dc19 100644 --- a/test/test_tree_2d_advection.jl +++ b/test/test_tree_2d_advection.jl @@ -247,6 +247,16 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_advection_variable_swirling_flow.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_variable_swirling_flow.jl"), + l2=[2.90963554e-01], + linf=[1.31858729e+00]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + # Coverage test for all initial conditions @testset "Linear scalar advection: Tests for initial conditions" begin # Linear scalar advection