Coverage for pySDC/implementations/problem_classes/nonlinear_ODE_1.py: 95%
41 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
3from pySDC.core.Errors import ProblemError
4from pySDC.core.Problem import ptype
5from pySDC.implementations.datatype_classes.mesh import mesh
8# noinspection PyUnusedLocal
9class nonlinear_ODE_1(ptype):
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