Feature: Mesh-IO and Save Solution support for DGMulti

[jmark]

Title says it all. It is still considered draft. Any comments, advice and contributions are welcome!

The pendant Trixi2Vtk PR: trixi-framework/Trixi2Vtk.jl#103

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[codecov]

Codecov Report

Attention: Patch coverage is 86.66667% with 20 lines in your changes missing coverage. Please review.

Project coverage is 96.77%. Comparing base (55d44d9) to head (2cbcc82).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/meshes/dgmulti_meshes.jl	65.38%	9 Missing ⚠️
src/callbacks_step/save_solution_dg.jl	84.00%	8 Missing ⚠️
src/meshes/mesh_io.jl	96.43%	2 Missing ⚠️
src/callbacks_step/save_solution.jl	88.89%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2259      +/-   ##
==========================================
- Coverage   96.80%   96.77%   -0.04%     
==========================================
  Files         495      495              
  Lines       40917    41055     +138     
==========================================
+ Hits        39608    39727     +119     
- Misses       1309     1328      +19     

Flag	Coverage Δ
unittests	96.77% <86.67%> (-0.04%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jlchan]

Thanks for putting this together @jmark! I'm still going through the related Trixi2Vtk PR. Is it correct to say that the main structural changes here are:

1. adding dg.basis to the mesh type,
2. making DGMultiMesh a mutable struct?

Some questions:

• why is making DGMultiMesh a mutable struct necessary?
• I saw that RefElemData is used to interpolate data in https://github.com/trixi-framework/Trixi2Vtk.jl/pull/103/files#. Is it necessary elsewhere? I am wondering if we can avoid adding RefElemData to the mesh type by modifying the interface for Trixi2Vtk to pass in a dg::DGMulti object

[jmark]

Thanks for putting this together @jmark! I'm still going through the related Trixi2Vtk PR. Is it correct to say that the main structural changes here are:

1. adding `dg.basis` to the mesh type,

Correct!

2. making `DGMultiMesh` a mutable struct?

Correct!

Some questions:

* why is making `DGMultiMesh` a mutable struct necessary?

Because the solution file interface of Trixi.jl modifies the fields

    current_filename :: String
    unsaved_changes  :: Bool

during the lifetime resp. after initialization of the mesh object.

* I saw that `RefElemData` is used to interpolate data in https://github.com/trixi-framework/Trixi2Vtk.jl/pull/103/files#. Is it necessary elsewhere? I am wondering if we can avoid adding `RefElemData` to the mesh type by modifying the interface for Trixi2Vtk to pass in a `dg::DGMulti` object

Good point! I agree with your proposal.

[jlchan]

Because the solution file interface of Trixi.jl modifies the fields

current_filename :: String
unsaved_changes  :: Bool

Gotcha! That makes sense to me.

[tristanmontoya]

• I saw that RefElemData is used to interpolate data in https://github.com/trixi-framework/Trixi2Vtk.jl/pull/103/files#. Is it necessary elsewhere? I am wondering if we can avoid adding RefElemData to the mesh type by modifying the interface for Trixi2Vtk to pass in a dg::DGMulti object

This doesn't seem straightforward to me, unless we make Trixi.jl always save the solver to file alongside mesh.h5 so that Trixi2Vtk.jl can load both. For solvers like DGSEM, we get away with loading only the mesh with Trixi.load_mesh, because it is simply assumed that the data is either equispaced or on LGL nodes.DGMulti is too general for that to work, so we need to access the RefElemData somehow. In my view, the simplest way to do this is to keep the field rd part of DGMultiMesh. If, going forward, we want to save and load a solver.h5 file instead, that could be an option as well, but we'd have to consider whether to do that just for DGMulti or for all Trixi.jl solvers. What do you think, @jlchan?

[jlchan]

This doesn't seem straightforward to me, unless we make Trixi.jl always save the solver to file alongside mesh.h5 so that Trixi2Vtk.jl can load both. For solvers like DGSEM, we get away with loading only the mesh with Trixi.load_mesh, because it is simply assumed that the data is either equispaced or on LGL nodes.DGMulti is too general for that to work, so we need to access the RefElemData somehow. In my view, the simplest way to do this is to keep the field rd part of DGMultiMesh. If, going forward, we want to save and load a solver.h5 file instead, that could be an option as well, but we'd have to consider whether to do that just for DGMulti or for all Trixi.jl solvers. What do you think, @jlchan?

