# Coverage for pySDC/implementations/problem_classes/nonlinear_ODE_1.py: 95%

## 41 statements

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

1import numpy as np

3from pySDC.core.errors import ProblemError

4from pySDC.core.problem import Problem

5from pySDC.implementations.datatype_classes.mesh import mesh

8# noinspection PyUnusedLocal

9class nonlinear_ODE_1(Problem):

10 r"""

11 This class implements a simple nonlinear ODE with a singularity in the derivative, taken from

12 https://www.osti.gov/servlets/purl/6111421 (Problem E-4). For :math:0 \leq t \leq 5, the problem is

13 given by

15 .. math::

16 \frac{du(t)}{dt} = \sqrt{1 - u(t)}

18 with initial condition :math:u(0) = 0. The exact solution is

20 .. math::

21 u(t) = t - \frac{t^2}{4}.

23 Parameters

24 ----------

25 u0 : float, optional

26 Initial condition.

27 newton_maxiter : float, optional

28 Maximum number of iterations for Newton's method.

29 newton_tol : float, optional

30 Tolerance for Newton's method to terminate.

31 stop_at_nan : bool, optional

32 Indicates that Newton solver has to stop if nan values arise.

34 Attributes

35 ----------

36 newton_itercount : int

37 Counts the Newton iterations.

38 newton_ncalls : int

39 Counts calls of Newton method.

40 """

42 dtype_u = mesh

43 dtype_f = mesh

45 def __init__(self, u0=0.0, newton_maxiter=200, newton_tol=5e-11, stop_at_nan=True):

46 nvars = 1

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

48 self._makeAttributeAndRegister(

49 'u0', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True

50 )

52 self.newton_itercount = 0

53 self.newton_ncalls = 0

55 def u_exact(self, t):

56 r"""

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

59 Parameters

60 ----------

61 t : float

62 Time of the exact solution.

64 Returns

65 -------

66 me : dtype_u

67 The exact solution.

68 """

70 me = self.dtype_u(self.init)

71 me[:] = t - t**2 / 4

72 return me

74 def eval_f(self, u, t):

75 """

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

78 Parameters

79 ----------

80 u : dtype_u

81 Current values of the numerical solution.

82 t : float

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

85 Returns

86 -------

87 f : dtype_f

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

89 """

91 f = self.dtype_f(self.init)

92 f[:] = np.sqrt(1 - u)

93 return f

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

96 """

97 Simple Newton solver for the nonlinear equation

99 Parameters

100 ----------

101 rhs : dtype_f

102 Right-hand side for the nonlinear system.

103 dt : float

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

105 u0 : dtype_u

106 Initial guess for the iterative solver.

107 t : float

108 Current time (e.g. for time-dependent BCs).

110 Returns

111 -------

112 u : dtype_u

113 The solution as mesh.

114 """

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

116 u = self.dtype_u(u0)

118 # start newton iteration

119 n = 0

120 res = 99

121 while n < self.newton_maxiter:

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

123 g = u - dt * np.sqrt(1 - u) - rhs

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

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

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

128 break

130 # assemble dg/du

131 dg = 1 - (-dt) / (2 * np.sqrt(1 - u))

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

133 u -= 1.0 / dg * g

135 # increase iteration count

136 n += 1

138 self.newton_itercount += 1

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

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

142 elif np.isnan(res):

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

145 if n == self.newton_maxiter:

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

148 self.newton_ncalls += 1

150 return u