Add Navier–Stokes flow-around-cylinder example using projection method#630
Add Navier–Stokes flow-around-cylinder example using projection method#630Ouardghi wants to merge 2 commits intoParallel-in-Time:masterfrom
Conversation
There was a problem hiding this comment.
Pull request overview
Adds a new 2D incompressible Navier–Stokes (flow around a cylinder) example for the StroemungsRaum project using a projection (fractional-step) method in a FEniCS-based problem class, plus a custom IMEX mass-matrix sweeper, a run script, and tests.
Changes:
- Introduce
fenics_NSE_2D_massproblem class implementing the projection-based NSE solve and RHS splitting. - Add
imex_1st_order_mass_NSEsweeper adapting the IMEX mass-matrix sweeper for NSE (including custom residual handling). - Add a runnable setup/run script and a FEniCS-marked pytest suite validating
solve_system,eval_f, and a benchmark-based drag/lift check.
Reviewed changes
Copilot reviewed 4 out of 5 changed files in this pull request and generated 8 comments.
| File | Description |
|---|---|
pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_FEniCS.py |
New FEniCS NSE problem class with projection method, BCs, and drag/lift measure. |
pySDC/projects/StroemungsRaum/sweepers/imex_1st_order_mass_NSE.py |
New IMEX mass sweeper variant calling the NSE projection solve and custom residual logic. |
pySDC/projects/StroemungsRaum/run_Navier_Stokes_equations_FEniCS.py |
New setup + execution helper for running the cylinder benchmark. |
pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_IMEX_FEniCS.py |
New tests for the NSE system solve, RHS evaluation, and benchmark coefficients. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
You can also share your feedback on Copilot code review. Take the survey.
| def compute_residual(self, stage=None): | ||
| """ | ||
| Computation of the residual using the collocation matrix Q | ||
|
|
||
| Args: | ||
| stage (str): The current stage of the step the level belongs to | ||
| """ | ||
|
|
||
| # get current level and problem description | ||
| L = self.level | ||
|
|
||
| # Check if we want to skip the residual computation to gain performance | ||
| # Keep in mind that skipping any residual computation is likely to give incorrect outputs of the residual! | ||
| if stage in self.params.skip_residual_computation: | ||
| L.status.residual = 0.0 if L.status.residual is None else L.status.residual | ||
| return None | ||
|
|
||
| # compute the residual for each node | ||
| # build QF(u) | ||
| res_norm = [] | ||
| res = [0] * (self.coll.num_nodes + 1) | ||
| for m in range(self.coll.num_nodes): | ||
|
|
||
| # compute the residual at node m, using the incremental criterion | ||
| if L.uold[m + 1] is not None: | ||
| res[m] = L.u[m + 1] - L.uold[m + 1] | ||
| else: | ||
| res[m] = L.u[m + 1] | ||
|
|
||
| # Due to different boundary conditions we might have to fix the residual | ||
| if L.prob.fix_bc_for_residual: | ||
| L.prob.fix_residual(res[m]) | ||
| # use abs function from data type here | ||
| res_norm.append(abs(res[m])) | ||
|
|
||
| # find maximal residual over the nodes | ||
| L.status.residual = max(res_norm) | ||
|
|
||
| if L.time == 3.1250000000e-04 and L.status.residual == 0.0: | ||
| L.status.residual = 1.0 | ||
|
|
There was a problem hiding this comment.
I don't like this magic number either..
There was a problem hiding this comment.
You are right, I should use the initial time t0 instead of a hard-coded value.
I introduced this check because zero is a valid solution to the problem. Since the residual is evaluated first, it becomes zero, and the simulation terminates immediately. To prevent this, I modify the residual only at the first time step, and only if it is zero, so that the iteration can proceed:
if L.time == L.prob.t0 and L.status.residual == 0.0:
L.status.residual = 1.0
There was a problem hiding this comment.
We fixed this issue here by always doing one SDC iteration, no? Is this really still needed?
In this file you can also see that you should pass level_params['e_tol'] = x to get an increment tolerance of x rather than assign the increment to the residual and use residual tolerance.
pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_IMEX_FEniCS.py
Outdated
Show resolved
Hide resolved
| class fenics_NSE_2D_mass(Problem): | ||
| r""" | ||
| Example implementing a forced two-dimensional incompressible Navier–Stokes problem for the DFG | ||
| benchmark flow around cyilder using FEniCS (dolfin). |
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
brownbaerchen
left a comment
brownbaerchen
left a comment
There was a problem hiding this comment.
Thanks for the PR, much cleaner than previous ones! I didn't really check the maths, but I trust you with it.
| """ | ||
| Custom sweeper class, implements Sweeper.py | ||
|
|
||
| First-order IMEX sweeper using implicit/explicit Euler as base integrator, with mass or weighting matrix |
There was a problem hiding this comment.
| First-order IMEX sweeper using implicit/explicit Euler as base integrator, with mass or weighting matrix | |
| First-order IMEX sweeper with mass or weighting matrix including a project step for Navier-Stokes equation |
Not sure if this is the best possible docstring, but the preconditioner is a parameter and should not be mentioned here.
| u = self.dtype_u(u0) | ||
| p = df.Function(self.Q) |
There was a problem hiding this comment.
This problem could probably benefit from its own datatype that gathers u and p into a single datatype. It's not urgently needed and since this project is at its end, maybe don't do it? What do you think @pancetta?
| # compute the residual at node m, using the incremental criterion | ||
| if L.uold[m + 1] is not None: | ||
| res[m] = L.u[m + 1] - L.uold[m + 1] | ||
| else: | ||
| res[m] = L.u[m + 1] |
There was a problem hiding this comment.
What's going on here? If you want to use an incremental tolerance for stopping the iteration, that is done elsewhere. The residual should be the residual and nothing else.
There was a problem hiding this comment.
It is hard to use the residual here, as the residual corresponds only to the first step of the projection method, where SDC is applied. However, the solver actually consists of multiple coupled substeps: after computing the tentative velocity, a pressure correction is performed, and the velocity is projected to enforce incompressibility.
I thought the easiest way to handle this was to modify the residual. Is there a simpler or alternative way to achieve this?
There was a problem hiding this comment.
As mentioned above, we have infrastructure for using increment tolerance, you can simply pass residual tolerance and no residual tolerance.
This PR adds a new example for the incompressible Navier–Stokes equations using a projection method for the classical flow around a cylinder benchmark.
The goal of this contribution is to provide a working example of a projection-based Navier–Stokes solver within the pySDC framework.
This PR includes
This example can serve as a reference for implementing incompressible flow solvers and for experimenting with projection-based approaches in pySDC.