SDIRK scheme implementation

[etiennemarin]

Description

This pull request implements the SDIRK (Singly Diagonally Implicit Runge-Kutta) time integration method for the unsteady Navier–Stokes equations.

The current implementation is quite brute-force and would surely benefit from a proper refactoring to improve efficiency and readability. The idea is to also use the SDIRK scheme for other applications beyond just lethe-fluid.

The code is structured to match as closely as possible the existing BDF structure. Here's how it works:

• At each time iteration, we loop over the stages. For "traditional" schemes (mainly BDF), the loop is straightforward: there is only one stage.
• At each stage, for the SDIRK scheme, we compute the vectors kᵢ.
• However, the structure of Lethe is designed to solve the nonlinear system on the velocity. So instead of directly solving kᵢ = F(tᵢ, uᵢ), we solve the equivalent system:

(uᵢ* - uₙ) / (h aᵢᵢ) - ∑_{j=0}^{i-1} (aᵢⱼ / aᵢᵢ) kⱼ = F(tᵢ, uᵢ)

• Once we have computed the $\mathbf{u}_i^*$ vector, we store the associated $\mathbf{k}_i$ vector in order to compute the next stages and the final sum.
• At the end of the loop, the solution is computed as:

uₙ₊₁ = uₙ + h ∑ bᵢ kᵢ

The main changes are implemented in the solvers/navier_stokes_base.cc file, within the iterate() function.
The Butcher table is stored in scratch-data because it's required during matrix assembly.

---

Testing

The code has been tested on the lethe-fluid/application-test Taylor–Green vortex (GLS formulation) and yields results close to the BDF2 scheme.
The test outputs, at each time step, the enstrophy, the kinetic energy, and the L2 error.

The goal is now to test this new scheme on multiple examples and with different Butcher tables.

---

Documentation

• The Butcher table is set in the scratch-data.
• The variable sum_over_previous_stages is now required for matrix assembly.

Miscellaneous (will be removed when merged)

Checklist (will be removed when merged)

See this page for more information about the pull request process.

Code related list:

• All in-code documentation related to this PR is up to date (Doxygen format)
• Copyright headers are present and up to date
• Lethe documentation is up to date
• New feature has unit test(s) (preferred) or application test(s), and restart files are in the generator folder
• The branch is rebased onto master
• Changelog (CHANGELOG.md) is up to date
• Code is indented with indent-all and .prm files (examples and tests) with prm-indent

Pull request related list:

• Labels are applied
• There are at least 2 reviewers (or 1 if small feature) excluding the responsible for the merge
• If this PR closes an issue or is related to a project, it is linked in the "Projects" or "Development" section
• If the fix is temporary, an issue is opened
• The PR description is cleaned and ready for merge

[blaisb]

Can you add the documentation? Furthermore, have you tested if you recoved the 2nd and 3rd order convergence rate for both SDIRK implementations?

[blaisb]

First round of comment.

[blaisb]

So you added a new bdf2 test case? Is it to have some sort of stored comparison between bdf2,sdirk22 and sdirk33?
In that case, I would suggest you do a case with a larger time-step so that the accuracy in the case be limited by the time-step and not the mesh resolution. Do. like a time-step of 0.2 maybe

[etiennemarin]

Yes that was my idea since the BDF2 is broadly used in Lethe. I found it more relevant than BDF1. I will do the simulation again with the larger time step

[blaisb]

needs a brief and a description of the return. PLease be also a bit more verbose. What does i mean here?

[etiennemarin]

Actually, since we are calculating the sum of a_ij * k_j in the loop in the iterate() function, I think we don't need these two functions. In the assembler, we only need to have the diagonal term, so we don't need to know at which stage we are. In the loop, we already know the stage. These two functions can maybe be usefull for more general implicit schemes, so I don't know if you want me to keep them. But for SDIRK, it is not necessary

[blaisb]

If you don't use them, then you can remove them for now i presume.

[blaisb]

Same as above comment.

[etiennemarin]

See the comment just above about the need (or not) to keep these two functions in the code

[blaisb]

Ok that's a lot of vectors... You really don't need that many! Please try to clean this to have as little as possible.

[blaisb]

??? You don't need this

[etiennemarin]

I need the Butcher's table in the iterate() function to build the sums. The table is in the scratch data. But I also need this table for the matrix assembly. Should I store the table in an other location ?

[blaisb]

You can recalculate the Butcher table when. you need it there. Just remake it and not store it as part of the class. It will be simpler. There,s no point in having a scratch data live just for a table. Just reget the table. Otherwise if you want to store it, store it in simulation control where the time related stuff is stored.

[blaisb]

Why don't you just do instead

Suggested change

	temp_sum_over_previous_stages.add(1.0, tmp);
	temp_sum_over_previous_stages.add(stage_data.a_ij[p] / a_ii, previous_k_j_solutions[p]);

[etiennemarin]

Yes, I don't need to have all these lines, your solution works. But we still need to have a non-locally relevant vector in the .add operation, and previous_k_j_solutions[p] is locally-relevant

[blaisb]

I understand, but why create a new locally_owned vector? You can use local evaluation point to store your temporary solution, this will remove one more vector.
Try to minimize the amount of local vectors you generate.

[blaisb]

Suggested change

	previous_k_j_solutions_p);
	previous_k_j_solutions[p]);

[etiennemarin]

In the operations (addition, multiplication, etc.) we need a non-locally relevant dofs vector, I think

[blaisb]

@etiennemarin Once you have done the convergence study, can you ask for two other reviewer outside of me? You can use the lethe channel in the slack to ask people. Don't forget to check the "check list" and ensure that you can really check all tasks of the checklist. It is a reminder of what you need to do when you make a PR.