Coverage for pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py: 100%
112 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 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_heat(Problem):
12 r"""
13 Example implementing the forced one-dimensional heat equation with Dirichlet boundary conditions
15 .. math::
16 \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f
18 for :math:`x \in \Omega:=[0,1]`, where the forcing term :math:`f` is defined by
20 .. math::
21 f(x, t) = -\sin(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
23 For initial conditions with constant c and
25 .. math::
26 u(x, 0) = \sin(\pi x) + c
28 the exact solution of the problem is given by
30 .. math::
31 u(x, t) = \sin(\pi x)\cos(t) + c.
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 = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.
39 The part containing the forcing term is treated explicitly, where it is interpolated in the function space.
40 The other part will be treated in an implicit way.
42 Parameters
43 ----------
44 c_nvars : int, optional
45 Spatial resolution, i.e., numbers of degrees of freedom in space.
46 t0 : float, optional
47 Starting time.
48 family : str, optional
49 Indicates the family of elements used to create the function space
50 for the trail and test functions. The default is ``'CG'``, which are the class
51 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
52 order : int, optional
53 Defines the order of the elements in the function space.
54 refinements : int, optional
55 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
56 nu : float, optional
57 Diffusion coefficient :math:`\nu`.
58 c: float, optional
59 Constant for the Dirichlet boundary condition :math: `c`
61 Attributes
62 ----------
63 V : FunctionSpace
64 Defines the function space of the trial and test functions.
65 M : scalar, vector, matrix or higher rank tensor
66 Denotes the expression :math:`\int_\Omega u_t v\,dx`.
67 K : scalar, vector, matrix or higher rank tensor
68 Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`.
69 g : Expression
70 The forcing term :math:`f` in the heat equation.
71 bc : DirichletBC
72 Denotes the Dirichlet boundary conditions.
74 References
75 ----------
76 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
77 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
78 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
79 Wells and others. Springer (2012).
80 """
82 dtype_u = fenics_mesh
83 dtype_f = rhs_fenics_mesh
85 def __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):
86 """Initialization routine"""
88 # define the Dirichlet boundary
89 def Boundary(x, on_boundary):
90 return on_boundary
92 # set logger level for FFC and dolfin
93 logging.getLogger('FFC').setLevel(logging.WARNING)
94 logging.getLogger('UFL').setLevel(logging.WARNING)
96 # set solver and form parameters
97 df.parameters["form_compiler"]["optimize"] = True
98 df.parameters["form_compiler"]["cpp_optimize"] = True
99 df.parameters['allow_extrapolation'] = True
101 # set mesh and refinement (for multilevel)
102 mesh = df.UnitIntervalMesh(c_nvars)
103 for _ in range(refinements):
104 mesh = df.refine(mesh)
106 # define function space for future reference
107 self.V = df.FunctionSpace(mesh, family, order)
108 tmp = df.Function(self.V)
109 print('DoFs on this level:', len(tmp.vector()[:]))
111 # invoke super init, passing number of dofs, dtype_u and dtype_f
112 super(fenics_heat, self).__init__(self.V)
113 self._makeAttributeAndRegister(
114 'c_nvars', 't0', 'family', 'order', 'refinements', 'nu', 'c', localVars=locals(), readOnly=True
115 )
117 # Stiffness term (Laplace)
118 u = df.TrialFunction(self.V)
119 v = df.TestFunction(self.V)
120 a_K = -1.0 * df.inner(df.nabla_grad(u), self.nu * df.nabla_grad(v)) * df.dx
122 # Mass term
123 a_M = u * v * df.dx
125 self.M = df.assemble(a_M)
126 self.K = df.assemble(a_K)
128 # set boundary values
129 self.bc = df.DirichletBC(self.V, df.Constant(c), Boundary)
130 self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)
132 # set forcing term as expression
133 self.g = df.Expression(
134 '-sin(a*x[0]) * (sin(t) - b*a*a*cos(t))',
135 a=np.pi,
136 b=self.nu,
137 t=self.t0,
138 degree=self.order,
139 )
141 def solve_system(self, rhs, factor, u0, t):
142 r"""
143 Dolfin's linear solver for :math:`(M - factor \cdot A) \vec{u} = \vec{rhs}`.
145 Parameters
146 ----------
147 rhs : dtype_f
148 Right-hand side for the nonlinear system.
149 factor : float
150 Abbrev. for the node-to-node stepsize (or any other factor required).
151 u0 : dtype_u
152 Initial guess for the iterative solver (not used here so far).
153 t : float
154 Current time.
156 Returns
157 -------
158 u : dtype_u
159 Solution.
160 """
162 b = self.apply_mass_matrix(rhs)
164 u = self.dtype_u(u0)
165 T = self.M - factor * self.K
166 self.bc.apply(T, b.values.vector())
167 df.solve(T, 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 (not used here).
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 self.g.t = t
189 fexpl = self.dtype_u(df.interpolate(self.g, self.V))
191 return fexpl
193 def __eval_fimpl(self, u, t):
194 """
195 Helper routine to evaluate the implicit part of the right-hand side.
197 Parameters
198 ----------
199 u : dtype_u
200 Current values of the numerical solution.
201 t : float
202 Current time at which the numerical solution is computed.
204 Returns
205 -------
206 fimpl : dtype_u
207 Explicit part of the right-hand side.
208 """
210 tmp = self.dtype_u(self.V)
211 self.K.mult(u.values.vector(), tmp.values.vector())
212 fimpl = self.__invert_mass_matrix(tmp)
214 return fimpl
216 def eval_f(self, u, t):
217 """
218 Routine to evaluate both parts of the right-hand side of the problem.
220 Parameters
221 ----------
222 u : dtype_u
223 Current values of the numerical solution.
224 t : float
225 Current time at which the numerical solution is computed.
227 Returns
228 -------
229 f : dtype_f
230 The right-hand side divided into two parts.
231 """
233 f = self.dtype_f(self.V)
234 f.impl = self.__eval_fimpl(u, t)
235 f.expl = self.__eval_fexpl(u, t)
236 return f
238 def apply_mass_matrix(self, u):
239 r"""
240 Routine to apply mass matrix.
242 Parameters
243 ----------
244 u : dtype_u
245 Current values of the numerical solution.
247 Returns
248 -------
249 me : dtype_u
250 The product :math:`M \vec{u}`.
251 """
253 me = self.dtype_u(self.V)
254 self.M.mult(u.values.vector(), me.values.vector())
256 return me
258 def __invert_mass_matrix(self, u):
259 r"""
260 Helper routine to invert mass matrix.
262 Parameters
263 ----------
264 u : dtype_u
265 Current values of the numerical solution.
267 Returns
268 -------
269 me : dtype_u
270 The product :math:`M^{-1} \vec{u}`.
271 """
273 me = self.dtype_u(self.V)
275 b = self.dtype_u(u)
276 M = self.M
277 self.bc_hom.apply(M, b.values.vector())
279 df.solve(M, me.values.vector(), b.values.vector())
280 return me
282 def u_exact(self, t):
283 r"""
284 Routine to compute the exact solution at time :math:`t`.
286 Parameters
287 ----------
288 t : float
289 Time of the exact solution.
291 Returns
292 -------
293 me : dtype_u
294 Exact solution.
295 """
297 u0 = df.Expression('sin(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order)
298 me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)
300 return me
303# noinspection PyUnusedLocal
304class fenics_heat_mass(fenics_heat):
305 r"""
306 Example implementing the forced one-dimensional heat equation with Dirichlet boundary conditions
308 .. math::
309 \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f
311 for :math:`x \in \Omega:=[0,1]`, where the forcing term :math:`f` is defined by
313 .. math::
314 f(x, t) = -\sin(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
316 For initial conditions with constant c and
318 .. math::
319 u(x, 0) = \sin(\pi x) + c
321 the exact solution of the problem is given by
323 .. math::
324 u(x, t) = \sin(\pi x)\cos(t) + c.
326 In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem
327 is reformulated to the *weak formulation*
329 .. math:
330 \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.
332 The forcing term is treated explicitly, and is expressed via the mass matrix resulting from the left-hand side term
333 :math:`\int_\Omega u_t v\,dx`, and the other part will be treated in an implicit way.
335 Parameters
336 ----------
337 c_nvars : int, optional
338 Spatial resolution, i.e., numbers of degrees of freedom in space.
339 t0 : float, optional
340 Starting time.
341 family : str, optional
342 Indicates the family of elements used to create the function space
343 for the trail and test functions. The default is ``'CG'``, which are the class
344 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
345 order : int, optional
346 Defines the order of the elements in the function space.
347 refinements : int, optional
348 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
349 nu : float, optional
350 Diffusion coefficient :math:`\nu`.
352 Attributes
353 ----------
354 V : FunctionSpace
355 Defines the function space of the trial and test functions.
356 M : scalar, vector, matrix or higher rank tensor
357 Denotes the expression :math:`\int_\Omega u_t v\,dx`.
358 K : scalar, vector, matrix or higher rank tensor
359 Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`.
360 g : Expression
361 The forcing term :math:`f` in the heat equation.
362 bc : DirichletBC
363 Denotes the Dirichlet boundary conditions.
364 bc_hom : DirichletBC
365 Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual
366 fix_bc_for_residual: boolean
367 flag to indicate that the residual requires special treatment due to boundary conditions
369 References
370 ----------
371 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
372 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
373 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
374 Wells and others. Springer (2012).
375 """
377 def __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):
378 """Initialization routine"""
380 super().__init__(c_nvars, t0, family, order, refinements, nu, c)
382 self.fix_bc_for_residual = True
384 def solve_system(self, rhs, factor, u0, t):
385 r"""
386 Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
388 Parameters
389 ----------
390 rhs : dtype_f
391 Right-hand side for the nonlinear system.
392 factor : float
393 Abbrev. for the node-to-node stepsize (or any other factor required).
394 u0 : dtype_u
395 Initial guess for the iterative solver (not used here so far).
396 t : float
397 Current time.
399 Returns
400 -------
401 u : dtype_u
402 Solution.
403 """
405 u = self.dtype_u(u0)
406 T = self.M - factor * self.K
407 b = self.dtype_u(rhs)
409 self.bc.apply(T, b.values.vector())
411 df.solve(T, u.values.vector(), b.values.vector())
413 return u
415 def eval_f(self, u, t):
416 """
417 Routine to evaluate both parts of the right-hand side.
419 Parameters
420 ----------
421 u : dtype_u
422 Current values of the numerical solution.
423 t : float
424 Current time at which the numerical solution is computed.
426 Returns
427 -------
428 f : dtype_f
429 The right-hand side divided into two parts.
430 """
432 f = self.dtype_f(self.V)
434 self.K.mult(u.values.vector(), f.impl.values.vector())
436 self.g.t = t
437 f.expl = self.dtype_u(df.interpolate(self.g, self.V))
438 f.expl = self.apply_mass_matrix(f.expl)
440 return f
442 def fix_residual(self, res):
443 """
444 Applies homogeneous Dirichlet boundary conditions to the residual
446 Parameters
447 ----------
448 res : dtype_u
449 Residual
450 """
451 self.bc_hom.apply(res.values.vector())
452 return None
455# noinspection PyUnusedLocal
456class fenics_heat_mass_timebc(fenics_heat_mass):
457 r"""
458 Example implementing the forced one-dimensional heat equation with time-dependent Dirichlet boundary conditions
460 .. math::
461 \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f
463 for :math:`x \in \Omega:=[0,1]`, where the forcing term :math:`f` is defined by
465 .. math::
466 f(x, t) = -\cos(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
468 and the boundary conditions are given by
470 .. math::
471 u(x, t) = \cos(\pi x)\cos(t).
473 The exact solution of the problem is given by
475 .. math::
476 u(x, t) = \cos(\pi x)\cos(t) + c.
478 In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem
479 is reformulated to the *weak formulation*
481 .. math:
482 \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.
484 The forcing term is treated explicitly, and is expressed via the mass matrix resulting from the left-hand side term
485 :math:`\int_\Omega u_t v\,dx`, and the other part will be treated in an implicit way.
487 Parameters
488 ----------
489 c_nvars : int, optional
490 Spatial resolution, i.e., numbers of degrees of freedom in space.
491 t0 : float, optional
492 Starting time.
493 family : str, optional
494 Indicates the family of elements used to create the function space
495 for the trail and test functions. The default is ``'CG'``, which are the class
496 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
497 order : int, optional
498 Defines the order of the elements in the function space.
499 refinements : int, optional
500 Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
501 nu : float, optional
502 Diffusion coefficient :math:`\nu`.
504 Attributes
505 ----------
506 V : FunctionSpace
507 Defines the function space of the trial and test functions.
508 M : scalar, vector, matrix or higher rank tensor
509 Denotes the expression :math:`\int_\Omega u_t v\,dx`.
510 K : scalar, vector, matrix or higher rank tensor
511 Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`.
512 g : Expression
513 The forcing term :math:`f` in the heat equation.
514 bc : DirichletBC
515 Denotes the time-dependent Dirichlet boundary conditions.
516 bc_hom : DirichletBC
517 Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual
518 fix_bc_for_residual: boolean
519 flag to indicate that the residual requires special treatment due to boundary conditions
521 References
522 ----------
523 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
524 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
525 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
526 Wells and others. Springer (2012).
527 """
529 def __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):
530 """Initialization routine"""
532 # define the Dirichlet boundary
533 def Boundary(x, on_boundary):
534 return on_boundary
536 super().__init__(c_nvars, t0, family, order, refinements, nu, c)
538 self.u_D = df.Expression('cos(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t0, degree=self.order)
539 self.bc = df.DirichletBC(self.V, self.u_D, Boundary)
540 self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)
542 # set forcing term as expression
543 self.g = df.Expression(
544 '-cos(a*x[0]) * (sin(t) - b*a*a*cos(t))',
545 a=np.pi,
546 b=self.nu,
547 t=self.t0,
548 degree=self.order,
549 )
551 def solve_system(self, rhs, factor, u0, t):
552 r"""
553 Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
555 Parameters
556 ----------
557 rhs : dtype_f
558 Right-hand side for the nonlinear system.
559 factor : float
560 Abbrev. for the node-to-node stepsize (or any other factor required).
561 u0 : dtype_u
562 Initial guess for the iterative solver (not used here so far).
563 t : float
564 Current time.
566 Returns
567 -------
568 u : dtype_u
569 Solution.
570 """
572 u = self.dtype_u(u0)
573 T = self.M - factor * self.K
574 b = self.dtype_u(rhs)
576 self.u_D.t = t
578 self.bc.apply(T, b.values.vector())
579 self.bc.apply(b.values.vector())
581 df.solve(T, u.values.vector(), b.values.vector())
583 return u
585 def u_exact(self, t):
586 r"""
587 Routine to compute the exact solution at time :math:`t`.
589 Parameters
590 ----------
591 t : float
592 Time of the exact solution.
594 Returns
595 -------
596 me : dtype_u
597 Exact solution.
598 """
600 u0 = df.Expression('cos(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order)
601 me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)
603 return me