Coverage for pySDC/core/problem.py: 96%
47 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 13:41 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Description
5-----------
7Module containing the base Problem class for pySDC
8"""
10import logging
11from typing import Any, Dict, Optional, Type, Callable
13from pySDC.core.common import RegisterParams
16class WorkCounter(object):
17 """
18 Utility class for counting iterations.
20 Contains one attribute `niter` initialized to zero during
21 instantiation, which can be incremented by calling object as
22 a function, e.g
24 >>> count = WorkCounter() # => niter = 0
25 >>> count() # => niter = 1
26 >>> count() # => niter = 2
27 """
29 def __init__(self) -> None:
30 """Start counting at zero."""
31 self.niter: int = 0
33 def __call__(self, *args: Any, **kwargs: Any) -> None:
34 # *args and **kwargs are necessary for gmres
35 self.niter += 1
37 def decrement(self) -> None:
38 """Decrease `niter` by one, e.g. to take back a count."""
39 self.niter -= 1
41 def __str__(self) -> str:
42 return f'{self.niter}'
45class Problem(RegisterParams):
46 r"""
47 Prototype class for problems, just defines the attributes essential to get started.
49 Parameters
50 ----------
51 init : list of args
52 Argument(s) used to initialize data types.
54 Attributes
55 ----------
56 logger: logging.Logger
57 custom logger for problem-related logging.
58 dtype_u : type
59 Variable data type. Should generate a data variable using dtype_u(init).
60 dtype_f : type
61 RHS data type. Should generate a data variable using dtype_f(init).
63 Notes
64 -----
65 **The correction-solve contract.** The deferred-correction ("delta-form") sweepers in
66 :mod:`pySDC.implementations.sweeper_classes.delta_form` hand the implicit solve a *correction*
67 equation instead of the state's. That is what lets the solver -- typically the expensive, often
68 external part -- work in, and return, a much lower precision while the iteration still
69 converges to full precision. Both sides of that call are specified here; nothing in it is
70 required, and a problem that ignores it keeps working.
72 The call is made per collocation node and sweep, by one of two routes:
74 - linear or affine implicit operator: the stock :meth:`solve_system`, with the sweeper
75 parameter ``linear_implicit=True``;
76 - nonlinear implicit operator: ``solve_system_delta(rhs, factor, base, f_base, t)``, which
77 solves :math:`\delta - factor\,[f(base+\delta) - f(base)] = rhs` and returns :math:`\delta`.
78 This is the contract for an external nonlinear solver.
80 Without either, the sweeper substitutes :math:`y = u^k_m + \delta_m` and calls
81 :meth:`solve_system`. That is exact, but the solver then sees an :math:`\mathcal{O}(1)`
82 unknown, and none of what follows holds.
84 *The sweeper guarantees the solver*
86 1. the unknown is a correction, the right-hand side is that of the correction equation, and it
87 shrinks with the iteration like the residual; the initial guess is zero;
88 2. on the linear route, a right-hand side of unit size: it is divided by its maximum before
89 the call and the result multiplied by it after, which is exact and keeps any format clear of
90 underflow. On the nonlinear route the scaled unknown would enter :math:`f`, so there the
91 solver scales internally if its format needs it;
92 3. ``base`` and ``f_base`` at the level's own precision, read only, and ``factor`` a plain
93 scalar;
94 4. that everything outside the solve -- the residual, the update :math:`u \leftarrow u +
95 \delta`, and the increment :math:`f(u+\delta) - f(u)` -- stays at the level's precision, so
96 every sweep is a step of iterative refinement around the solve.
98 *The solver guarantees the sweeper*
100 1. it returns the correction, not the state, in any precision: the sweeper widens it to the
101 level's before using it;
102 2. a relative accuracy :math:`\eta`, :math:`\|\tilde\delta - \delta\| \le \eta\|\delta\|`, below
103 what the method needs. How it gets there -- its internal precision, iterative refinement, a
104 Krylov tolerance, multigrid cycles -- is its own business;
105 3. relative tolerances only, and none below what its own precision can deliver: an absolute
106 tolerance caps the attainable accuracy, and an unreachable one costs iterations or, in a
107 Krylov method, breaks it down;
108 4. no boundary data or affine terms re-imposed on the correction, whose boundary conditions are
109 homogeneous.
111 Then the iteration converges to the level's precision, at an iteration cost set by
112 :math:`\eta` alone. Measured on the heat equation: SDC needs :math:`\eta` of about 1e-4, an
113 MLSDC fine level about 1e-6 and an MLSDC coarse level about 1e-2 to lose at most one iteration.
115 One more optional method serves the contract. ``eval_f_increment(base, delta, t)`` returns
116 :math:`f(base+\delta) - f(base)` with the same splitting as :meth:`eval_f`, expanded so that
117 every term carries an explicit factor :math:`\delta`. Without it the sweeper subtracts two
118 stored right-hand sides, a cancellation that carries the operator norm: on the FEniCS heat
119 equation, where :math:`|M^{-1}Ku|` is of order :math:`10^5`, it is already 7.7e-11 off in
120 double precision, and below double precision it binds before the solver does.
121 """
123 logger: logging.Logger = logging.getLogger('problem')
124 dtype_u: Optional[Type[Any]] = None
125 dtype_f: Optional[Type[Any]] = None
127 def __init__(self, init: Any) -> None:
128 """Store `init`, see the class docstring, and start without work counters."""
129 self.work_counters: Dict[str, WorkCounter] = {} # Dictionary to store WorkCounter objects
130 self.init: Any = init # Initialization parameter to instantiate data types
132 @property
133 def u_init(self) -> Any:
134 """Generate a data variable for u"""
135 return self.dtype_u(self.init)
137 @property
138 def f_init(self) -> Any:
139 """Generate a data variable for RHS"""
140 return self.dtype_f(self.init)
142 @classmethod
143 def get_default_sweeper_class(cls) -> Type[Any]:
144 """
145 Get the sweeper class to use with this problem if none is given, e.g. in `pySDC.helpers.setup_helper`.
147 Returns
148 -------
149 type
150 Sweeper class.
152 Raises
153 ------
154 NotImplementedError
155 If the problem class does not define a default.
156 """
157 raise NotImplementedError(f'No default sweeper class implemented for {cls} problem!')
159 def setUpFieldsIO(self) -> None:
160 """
161 Set up FieldsIO for MPI with the space decomposition of this problem
162 """
163 pass
165 def getOutputFile(self, fileName: str) -> Any:
166 """
167 Create the file the `LogToFile` hook writes the solution to. Implemented by problems that support file output.
169 Parameters
170 ----------
171 fileName : str
172 Name of the file.
174 Returns
175 -------
176 pySDC.helpers.fieldsIO.FieldsIO
177 File object with the header set and initialized.
178 """
179 raise NotImplementedError(f'No output implemented file for {type(self).__name__}')
181 def processSolutionForOutput(self, u: Any) -> Any:
182 """
183 Convert a solution to the data written to the file from `getOutputFile`.
185 Parameters
186 ----------
187 u : dtype_u
188 Solution.
190 Returns
191 -------
192 np.ndarray
193 Data to write, `u` itself by default.
194 """
195 return u
197 def eval_f(self, u: Any, t: float) -> Any:
198 """
199 Abstract interface to RHS computation of the ODE
201 Parameters
202 ----------
203 u : dtype_u
204 Current values.
205 t : float
206 Current time.
208 Returns
209 -------
210 f : dtype_f
211 The RHS values.
212 """
213 raise NotImplementedError('ERROR: problem has to implement eval_f(self, u, t)')
215 def apply_mass_matrix(self, u: Any) -> Any:
216 r"""
217 Apply the problem's mass matrix, :math:`M \vec{u}`.
219 There is deliberately no default implementation. The mass-matrix sweepers and transfers
220 (``generic_implicit_mass``, ``imex_1st_order_mass``, ``base_transfer_mass``) build ``u0``,
221 the FAS ``tau`` and the residual out of this, so a problem that does not implement it would
222 silently compute the wrong thing instead of failing. The identity default that used to sit
223 here did exactly that: it produced ``nan`` several sweeps later, with nothing pointing back
224 to the missing method.
225 """
226 raise NotImplementedError(
227 f'{type(self).__name__} is being used in a mass-matrix formulation but does not '
228 'implement apply_mass_matrix'
229 )
231 def generate_scipy_reference_solution(
232 self, eval_rhs: Callable, t: float, u_init: Optional[Any] = None, t_init: Optional[float] = None, **kwargs: Any
233 ) -> Any:
234 """
235 Compute a reference solution using `scipy.solve_ivp` with very small tolerances.
236 Keep in mind that scipy needs the solution to be a one dimensional array. If you are solving something higher
237 dimensional, you need to make sure the function `eval_rhs` takes a flattened one-dimensional version as an input
238 and output, but reshapes to whatever the problem needs for evaluation.
240 The keyword arguments will be passed to `scipy.solve_ivp`. You should consider passing `method='BDF'` for stiff
241 problems and to accelerate that you can pass a function that evaluates the Jacobian with arguments `jac(t, u)`
242 as `jac=jac`.
244 Args:
245 eval_rhs (function): Function evaluate the full right hand side. Must have signature `eval_rhs(float: t, numpy.1darray: u)`
246 t (float): current time
247 u_init (pySDC.implementations.problem_classes.Lorenz.dtype_u): initial conditions for getting the exact solution
248 t_init (float): the starting time
250 Returns:
251 numpy.ndarray: Reference solution
252 """
253 import numpy as np
254 from scipy.integrate import solve_ivp
256 kwargs = {
257 'atol': 100 * np.finfo(float).eps,
258 'rtol': 100 * np.finfo(float).eps,
259 **kwargs,
260 }
261 u_init = self.u_exact(t=0) if u_init is None else u_init * 1.0
262 t_init = 0 if t_init is None else t_init
264 u_shape = u_init.shape
265 return solve_ivp(eval_rhs, (t_init, t), u_init.flatten(), **kwargs).y[:, -1].reshape(u_shape)
267 def get_fig(self) -> Any:
268 """
269 Get a figure suitable to plot the solution of this problem
271 Returns
272 -------
273 self.fig : matplotlib.pyplot.figure.Figure
274 """
275 raise NotImplementedError
277 def plot(self, u: Any, t: Optional[float] = None, fig: Optional[Any] = None) -> None:
278 r"""
279 Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
281 Parameters
282 ----------
283 u : dtype_u
284 Solution to be plotted
285 t : float
286 Time to display at the top of the figure
287 fig : matplotlib.pyplot.figure.Figure
288 Figure with the correct structure
290 Returns
291 -------
292 None
293 """
294 raise NotImplementedError
296 def solve_system(self, rhs: Any, dt: float, u0: Any, t: float) -> Any:
297 """
298 Perform an Euler step.
300 Args:
301 rhs: Right hand side for the Euler step
302 dt (float): Step size for the Euler step
303 u0: Initial guess
304 t (float): Current time
306 Returns:
307 solution to the Euler step
308 """
309 raise NotImplementedError
311 def use_restricted_boundary_values(self, u: list, t: list) -> None:
312 """
313 Called by the FAS restriction on the coarse problem with the restricted fine solution `u` at the coarse
314 collocation nodes `t` (one entry per node held by this rank). Does nothing by default.
316 A problem that imposes boundary data in `solve_system` should impose the boundary values of `u` in its solves
317 at these times instead, until the next call. Then the coarse correction has homogeneous boundary conditions, as
318 in multigrid, and the restricted converged fine solution is a fixed point of the coarse sweep. The tau
319 correction cannot do this, because the boundary rows of the right hand side are replaced by the data. With
320 its own data at the coarse nodes, MLSDC and PFASST converge to a different solution than SDC whenever the data
321 depends on time: the restricted boundary values interpolate the data at the fine nodes, which is not the data
322 at the coarse nodes.
324 Args:
325 u (list): Restricted solution at the coarse nodes
326 t (list): Times of the coarse nodes
327 """
328 pass
330 def solve_jacobian(
331 self, rhs: Any, dt: float, u: Optional[Any] = None, u0: Optional[Any] = None, t: float = 0, **kwargs: Any
332 ) -> Any:
333 """
334 Solve the Jacobian for an Euler step, linearized around u.
335 This defaults to an Euler step to accommodate linear problems.
337 Args:
338 rhs: Right hand side for the Euler step
339 dt (float): Step size for the Euler step
340 u: Solution to linearize around
341 u0: Initial guess
342 t (float): Current time
344 Returns:
345 Solution
346 """
347 return self.solve_system(rhs, dt, u0=u, t=t, **kwargs)