Coverage for pySDC/implementations/problem_classes/GrayScott_1D_FEniCS_implicit.py: 0%
85 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import logging
2import random
4import dolfin as df
5import numpy as np
7from pySDC.core.Errors import ParameterError
8from pySDC.core.Problem import ptype
9from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh
12# noinspection PyUnusedLocal
13class fenics_grayscott(ptype):
14 r"""
15 The Gray-Scott system [1]_ describes a reaction-diffusion process of two substances :math:`u` and :math:`v`,
16 where they diffuse over time. During the reaction :math:`u` is used up with overall decay rate :math:`B`,
17 whereas :math:`v` is produced with feed rate :math:`A`. :math:`D_u,\, D_v` are the diffusion rates for
18 :math:`u,\, v`. This process is described by the one-dimensional model using Dirichlet boundary conditions
20 .. math::
21 \frac{\partial u}{\partial t} = D_u \Delta u - u v^2 + A (1 - u),
23 .. math::
24 \frac{\partial v}{\partial t} = D_v \Delta v + u v^2 - B u
26 for :math:`x \in \Omega:=[0, 100]`. The *weak formulation* of the problem can be obtained by multiplying the
27 system with a test function :math:`q`:
29 .. math::
30 \int_\Omega u_t q\,dx = \int_\Omega D_u \Delta u q - u v^2 q + A (1 - u) q\,dx,
32 .. math::
33 \int_\Omega v_t q\,dx = \int_\Omega D_v \Delta v q + u v^2 q - B u q\,dx,
35 The spatial solve of the weak formulation is realized by ``FEniCS`` [2]_.
37 Parameters
38 ----------
39 c_nvars : int, optional
40 Spatial resolution, i.e., number of degrees of freedom in space.
41 t0 : float, optional
42 Starting time.
43 family : str, optional
44 Indicates the family of elements used to create the function space
45 for the trail and test functions. The default is ``'CG'``, which are the class
46 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [3]_.
47 order : int, optional
48 Defines the order of the elements in the function space.
49 refinements : list or tuple, optional
50 Defines the refinement for the spatial grid. Needs to be a list or tuple, e.g.
51 ``refinements=[2, 2]`` or ``refinements=(2, 2)``.
52 Du : float, optional
53 Diffusion rate for :math:`u`.
54 Dv: float, optional
55 Diffusion rate for :math:`v`.
56 A : float, optional
57 Feed rate for :math:`v`.
58 B : float, optional
59 Overall decay rate for :math:`u`.
61 Attributes
62 ----------
63 V : FunctionSpace
64 Defines the function space of the trial and test functions.
65 w : Function
66 Function for the right-hand side.
67 w1 : Function
68 Split of w, part 1.
69 w2 : Function
70 Split of w, part 2.
71 F1 : scalar, vector, matrix or higher rank tensor
72 Weak form of right-hand side, first part.
73 F2 : scalar, vector, matrix or higher rank tensor
74 Weak form of right-hand side, second part.
75 F : scalar, vector, matrix or higher rank tensor
76 Weak form of full right-hand side.
77 M : matrix
78 Full mass matrix for both parts.
80 References
81 ----------
82 .. [1] Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms
83 of multistability. P. Gray, S. K. Scott. Chem. Eng. Sci. 38, 1 (1983).
84 .. [2] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
85 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
86 .. [3] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
87 Wells and others. Springer (2012).
88 """
90 dtype_u = fenics_mesh
91 dtype_f = fenics_mesh
93 def __init__(self, c_nvars=256, t0=0.0, family='CG', order=4, refinements=None, Du=1.0, Dv=0.01, A=0.09, B=0.086):
94 """Initialization routine"""
96 if refinements is None:
97 refinements = [1, 0]
99 # define the Dirichlet boundary
100 def Boundary(x, on_boundary):
101 return on_boundary
103 # set logger level for FFC and dolfin
104 df.set_log_level(df.WARNING)
105 logging.getLogger('FFC').setLevel(logging.WARNING)
107 # set solver and form parameters
108 df.parameters["form_compiler"]["optimize"] = True
109 df.parameters["form_compiler"]["cpp_optimize"] = True
111 # set mesh and refinement (for multilevel)
112 mesh = df.IntervalMesh(c_nvars, 0, 100)
113 for _ in range(refinements):
114 mesh = df.refine(mesh)
116 # define function space for future reference
117 V = df.FunctionSpace(mesh, family, order)
118 self.V = V * V
120 # invoke super init, passing number of dofs
121 super(fenics_grayscott).__init__(V)
122 self._makeAttributeAndRegister(
123 'c_nvars', 't0', 'family', 'order', 'refinements', 'Du', 'Dv', 'A', 'B', localVars=locals(), readOnly=True
124 )
125 # rhs in weak form
126 self.w = df.Function(self.V)
127 q1, q2 = df.TestFunctions(self.V)
129 self.w1, self.w2 = df.split(self.w)
131 self.F1 = (
132 -self.Du * df.inner(df.nabla_grad(self.w1), df.nabla_grad(q1))
133 - self.w1 * (self.w2**2) * q1
134 + self.A * (1 - self.w1) * q1
135 ) * df.dx
136 self.F2 = (
137 -self.Dv * df.inner(df.nabla_grad(self.w2), df.nabla_grad(q2))
138 + self.w1 * (self.w2**2) * q2
139 - self.B * self.w2 * q2
140 ) * df.dx
141 self.F = self.F1 + self.F2
143 # mass matrix
144 u1, u2 = df.TrialFunctions(self.V)
145 a_M = u1 * q1 * df.dx
146 M1 = df.assemble(a_M)
147 a_M = u2 * q2 * df.dx
148 M2 = df.assemble(a_M)
149 self.M = M1 + M2
151 def __invert_mass_matrix(self, u):
152 r"""
153 Helper routine to invert mass matrix.
155 Parameters
156 ----------
157 u : dtype_u
158 Current values of the numerical solution.
160 Returns
161 -------
162 me : dtype_u
163 The product :math:`M^{-1} \vec{u}`.
164 """
166 me = self.dtype_u(self.V)
168 A = 1.0 * self.M
169 b = self.dtype_u(u)
171 df.solve(A, me.values.vector(), b.values.vector())
173 return me
175 def solve_system(self, rhs, factor, u0, t):
176 r"""
177 Dolfin's linear solver for :math:`(M - factor \cdot A) \vec{u} = \vec{rhs}`.
179 Parameters
180 ----------
181 rhs : dtype_f
182 Right-hand side for the nonlinear system.
183 factor : float
184 Abbrev. for the node-to-node stepsize (or any other factor required).
185 u0 : dtype_u
186 Initial guess for the iterative solver (not used here so far).
187 t : float
188 Current time.
190 Returns
191 -------
192 me : dtype_u
193 Solution as mesh.
194 """
196 sol = self.dtype_u(self.V)
198 self.w.assign(sol.values)
200 # fixme: is this really necessary to do each time?
201 q1, q2 = df.TestFunctions(self.V)
202 w1, w2 = df.split(self.w)
203 r1, r2 = df.split(rhs.values)
204 F1 = w1 * q1 * df.dx - factor * self.F1 - r1 * q1 * df.dx
205 F2 = w2 * q2 * df.dx - factor * self.F2 - r2 * q2 * df.dx
206 F = F1 + F2
207 du = df.TrialFunction(self.V)
208 J = df.derivative(F, self.w, du)
210 problem = df.NonlinearVariationalProblem(F, self.w, [], J)
211 solver = df.NonlinearVariationalSolver(problem)
213 prm = solver.parameters
214 prm['newton_solver']['absolute_tolerance'] = 1e-09
215 prm['newton_solver']['relative_tolerance'] = 1e-08
216 prm['newton_solver']['maximum_iterations'] = 100
217 prm['newton_solver']['relaxation_parameter'] = 1.0
219 solver.solve()
221 sol.values.assign(self.w)
223 return sol
225 def eval_f(self, u, t):
226 """
227 Routine to evaluate both parts of the right-hand side.
229 Parameters
230 ----------
231 u : dtype_u
232 Current values of the numerical solution.
233 t : float
234 Current time at which the numerical solution is computed.
236 Returns
237 -------
238 f : dtype_f
239 The right-hand side divided into two parts.
240 """
242 f = self.dtype_f(self.V)
244 self.w.assign(u.values)
245 f.values = df.Function(self.V, df.assemble(self.F))
247 f = self.__invert_mass_matrix(f)
249 return f
251 def u_exact(self, t):
252 r"""
253 Routine to compute the exact solution at time :math:`t`.
255 Parameters
256 ----------
257 t : float
258 Time of the exact solution.
260 Returns
261 -------
262 me : dtype_u
263 Exact solution (only at :math:`t_0 = 0.0`).
264 """
266 class InitialConditions(df.Expression):
267 def __init__(self):
268 # fixme: why do we need this?
269 random.seed(2)
270 pass
272 def eval(self, values, x):
273 values[0] = 1 - 0.5 * np.power(np.sin(np.pi * x[0] / 100), 100)
274 values[1] = 0.25 * np.power(np.sin(np.pi * x[0] / 100), 100)
276 def value_shape(self):
277 return (2,)
279 assert t == 0, 'ERROR: u_exact only valid for t=0'
281 uinit = InitialConditions()
283 me = self.dtype_u(self.V)
284 me.values = df.interpolate(uinit, self.V)
286 return me