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