Coverage for pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py: 0%
108 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +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 This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions
14 in :math:`[0, 1]^2`
16 .. math::
17 \frac{\partial w}{\partial t} = \nu \Delta w
19 for some parameter :math:`\nu`. In this class the problem is implemented that the spatial part is solved
20 using ``FEniCS`` [1]_. Hence, the problem is reformulated to the *weak formulation*
22 .. math::
23 \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx
25 This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.
26 The mass matrix needs inversion for this type of problem class, see the derived one for the mass-matrix version without inversion.
28 Parameters
29 ----------
30 c_nvars : List of int tuple, optional
31 Spatial resolution, i.e., numbers of degrees of freedom in space, e.g. ``c_nvars=[(128, 128)]``.
32 family : str, optional
33 Indicates the family of elements used to create the function space
34 for the trail and test functions. The default is ``'CG'``, which are the class
35 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
36 order : int, optional
37 Defines the order of the elements in the function space.
38 refinements : int, optional
39 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
40 nu : float, optional
41 Diffusion coefficient :math:`\nu`.
42 rho : int, optional
43 Problem parameter.
44 delta : float, optional
45 Problem parameter.
47 Attributes
48 ----------
49 V : FunctionSpace
50 Defines the function space of the trial and test functions.
51 M : scalar, vector, matrix or higher rank tensor
52 Mass matrix for FENiCS.
53 K : scalar, vector, matrix or higher rank tensor
54 Stiffness matrix including diffusion coefficient (and correct sign).
56 References
57 ----------
58 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
59 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
60 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
61 Wells and others. Springer (2012).
62 """
64 dtype_u = fenics_mesh
65 dtype_f = rhs_fenics_mesh
67 def __init__(self, c_nvars=None, family='CG', order=4, refinements=None, nu=0.01, rho=50, delta=0.05):
68 """
69 Initialization routine
71 Args:
72 problem_params (dict): custom parameters for the example
73 dtype_u: FEniCS mesh data type (will be passed to parent class)
74 dtype_f: FEniCS mesh data data type with implicit and explicit parts (will be passed to parent class)
75 """
77 if c_nvars is None:
78 c_nvars = [(32, 32)]
80 if refinements is None:
81 refinements = 1
83 # Subdomain for Periodic boundary condition
84 class PeriodicBoundary(df.SubDomain):
85 # Left boundary is "target domain" G
86 def inside(self, x, on_boundary):
87 # return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
88 return bool(
89 (df.near(x[0], 0) or df.near(x[1], 0))
90 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))))
91 and on_boundary
92 )
94 def map(self, x, y):
95 if df.near(x[0], 1) and df.near(x[1], 1):
96 y[0] = x[0] - 1.0
97 y[1] = x[1] - 1.0
98 elif df.near(x[0], 1):
99 y[0] = x[0] - 1.0
100 y[1] = x[1]
101 else: # near(x[1], 1)
102 y[0] = x[0]
103 y[1] = x[1] - 1.0
105 # set logger level for FFC and dolfin
106 logging.getLogger('FFC').setLevel(logging.WARNING)
107 logging.getLogger('UFL').setLevel(logging.WARNING)
109 # set solver and form parameters
110 df.parameters["form_compiler"]["optimize"] = True
111 df.parameters["form_compiler"]["cpp_optimize"] = True
113 # set mesh and refinement (for multilevel)
114 mesh = df.UnitSquareMesh(c_nvars[0], c_nvars[1])
115 for _ in range(refinements):
116 mesh = df.refine(mesh)
118 self.mesh = df.Mesh(mesh)
120 # define function space for future reference
121 self.V = df.FunctionSpace(mesh, family, order, constrained_domain=PeriodicBoundary())
122 tmp = df.Function(self.V)
123 print('DoFs on this level:', len(tmp.vector()[:]))
124 self.fix_bc_for_residual = False
126 # invoke super init, passing number of dofs, dtype_u and dtype_f
127 super(fenics_vortex_2d, self).__init__(self.V)
128 self._makeAttributeAndRegister(
129 'c_nvars', 'family', 'order', 'refinements', 'nu', 'rho', 'delta', localVars=locals(), readOnly=True
130 )
132 w = df.TrialFunction(self.V)
133 v = df.TestFunction(self.V)
135 # Stiffness term (diffusion)
136 a_K = df.inner(df.nabla_grad(w), df.nabla_grad(v)) * df.dx
138 # Mass term
139 a_M = w * v * df.dx
141 self.M = df.assemble(a_M)
142 self.K = df.assemble(a_K)
144 def solve_system(self, rhs, factor, u0, t):
145 r"""
146 Dolfin's linear solver for :math:`(M - factor \cdot A)\vec{u} = \vec{rhs}`.
148 Parameters
149 ----------
150 rhs : dtype_f
151 Right-hand side for the nonlinear system.
152 factor : float
153 Abbrev. for the node-to-node stepsize (or any other factor required).
154 u0 : dtype_u
155 Initial guess for the iterative solver (not used here so far).
156 t : float
157 Current time.
159 Returns
160 -------
161 u : dtype_u
162 The solution as mesh.
163 """
165 A = self.M + self.nu * factor * self.K
166 b = self.apply_mass_matrix(rhs)
167 u = self.dtype_u(u0)
168 df.solve(A, u.values.vector(), b.values.vector())
170 return u
172 def __eval_fexpl(self, u, t):
173 """
174 Helper routine to evaluate the explicit part of the right-hand side.
176 Parameters
177 ----------
178 u : dtype_u
179 Current values of the numerical solution.
180 t : float
181 Current time at which the numerical solution is computed.
183 Returns
184 -------
185 fexpl : dtype_u
186 Explicit part of the right-hand side.
187 """
189 b = self.apply_mass_matrix(u)
190 psi = self.dtype_u(self.V)
191 df.solve(self.K, psi.values.vector(), b.values.vector())
193 fexpl = self.dtype_u(self.V)
194 fexpl.values = df.project(
195 df.Dx(psi.values, 1) * df.Dx(u.values, 0) - df.Dx(psi.values, 0) * df.Dx(u.values, 1), self.V
196 )
198 return fexpl
200 def __eval_fimpl(self, u, t):
201 """
202 Helper routine to evaluate the implicit part of the right-hand side.
204 Parameters
205 ----------
206 u : dtype_u
207 Current values of the numerical solution.
208 t : float
209 Current time at which the numerical solution is computed.
211 Returns
212 -------
213 fimpl : dtype_u
214 Implicit part of the right-hand side.
215 """
217 A = -self.nu * self.K
218 fimpl = self.dtype_u(self.V)
219 A.mult(u.values.vector(), fimpl.values.vector())
220 fimpl = self.__invert_mass_matrix(fimpl)
222 return fimpl
224 def eval_f(self, u, t):
225 """
226 Routine to evaluate both parts of the right-hand side.
228 Parameters
229 ----------
230 u : dtype_u
231 Current values of the numerical solution.
232 t : float
233 Current time at which the numerical solution is computed.
235 Returns
236 -------
237 f : dtype_f
238 The right-hand side divided into two parts.
239 """
241 f = self.dtype_f(self.V)
242 f.impl = self.__eval_fimpl(u, t)
243 f.expl = self.__eval_fexpl(u, t)
244 return f
246 def apply_mass_matrix(self, u):
247 r"""
248 Routine to apply mass matrix.
250 Parameters
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 This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions
330 in :math:`[0, 1]^2`
332 .. math::
333 \frac{\partial w}{\partial t} = \nu \Delta w
335 for some parameter :math:`\nu`. In this class the problem is implemented that the spatial part is solved
336 using ``FEniCS`` [1]_. Hence, the problem is reformulated to the *weak formulation*
338 .. math::
339 \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx
341 This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.
342 No mass matrix inversion is needed here, i.e. using this problem class requires the imex_1st_order_mass sweeper.
344 Parameters
345 ----------
346 c_nvars : List of int tuple, optional
347 Spatial resolution, i.e., numbers of degrees of freedom in space, e.g. ``c_nvars=[(128, 128)]``.
348 family : str, optional
349 Indicates the family of elements used to create the function space
350 for the trail and test functions. The default is ``'CG'``, which are the class
351 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
352 order : int, optional
353 Defines the order of the elements in the function space.
354 refinements : int, optional
355 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
356 nu : float, optional
357 Diffusion coefficient :math:`\nu`.
358 rho : int, optional
359 Problem parameter.
360 delta : float, optional
361 Problem parameter.
363 Attributes
364 ----------
365 V : FunctionSpace
366 Defines the function space of the trial and test functions.
367 M : scalar, vector, matrix or higher rank tensor
368 Mass matrix for FENiCS.
369 K : scalar, vector, matrix or higher rank tensor
370 Stiffness matrix including diffusion coefficient (and correct sign).
372 References
373 ----------
374 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
375 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
376 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
377 Wells and others. Springer (2012).
378 """
380 def solve_system(self, rhs, factor, u0, t):
381 r"""
382 Dolfin's linear solver for :math:`(M - factor \cdot A)\vec{u} = \vec{rhs}`.
384 Parameters
385 ----------
386 rhs : dtype_f
387 Right-hand side for the nonlinear system.
388 factor : float
389 Abbrev. for the node-to-node stepsize (or any other factor required).
390 u0 : dtype_u
391 Initial guess for the iterative solver (not used here so far).
392 t : float
393 Current time.
395 Returns
396 -------
397 u : dtype_u
398 The solution as mesh.
399 """
401 A = self.M + self.nu * factor * self.K
402 b = self.dtype_u(rhs)
403 u = self.dtype_u(u0)
404 df.solve(A, u.values.vector(), b.values.vector())
406 return u
408 def __eval_fexpl(self, u, t):
409 """
410 Helper routine to evaluate the explicit part of the right-hand side.
412 Parameters
413 ----------
414 u : 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 fexpl : dtype_u
422 Explicit part of the right-hand side.
423 """
425 b = self.apply_mass_matrix(u)
426 psi = self.dtype_u(self.V)
427 df.solve(self.K, psi.values.vector(), b.values.vector())
429 fexpl = self.dtype_u(self.V)
430 fexpl.values = df.project(
431 df.Dx(psi.values, 1) * df.Dx(u.values, 0) - df.Dx(psi.values, 0) * df.Dx(u.values, 1), self.V
432 )
433 fexpl = self.apply_mass_matrix(fexpl)
435 return fexpl
437 def __eval_fimpl(self, u, t):
438 """
439 Helper routine to evaluate the implicit part of the right-hand side.
441 Parameters
442 ----------
443 u : dtype_u
444 Current values of the numerical solution.
445 t : float
446 Current time at which the numerical solution is computed.
448 Returns
449 -------
450 fimpl : dtype_u
451 Implicit part of the right-hand side.
452 """
454 A = -self.nu * self.K
455 fimpl = self.dtype_u(self.V)
456 A.mult(u.values.vector(), fimpl.values.vector())
458 return fimpl
460 def eval_f(self, u, t):
461 """
462 Routine to evaluate both parts of the right-hand side.
464 Note: Need to add this here, because otherwise the parent class will call the "local" functions __eval_*
465 and not the ones of the child class.
467 Parameters
468 ----------
469 u : dtype_u
470 Current values of the numerical solution.
471 t : float
472 Current time at which the numerical solution is computed.
474 Returns
475 -------
476 f : dtype_f
477 The right-hand side divided into two parts.
478 """
480 f = self.dtype_f(self.V)
481 f.impl = self.__eval_fimpl(u, t)
482 f.expl = self.__eval_fexpl(u, t)
483 return f