First, sorry for the slow reply - just got back from a conference.

I really don't like adding RefElemData as a field to the mesh since it would make basis::RefElemData now a duplicated and redundant field.

I don't understand why other meshes don't take the basis into account. Aren't there finite difference SBP solvers which use a non-LGL node set? I guess I'm not as familiar with the save solution and i/o interface - is there no way to wrap a DGMultiMesh and DGMulti basis into a struct, write that single instance to hdf5, then load that?

[tristanmontoya]

I don't understand why other meshes don't take the basis into account. Aren't there finite difference SBP solvers which use a non-LGL node set?

By default, it's assumed that everything is LGL unless one specifies data_is_uniform = true in Trixi2Vtk.

I guess I'm not as familiar with the save solution and i/o interface - is there no way to wrap a DGMultiMesh and DGMulti basis into a struct, write that single instance to hdf5, then load that?

I think it is possible - I do something similar in StableSpectralElements.jl with JLD2 - but saving custom Julia structs doesn't seem like the favoured approach in Trixi.jl, which is understandable for backward compatibility reasons. I actually misunderstood what @jmark implemented, as I assumed he was directly saving a DGMultiMesh object (including both MeshData and RefElemData) to file, but what he is instead doing is saving the element type and polynomial degree to file and reconstructing a RefElemData based on that when Trixi.load_mesh is called from Trixi2Vtk.

This actually makes it easy to avoid duplicating RefElemData as it can just be reconstructed separately from the mesh object, but the approach is also restrictive as one has to assume that solution is stored on whatever StartUpDG's default interpolation nodes are for that reference element type. I suppose this could be generalized by additionally storing rd.VDM (or the reference node positions, if we restrict ourselves to storing the solution as nodal values and not modal coefficients) in mesh.h5.

[jlchan]

I guess I'm not as familiar with the save solution and i/o interface - is there no way to wrap a DGMultiMesh and DGMulti basis into a struct, write that single instance to hdf5, then load that?

I think it is possible - I do something similar in StableSpectralElements.jl with JLD2 - but saving custom Julia structs doesn't seem like the favoured approach in Trixi.jl, which is understandable for backward compatibility reasons. I actually misunderstood what @jmark implemented, as I assumed he was directly saving a DGMultiMesh object (including both MeshData and RefElemData) to file, but what he is instead doing is saving the element type and polynomial degree to file and reconstructing a RefElemData based on that when Trixi.load_mesh is called from Trixi2Vtk.
This actually makes it easy to avoid duplicating RefElemData as it can just be reconstructed separately from the mesh object, but the approach is also restrictive as one has to assume that solution is stored on whatever StartUpDG's default interpolation nodes are for that reference element type. I suppose this could be generalized by additionally storing rd.VDM in mesh.h5.

I see - but RefElemData is still being added to MeshData in this PR, which I'm guessing is because save_mesh_file only takes the mesh as input?

I'd also note is that this approach might not work if using approximation_type = SBP() on triangles. We'd also have to save the rd.approximation_type to reconstruct general RefElemData.

[jlchan]

I'll consider adding RefElemData to DGMultiMesh, but I think it's pretty hacky and would like to avoid it if possible.

To clarify, is save_mesh_file the problematic function here that requires adding RefElemData to DGMultiMesh? If so, my preference would be to specialize further up in the call chain; for example,

Trixi.jl/src/callbacks_step/save_solution.jl

Lines 152 to 166 in c36369f

	function save_mesh(semi::AbstractSemidiscretization, output_directory, timestep = 0)
	mesh, _, _, _ = mesh_equations_solver_cache(semi)
	if mesh.unsaved_changes
	# We only append the time step number to the mesh file name if it has
	# changed during the simulation due to AMR. We do not append it for
	# the first time step.
	if timestep == 0
	mesh.current_filename = save_mesh_file(mesh, output_directory)
	else
	mesh.current_filename = save_mesh_file(mesh, output_directory, timestep)
	end
	mesh.unsaved_changes = false
	end
	end

could be specialized for DGMulti semi-discretizations, with both the basis and mesh being passed in to save_mesh_file.

@tristanmontoya @jmark thoughts?

[tristanmontoya]

@tristanmontoya @jmark thoughts?

Yes, I think that's the main thing that would need to be changed. Even with RefElemData part of DGMultiMesh as in this PR, we still are not saving the approximation type, and just assuming it is Polynomial. What I would propose as a first step is to update this PR to remove RefElemData from DGMultiMesh, and add the specialization on the semidiscretization as you suggested. I could do this and then update the Trixi2Vtk PR accordingly.

Wouldn't it still be possible to visualize solutions from SBP approximations as long as we don't reinterpolate?

[jlchan]

As discussed on Slack:

[make] a special case to allow both tuple and integer polynomial degrees for wedges.

Adding a test for the case when using Wedge with N::Int would also be great.

Thanks!

[jlchan]

Thanks for differentiating between Wedges and TensorProductWedges! This looks good to me. Can you open an issue on support for non-Polynomial approximation types?

[tristanmontoya]

Thanks for differentiating between Wedges and TensorProductWedges! This looks good to me. Can you open an issue on support for non-Polynomial approximation types?

Sure. I am still trying to figure out why the ubuntu CI fails for the prism test with scalar polynomial degree (on my Mac get linf = 0.0006597931027270132 while the ubuntu CI gets linf = 0.0006597931292229298). Not sure what to do there.

[jlchan]

Does changing the timestepper make the test results more consistent? Can you link the elixir that is failing?

[jlchan]

Looks like tests pass now - guess it was an issue of tolerances?

[tristanmontoya]

Hi @jlchan, the elixir that failed was https://github.com/trixi-framework/Trixi.jl/blob/dgmulti-mesh-io-support/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl, but only when I use a scalar polydeg=3 instead of the tensor-product polynomial degree (3, 4). Relaxing the atol to 1e-10 makes the tests pass.

I can get the difference in linf error values between Linux and Mac to be a bit below 1e-11 with a smaller or fixed time step, but beyond that doesn't seem to be possible, so I don't think it's an issue with OrdinaryDiffEq.jl. Altogether, I don't like having to relax the tolerance for this specific test without understanding why it's needed, but I am also not sure how tractable the issue is. Let me know if you're ok with merging this PR with the relaxed tolerance on this test.

[jlchan]

Hi @jlchan, the elixir that failed was https://github.com/trixi-framework/Trixi.jl/blob/dgmulti-mesh-io-support/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl, but only when I use a scalar polydeg=3 instead of the tensor-product polynomial degree (3, 4). Relaxing the atol to 1e-10 makes the tests pass.

I can get the difference in linf error values between Linux and Mac to be a bit below 1e-11 with a smaller or fixed time step, but beyond that doesn't seem to be possible, so I don't think it's an issue with OrdinaryDiffEq.jl. Altogether, I don't like having to relax the tolerance for this specific test without understanding why it's needed, but I am also not sure how tractable the issue is. Let me know if you're ok with merging this PR with the relaxed tolerance on this test.

Gotcha. That could be a conditioning issue since the polydeg=3 case directly does linear algebra with larger 3D matrices, so their could be some roundoff stuff. polydeg=3 seems too small to observe these issues, however.

One other check - does lowering the CFL make a difference?

[tristanmontoya]

Whether using a (much) smaller CFL, a different time-marching method, or turning off adaptive time stepping, linf bottoms out between around 1e-12 and 1e-11 and thus doesn't meet Trixi.jl's default test tolerances.

[tristanmontoya]

If you'd rather not merge until the tolerance issue is resolved, I could go through some of the DGMulti spatial discretization code to see where the larger floating-point errors arise. I wouldn't have much time to look at it this week though, so it may have to wait.

[jlchan]

If you'd rather not merge until the tolerance issue is resolved, I could go through some of the DGMulti spatial discretization code to see where the larger floating-point errors arise. I wouldn't have much time to look at it this week though, so it may have to wait.

Thanks, but I think it's ok to merge this with the relaxed tolerance for now and make an issue describing the CI failure under a tighter tolerance. Can you do file that issue?

[tristanmontoya]

@jlchan Sounds good to me - I have opened the issues you suggested in #2430 and #2431.

[jlchan]

Thanks for pushing this through @tristanmontoya!