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

57 statements  

« prev     ^ index     » next       coverage.py v7.6.9, created at 2024-12-20 14:51 +0000

1import numpy as np 

2import scipy.sparse as sp 

3from scipy.sparse.linalg import spsolve 

4 

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 

9 

10 

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]_ 

17 

18 .. math:: 

19 \frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu) 

20 

21 with initial condition 

22 

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} 

27 

28 for :math:`x \in \mathbb{R}`. For 

29 

30 .. math:: 

31 \delta = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2}, 

32 

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], 

37 

38 the exact solution is 

39 

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}. 

44 

45 Spatial discretization is done by centered finite differences. 

46 

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. 

63 

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. 

70 

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 """ 

76 

77 dtype_u = mesh 

78 dtype_f = mesh 

79 

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""" 

84 

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') 

88 

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 ) 

102 

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 ) 

114 

115 # noinspection PyTypeChecker 

116 def solve_system(self, rhs, factor, u0, t): 

117 """ 

118 Simple Newton solver. 

119 

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). 

130 

131 Returns 

132 ------- 

133 u : dtype_u 

134 Solution. 

135 """ 

136 

137 u = self.dtype_u(u0) 

138 

139 nu = self.nu 

140 lambda0 = self.lambda0 

141 

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) 

147 

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 

155 

156 # if g is close to 0, then we are done 

157 res = np.linalg.norm(g, np.inf) 

158 

159 if res < self.newton_tol: 

160 break 

161 

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 ) 

166 

167 # newton update: u1 = u0 - g/dg 

168 u -= spsolve(dg, g) 

169 

170 # increase iteration count 

171 n += 1 

172 

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) 

177 

178 if n == self.newton_maxiter: 

179 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res)) 

180 

181 return u 

182 

183 def eval_f(self, u, t): 

184 """ 

185 Routine to evaluate the right-hand side of the problem. 

186 

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. 

193 

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 ) 

208 

209 uext = np.concatenate(([ul], u, [ur])) 

210 

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 

214 

215 def u_exact(self, t): 

216 r""" 

217 Routine to compute the exact solution at time :math:`t`. 

218 

219 Parameters 

220 ---------- 

221 t : float 

222 Time of the exact solution. 

223 

224 Returns 

225 ------- 

226 me : dtype_u 

227 Exact solution. 

228 """ 

229 

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)]) 

232 

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