Skip to content

Add auxiliary variables, part 1#3008

Open
benegee wants to merge 100 commits into
mainfrom
bg/feature-aux-vars-merge
Open

Add auxiliary variables, part 1#3008
benegee wants to merge 100 commits into
mainfrom
bg/feature-aux-vars-merge

Conversation

@benegee
Copy link
Copy Markdown
Contributor

@benegee benegee commented May 12, 2026

Continuing #2348

As a first step, auxiliary variables are implemented for TreeMesh 2D. AcousticPerturbationEquations2DAuxVars demonstrates how auxiliary variables can be used instead of augmenting the solution state vector, as previously done, e.g., in AcousticPerturbationEquations2D.

Some comments:

  • Auxiliary variables could potentially be used everywhere in Trixi.jl, which makes infrastructural changes necessary almost everywhere. It usually boils down to querying have_aux_node_vars(equations) and using the result for multiple dispatch.
  • Specializations for have_aux_node_vars::True are added only for the 2D TreeMesh implementation in this PR. The code changes are trivial, usually just adding get_aux_node_vars or get_aux_surface_node_vars where get_node_vars or get_surface_node_vars were already called, and passing the resulting aux vars down, e.g. to flux computations.
  • A less trivial part is the setup of the aux vars in solvers/dgsem_tree/containers_2d.jl, although it mostly follows corresponding prolong2* functions.
  • A non-trivial part is the setup and usage of MPI interfaces, where the aux vars are directly evaluated on both sides and not communicated. This assumes the aux vars field is steady, never discontinuous.
  • A truly non-trivial part is the setup and use of mortars. So far, aux vars are always evaluated for the small mortar elements and these values are used on both sides of the mortar, unlike the usual variables, which also exist on the large elements and need to be interpolated to the small elements. Similar to the discussion in Add auxiliary variables #2348 (comment), an alternative interpolating version could be added.

Supersedes #2907

benegee added 30 commits April 3, 2025 17:45
- add container for node and surface auxiliary variables
- add initializer to SemidiscretizationHyperbolic
- dispatch on have_auxiliary_node_vars for 2D TreeMesh rhs!
to demonstrate how auxiliary variables can be used instead of augmenting
the solution state vector
ŕesults have to be identical to the corresponing AcousticPerturbationEquations2D results
@benegee benegee marked this pull request as ready for review May 15, 2026 09:07
@benegee
Copy link
Copy Markdown
Contributor Author

benegee commented May 15, 2026

All conflicts resolved. This is now ready for review.

Downstream tests are expected to fail because TrixiShallowWater.jl and TrixiAtmo.jl overwrite internal Trixi.jl functions, e.g. https://github.com/trixi-framework/TrixiShallowWater.jl/blob/14552f2fdd0df83dbc094d2700c9e37dfaeab037/src/solvers/dgsem_p4est/dg_2d.jl#L314 would now need an additional have_aux_node_vars argument.

Copy link
Copy Markdown
Contributor

@MarcoArtiano MarcoArtiano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Benedict! That looks so interesting and great to have in Trixi. I have not made a review, but just left a few comments to understand better the infrastructure. I'm not actually sure if there are some problems or issues that could increase the complexity of the implementation.

###############################################################################
# semidiscretization of the acoustic perturbation equations

function auxiliary_variables_mean_values(x, equations)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it make sense (or be interesting) to compute first the initial condition and then pass also the initial condition u0 as a parameter to the auxiliary field function?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see this would be convenient when e.g. setting up perturbations of a steady background state.

However, u0 is computed using compute_coefficients(t, semi), which takes a Semidiscretization object, right?

Comment thread examples/tree_2d_dgsem/elixir_acoustics_convergence_auxvars.jl
@DanielDoehring DanielDoehring added the enhancement New feature or request label May 17, 2026
@DanielDoehring
Copy link
Copy Markdown
Member

I suggest this splitting into multiple PRs, starting from 1D TreeMesh?

Copy link
Copy Markdown
Member

@ranocha ranocha left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for your efforts. How do you see the chances to allow discontinuities at interfaces?

Comment on lines +8 to +10
# constant auxiliary variables (mean state)
return global_mean_vars(equations)
end
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does it really make sense to introduce this set of equations if the auxiliary variables are only given by some constants that are also stored in the equations? Could it be better to use something that appears to be less artificial, e.g., a 2D linear advection equation with variable velocity field? If it is required to be divergence free, a reasonable discretization can be a split form (but we can test the weak form as well, of course).

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is indeed not really useful; merely a proof of concept that can be compared to something in main as a sanity check.

Non trivial velocity fields would be more exciting. I like the swirling flow by LeVeque:
swirling

Although in this form it requires time dependent aux vars, which I added quick and dirty here

# Recompute aux vars
@trixi_timeit_ext backend timer() "recompute aux vars" begin
recompute_aux_vars!(mesh, have_aux_node_vars(equations), equations, dg, cache,
t)

One could of course start with a time constant variant of it.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice example! Something like this with constant-in-time (but variable-in-space) velocity would already be nice to ensure that the non-constant case works as expected.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done in examples/tree_2d_dgsem/elixir_advection_variable_swirling_flow.jl

Comment thread src/callbacks_step/analysis_dg1d.jl Outdated
Comment thread src/callbacks_step/analysis_dg1d.jl Outdated
Comment on lines 619 to 626
function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
return analyze(quantity, du, u, t, mesh, equations, solver, cache)
return analyze(quantity, du, u, t, mesh, have_aux_node_vars(equations), equations,
solver, cache)
end
function analyze(quantity, du, u, t, mesh, equations, solver, cache)
function analyze(quantity, du, u, t, mesh, have_aux_node_vars, equations, solver, cache)
return integrate(quantity, u, mesh, equations, solver, cache, normalize = true)
end
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we acutally need the dispatch at this high level? It introduces some complexity here where I do not see the benefits of it. Would it be sufficient to extract the trait when it is required, e.g., for analyze(::typeof(entropy_timederivative), ... later?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure.

I need the aux vars here:

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::DG, 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)
dot(cons2entropy(u_node, aux_node, equations), du_node)
end
end

