# Coverage for pySDC/implementations/problem_classes/GeneralizedFisher_1D_FD_implicit.py: 96%

## 57 statements

, created at 2024-09-09 14:59 +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 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(),

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