Coverage for pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py: 100%
136 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1import logging
3import dolfin as df
4import numpy as np
6from pySDC.core.problem import Problem
7from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh
8from pySDC.projects.StroemungsRaum.problem_classes.newton_step import NewtonStep
11class _PeriodicX(df.SubDomain):
12 """
13 Identifies the right boundary :math:`x = 0.5` with the left boundary :math:`x = -0.5`.
15 Passed to ``FunctionSpace`` as ``constrained_domain``, so periodicity is handled by
16 the dof map itself and the periodic dofs never enter the linear system.
17 """
19 def inside(self, x, on_boundary):
20 return bool(df.near(x[0], -0.5) and on_boundary)
22 def map(self, x, y):
23 y[0] = x[0] - 1.0
24 y[1] = x[1]
27class fenics_NSE_2D_TaylorGreen(Problem):
28 r"""
29 Forced two-dimensional incompressible Navier-Stokes equations on :math:`\Omega = [-0.5, 0.5]^2`,
30 set up to expose the order reduction caused by time-dependent Dirichlet boundary conditions.
32 .. math::
33 \frac{\partial u}{\partial t} = - u \cdot \nabla u + \nu \Delta u - \nabla p + g,
34 \qquad \nabla \cdot u = 0
36 The forcing :math:`g` is manufactured from the analytical solution
38 .. math::
39 u(x, y, t) &= 1 - e^{-8\pi^2\nu t}\sin(2\pi(x - t))\sin(\pi y)\cos(\pi y) \\
40 v(x, y, t) &= - e^{-8\pi^2\nu t}\cos(2\pi(x - t))\cos^2(\pi y) \\
41 p(x, y, t) &= 1 + \frac{4}{17}e^{-16\pi^2\nu t}\cos(4\pi(x - t))\cos(\pi y)
43 which is divergence free and exactly one-periodic in :math:`x`. On :math:`y = \pm 0.5` it
44 collapses to the constants :math:`u = (1, 0)`, :math:`p = 1`, so the top and bottom boundary
45 data is time-independent.
47 Because the solution is genuinely periodic in :math:`x`, the *same* exact solution satisfies
48 both variants selected by ``periodic``:
50 - ``periodic=False``: time-dependent Dirichlet conditions on :math:`x = \pm 0.5`,
51 - ``periodic=True``: periodic conditions on :math:`x = \pm 0.5`.
53 The only difference between the two runs is therefore the presence of time-dependent boundary
54 data, which is what isolates the order reduction.
56 On :math:`x = \pm 0.5` the *pressure* is prescribed from the exact solution as well. That is
57 not neutral: it acts as a partial lifting of the algebraic constraint, and it lifts the
58 observed pressure order from :math:`M` to :math:`M+1`. Dropping it gives order :math:`M`
59 and a roughly 20 times larger error, so the gap measured here understates what a setup
60 without prescribed boundary pressure would show.
62 Setting ``differentiated_bc`` imposes the time-dependent data in differentiated form and
63 recovers most of the lost order, following the remedy explored for a time-dependent
64 *constraint* in pull request #641. It requires the ``generic_implicit_mass_diffbc`` sweeper.
65 See :meth:`prepare_step` for the construction and its measured effect.
67 Note that the number of collocation nodes decides whether anything can be seen at all.
68 RADAU-RIGHT with :math:`M` nodes has design order :math:`2M-1` and falls back to the stiff
69 order :math:`M+1` in the presence of time-dependent boundary data, so the gap on offer is
70 :math:`M-2`: **zero for M = 2**, where both are 3. Use :math:`M \geq 4`; at :math:`M = 4`
71 the measured pressure orders are 7 (periodic) against 5 (Dirichlet).
73 The problem is discretized in space with Taylor-Hood elements on a mixed velocity-pressure
74 space and solved monolithically, so the semi-discrete system is the differential-algebraic
75 system :math:`M \dot{w} = f(w, t)` with the singular mass matrix :math:`M = \mathrm{diag}(M_v, 0)`.
76 It therefore requires ``generic_implicit_mass`` as sweeper, which applies :math:`M` where needed
77 instead of inverting it.
79 Parameters
80 ----------
81 nelems : int, optional
82 Number of elements per spatial direction.
83 t0 : float, optional
84 Starting time.
85 order : int, optional
86 Polynomial degree of the velocity space; the pressure space uses ``order - 1``.
87 nu : float, optional
88 Kinematic viscosity :math:`\nu`.
89 periodic : bool, optional
90 Use periodic instead of time-dependent Dirichlet conditions on :math:`x = \pm 0.5`.
91 differentiated_bc : bool, optional
92 Impose the time-dependent boundary data in differentiated form; needs ``periodic=False``
93 and the ``generic_implicit_mass_diffbc`` sweeper.
94 Sol_tol : float, optional
95 Absolute tolerance for the Newton solver.
97 Attributes
98 ----------
99 V : FunctionSpace
100 Velocity space.
101 Q : FunctionSpace
102 Pressure space.
103 W : FunctionSpace
104 Mixed velocity-pressure space.
105 M : Matrix
106 The velocity mass matrix :math:`\mathrm{diag}(M_v, 0)` on the mixed space.
107 g : Expression
108 Manufactured forcing term.
109 bc : list of DirichletBC
110 Dirichlet boundary conditions, time-dependent unless ``periodic``.
111 bc_hom : list of DirichletBC
112 Homogeneous conditions on the Dirichlet part of the boundary only, used to fix the residual.
113 fix_bc_for_residual : bool
114 Flag indicating that the residual requires special treatment due to boundary conditions.
116 References
117 ----------
118 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
119 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
120 """
122 dtype_u = fenics_mesh
123 dtype_f = fenics_mesh
125 df.set_log_active(False)
127 def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, differentiated_bc=False, Sol_tol=1e-10):
129 # set logger level for FFC and dolfin
130 logging.getLogger('FFC').setLevel(logging.WARNING)
131 logging.getLogger('UFL').setLevel(logging.WARNING)
133 # set solver and form parameters
134 df.parameters["form_compiler"]["optimize"] = True
135 df.parameters["form_compiler"]["cpp_optimize"] = True
137 mesh = df.RectangleMesh(df.Point(-0.5, -0.5), df.Point(0.5, 0.5), nelems, nelems)
139 # define function spaces (Taylor-Hood); periodicity is baked into the dof map
140 P2 = df.VectorElement("P", mesh.ufl_cell(), order)
141 P1 = df.FiniteElement("P", mesh.ufl_cell(), order - 1)
142 constraint = _PeriodicX() if periodic else None
143 self.W = df.FunctionSpace(mesh, df.MixedElement([P2, P1]), constrained_domain=constraint)
144 self.V = df.FunctionSpace(mesh, P2, constrained_domain=constraint)
145 self.Q = df.FunctionSpace(mesh, P1, constrained_domain=constraint)
147 super().__init__(self.W)
148 self._makeAttributeAndRegister(
149 'nelems',
150 't0',
151 'order',
152 'nu',
153 'periodic',
154 'differentiated_bc',
155 'Sol_tol',
156 localVars=locals(),
157 readOnly=True,
158 )
160 self.logger.debug('DoFs on this level: %d', self.W.dim())
162 # trial and test functions on the mixed space
163 self.u, self.p = df.TrialFunctions(self.W)
164 self.v, self.q = df.TestFunctions(self.W)
166 # velocity mass matrix on the mixed space, i.e. diag(M_v, 0)
167 self.M = df.assemble(df.inner(self.u, self.v) * df.dx)
169 # manufactured solution and the forcing term derived from it
170 self.u_ex = df.Expression(
171 (
172 '1.0 - exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])',
173 '-exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])',
174 ),
175 pi=np.pi,
176 nu=nu,
177 t=t0,
178 degree=order + 2,
179 )
180 self.p_ex = df.Expression(
181 '1.0 + (4.0/17.0)*exp(-16*pi*pi*nu*t)*cos(4*pi*(x[0] - t))*cos(pi*x[1])',
182 pi=np.pi,
183 nu=nu,
184 t=t0,
185 degree=order + 2,
186 )
187 self.g = df.Expression(
188 (
189 'pi/34.0*exp(-16*pi*pi*nu*t)*sin(4*pi*(t - x[0]))*cos(pi*x[1])*(32.0 - 17.0*cos(pi*x[1]))',
190 '2*pi*pi*nu*exp(-8*pi*pi*nu*t)*cos(2*pi*(t - x[0]))'
191 ' - pi*exp(-16*pi*pi*nu*t)*sin(pi*x[1])'
192 '*(2.0*pow(cos(pi*x[1]), 3) + 4.0/17.0*cos(4*pi*(t - x[0])))',
193 ),
194 pi=np.pi,
195 nu=nu,
196 t=t0,
197 degree=order + 2,
198 )
200 # on y = +-0.5 the exact solution is constant in space and time
201 top_bottom = 'near(x[1], -0.5) || near(x[1], 0.5)'
202 self.left_right = 'near(x[0], -0.5) || near(x[0], 0.5)'
203 self.bc_fixed = [
204 df.DirichletBC(self.W.sub(0), df.Constant((1.0, 0.0)), top_bottom),
205 df.DirichletBC(self.W.sub(1), df.Constant(1.0), top_bottom),
206 ]
207 self.bc = list(self.bc_fixed)
208 if not periodic:
209 self.bc += [
210 df.DirichletBC(self.W.sub(0), self.u_ex, self.left_right),
211 df.DirichletBC(self.W.sub(1), self.p_ex, self.left_right),
212 ]
214 # boundary conditions per collocation node, filled in by prepare_step
215 self._node_times = None
216 self._node_bcs = None
217 if differentiated_bc:
218 if periodic:
219 raise ValueError('differentiated_bc has no effect without time-dependent boundary data')
220 self.u_dot, self.p_dot = self._boundary_derivatives(nu, order, t0)
222 # the residual is meaningless where the solution is prescribed, but only there: with
223 # periodicity the dofs on x = +-0.5 are unknowns and their residual has to be kept
224 dirichlet = top_bottom if periodic else 'on_boundary'
225 self.bc_hom = [
226 df.DirichletBC(self.W.sub(0), df.Constant((0.0, 0.0)), dirichlet),
227 df.DirichletBC(self.W.sub(1), df.Constant(0.0), dirichlet),
228 ]
229 self.fix_bc_for_residual = True
231 # residual form for a single node-to-node step, assembled once; `factor` and the
232 # boundary/forcing expressions carry the time dependence
233 self.factor = df.Constant(0.0)
234 self.w = df.Function(self.W)
235 u, p = df.split(self.w)
237 F = df.dot(u, self.v) * df.dx
238 F += self.factor * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx
239 F += self.factor * self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx
240 F -= self.factor * df.dot(p, df.div(self.v)) * df.dx
241 F -= self.factor * df.dot(self.g, self.v) * df.dx
242 F -= self.factor * df.dot(df.div(u), self.q) * df.dx
244 self.step = NewtonStep(F, df.derivative(F, self.w))
245 self.newton = df.NewtonSolver()
246 self.newton.parameters['absolute_tolerance'] = Sol_tol
248 @staticmethod
249 def _boundary_derivatives(nu, order, t0):
250 r"""
251 Time derivatives of the boundary data, needed to impose it in differentiated form.
253 Returns
254 -------
255 u_dot, p_dot : Expression
256 :math:`\partial_t u` and :math:`\partial_t p` of the manufactured solution.
257 """
258 kwargs = dict(pi=np.pi, nu=nu, t=t0, degree=order + 2)
259 u_dot = df.Expression(
260 (
261 '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])'
262 ' + 2*pi*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])',
263 '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])'
264 ' - 2*pi*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])',
265 ),
266 **kwargs,
267 )
268 p_dot = df.Expression(
269 '(4.0/17.0)*cos(pi*x[1])*('
270 '-16*pi*pi*nu*exp(-16*pi*pi*nu*t)*cos(4*pi*(x[0] - t))'
271 ' + 4*pi*exp(-16*pi*pi*nu*t)*sin(4*pi*(x[0] - t)))',
272 **kwargs,
273 )
274 return u_dot, p_dot
276 def prepare_step(self, t0, dt, coll):
277 r"""
278 Build the differentiated boundary conditions for every collocation node of a step.
280 Rather than evaluating the boundary data pointwise at the node, :math:`u_B(\tau_m) =
281 g(\tau_m)`, the condition is imposed on the *derivative* and the stage value recovered
282 by the collocation quadrature,
284 .. math::
285 u_B(\tau_m) = g(t_0) + \Delta t \sum_j Q_{mj}\, \dot{g}(\tau_j).
287 The two differ by the quadrature error :math:`O(\Delta t^{M+1})`, but the second is
288 consistent with the collocation polynomial instead of pointwise exact, which is what
289 recovers the order lost to time-dependent boundary data.
291 Measured at :math:`M = 4`, ``nelems=24``, ``nu=0.1``, orders and errors in the pressure
292 taken from consecutive step sizes:
294 ========================= ========== =====================
295 boundary condition order error at ``dt = 0.1``
296 ========================= ========== =====================
297 periodic (best possible) 6.32 7.4e-08
298 pointwise 5.74 9.7e-07
299 differentiated 6.30 1.3e-07
300 ========================= ========== =====================
302 The remaining factor of 1.8 against the periodic case is a constant, not a rate. Note
303 that the observed orders here are pre-asymptotic -- the periodic reference does not
304 reach its design order 7 either -- so these numbers show that the remedy works, not
305 that it restores exactly :math:`2M-1`.
307 Called once per step by :class:`generic_implicit_mass_diffbc`; ``solve_system`` then
308 picks the condition belonging to the node it is asked to solve at.
310 Parameters
311 ----------
312 t0 : float
313 Left end of the step.
314 dt : float
315 Step size.
316 coll : pySDC.core.collocation.CollBase
317 Collocation rule of the sweeper, supplying the nodes and the matrix Q.
318 """
319 if not self.differentiated_bc:
320 raise RuntimeError(
321 'prepare_step builds the differentiated boundary conditions, which this problem '
322 'was not set up for; use generic_implicit_mass or pass differentiated_bc=True'
323 )
325 M = coll.num_nodes
326 Q = coll.Qmat[1:, 1:]
327 self._node_times = t0 + dt * np.asarray(coll.nodes)
329 u_rate, p_rate = [], []
330 for j in range(M):
331 self.u_dot.t = self._node_times[j]
332 self.p_dot.t = self._node_times[j]
333 u_rate.append(df.interpolate(self.u_dot, self.V))
334 p_rate.append(df.interpolate(self.p_dot, self.Q))
336 self.u_ex.t = t0
337 self.p_ex.t = t0
338 u_base = df.interpolate(self.u_ex, self.V)
339 p_base = df.interpolate(self.p_ex, self.Q)
341 self._node_bcs = []
342 for m in range(M):
343 gu, gp = df.Function(self.V), df.Function(self.Q)
344 gu.assign(u_base)
345 gp.assign(p_base)
346 for j in range(M):
347 gu.vector().axpy(dt * Q[m, j], u_rate[j].vector())
348 gp.vector().axpy(dt * Q[m, j], p_rate[j].vector())
349 self._node_bcs.append(
350 self.bc_fixed
351 + [
352 df.DirichletBC(self.W.sub(0), gu, self.left_right),
353 df.DirichletBC(self.W.sub(1), gp, self.left_right),
354 ]
355 )
357 def solve_system(self, rhs, factor, u0, t):
358 r"""
359 Newton solver for :math:`M w + factor \cdot N(w, t) = rhs`, where :math:`N` collects the
360 convective, viscous, pressure and divergence terms and ``rhs`` is the mass-weighted
361 right-hand side assembled by the sweeper.
363 Parameters
364 ----------
365 rhs : dtype_f
366 Right-hand side for the nonlinear system.
367 factor : float
368 Abbrev. for the node-to-node stepsize (or any other factor required).
369 u0 : dtype_u
370 Initial guess for the iterative solver.
371 t : float
372 Current time.
374 Returns
375 -------
376 w : dtype_u
377 Solution.
378 """
379 self.factor.assign(factor)
380 self.u_ex.t = t
381 self.p_ex.t = t
382 self.g.t = t
384 if self.differentiated_bc:
385 if self._node_bcs is None:
386 raise RuntimeError(
387 'differentiated_bc requires the generic_implicit_mass_diffbc sweeper, '
388 'which calls prepare_step once per step'
389 )
390 node = np.flatnonzero(self._node_times == t)
391 if node.size != 1:
392 raise RuntimeError(
393 f'no collocation node of the prepared step is at t = {t}; the prepared step '
394 f'covers {self._node_times}'
395 )
396 self.bc = self._node_bcs[node[0]]
398 self.w.vector()[:] = u0.values.vector()[:]
399 self.step.rhs = rhs.values.vector()
400 self.step.bcs = self.bc
401 self.newton.solve(self.step, self.w.vector())
403 me = self.dtype_u(self.W)
404 me.values.vector()[:] = self.w.vector()[:]
405 return me
407 def eval_f(self, w, t):
408 r"""
409 Routine to evaluate the right-hand side of the problem in weak form, i.e. *without*
410 applying :math:`M^{-1}`.
412 Parameters
413 ----------
414 w : dtype_u
415 Current values of the numerical solution.
416 t : float
417 Current time at which the numerical solution is computed.
419 Returns
420 -------
421 f : dtype_f
422 The right-hand side.
423 """
424 u, p = df.split(w.values)
425 self.g.t = t
427 F = -df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx
428 F -= self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx
429 F += df.dot(p, df.div(self.v)) * df.dx
430 F += df.dot(self.g, self.v) * df.dx
431 F += df.dot(df.div(u), self.q) * df.dx
433 f = self.dtype_f(self.W)
434 df.assemble(F, tensor=f.values.vector())
435 return f
437 def apply_mass_matrix(self, w):
438 r"""
439 Routine to apply the velocity mass matrix.
441 Parameters
442 ----------
443 w : dtype_u
444 Current values of the numerical solution.
446 Returns
447 -------
448 me : dtype_u
449 The product :math:`M \vec{w}`.
450 """
451 me = self.dtype_u(self.W)
452 self.M.mult(w.values.vector(), me.values.vector())
453 return me
455 def u_exact(self, t):
456 r"""
457 Routine to compute the exact solution at time :math:`t`.
459 Parameters
460 ----------
461 t : float
462 Time of the exact solution.
464 Returns
465 -------
466 me : dtype_u
467 Exact solution.
468 """
469 self.u_ex.t = t
470 self.p_ex.t = t
472 me = self.dtype_u(self.W)
473 df.assign(me.values.sub(0), df.interpolate(self.u_ex, self.V))
474 df.assign(me.values.sub(1), df.interpolate(self.p_ex, self.Q))
475 return me
477 def fix_residual(self, res):
478 """
479 Applies homogeneous Dirichlet boundary conditions to the residual, on the Dirichlet part
480 of the boundary only.
482 Parameters
483 ----------
484 res : dtype_u
485 Residual.
486 """
487 for bc in self.bc_hom:
488 bc.apply(res.values.vector())
489 return None