Would I need to specialize integrate_via_indices and the anonymous function (defined by the do syntax). So far, I was not successful.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can't you introduce the additional dispatch only for this method of analyze, e.g., something along the lines of

 function analyze(::typeof(entropy_timederivative), du, u, t, 
                  mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, 
                              UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, 
                  equations, dg::DG, cache)
    # 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

and then the specialization you wrote above?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, obviously... Thanks for the hint!

Comment thread src/callbacks_step/stepsize.jl Outdated
Comment thread src/solvers/dgsem_tree/containers_2d.jl
Comment thread src/solvers/dgsem_tree/dg_1d.jl Outdated
Comment thread src/solvers/dgsem_tree/dg_2d.jl
Comment thread NEWS.md Outdated
Comment thread NEWS.md Outdated
Copy link
Copy Markdown
Contributor Author

@benegee benegee left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks a lot for your feedback! I could easily address most of the comments.

How do you see the chances to allow discontinuities at interfaces?

This is something to discuss!

I think this would be possible as well.

  • We can still define the interface as we see fit. This is the current call:
    aux_field(x_local, t, equations),

    We could additionally pass the element. Or how would we usually "detect" an interface. How do we do it with initial conditions?
  • So far I used the assumption the aux vars be continuous to spare some MPI communication at the interfaces, and to use non-interpolated data at mortar interfaces (comment in containers_2d), but this could be adapted appropriately.

Comment on lines +8 to +10
# constant auxiliary variables (mean state)
return global_mean_vars(equations)
end
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is indeed not really useful; merely a proof of concept that can be compared to something in main as a sanity check.

Non trivial velocity fields would be more exciting. I like the swirling flow by LeVeque:
swirling

Although in this form it requires time dependent aux vars, which I added quick and dirty here

# Recompute aux vars
@trixi_timeit_ext backend timer() "recompute aux vars" begin
recompute_aux_vars!(mesh, have_aux_node_vars(equations), equations, dg, cache,
t)

One could of course start with a time constant variant of it.

Comment on lines 619 to 626
function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
return analyze(quantity, du, u, t, mesh, equations, solver, cache)
return analyze(quantity, du, u, t, mesh, have_aux_node_vars(equations), equations,
solver, cache)
end
function analyze(quantity, du, u, t, mesh, equations, solver, cache)
function analyze(quantity, du, u, t, mesh, have_aux_node_vars, equations, solver, cache)
return integrate(quantity, u, mesh, equations, solver, cache, normalize = true)
end
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure.

I need the aux vars here:

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::DG, 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)
dot(cons2entropy(u_node, aux_node, equations), du_node)
end
end

Would I need to specialize integrate_via_indices and the anonymous function (defined by the do syntax). So far, I was not successful.

Comment thread src/semidiscretization/semidiscretization_hyperbolic.jl Outdated
Comment thread src/solvers/dgsem_tree/dg_2d.jl
Comment thread examples/tree_2d_dgsem/elixir_advection_variable_swirling_flow.jl Outdated
Comment thread src/callbacks_step/analysis_dg2d.jl Outdated
Comment thread src/callbacks_step/analysis_dg2d.jl Outdated
Comment thread src/callbacks_step/analysis_dg2d.jl Outdated
Comment thread src/equations/linear_variable_scalar_advection_2d.jl Outdated
Comment thread src/equations/linear_variable_scalar_advection_2d.jl Outdated
Comment thread src/equations/linear_variable_scalar_advection_2d.jl Outdated
Comment thread src/semidiscretization/semidiscretization_hyperbolic.jl Outdated
Comment thread test/test_mpi_tree.jl Outdated
Comment thread test/test_tree_2d_advection.jl Outdated
Copy link
Copy Markdown
Member

@tristanmontoya tristanmontoya left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The core implementation looks good, although if we can remove the trait and just use zero-length vectors without overhead like you mentioned here, that would probably be cleaner, assuming there's negligible overhead. My suggestions are basically just to add clearer documentation in places (see comments).

Comment on lines 233 to +235
function analyze(::typeof(entropy_timederivative), du, u, t,
mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations,
dg::Union{DGSEM, FDSBP}, cache)
mesh::Union{TreeMesh{1}, StructuredMesh{1}},
equations, dg::Union{DGSEM, FDSBP}, cache)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are some lines in this PR that just seem to be formatting changes. Not a big deal but it makes the diff look a bit bigger than it actually is.

Comment on lines +409 to +410
Trait function determining whether `equations` need to access additional auxiliary
variables.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would give a more complete description of auxiliary variables here, or if you define it elsewhere, link it here. For example:

Suggested change
Trait function determining whether `equations` need to access additional auxiliary
variables.
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:
```
flux(u, aux_vars, orientation_or_normal_direction, equations)
```
The auxiliary variables are computed by passing an additional keyword argument
`aux_field` into [`SemidiscretizationHyperbolic`](@ref), where
`aux_field(x, equations)` returns an `SVector` of size `n_aux_node_vars(equations)`
at every point `x` in the spatial domain.
!!! warning "Experimental implementation"
The use of auxiliary variables is an experimental feature and may change in future
releases. Currently, only 2D `TreeMesh` is supported.

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.
NOTE: This is experimental!
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
NOTE: This is experimental!
!!! warning "Experimental implementation"
The use of auxiliary variables is an experimental feature and may change in future
releases. Currently, only 2D `TreeMesh` is supported.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants