Coverage for pySDC / projects / StroemungsRaum / problem_classes / NavierStokes_2D_monolithic_FEniCS.py: 100%
106 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-14 15:19 +0000
1import logging
2import os
4import dolfin as df
5import numpy as np
7from pySDC.core.problem import Problem
8from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh
11class fenics_NSE_2D_Monolithic(Problem):
12 r"""
13 Example implementing the forced two-dimensional incompressible Navier-Stokes equations with
14 time-dependent Dirichlet boundary conditions
16 .. math::
17 \frac{\partial u}{\partial t} = - u \cdot \nabla u + \nu \Delta u - \nabla p + f
18 0 = \nabla \cdot u
20 for :math:`x \in \Omega`, where the forcing term :math:`f` is defined by
22 .. math::
23 f(x, t) = (0, 0).
25 This implementation follows a monolithic approach, where velocity and pressure are
26 solved simultaneously in a coupled system using a mixed finite element formulation.
28 Boundary conditions are applied on subsets of the boundary:
29 - no-slip on channel walls and cylinder surface,
30 - a time-dependent inflow profile at the inlet,
31 - pressure condition at the outflow.
33 In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem
34 is reformulated to the *weak formulation*
36 .. math::
37 \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
38 \int_\Omega \nabla \cdot u q\,dx = 0
40 Parameters
41 ----------
42 t0 : float, optional
43 Starting time.
44 order : int, optional
45 Defines the order of the elements in the function space.
46 nu : float, optional
47 Diffusion coefficient :math:`\nu`.
48 Sol_tol : float, optional
49 Tolerance for the nonlinear solver.
51 Attributes
52 ----------
53 V : FunctionSpace
54 Defines the function space of the trial and test functions for velocity.
55 Q : FunctionSpace
56 Defines the function space of the trial and test functions for pressure.
57 W : FunctionSpace
58 Defines the mixed function space for the coupled velocity-pressure system.
59 M : scalar, vector, matrix or higher rank tensor
60 Denotes the expression :math:`\int_\Omega u_t v\,dx`.
61 Mf : scalar, vector, matrix or higher rank tensor
62 Denotes the expression :math:`\int_\Omega u v\,dx + \int_\Omega p q\,dx`.
63 g : Expression
64 The forcing term :math:`f` in the Navier-Stokes momentum equation.
65 bc : DirichletBC
66 Denotes the time-dependent Dirichlet boundary conditions.
67 bc_hom : DirichletBC
68 Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual
69 fix_bc_for_residual: boolean
70 flag to indicate that the residual requires special treatment due to boundary conditions
72 References
73 ----------
74 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
75 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
76 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
77 Wells and others. Springer (2012).
78 """
80 dtype_u = fenics_mesh
81 dtype_f = fenics_mesh
83 df.set_log_active(False)
85 def __init__(self, t0=0.0, order=2, nu=0.001, Sol_tol=1e-10):
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 spaces for future reference (Taylor-Hood)
101 P2 = df.VectorElement("P", mesh.ufl_cell(), order)
102 P1 = df.FiniteElement("P", mesh.ufl_cell(), order - 1)
103 TH = df.MixedElement([P2, P1])
104 self.W = df.FunctionSpace(mesh, TH)
105 self.V = df.FunctionSpace(mesh, P2)
106 self.Q = df.FunctionSpace(mesh, P1)
108 # print the number of DoFs for debugging purposes
109 tmp = df.Function(self.W)
110 self.logger.debug('DoFs on this level: %d', len(tmp.vector()[:]))
112 super().__init__(self.W)
113 self._makeAttributeAndRegister('t0', 'order', 'nu', 'Sol_tol', localVars=locals(), readOnly=True)
115 # Trial and test function for the Mixed FE space
116 self.u, self.p = df.TrialFunctions(self.W)
117 self.v, self.q = df.TestFunctions(self.W)
119 # velocity mass matrix
120 a_M = df.inner(self.u, self.v) * df.dx
121 self.M = df.assemble(a_M)
123 # full mass matrix
124 a_Mf = df.inner(self.u, self.v) * df.dx + df.inner(self.p, self.q) * df.dx
125 Mf = df.assemble(a_Mf)
127 # define the time-dependent inflow profile as an Expression
128 Uin = '4.0*1.5*sin(pi*t/8)*x[1]*(0.41 - x[1]) / pow(0.41, 2)'
129 self.u_in = df.Expression((Uin, '0'), pi=np.pi, t=t0, degree=self.order)
131 # define boundaries
132 inflow = 'near(x[0], 0)'
133 outflow = 'near(x[0], 2.2)'
134 walls = 'near(x[1], 0) || near(x[1], 0.41)'
135 cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'
137 # define boundary conditions
138 bc_in = df.DirichletBC(self.W.sub(0), self.u_in, inflow)
139 bc_out = df.DirichletBC(self.W.sub(1), 0, outflow)
140 bc_walls = df.DirichletBC(self.W.sub(0), (0, 0), walls)
141 bc_cylinder = df.DirichletBC(self.W.sub(0), (0, 0), cylinder)
142 self.bc = [bc_cylinder, bc_walls, bc_out, bc_in]
144 # homogeneous boundary conditions for fixing the residual
145 bc_hom_u = df.DirichletBC(self.W.sub(0), df.Constant((0, 0)), 'on_boundary')
146 bc_hom_p = df.DirichletBC(self.W.sub(1), df.Constant(0), 'on_boundary')
147 self.bc_hom = [bc_hom_u, bc_hom_p]
148 self.fix_bc_for_residual = True
150 # define measure for drag and lift computation
151 Cylinder = df.CompiledSubDomain(cylinder)
152 CylinderBoundary = df.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
153 Cylinder.mark(CylinderBoundary, 1)
154 self.dsc = df.Measure("ds", domain=mesh, subdomain_data=CylinderBoundary, subdomain_id=1)
156 # set forcing term as expression
157 self.g = df.Expression(('0', '0'), a=np.pi, b=self.nu, t=self.t0, degree=self.order)
159 # initialize XDMF files for velocity and pressure if needed
160 self.xdmffile_p = None
161 self.xdmffile_u = None
163 # set up linear solver for the inversion of the mass matrix
164 self.solver = df.LUSolver(Mf)
165 # self.solver.parameters['reuse_factorization'] = True
167 def solve_system(self, rhs, factor, u0, t):
168 r"""
169 Dolfin's nonlinear solver for :
170 .. math::
171 (u,v) + factor * (u \cdot \nabla u, v) + factor * \nu (\nabla u, \nabla v) - factor * (p, \nabla \cdot v) - factor * (g, v) - factor * (div(u), q) = (rhs_u, v) + (rhs_p, q)
173 Parameters
174 ----------
175 rhs : dtype_f
176 Right-hand side for the nonlinear system.
177 factor : float
178 Abbrev. for the node-to-node stepsize (or any other factor required).
179 u0 : dtype_u
180 Initial guess for the iterative solver (not used here so far).
181 t : float
182 Current time.
184 Returns
185 -------
186 w : dtype_u
187 Solution.
188 """
189 # introduce the coupled solution vector for velocity and pressure
190 w = self.dtype_u(u0)
191 u, p = df.split(w.values)
193 # get the SDC right-hand side
194 rhs = self.__invert_mass_matrix(rhs)
195 rhs_u, rhs_p = df.split(rhs.values)
197 # update time in boundary conditions
198 self.u_in.t = t
200 # get the forcing term
201 self.g.t = t
202 g = df.interpolate(self.g, self.V)
204 # build the variational form for the coupled system
205 F = df.dot(u, self.v) * df.dx
206 F += factor * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx
207 F += factor * self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx
208 F -= factor * df.dot(p, df.div(self.v)) * df.dx
209 F -= factor * df.dot(g, self.v) * df.dx
210 F -= factor * df.dot(df.div(u), self.q) * df.dx
211 F -= df.dot(rhs_u, self.v) * df.dx
212 F -= df.dot(rhs_p, self.q) * df.dx
214 # solve the nonlinear system using Newton's method
215 df.solve(F == 0, w.values, self.bc, solver_parameters={"newton_solver": {"absolute_tolerance": self.Sol_tol}})
217 return w
219 def eval_f(self, w, t):
220 """
221 Routine to evaluate both parts of the right-hand side of the problem.
223 Parameters
224 ----------
225 w : dtype_u
226 Current values of the numerical solution.
227 t : float
228 Current time at which the numerical solution is computed.
230 Returns
231 -------
232 f : dtype_f
233 The right-hand side.
234 """
236 f = self.dtype_f(self.W)
238 u, p = df.split(w.values)
240 # get the forcing term
241 self.g.t = t
242 g = self.dtype_f(df.interpolate(self.g, self.V))
244 F = -1.0 * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx
245 F -= self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx
246 F += df.dot(p, df.div(self.v)) * df.dx
247 F += df.dot(g.values, self.v) * df.dx
248 F += df.dot(df.div(u), self.q) * df.dx
250 df.assemble(F, tensor=f.values.vector())
252 return f
254 def apply_mass_matrix(self, w):
255 r"""
256 Routine to apply velocity mass matrix.
258 Parameters
259 ----------
260 w : dtype_u
261 Current values of the numerical solution.
263 Returns
264 -------
265 me : dtype_u
266 The product :math:`M \vec{w}`.
267 """
269 me = self.dtype_u(self.W)
270 self.M.mult(w.values.vector(), me.values.vector())
272 return me
274 def __invert_mass_matrix(self, w):
275 r"""
276 Helper routine to invert the full mass matrix Mf.
278 Parameters
279 ----------
280 w : dtype_u
281 Current values of the numerical solution.
283 Returns
284 -------
285 me : dtype_u
286 The product :math:`Mf^{-1} \vec{w}`.
287 """
289 me = self.dtype_u(self.W)
290 self.solver.solve(me.values.vector(), w.values.vector())
291 return me
293 def u_exact(self, t):
294 r"""
295 Routine to compute the exact solution at time :math:`t`.
297 Parameters
298 ----------
299 t : float
300 Time of the exact solution.
302 Returns
303 -------
304 me : dtype_u
305 Exact solution.
306 """
307 # define the exact solution
308 w = df.Function(self.W)
310 # assign the exact solution as an Expression
311 df.assign(w.sub(0), df.interpolate(df.Expression(('0.0', '0.0'), degree=self.order), self.V))
312 df.assign(w.sub(1), df.interpolate(df.Expression('0.0', degree=self.order - 1), self.Q))
314 # update time in Boundary conditions
315 self.u_in.t = t
316 [bc.apply(w.vector()) for bc in self.bc]
318 me = self.dtype_u(w, val=self.W)
320 return me
322 def fix_residual(self, res):
323 """
324 Applies homogeneous Dirichlet boundary conditions to the residual
326 Parameters
327 ----------
328 res : dtype_u
329 Residual
330 """
331 [bc.apply(res.values.vector()) for bc in self.bc_hom]
333 return None