Coverage for pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_FEniCS.py: 100%
94 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import logging
2import dolfin as df
3import numpy as np
4import os
6from pySDC.core.problem import Problem
7from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh
10class fenics_NSE_2D_mass(Problem):
11 r"""
12 Example implementing a forced two-dimensional incompressible Navier–Stokes problem for the DFG
13 benchmark flow around cylinder using FEniCS (dolfin).
15 The unknowns are the velocity field :math:`\mathbf{u}(x,y,t)` and the pressure field
16 :math:`p(x,y,t)` on a 2D domain :math:`\Omega` (here loaded from ``cylinder.xml``).
17 The model is
19 .. math::
20 \frac{\partial \mathbf{u}}{\partial t}
21 + (\mathbf{u}\cdot\nabla)\mathbf{u}
22 - \nu \Delta \mathbf{u}
23 + \nabla p
24 = \mathbf{f},
25 \qquad
26 \nabla\cdot \mathbf{u} = 0,
28 where :math:`\nu` is the kinematic viscosity and :math:`\mathbf{f}` is a forcing term.
30 This implementation uses a fractional-step /projection method:
31 a viscous velocity solve, a pressure Poisson solve enforcing incompressibility, and a velocity
32 correction using the pressure gradient.
34 Boundary conditions are applied on subsets of the boundary:
35 - no-slip on channel walls and cylinder surface,
36 - a time-dependent inflow profile at the inlet,
37 - pressure condition at the outflow.
39 Parameters
40 ----------
41 t0 : float, optional
42 Starting time.
43 family : str, optional
44 Finite element family for velocity and pressure spaces. Default is ``'CG'``.
45 order : int, optional
46 Polynomial degree for the velocity space. The pressure degree is chosen as ``order - 1``.
47 nu : float, optional
48 Kinematic viscosity :math:`\nu`.
50 Attributes
51 ----------
52 mesh : Mesh
53 Computational mesh loaded from ``cylinder.xml``.
54 V : VectorFunctionSpace
55 Function space for the velocity field.
56 Q : FunctionSpace
57 Function space for the pressure field.
58 M : Matrix
59 Assembled velocity mass matrix, corresponding to :math:`\int_\Omega \mathbf{u}\cdot\mathbf{v}\,dx`.
60 K : Matrix
61 Assembled viscous operator matrix, corresponding to
62 :math:`- \int_\Omega \nabla \mathbf{u} : (\nu \nabla \mathbf{v})\,dx`.
63 Mp : Matrix
64 Assembled pressure mass matrix, corresponding to :math:`\int_\Omega p\,q\,dx`.
65 g : Expression
66 Forcing term :math:`\mathbf{f}` (here set to zero).
67 bcu : list[DirichletBC]
68 Velocity Dirichlet boundary conditions (walls, inflow, cylinder).
69 bcp : list[DirichletBC]
70 Pressure boundary conditions (outflow).
71 bc_hom : DirichletBC
72 Homogeneous velocity boundary condition used to fix the residual when requested.
74 Notes
75 -----
76 The right-hand side evaluation splits into an implicit viscous part and an explicit part
77 containing forcing and nonlinear convection:
78 - implicit: :math:`K\,\mathbf{u}`
79 - explicit: :math:`M\,\mathbf{f} - M\,(\mathbf{u}\cdot\nabla)\mathbf{u}`
80 """
82 dtype_u = fenics_mesh
83 dtype_f = rhs_fenics_mesh
85 def __init__(self, t0=0.0, family='CG', order=2, nu=0.001):
87 # set logger level for FFC and dolfin
88 logging.getLogger('FFC').setLevel(logging.WARNING)
89 logging.getLogger('UFL').setLevel(logging.WARNING)
91 # set solver and form parameters
92 df.parameters["form_compiler"]["optimize"] = True
93 df.parameters["form_compiler"]["cpp_optimize"] = True
94 df.parameters['allow_extrapolation'] = True
96 # load mesh
97 path = f"{os.path.dirname(__file__)}/../meshs/"
98 mesh = df.Mesh(path + '/cylinder.xml')
100 # define function space for future reference
101 self.V = df.VectorFunctionSpace(mesh, family, order)
102 self.Q = df.FunctionSpace(mesh, family, order - 1)
104 tmp_u = df.Function(self.V)
105 tmp_p = df.Function(self.Q)
106 self.logger.debug('Velocity DoFs on this level:', len(tmp_u.vector()[:]))
107 self.logger.debug('Pressure DoFs on this level:', len(tmp_p.vector()[:]))
109 super().__init__(self.V)
110 self._makeAttributeAndRegister('t0', 'family', 'order', 'nu', localVars=locals(), readOnly=True)
112 # define trial and test functions
113 self.u = df.TrialFunction(self.V)
114 self.v = df.TestFunction(self.V)
115 self.p = df.TrialFunction(self.Q)
116 self.q = df.TestFunction(self.Q)
118 # initialize solution functions for pressure
119 self.pn = df.Function(self.Q)
121 # mass and stiffness terms
122 a_K = -1.0 * df.inner(df.nabla_grad(self.u), self.nu * df.nabla_grad(self.v)) * df.dx
123 a_M = df.inner(self.u, self.v) * df.dx
124 a_S = df.inner(df.nabla_grad(self.p), df.nabla_grad(self.q)) * df.dx
125 a_Mp = df.inner(self.p, self.q) * df.dx
127 # assemble the forms
128 self.M = df.assemble(a_M)
129 self.K = df.assemble(a_K)
130 self.S = df.assemble(a_S)
131 self.Mp = df.assemble(a_Mp)
133 # set inflow profile at the domain inlet
134 self.inflow_profile = df.Expression(
135 ('4.0*1.5*sin(pi*t/8)*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0'), pi=np.pi, t=t0, degree=self.order
136 )
138 # define boundaries
139 inflow = 'near(x[0], 0)'
140 outflow = 'near(x[0], 2.2)'
141 walls = 'near(x[1], 0) || near(x[1], 0.41)'
142 cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'
144 # define boundary conditions
145 bcu_noslip = df.DirichletBC(self.V, df.Constant((0, 0)), walls)
146 bcu_inflow = df.DirichletBC(self.V, self.inflow_profile, inflow)
147 bcu_cylinder = df.DirichletBC(self.V, df.Constant((0, 0)), cylinder)
148 bcp_outflow = df.DirichletBC(self.Q, df.Constant(0), outflow)
150 # set boundary values
151 self.bcu = [bcu_noslip, bcu_inflow, bcu_cylinder]
152 self.bcp = [bcp_outflow]
153 self.bc_hom = df.DirichletBC(self.V, df.Constant((0.0, 0.0)), 'on_boundary')
154 self.fix_bc_for_residual = True
156 # Define measure for drag and lift computation
157 Cylinder = df.CompiledSubDomain(cylinder)
158 CylinderBoundary = df.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
159 Cylinder.mark(CylinderBoundary, 1)
160 self.dsc = df.Measure("ds", domain=mesh, subdomain_data=CylinderBoundary, subdomain_id=1)
162 # set forcing term as expression
163 self.g = df.Expression(('0.0', '0.0'), t=self.t0, degree=self.order)
165 # initialize XDMF files for velocity and pressure if needed
166 self.xdmffile_p = None
167 self.xdmffile_u = None
169 def use_restricted_boundary_values(self, u, t):
170 """
171 Refuse multilevel runs: the solves impose the time-dependent boundary data, which at the coarse nodes is not
172 the restricted fine boundary values, so MLSDC and PFASST would not converge to the collocation solution. See
173 :meth:`pySDC.core.problem.Problem.use_restricted_boundary_values`.
174 """
175 raise NotImplementedError(f'{type(self).__name__} cannot impose restricted boundary values on coarse levels')
177 def solve_system(self, rhs, factor, u0, t, dtau):
178 r"""
179 Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
181 Parameters
182 ----------
183 rhs : dtype_f
184 Right-hand side for the nonlinear system.
185 factor : float
186 Factor for the linear system.
187 u0 : dtype_u
188 Initial guess for the iterative solver (not used here so far).
189 t : float
190 Current time.
191 dtau : float
192 Zero-to-node step size in the SDC sweep, used in the pressure correction step.
194 Returns
195 -------
196 u : dtype_u
197 Solution: velocity.
198 p : space function
199 Solution: pressure.
201 """
203 u = self.dtype_u(u0)
204 p = df.Function(self.Q)
205 T = self.M - factor * self.K
207 # solve for the intermediate velocity
208 self.inflow_profile.t = t
209 [bc.apply(T, rhs.values.vector()) for bc in self.bcu]
210 df.solve(T, u.values.vector(), rhs.values.vector())
212 # solve for the pressure correction
213 L2 = -(1 / dtau) * df.div(u.values) * self.q * df.dx
214 Rhs2 = df.assemble(L2)
215 [bc.apply(self.S, Rhs2) for bc in self.bcp]
216 df.solve(self.S, p.vector(), Rhs2)
218 # velocity correction
219 L3 = df.dot(u.values, self.v) * df.dx - dtau * df.dot(df.nabla_grad(p), self.v) * df.dx
220 Rhs3 = df.assemble(L3)
221 [bc.apply(self.M, Rhs3) for bc in self.bcu]
222 df.solve(self.M, u.values.vector(), Rhs3)
224 return u, p
226 def eval_f(self, u, t):
227 """
228 Routine to evaluate both parts of the right-hand side.
230 Parameters
231 ----------
232 u : dtype_u
233 Current values of the numerical solution.
234 t : float
235 Current time at which the numerical solution is computed.
237 Returns
238 -------
239 f : dtype_f
240 The right-hand side divided into two parts.
241 """
242 self.g.t = t
243 conv = -1.0 * df.dot(u.values, df.nabla_grad(u.values))
245 f = self.dtype_f(self.V)
247 self.K.mult(u.values.vector(), f.impl.values.vector())
248 f.expl = self.apply_mass_matrix(self.dtype_u(df.project(conv, self.V)))
249 f.expl += self.apply_mass_matrix(self.dtype_u(df.interpolate(self.g, self.V)))
251 return f
253 def apply_mass_matrix(self, u):
254 r"""
255 Routine to apply mass matrix.
257 Parameters
258 ----------
259 u : dtype_u
260 Current values of the numerical solution.
262 Returns
263 -------
264 me : dtype_u
265 The product :math:`M \vec{u}`.
266 """
268 me = self.dtype_u(self.V)
269 self.M.mult(u.values.vector(), me.values.vector())
271 return me
273 def u_exact(self, t):
274 r"""
275 Routine to compute the exact solution at time :math:`t`.
277 Parameters
278 ----------
279 t : float
280 Time of the exact solution.
282 Returns
283 -------
284 me : dtype_u
285 Exact solution.
286 """
288 u0 = df.Expression(('0.0', '0.0'), degree=self.order)
289 me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)
291 return me
293 def fix_residual(self, res):
294 """
295 Applies homogeneous Dirichlet boundary conditions to the residual
297 Parameters
298 ----------
299 res : dtype_u
300 Residual
301 """
302 self.bc_hom.apply(res.values.vector())
303 return None