Coverage for pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py: 79%
111 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 13:41 +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, rhs_fenics_mesh
10# noinspection PyUnusedLocal
11class fenics_vortex_2d(Problem):
12 r"""
13 2D periodic vorticity-velocity problem with FEniCS, IMEX (diffusion implicit), with the mass matrix inverted.
15 This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions
16 in :math:`[0, 1]^2`
18 .. math::
19 \frac{\partial w}{\partial t} = \nu \Delta w
21 for some parameter :math:`\nu`. In this class the problem is implemented that the spatial part is solved
22 using ``FEniCS`` [#]_. Hence, the problem is reformulated to the *weak formulation*
24 .. math::
25 \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx
27 This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.
28 The mass matrix needs inversion for this type of problem class, see the derived one for the mass-matrix version without inversion.
30 Parameters
31 ----------
32 c_nvars : tuple of int, optional
33 Number of cells of the coarse mesh in :math:`x` and :math:`y` direction, e.g. ``c_nvars=(128, 128)``.
34 family : str, optional
35 Indicates the family of elements used to create the function space
36 for the trail and test functions. The default is ``'CG'``, which are the class
37 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [#]_.
38 order : int, optional
39 Defines the order of the elements in the function space.
40 refinements : int, optional
41 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
42 nu : float, optional
43 Diffusion coefficient :math:`\nu`.
44 rho : int, optional
45 Steepness of the shear layers in the initial condition.
46 delta : float, optional
47 Amplitude of the perturbation in the initial condition.
49 Attributes
50 ----------
51 V : FunctionSpace
52 Defines the function space of the trial and test functions.
53 M : scalar, vector, matrix or higher rank tensor
54 Mass matrix for FENiCS.
55 K : scalar, vector, matrix or higher rank tensor
56 Stiffness matrix including diffusion coefficient (and correct sign).
57 dtype_u : fenics_mesh
58 FEniCS mesh data type.
59 dtype_f : rhs_fenics_mesh
60 FEniCS mesh data type with implicit and explicit parts.
62 References
63 ----------
64 .. [#] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
65 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
66 .. [#] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
67 Wells and others. Springer (2012).
68 """
70 dtype_u = fenics_mesh
71 dtype_f = rhs_fenics_mesh
73 def __init__(self, c_nvars=None, family='CG', order=4, refinements=None, nu=0.01, rho=50, delta=0.05):
74 """Initialization routine"""
76 if c_nvars is None:
77 c_nvars = (32, 32)
79 if refinements is None:
80 refinements = 1
82 # Subdomain for Periodic boundary condition
83 class PeriodicBoundary(df.SubDomain):
84 # Left boundary is "target domain" G
85 def inside(self, x, on_boundary):
86 # return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
87 return bool(
88 (df.near(x[0], 0) or df.near(x[1], 0))
89 and (not ((df.near(x[0], 0) and df.near(x[1], 1)) or (df.near(x[0], 1) and df.near(x[1], 0))))
90 and on_boundary
91 )
93 def map(self, x, y):
94 if df.near(x[0], 1) and df.near(x[1], 1):
95 y[0] = x[0] - 1.0
96 y[1] = x[1] - 1.0
97 elif df.near(x[0], 1):
98 y[0] = x[0] - 1.0
99 y[1] = x[1]
100 else: # near(x[1], 1)
101 y[0] = x[0]
102 y[1] = x[1] - 1.0
104 # set logger level for FFC and dolfin
105 logging.getLogger('FFC').setLevel(logging.WARNING)
106 logging.getLogger('UFL').setLevel(logging.WARNING)
108 # set solver and form parameters
109 df.parameters["form_compiler"]["optimize"] = True
110 df.parameters["form_compiler"]["cpp_optimize"] = True
112 # set mesh and refinement (for multilevel)
113 mesh = df.UnitSquareMesh(c_nvars[0], c_nvars[1])
114 for _ in range(refinements):
115 mesh = df.refine(mesh)
117 self.mesh = df.Mesh(mesh)
119 # define function space for future reference
120 self.V = df.FunctionSpace(mesh, family, order, constrained_domain=PeriodicBoundary())
121 tmp = df.Function(self.V)
122 print('DoFs on this level:', len(tmp.vector()[:]))
123 self.fix_bc_for_residual = False
125 # invoke super init, passing number of dofs, dtype_u and dtype_f
126 super(fenics_vortex_2d, self).__init__(self.V)
127 self._makeAttributeAndRegister(
128 'c_nvars', 'family', 'order', 'refinements', 'nu', 'rho', 'delta', localVars=locals(), readOnly=True
129 )
131 w = df.TrialFunction(self.V)
132 v = df.TestFunction(self.V)
134 # Stiffness term (diffusion)
135 a_K = df.inner(df.nabla_grad(w), df.nabla_grad(v)) * df.dx
137 # Mass term
138 a_M = w * v * df.dx
140 self.M = df.assemble(a_M)
141 self.K = df.assemble(a_K)
143 def solve_system(self, rhs, factor, u0, t):
144 r"""
145 Dolfin's linear solver for :math:`(M - factor \cdot A)\vec{u} = \vec{rhs}`.
147 Parameters
148 ----------
149 rhs : dtype_f
150 Right-hand side for the nonlinear system.
151 factor : float
152 Abbrev. for the node-to-node stepsize (or any other factor required).
153 u0 : dtype_u
154 Initial guess for the iterative solver (not used here so far).
155 t : float
156 Current time.
158 Returns
159 -------
160 u : dtype_u
161 The solution as mesh.
162 """
164 A = self.M + self.nu * factor * self.K
165 b = self.apply_mass_matrix(rhs)
166 u = self.dtype_u(u0)
167 df.solve(A, u.values.vector(), b.values.vector())
169 return u
171 def __eval_fexpl(self, u, t):
172 """
173 Helper routine to evaluate the explicit part of the right-hand side.
175 Parameters
176 ----------
177 u : dtype_u
178 Current values of the numerical solution.
179 t : float
180 Current time at which the numerical solution is computed.
182 Returns
183 -------
184 fexpl : dtype_u
185 Explicit part of the right-hand side.
186 """
188 b = self.apply_mass_matrix(u)
189 psi = self.dtype_u(self.V)
190 df.solve(self.K, psi.values.vector(), b.values.vector())
192 fexpl = self.dtype_u(self.V)
193 fexpl.values = df.project(
194 df.Dx(psi.values, 1) * df.Dx(u.values, 0) - df.Dx(psi.values, 0) * df.Dx(u.values, 1), self.V
195 )
197 return fexpl
199 def __eval_fimpl(self, u, t):
200 """
201 Helper routine to evaluate the implicit part of the right-hand side.
203 Parameters
204 ----------
205 u : dtype_u
206 Current values of the numerical solution.
207 t : float
208 Current time at which the numerical solution is computed.
210 Returns
211 -------
212 fimpl : dtype_u
213 Implicit part of the right-hand side.
214 """
216 A = -self.nu * self.K
217 fimpl = self.dtype_u(self.V)
218 A.mult(u.values.vector(), fimpl.values.vector())
219 fimpl = self.__invert_mass_matrix(fimpl)
221 return fimpl
223 def eval_f(self, u, t):
224 """
225 Routine to evaluate both parts of the right-hand side.
227 Parameters
228 ----------
229 u : dtype_u
230 Current values of the numerical solution.
231 t : float
232 Current time at which the numerical solution is computed.
234 Returns
235 -------
236 f : dtype_f
237 The right-hand side divided into two parts.
238 """
240 f = self.dtype_f(self.V)
241 f.impl = self.__eval_fimpl(u, t)
242 f.expl = self.__eval_fexpl(u, t)
243 return f
245 def apply_mass_matrix(self, u):
246 r"""
247 Routine to apply mass matrix.
249 Parameters
250 ----------
251 u : dtype_u
252 Current values of the numerical solution.
254 Returns
255 -------
256 me : dtype_u
257 The product :math:`M\vec{u}`.
258 """
260 me = self.dtype_u(self.V)
261 self.M.mult(u.values.vector(), me.values.vector())
263 return me
265 def __invert_mass_matrix(self, u):
266 """
267 Helper routine to invert mass matrix.
269 Parameters
270 ----------
271 u : dtype_u
272 Current values of the numerical solution.
274 Returns
275 -------
276 me : dtype_u
277 The product :math:`M^{-1} \vec{u}`.
278 """
280 me = self.dtype_u(self.V)
281 df.solve(self.M, me.values.vector(), u.values.vector())
282 return me
284 def u_exact(self, t):
285 r"""
286 Routine to compute the exact solution at time :math:`t`.
288 Parameters
289 ----------
290 t : float
291 Time of the exact solution.
293 Returns
294 -------
295 me : dtype_u
296 The exact solution.
297 """
298 assert t == 0, 'ERROR: u_exact only valid for t=0'
300 w = df.Expression(
301 'r*(1-pow(tanh(r*((0.75-4) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-4))),2)) - \
302 r*(1-pow(tanh(r*((0.75-3) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-3))),2)) - \
303 r*(1-pow(tanh(r*((0.75-2) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-2))),2)) - \
304 r*(1-pow(tanh(r*((0.75-1) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-1))),2)) - \
305 r*(1-pow(tanh(r*((0.75-0) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-0))),2)) - \
306 r*(1-pow(tanh(r*((0.75+1) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+1))),2)) - \
307 r*(1-pow(tanh(r*((0.75+2) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+2))),2)) - \
308 r*(1-pow(tanh(r*((0.75+3) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+3))),2)) - \
309 r*(1-pow(tanh(r*((0.75+4) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+4))),2)) - \
310 d*2*a*cos(2*a*(x[0]+0.25))',
311 d=self.delta,
312 r=self.rho,
313 a=np.pi,
314 degree=self.order,
315 )
317 me = self.dtype_u(self.V)
318 me.values = df.interpolate(w, self.V)
320 # df.plot(me.values)
321 # df.interactive()
322 # exit()
324 return me
327class fenics_vortex_2d_mass(fenics_vortex_2d):
328 r"""
329 2D periodic vorticity-velocity problem with FEniCS, IMEX (diffusion implicit), with the mass matrix applied.
331 This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions
332 in :math:`[0, 1]^2`
334 .. math::
335 \frac{\partial w}{\partial t} = \nu \Delta w
337 for some parameter :math:`\nu`. In this class the problem is implemented that the spatial part is solved
338 using ``FEniCS`` [#]_. Hence, the problem is reformulated to the *weak formulation*
340 .. math::
341 \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx
343 This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.
344 No mass matrix inversion is needed here, i.e. using this problem class requires the imex_1st_order_mass sweeper.
346 Parameters
347 ----------
348 c_nvars : tuple of int, optional
349 Number of cells of the coarse mesh in :math:`x` and :math:`y` direction, e.g. ``c_nvars=(128, 128)``.
350 family : str, optional
351 Indicates the family of elements used to create the function space
352 for the trail and test functions. The default is ``'CG'``, which are the class
353 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [#]_.
354 order : int, optional
355 Defines the order of the elements in the function space.
356 refinements : int, optional
357 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
358 nu : float, optional
359 Diffusion coefficient :math:`\nu`.
360 rho : int, optional
361 Steepness of the shear layers in the initial condition.
362 delta : float, optional
363 Amplitude of the perturbation in the initial condition.
365 Attributes
366 ----------
367 V : FunctionSpace
368 Defines the function space of the trial and test functions.
369 M : scalar, vector, matrix or higher rank tensor
370 Mass matrix for FENiCS.
371 K : scalar, vector, matrix or higher rank tensor
372 Stiffness matrix including diffusion coefficient (and correct sign).
373 dtype_u : fenics_mesh
374 FEniCS mesh data type.
375 dtype_f : rhs_fenics_mesh
376 FEniCS mesh data type with implicit and explicit parts.
378 References
379 ----------
380 .. [#] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
381 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
382 .. [#] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
383 Wells and others. Springer (2012).
384 """
386 def solve_system(self, rhs, factor, u0, t):
387 r"""
388 Dolfin's linear solver for :math:`(M - factor \cdot A)\vec{u} = \vec{rhs}`.
390 Parameters
391 ----------
392 rhs : dtype_f
393 Right-hand side for the nonlinear system.
394 factor : float
395 Abbrev. for the node-to-node stepsize (or any other factor required).
396 u0 : dtype_u
397 Initial guess for the iterative solver (not used here so far).
398 t : float
399 Current time.
401 Returns
402 -------
403 u : dtype_u
404 The solution as mesh.
405 """
407 A = self.M + self.nu * factor * self.K
408 b = self.dtype_u(rhs)
409 u = self.dtype_u(u0)
410 df.solve(A, u.values.vector(), b.values.vector())
412 return u
414 def __eval_fexpl(self, u, t):
415 """
416 Helper routine to evaluate the explicit part of the right-hand side.
418 Parameters
419 ----------
420 u : dtype_u
421 Current values of the numerical solution.
422 t : float
423 Current time at which the numerical solution is computed.
425 Returns
426 -------
427 fexpl : dtype_u
428 Explicit part of the right-hand side.
429 """
431 b = self.apply_mass_matrix(u)
432 psi = self.dtype_u(self.V)
433 df.solve(self.K, psi.values.vector(), b.values.vector())
435 # Assemble the load vector of the advection term directly. The previous version projected the
436 # expression onto V and multiplied the result by M again, which is M M^-1 = identity: it paid
437 # for a mass solve per right-hand side evaluation and produced the same vector. It is also the
438 # form FFC miscompiles under PFASST, emitting a duplicate J_c1 declaration.
439 # u may have been handed over from another step, and fenics_mesh copies keep the function
440 # space they were created on. Under PFASST every step owns a separate (but identical) mesh, so
441 # a form built straight from u.values would mix two meshes and UFL refuses to order them.
442 # Copying the coefficients onto this problem's own space keeps the form single-domain.
443 w = df.Function(self.V)
444 w.vector()[:] = u.values.vector()[:]
446 v = df.TestFunction(self.V)
447 expr = df.Dx(psi.values, 1) * df.Dx(w, 0) - df.Dx(psi.values, 0) * df.Dx(w, 1)
449 fexpl = self.dtype_u(self.V)
450 fexpl.values = df.Function(self.V, df.assemble(expr * v * df.dx))
452 return fexpl
454 def __eval_fimpl(self, u, t):
455 """
456 Helper routine to evaluate the implicit part of the right-hand side.
458 Parameters
459 ----------
460 u : dtype_u
461 Current values of the numerical solution.
462 t : float
463 Current time at which the numerical solution is computed.
465 Returns
466 -------
467 fimpl : dtype_u
468 Implicit part of the right-hand side.
469 """
471 A = -self.nu * self.K
472 fimpl = self.dtype_u(self.V)
473 A.mult(u.values.vector(), fimpl.values.vector())
475 return fimpl
477 def eval_f(self, u, t):
478 """
479 Routine to evaluate both parts of the right-hand side.
481 Note: Need to add this here, because otherwise the parent class will call the "local" functions __eval_*
482 and not the ones of the child class.
484 Parameters
485 ----------
486 u : dtype_u
487 Current values of the numerical solution.
488 t : float
489 Current time at which the numerical solution is computed.
491 Returns
492 -------
493 f : dtype_f
494 The right-hand side divided into two parts.
495 """
497 f = self.dtype_f(self.V)
498 f.impl = self.__eval_fimpl(u, t)
499 f.expl = self.__eval_fexpl(u, t)
500 return f