Coverage for pySDC/implementations/problem_classes/GeneralizedFisher_1D_FD_implicit.py: 96%
57 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
2import scipy.sparse as sp
3from scipy.sparse.linalg import spsolve
5from pySDC.core.Errors import ParameterError, ProblemError
6from pySDC.core.Problem import ptype
7from pySDC.helpers import problem_helper
8from pySDC.implementations.datatype_classes.mesh import mesh
11# noinspection PyUnusedLocal
12class generalized_fisher(ptype):
13 r"""
14 The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
15 be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
16 problem [1]_
18 .. math::
19 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
21 with initial condition
23 .. math::
24 u(x, 0) = \left[
25 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\delta x\right)
26 \right]^{2 / \nu}
28 for :math:`x \in \mathbb{R}`. For
30 .. math::
31 \delta = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
33 .. math::
34 \lambda_1 = \frac{\lambda_0}{2} \left[
35 \left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
36 \right],
38 the exact solution is
40 .. math::
41 u(x, t) = \left(
42 1 + \left(2^{\nu / 2} - 1\right) \exp\left(-\frac{\nu}{2}\delta (x + 2 \lambda_1 t)\right)
43 \right)^{-2 / n}.
45 Spatial discretization is done by centered finite differences.
47 Parameters
48 ----------
49 nvars : int, optional
50 Spatial resolution, i.e., number of degrees of freedom in space.
51 nu : float, optional
52 Problem parameter :math:`\nu`.
53 lambda0 : float, optional
54 Problem parameter :math:`\lambda_0`.
55 newton_maxiter : int, optional
56 Maximum number of Newton iterations to solve the nonlinear system.
57 newton_tol : float, optional
58 Tolerance for Newton to terminate.
59 interval : tuple, optional
60 Defines the spatial domain.
61 stop_at_nan : bool, optional
62 Indicates if the nonlinear solver should stop if ``nan`` values arise.
64 Attributes
65 ----------
66 A : sparse matrix (CSC)
67 Second-order FD discretization of the 1D Laplace operator.
68 dx : float
69 Distance between two spatial nodes.
71 References
72 ----------
73 .. [1] Z. Feng. Traveling wave behavior for a generalized fisher equation. Chaos Solitons Fractals 38(2),
74 481–488 (2008).
75 """
77 dtype_u = mesh
78 dtype_f = mesh
80 def __init__(
81 self, nvars=127, nu=1.0, lambda0=2.0, newton_maxiter=100, newton_tol=1e-12, interval=(-5, 5), stop_at_nan=True
82 ):
83 """Initialization routine"""
85 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
86 if (nvars + 1) % 2 != 0:
87 raise ProblemError('setup requires nvars = 2^p - 1')
89 # invoke super init, passing number of dofs, dtype_u and dtype_f
90 super().__init__((nvars, None, np.dtype('float64')))
91 self._makeAttributeAndRegister(
92 'nvars',
93 'nu',
94 'lambda0',
95 'newton_maxiter',
96 'newton_tol',
97 'interval',
98 'stop_at_nan',
99 localVars=locals(),
100 readOnly=True,
101 )
103 # compute dx and get discretization matrix A
104 self.dx = (self.interval[1] - self.interval[0]) / (self.nvars + 1)
105 self.A, _ = problem_helper.get_finite_difference_matrix(
106 derivative=2,
107 order=2,
108 stencil_type='center',
109 dx=self.dx,
110 size=self.nvars + 2,
111 dim=1,
112 bc='dirichlet-zero',
113 )
115 # noinspection PyTypeChecker
116 def solve_system(self, rhs, factor, u0, t):
117 """
118 Simple Newton solver.
120 Parameters
121 ----------
122 rhs : dtype_f
123 Right-hand side for the nonlinear system.
124 factor : float
125 Abbrev. for the node-to-node stepsize (or any other factor required).
126 u0 : dtype_u
127 Initial guess for the iterative solver.
128 t : float
129 urrent time (required here for the BC).
131 Returns
132 -------
133 u : dtype_u
134 Solution.
135 """
137 u = self.dtype_u(u0)
139 nu = self.nu
140 lambda0 = self.lambda0
142 # set up boundary values to embed inner points
143 lam1 = lambda0 / 2.0 * ((nu / 2.0 + 1) ** 0.5 + (nu / 2.0 + 1) ** (-0.5))
144 sig1 = lam1 - np.sqrt(lam1**2 - lambda0**2)
145 ul = (1 + (2 ** (nu / 2.0) - 1) * np.exp(-nu / 2.0 * sig1 * (self.interval[0] + 2 * lam1 * t))) ** (-2.0 / nu)
146 ur = (1 + (2 ** (nu / 2.0) - 1) * np.exp(-nu / 2.0 * sig1 * (self.interval[1] + 2 * lam1 * t))) ** (-2.0 / nu)
148 # start newton iteration
149 n = 0
150 res = 99
151 while n < self.newton_maxiter:
152 # form the function g with g(u) = 0
153 uext = np.concatenate(([ul], u, [ur]))
154 g = u - factor * (self.A.dot(uext)[1:-1] + lambda0**2 * u * (1 - u**nu)) - rhs
156 # if g is close to 0, then we are done
157 res = np.linalg.norm(g, np.inf)
159 if res < self.newton_tol:
160 break
162 # assemble dg
163 dg = sp.eye(self.nvars) - factor * (
164 self.A[1:-1, 1:-1] + sp.diags(lambda0**2 - lambda0**2 * (nu + 1) * u**nu, offsets=0)
165 )
167 # newton update: u1 = u0 - g/dg
168 u -= spsolve(dg, g)
170 # increase iteration count
171 n += 1
173 if np.isnan(res) and self.stop_at_nan:
174 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
175 elif np.isnan(res):
176 self.logger.warning('Newton got nan after %i iterations...' % n)
178 if n == self.newton_maxiter:
179 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
181 return u
183 def eval_f(self, u, t):
184 """
185 Routine to evaluate the right-hand side of the problem.
187 Parameters
188 ----------
189 u : dtype_u
190 Current values of the numerical solution.
191 t : float
192 Current time of the numerical solution is computed.
194 Returns
195 -------
196 f : dtype_f
197 The right-hand side of the problem.
198 """
199 # set up boundary values to embed inner points
200 lam1 = self.lambda0 / 2.0 * ((self.nu / 2.0 + 1) ** 0.5 + (self.nu / 2.0 + 1) ** (-0.5))
201 sig1 = lam1 - np.sqrt(lam1**2 - self.lambda0**2)
202 ul = (1 + (2 ** (self.nu / 2.0) - 1) * np.exp(-self.nu / 2.0 * sig1 * (self.interval[0] + 2 * lam1 * t))) ** (
203 -2 / self.nu
204 )
205 ur = (1 + (2 ** (self.nu / 2.0) - 1) * np.exp(-self.nu / 2.0 * sig1 * (self.interval[1] + 2 * lam1 * t))) ** (
206 -2 / self.nu
207 )
209 uext = np.concatenate(([ul], u, [ur]))
211 f = self.dtype_f(self.init)
212 f[:] = self.A.dot(uext)[1:-1] + self.lambda0**2 * u * (1 - u**self.nu)
213 return f
215 def u_exact(self, t):
216 r"""
217 Routine to compute the exact solution at time :math:`t`.
219 Parameters
220 ----------
221 t : float
222 Time of the exact solution.
224 Returns
225 -------
226 me : dtype_u
227 Exact solution.
228 """
230 me = self.dtype_u(self.init)
231 xvalues = np.array([(i + 1 - (self.nvars + 1) / 2) * self.dx for i in range(self.nvars)])
233 lam1 = self.lambda0 / 2.0 * ((self.nu / 2.0 + 1) ** 0.5 + (self.nu / 2.0 + 1) ** (-0.5))
234 sig1 = lam1 - np.sqrt(lam1**2 - self.lambda0**2)
235 me[:] = (1 + (2 ** (self.nu / 2.0) - 1) * np.exp(-self.nu / 2.0 * sig1 * (xvalues + 2 * lam1 * t))) ** (
236 -2.0 / self.nu
237 )
238 return me