Add a monolithic incompressible Navier-Stokes example and a fully implicit mass sweeper#642
Add a monolithic incompressible Navier-Stokes example and a fully implicit mass sweeper#642Ouardghi wants to merge 1 commit intoParallel-in-Time:masterfrom
Conversation
…licit mass sweeper
|
Thanks! Is this the Taylor-Green example that shows order reduction in the presence of time-dependent BCs? If not, could you implement this, too (maybe as a separate PR once this here is merged)? See #640 |
There was a problem hiding this comment.
Pull request overview
This PR adds a monolithic (mixed velocity–pressure) incompressible Navier–Stokes FEniCS problem together with a new fully-implicit sweeper that supports mass-matrix formulations, plus a runnable example script and regression tests.
Changes:
- Added
fenics_NSE_2D_Monolithicimplementing a coupled (monolithic) Navier–Stokes discretization in a mixed FE space. - Added
generic_implicit_mass, extendinggeneric_implicitto handle problems formulated with a mass matrix. - Added a runnable monolithic NSE driver script and a dedicated FEniCS test module covering
solve_system,eval_f, and a short benchmark step.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 8 comments.
| File | Description |
|---|---|
pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py |
New monolithic mixed FE Navier–Stokes problem with mass-matrix helpers and drag/lift boundary measure. |
pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py |
New fully implicit sweeper variant that incorporates a mass matrix into node updates/residual computation. |
pySDC/projects/StroemungsRaum/run_Navier_Stokes_equations_monolithic_FEniCS.py |
New example script to run the monolithic NSE setup with the new sweeper. |
pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_monolithic_FEniCS.py |
New FEniCS tests for the monolithic problem + example setup. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| sweeper_params = dict() | ||
| sweeper_params['quad_type'] = 'RADAU-RIGHT' | ||
| sweeper_params['num_nodes'] = 2 | ||
| sweeper_params['QI'] = ['LU'] |
There was a problem hiding this comment.
sweeper_params['QI'] is set to a list (['LU']), but generic_implicit/Sweeper.get_Qdelta_implicit expects a string key (e.g. 'LU'). Passing a list will raise a TypeError: unhashable type: 'list' when building the QΔ generator, which will also break test_problem_class() since it imports setup().
Set this to 'LU' (string), and similarly keep QE/QI/num_nodes scalar values in this script unless you explicitly handle parameter sweeps elsewhere.
| sweeper_params['QI'] = ['LU'] | |
| sweeper_params['QI'] = 'LU' |
| from pathlib import Path | ||
| import json | ||
|
|
There was a problem hiding this comment.
Path and json are imported but never used in this script. This adds noise and can trigger linting failures in environments that enforce unused-import checks.
Remove the unused imports, or add the intended postprocessing code that uses them (similar to run_Navier_Stokes_equations_FEniCS.py).
| from pathlib import Path | |
| import json |
| # implicit solve with prefactor stemming from the diagonal of Qd | ||
| alpha = L.dt * self.QI[m + 1, m + 1] | ||
| if alpha == 0: | ||
| L.u[m + 1] = rhs |
There was a problem hiding this comment.
For mass-matrix problems, the alpha == 0 shortcut is not equivalent to the non-mass case: the stage equation becomes M u = rhs (so you need u = M^{-1} rhs). Assigning L.u[m+1] = rhs will store the mass-weighted vector as the solution, which is inconsistent with the rest of the sweeper (and will be wrong for collocation rules/QΔ choices that yield alpha=0, e.g. first node at 0).
Instead of assigning directly, call P.solve_system(rhs, 0.0, ...) (letting the problem handle the mass inversion), or explicitly apply the appropriate mass-matrix inverse before storing L.u[m+1].
| L.u[m + 1] = rhs | |
| L.u[m + 1] = P.solve_system(rhs, 0.0, L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) |
| class generic_implicit_mass(generic_implicit): | ||
|
|
||
| def update_nodes(self): | ||
| """ | ||
| This sweeper extends the generic_implicit sweeper from the implementations | ||
| package for a monolithic discretization of the incompressible Navier-Stokes | ||
| equations with a mass matrix. It updates the solution and right-hand side | ||
| values at all collocation nodes during a single sweep. | ||
|
|
||
| Returns: | ||
| None | ||
| """ |
There was a problem hiding this comment.
generic_implicit_mass modifies the sweep to use a mass matrix, but it does not override compute_end_point(). If a collocation rule is used where right_is_node is false (or do_coll_update is forced true), the inherited endpoint update will perform an explicit Picard-style update without applying M^{-1}, which is not valid for M u' = f(u) when f is assembled in the mass-weighted form.
Consider overriding compute_end_point() to either (a) require right_is_node and copy the last node (as in imex_1st_order_mass), or (b) implement a correct mass-aware endpoint update.
| r""" | ||
| Example implementing the forced two-dimensional incompressible Navier-Stokes equations with | ||
| time-dependent Dirichlet boundary conditions | ||
|
|
||
| .. math:: | ||
| \frac{\partial u}{\partial t} = - u \cdot \nabla u + \nu \nabla u - \nabla p + f | ||
| 0 = \nabla \cdot u | ||
|
|
||
| for :math:`x \in \Omega`, where the forcing term :math:`f` is defined by | ||
|
|
||
| .. math:: | ||
| f(x, t) = (0, 0). | ||
|
|
||
| This implementation follows a monolithic approach, where velocity and pressure are | ||
| solved simultaneously in a coupled system using a mixed finite element formulation. | ||
|
|
||
| Boundary conditions are applied on subsets of the boundary: | ||
| - no-slip on channel walls and cylinder surface, | ||
| - a time-dependent inflow profile at the inlet, | ||
| - pressure condition at the outflow. | ||
|
|
||
| In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem | ||
| is reformulated to the *weak formulation* | ||
|
|
||
| .. math: | ||
| \int_\Omega u_t v\,dx = - \int_\Omega u \cdot \nabla u v\,dx - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega p \nabla \cdot v\,dx + \int_\Omega f v\,dx | ||
| \int_\Omega \nabla \cdot u q\,dx = 0 |
There was a problem hiding this comment.
The docstring has a few mathematical / Sphinx issues that make it misleading or render incorrectly:
- The strong form uses
\nu \nabla ubut should be\nu \Delta u(and typically-\nu\Delta udepending on sign convention). - The “weak formulation” block starts with
.. math:instead of.. math::. - Attribute description says
gis the forcing term “in the heat equation”, which is unrelated here.
Please fix these so the documentation matches the implemented monolithic Navier–Stokes formulation and builds cleanly in Sphinx.
|
|
||
| # print the number of DoFs for debugging purposes | ||
| tmp = df.Function(self.W) | ||
| self.logger.debug('DoFs on this level:', len(tmp.vector()[:])) |
There was a problem hiding this comment.
This debug log call uses the logging API in a way that can raise a TypeError if the logger is set to DEBUG: the extra argument is treated as %-formatting args, but the message string has no placeholder.
Use self.logger.debug('DoFs on this level: %d', n_dofs) or an f-string so enabling debug logging doesn’t break initialization.
| self.logger.debug('DoFs on this level:', len(tmp.vector()[:])) | |
| self.logger.debug('DoFs on this level: %d', len(tmp.vector()[:])) |
| """ | ||
|
|
||
| me = self.dtype_u(self.W) | ||
| df.solve(self.Mf, me.values.vector(), w.values.vector()) |
There was a problem hiding this comment.
__invert_mass_matrix calls df.solve(self.Mf, ...) each time solve_system runs. Since Mf is constant, repeatedly refactorizing/setting up the linear solve can dominate runtime (this is called once per node per sweep).
Consider creating and caching a df.LUSolver/df.KrylovSolver for Mf in __init__ (or caching the factorization) and reusing it here to avoid repeated setup costs.
| df.solve(self.Mf, me.values.vector(), w.values.vector()) | |
| if not hasattr(self, '_Mf_solver'): | |
| self._Mf_solver = df.LUSolver(self.Mf) | |
| self._Mf_solver.solve(me.values.vector(), w.values.vector()) |
| F += df.dot(g.values, self.v) * df.dx | ||
| F += df.dot(df.div(u), self.q) * df.dx | ||
|
|
||
| f.values.vector()[:] = df.assemble(F) |
There was a problem hiding this comment.
eval_f assembles into a temporary vector and then assigns into f.values.vector()[:]. This forces an allocation every call and can be a significant cost in SDC sweeps.
Prefer assembling directly into the existing vector (e.g. via df.assemble(F, tensor=f.values.vector()) or an equivalent in-place assemble) to avoid repeated allocations and improve performance.
| f.values.vector()[:] = df.assemble(F) | |
| df.assemble(F, tensor=f.values.vector()) |
| u, p = df.split(w.values) | ||
|
|
||
| # get the SDC right-hand side | ||
| rhs = self.__invert_mass_matrix(rhs) |
There was a problem hiding this comment.
Do you really need to invert the mass matrix? Can't you rewrite all of this multiplying the whole thing with the mass matrix?
This PR adds a new example for solving the incompressible Navier-Stokes equations using a monolithic formulation.
The PR introduces a new problem class together with a new sweeper,
generic_implicit_mass.This sweeper extends the existinggeneric_implicitsweeper to support fully implicit problems with a mass matrix.The existing mass-based sweepers are designed for IMEX methods, so they cannot be used directly for a fully implicit monolithic Navier-Stokes discretization.
Main changes:
generic_implicit_mass,a fully implicit sweeper with mass matrix support.