Coverage for pySDC/implementations/problem_classes/GrayScott_1D_FEniCS_implicit.py: 100%
116 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import logging
2import random
4import dolfin as df
5import numpy as np
7from pySDC.core.problem import Problem
8from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh
11# noinspection PyUnusedLocal
12class fenics_grayscott(Problem):
13 r"""
14 1D Gray-Scott reaction-diffusion system with FEniCS finite elements, fully implicit with FEniCS' Newton solver.
16 The Gray-Scott system [1]_ describes a reaction-diffusion process of two substances :math:`u` and :math:`v`,
17 where they diffuse over time. During the reaction :math:`u` is used up with overall decay rate :math:`B`,
18 whereas :math:`v` is produced with feed rate :math:`A`. :math:`D_u,\, D_v` are the diffusion rates for
19 :math:`u,\, v`. This process is described by the one-dimensional model with homogeneous Neumann boundary conditions
21 .. math::
22 \frac{\partial u}{\partial t} = D_u \Delta u - u v^2 + A (1 - u),
24 .. math::
25 \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B u
27 for :math:`x \in \Omega:=[0, 100]`. The *weak formulation* of the problem can be obtained by multiplying the
28 system with a test function :math:`q`:
30 .. math::
31 \int_\Omega u_t q\,dx = \int_\Omega D_u \Delta u q - u v^2 q + A (1 - u) q\,dx,
33 .. math::
34 \int_\Omega v_t q\,dx = \int_\Omega D_v \Delta v q + u v^2 q - B u q\,dx,
36 The Laplacians are integrated by parts, and no Dirichlet conditions are imposed, so the boundary terms drop out and
37 the boundaries are zero-flux. The spatial solve of the weak formulation is realized by ``FEniCS`` [2]_.
39 Parameters
40 ----------
41 c_nvars : int, optional
42 Spatial resolution, i.e., number of degrees of freedom in space.
43 t0 : float, optional
44 Starting time.
45 family : str, optional
46 Indicates the family of elements used to create the function space
47 for the trail and test functions. The default is ``'CG'``, which are the class
48 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [3]_.
49 order : int, optional
50 Defines the order of the elements in the function space.
51 refinements : int, optional
52 How many times the mesh of ``c_nvars`` cells is refined uniformly, halving the cell size each time:
53 ``refinements=2`` gives ``4 * c_nvars`` cells. For a multilevel hierarchy, pass one value per level in the
54 description, e.g. ``[1, 0]``, which the controller splits across the levels. The default is ``0``.
55 Du : float, optional
56 Diffusion rate for :math:`u`.
57 Dv: float, optional
58 Diffusion rate for :math:`v`.
59 A : float, optional
60 Feed rate for :math:`v`.
61 B : float, optional
62 Overall decay rate for :math:`u`.
63 newton_tol : float, optional
64 Absolute tolerance of the node-local Newton solve. It has to be tighter than the SDC residual
65 tolerance, otherwise the node-local solve, not the SDC iteration, sets the accuracy floor.
66 newton_rtol : float, optional
67 Relative tolerance of the node-local Newton solve. Dolfin stops on whichever bar is met first,
68 so this one has to be lowered alongside ``newton_tol`` to actually tighten the solve.
69 newton_maxiter : int, optional
70 Maximum number of node-local Newton iterations.
72 Attributes
73 ----------
74 V : FunctionSpace
75 Defines the function space of the trial and test functions.
76 w : Function
77 Function for the right-hand side.
78 w1 : Function
79 Split of w, part 1.
80 w2 : Function
81 Split of w, part 2.
82 F1 : scalar, vector, matrix or higher rank tensor
83 Weak form of right-hand side, first part.
84 F2 : scalar, vector, matrix or higher rank tensor
85 Weak form of right-hand side, second part.
86 F : scalar, vector, matrix or higher rank tensor
87 Weak form of full right-hand side.
88 M : matrix
89 Full mass matrix for both parts.
91 References
92 ----------
93 .. [1] Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms
94 of multistability. P. Gray, S. K. Scott. Chem. Eng. Sci. 38, 1 (1983).
95 .. [2] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
96 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
97 .. [3] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
98 Wells and others. Springer (2012).
99 """
101 dtype_u = fenics_mesh
102 dtype_f = fenics_mesh
104 def __init__(
105 self,
106 c_nvars=256,
107 t0=0.0,
108 family='CG',
109 order=4,
110 refinements=0,
111 Du=1.0,
112 Dv=0.01,
113 A=0.09,
114 B=0.086,
115 newton_tol=1e-9,
116 newton_rtol=1e-8,
117 newton_maxiter=100,
118 ):
119 """Initialization routine"""
121 # set logger level for FFC and dolfin
122 warning_level = getattr(df, 'WARNING', None)
123 if warning_level is None and hasattr(df, 'LogLevel'):
124 warning_level = df.LogLevel.WARNING
125 if warning_level is not None:
126 df.set_log_level(warning_level)
127 logging.getLogger('FFC').setLevel(logging.WARNING)
129 # set solver and form parameters
130 df.parameters["form_compiler"]["optimize"] = True
131 df.parameters["form_compiler"]["cpp_optimize"] = True
133 # set mesh and refinement (for multilevel)
134 mesh = df.IntervalMesh(c_nvars, 0, 100)
135 for _ in range(refinements):
136 mesh = df.refine(mesh)
138 # define mixed function space for future reference. `V * V` was removed in DOLFIN 2019.1,
139 # so the mixed space is built from a MixedElement instead.
140 element = df.FiniteElement(family, mesh.ufl_cell(), order)
141 self.V = df.FunctionSpace(mesh, df.MixedElement([element, element]))
143 super().__init__(self.V)
144 self._makeAttributeAndRegister(
145 'c_nvars', 't0', 'family', 'order', 'refinements', 'Du', 'Dv', 'A', 'B', localVars=locals(), readOnly=True
146 )
147 self._makeAttributeAndRegister(
148 'newton_tol', 'newton_rtol', 'newton_maxiter', localVars=locals(), readOnly=False
149 )
150 # rhs in weak form
151 self.w = df.Function(self.V)
152 q1, q2 = df.TestFunctions(self.V)
154 self.w1, self.w2 = df.split(self.w)
156 self.F1 = (
157 -self.Du * df.inner(df.nabla_grad(self.w1), df.nabla_grad(q1))
158 - self.w1 * (self.w2**2) * q1
159 + self.A * (1 - self.w1) * q1
160 ) * df.dx
161 self.F2 = (
162 -self.Dv * df.inner(df.nabla_grad(self.w2), df.nabla_grad(q2))
163 + self.w1 * (self.w2**2) * q2
164 - self.B * self.w2 * q2
165 ) * df.dx
166 self.F = self.F1 + self.F2
168 # mass matrix
169 u1, u2 = df.TrialFunctions(self.V)
170 a_M = u1 * q1 * df.dx
171 M1 = df.assemble(a_M)
172 a_M = u2 * q2 * df.dx
173 M2 = df.assemble(a_M)
174 self.M = M1 + M2
176 def apply_mass_matrix(self, u):
177 r"""
178 Apply the mass matrix, :math:`M \vec{u}`.
180 Required by the mass-matrix sweepers and transfers; ``Problem`` has no default.
182 Parameters
183 ----------
184 u : dtype_u
185 Current values of the numerical solution.
187 Returns
188 -------
189 me : dtype_u
190 The product :math:`M \vec{u}`.
191 """
193 me = self.dtype_u(self.V)
194 self.M.mult(u.values.vector(), me.values.vector())
196 return me
198 def __invert_mass_matrix(self, u):
199 r"""
200 Helper routine to invert mass matrix.
202 Parameters
203 ----------
204 u : dtype_u
205 Current values of the numerical solution.
207 Returns
208 -------
209 me : dtype_u
210 The product :math:`M^{-1} \vec{u}`.
211 """
213 me = self.dtype_u(self.V)
215 A = 1.0 * self.M
216 b = self.dtype_u(u)
218 df.solve(A, me.values.vector(), b.values.vector())
220 return me
222 def solve_system(self, rhs, factor, u0, t):
223 r"""
224 Dolfin's linear solver for :math:`(M - factor \cdot A) \vec{u} = \vec{rhs}`.
226 Parameters
227 ----------
228 rhs : dtype_f
229 Right-hand side for the nonlinear system.
230 factor : float
231 Abbrev. for the node-to-node stepsize (or any other factor required).
232 u0 : dtype_u
233 Initial guess for the iterative solver (not used here so far).
234 t : float
235 Current time.
237 Returns
238 -------
239 me : dtype_u
240 Solution as mesh.
241 """
243 sol = self.dtype_u(self.V)
245 self.w.assign(sol.values)
247 # fixme: is this really necessary to do each time?
248 q1, q2 = df.TestFunctions(self.V)
249 w1, w2 = df.split(self.w)
250 r1, r2 = df.split(rhs.values)
251 F1 = w1 * q1 * df.dx - factor * self.F1 - r1 * q1 * df.dx
252 F2 = w2 * q2 * df.dx - factor * self.F2 - r2 * q2 * df.dx
253 F = F1 + F2
254 du = df.TrialFunction(self.V)
255 J = df.derivative(F, self.w, du)
257 problem = df.NonlinearVariationalProblem(F, self.w, [], J)
258 solver = df.NonlinearVariationalSolver(problem)
260 prm = solver.parameters
261 prm['newton_solver']['absolute_tolerance'] = self.newton_tol
262 prm['newton_solver']['relative_tolerance'] = self.newton_rtol
263 prm['newton_solver']['maximum_iterations'] = self.newton_maxiter
264 prm['newton_solver']['relaxation_parameter'] = 1.0
266 solver.solve()
268 sol.values.assign(self.w)
270 return sol
272 def eval_f(self, u, t):
273 """
274 Routine to evaluate both parts of the right-hand side.
276 Parameters
277 ----------
278 u : dtype_u
279 Current values of the numerical solution.
280 t : float
281 Current time at which the numerical solution is computed.
283 Returns
284 -------
285 f : dtype_f
286 The right-hand side divided into two parts.
287 """
289 f = self.dtype_f(self.V)
291 self.w.assign(u.values)
292 f.values = df.Function(self.V, df.assemble(self.F))
294 f = self.__invert_mass_matrix(f)
296 return f
298 def u_exact(self, t):
299 r"""
300 Routine to compute the exact solution at time :math:`t`.
302 Parameters
303 ----------
304 t : float
305 Time of the exact solution.
307 Returns
308 -------
309 me : dtype_u
310 Exact solution (only at :math:`t_0 = 0.0`).
311 """
313 # subclassing df.Expression was removed in DOLFIN 2018.1; UserExpression is the
314 # replacement and requires the base initialiser to run.
315 class InitialConditions(df.UserExpression):
316 def __init__(self, **kwargs):
317 random.seed(2)
318 super().__init__(**kwargs)
320 def eval(self, values, x):
321 values[0] = 1 - 0.5 * np.power(np.sin(np.pi * x[0] / 100), 100)
322 values[1] = 0.25 * np.power(np.sin(np.pi * x[0] / 100), 100)
324 def value_shape(self):
325 return (2,)
327 assert t == 0, 'ERROR: u_exact only valid for t=0'
329 uinit = InitialConditions(degree=self.order)
331 me = self.dtype_u(self.V)
332 me.values = df.interpolate(uinit, self.V)
334 return me
337class fenics_grayscott_mass(fenics_grayscott):
338 r"""
339 Gray-Scott in mass-matrix form, :math:`M \vec{u}' = F(\vec{u})`, with no mass inversion.
341 ``eval_f`` returns the assembled weak-form residual (a load vector) rather than
342 :math:`M^{-1} F`, and ``solve_system`` takes an rhs that is already in the dual space. Use
343 with :class:`generic_implicit_mass` and, for MLSDC, :class:`base_transfer_mass`.
345 The Newton loop is written out rather than handed to ``NonlinearVariationalSolver`` because
346 the rhs is a vector, not a form, so it cannot be folded into the variational problem.
347 """
349 def eval_f(self, u, t):
350 r"""
351 Evaluate the right-hand side in weak (dual) form, without inverting the mass matrix.
353 Parameters
354 ----------
355 u : dtype_u
356 Current values of the numerical solution.
357 t : float
358 Current time.
360 Returns
361 -------
362 f : dtype_f
363 The assembled load vector :math:`F(\vec{u})`.
364 """
366 f = self.dtype_f(self.V)
367 self.w.assign(u.values)
368 f.values = df.Function(self.V, df.assemble(self.F))
370 return f
372 def solve_system(self, rhs, factor, u0, t):
373 r"""
374 Solve :math:`M \vec{u} - factor \cdot F(\vec{u}) = \vec{rhs}` by Newton, rhs already dual.
376 Parameters
377 ----------
378 rhs : dtype_f
379 Right-hand side, in the dual space.
380 factor : float
381 Node-to-node stepsize.
382 u0 : dtype_u
383 Initial guess (unused; the loop starts from zero, matching the mass-inverse variant so
384 that the two are directly comparable).
385 t : float
386 Current time.
388 Returns
389 -------
390 sol : dtype_u
391 Solution.
392 """
394 sol = self.dtype_u(self.V)
395 self.w.assign(sol.values)
397 du = df.TrialFunction(self.V)
398 Jform = df.derivative(self.F, self.w, du)
400 res = df.Function(self.V)
401 delta = df.Function(self.V)
402 res0 = None
404 for _ in range(self.newton_maxiter):
405 self.M.mult(self.w.vector(), res.vector())
406 res.vector().axpy(-factor, df.assemble(self.F))
407 res.vector().axpy(-1.0, rhs.values.vector())
409 norm = res.vector().norm('l2')
410 res0 = norm if res0 is None else res0
411 if norm < self.newton_tol or norm < self.newton_rtol * res0:
412 break
414 df.solve(self.M - factor * df.assemble(Jform), delta.vector(), res.vector())
415 self.w.vector().axpy(-1.0, delta.vector())
417 sol.values.assign(self.w)
419 return sol