# Coverage for pySDC/implementations/problem_classes/GrayScott_1D_FEniCS_implicit.py: 0%

## 85 statements

, created at 2024-09-09 14:59 +0000

1import logging

2import random

4import dolfin as df

5import numpy as np

7from pySDC.core.errors import ParameterError

8from pySDC.core.problem import Problem

9from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh

12# noinspection PyUnusedLocal

13class fenics_grayscott(Problem):

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 = (

133 - self.w1 * (self.w2**2) * q1

134 + self.A * (1 - self.w1) * q1

135 ) * df.dx

136 self.F2 = (

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