# Coverage for pySDC/implementations/problem_classes/odeScalar.py: 100%

## 53 statements

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

1#!/usr/bin/env python3

2# -*- coding: utf-8 -*-

3"""

4Implementation of scalar test problem ODEs.

7Reference :

9Van der Houwen, P. J., & Sommeijer, B. P. (1991). Iterated Runge–Kutta methods

10on parallel computers. SIAM journal on scientific and statistical computing,

1112(5), 1000-1028.

12"""

13import numpy as np

15from pySDC.core.errors import ProblemError

16from pySDC.core.problem import Problem, WorkCounter

17from pySDC.implementations.datatype_classes.mesh import mesh

20class ProtheroRobinson(Problem):

21 r"""

22 Implement the Prothero-Robinson problem:

24 .. math::

25 \frac{du}{dt} = -\frac{u-g(t)}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0).,

27 with :math:\epsilon a stiffness parameter, that makes the problem more stiff

28 the smaller it is (usual taken value is :math:\epsilon=1e^{-3}).

29 Exact solution is given by :math:u(t)=g(t), and this implementation uses

30 :math:g(t)=\cos(t).

32 Implement also the non-linear form of this problem:

34 .. math::

35 \frac{du}{dt} = -\frac{u^3-g(t)^3}{\epsilon} + \frac{dg}{dt}, \quad u(0) = g(0).

37 To use an other exact solution, one just have to derivate this class

38 and overload the g and dg methods. For instance,

39 to use :math:g(t)=e^{-0.2*t}, define and use the following class:

41 >>> class MyProtheroRobinson(ProtheroRobinson):

42 >>>

43 >>> def g(self, t):

44 >>> return np.exp(-0.2 * t)

45 >>>

46 >>> def dg(self, t):

47 >>> return (-0.2) * np.exp(-0.2 * t)

49 Parameters

50 ----------

51 epsilon : float, optional

52 Stiffness parameter. The default is 1e-3.

53 nonLinear : bool, optional

54 Wether or not to use the non-linear form of the problem. The default is False.

55 newton_maxiter : int, optional

56 Maximum number of Newton iteration in solve_system. The default is 200.

57 newton_tol : float, optional

58 Residuum tolerance for Newton iteration in solve_system. The default is 5e-11.

59 stop_at_nan : bool, optional

60 Wheter to stop or not solve_system when getting NAN. The default is True.

62 Reference

63 ---------

64 A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving

65 stiff systems of ordinary differential equations, Mathematics of Computation, 28 (1974),

66 pp. 145–162.

67 """

69 dtype_u = mesh

70 dtype_f = mesh

72 def __init__(self, epsilon=1e-3, nonLinear=False, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):

73 nvars = 1

74 super().__init__((nvars, None, np.dtype('float64')))

76 self.f = self.f_NONLIN if nonLinear else self.f_LIN

77 self.jac = self.jac_NONLIN if nonLinear else self.jac_LIN

78 self._makeAttributeAndRegister(

79 'epsilon', 'nonLinear', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True

80 )

81 self.work_counters['newton'] = WorkCounter()

82 self.work_counters['rhs'] = WorkCounter()

84 # -------------------------------------------------------------------------

85 # g function (analytical solution), and its first derivative

86 # -------------------------------------------------------------------------

87 def g(self, t):

88 return np.cos(t)

90 def dg(self, t):

91 return -np.sin(t)

93 # -------------------------------------------------------------------------

94 # f(u,t) and Jacobian functions

95 # -------------------------------------------------------------------------

96 def f(self, u, t):

97 raise NotImplementedError()

99 def f_LIN(self, u, t):

100 return -self.epsilon ** (-1) * (u - self.g(t)) + self.dg(t)

102 def f_NONLIN(self, u, t):

103 return -self.epsilon ** (-1) * (u**3 - self.g(t) ** 3) + self.dg(t)

105 def jac(self, u, t):

106 raise NotImplementedError()

108 def jac_LIN(self, u, t):

109 return -self.epsilon ** (-1)

111 def jac_NONLIN(self, u, t):

112 return -self.epsilon ** (-1) * 3 * u**2

114 # -------------------------------------------------------------------------

115 # pySDC required methods

116 # -------------------------------------------------------------------------

117 def u_exact(self, t, u_init=None, t_init=None):

118 r"""

119 Routine to return initial conditions or exact solution.

121 Parameters

122 ----------

123 t : float

124 Time at which the exact solution is computed.

125 u_init : dtype_u

126 Initial conditions for getting the exact solution.

127 t_init : float

128 The starting time.

130 Returns

131 -------

132 u : dtype_u

133 The exact solution.

134 """

135 u = self.dtype_u(self.init)

136 u[:] = self.g(t)

137 return u

139 def eval_f(self, u, t):

140 """

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

143 Parameters

144 ----------

145 u : dtype_u

146 Current values of the numerical solution.

147 t : float

148 Current time of the numerical solution is computed (not used here).

150 Returns

151 -------

152 f : dtype_f

153 The right-hand side of the problem (one component).

154 """

156 f = self.dtype_f(self.init)

157 f[:] = self.f(u, t)

158 self.work_counters['rhs']()

159 return f

161 def solve_system(self, rhs, dt, u0, t):

162 """

163 Simple Newton solver for the nonlinear equation

165 Parameters

166 ----------

167 rhs : dtype_f

168 Right-hand side for the nonlinear system.

169 dt : float

170 Abbrev. for the node-to-node stepsize (or any other factor required).

171 u0 : dtype_u

172 Initial guess for the iterative solver.

173 t : float

174 Time of the updated solution (e.g. for time-dependent BCs).

176 Returns

177 -------

178 u : dtype_u

179 The solution as mesh.

180 """

181 # create new mesh object from u0 and set initial values for iteration

182 u = self.dtype_u(u0)

184 # start newton iteration

185 n, res = 0, np.inf

186 while n < self.newton_maxiter:

187 # form the function g with g(u) = 0

188 g = u - dt * self.f(u, t) - rhs

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

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

192 if res < self.newton_tol or np.isnan(res):

193 break

195 # assemble dg/du

196 dg = 1 - dt * self.jac(u, t)

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

199 u -= dg ** (-1) * g

201 # increase iteration count and work counter

202 n += 1

203 self.work_counters['newton']()

205 if np.isnan(res) and self.stop_at_nan:

206 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)

207 elif np.isnan(res): # pragma: no cover

208 self.logger.warning('Newton got nan after %i iterations...' % n)

210 if n == self.newton_maxiter:

211 raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))

213 return u