Coverage for pySDC/implementations/problem_classes/nonlinear_ODE_1.py: 100%
41 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +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 """Initialization routine"""
47 nvars = 1
48 super().__init__((nvars, None, np.dtype('float64')))
49 self._makeAttributeAndRegister(
50 'u0', 'newton_maxiter', 'newton_tol', 'stop_at_nan', localVars=locals(), readOnly=True
51 )
53 self.newton_itercount = 0
54 self.newton_ncalls = 0
56 def u_exact(self, t):
57 r"""
58 Routine to compute the exact solution at time :math:`t`.
60 Parameters
61 ----------
62 t : float
63 Time of the exact solution.
65 Returns
66 -------
67 me : dtype_u
68 The exact solution.
69 """
71 me = self.dtype_u(self.init)
72 me[:] = t - t**2 / 4
73 return me
75 def eval_f(self, u, t):
76 """
77 Routine to evaluate the right-hand side of the problem.
79 Parameters
80 ----------
81 u : dtype_u
82 Current values of the numerical solution.
83 t : float
84 Current time of the numerical solution is computed (not used here).
86 Returns
87 -------
88 f : dtype_f
89 The right-hand side of the problem (one component).
90 """
92 f = self.dtype_f(self.init)
93 f[:] = np.sqrt(1 - u)
94 return f
96 def solve_system(self, rhs, dt, u0, t):
97 """
98 Simple Newton solver for the nonlinear equation
100 Parameters
101 ----------
102 rhs : dtype_f
103 Right-hand side for the nonlinear system.
104 dt : float
105 Abbrev. for the node-to-node stepsize (or any other factor required).
106 u0 : dtype_u
107 Initial guess for the iterative solver.
108 t : float
109 Current time (e.g. for time-dependent BCs).
111 Returns
112 -------
113 u : dtype_u
114 The solution as mesh.
115 """
116 # create new mesh object from u0 and set initial values for iteration
117 u = self.dtype_u(u0)
119 # start newton iteration
120 n = 0
121 res = 99
122 while n < self.newton_maxiter:
123 # form the function g with g(u) = 0
124 g = u - dt * np.sqrt(1 - u) - rhs
126 # if g is close to 0, then we are done
127 res = np.linalg.norm(g, np.inf)
128 if res < self.newton_tol or np.isnan(res):
129 break
131 # assemble dg/du
132 dg = 1 - (-dt) / (2 * np.sqrt(1 - u))
133 # newton update: u1 = u0 - g/dg
134 u -= 1.0 / dg * g
136 # increase iteration count
137 n += 1
139 self.newton_itercount += 1
141 if np.isnan(res) and self.stop_at_nan:
142 raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
143 elif np.isnan(res):
144 self.logger.warning('Newton got nan after %i iterations...' % n)
146 if n == self.newton_maxiter:
147 self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
149 self.newton_ncalls += 1
151 return u