Coverage for pySDC/implementations/sweeper_classes/Multistep.py: 93%
91 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1from pySDC.core.sweeper import Sweeper, _Pars
2from pySDC.core.level import Level
5class Cache(object):
6 """
7 Class for managing solutions and right hand side evaluations of previous steps for the MultiStep "sweeper".
9 Attributes:
10 - u (list): Contains solution from previous steps
11 - f (list): Contains right hand side evaluations from previous steps
12 - t (list): Contains time of previous steps
13 """
15 def __init__(self, num_steps):
16 """
17 Initialization routing
19 Args:
20 num_steps (int): Number of entries for the cache
21 """
22 self.num_steps = num_steps
23 self.u = [None] * num_steps
24 self.f = [None] * num_steps
25 self.t = [None] * num_steps
27 def update(self, t, u, f):
28 """
29 Add a new value to the cache and remove the oldest one.
31 Args:
32 t (float): Time of new step
33 u (dtype_u): Solution of new step
34 f (dtype_f): Right hand side evaluation at new step
35 """
36 self.u[:-1] = self.u[1:]
37 self.f[:-1] = self.f[1:]
38 self.t[:-1] = self.t[1:]
40 self.u[-1] = u
41 self.f[-1] = f
42 self.t[-1] = t
44 def __str__(self):
45 """
46 Print the contents of the cache for debugging purposes.
47 """
48 string = ''
49 for i in range(self.num_steps):
50 string = f'{string} t={self.t[i]}: u={self.u[i]}, f={self.f[i]}'
52 return string
55class MultiStep(Sweeper):
56 """
57 Base class for linear multistep methods, given by alpha and beta coefficients and a cache of previous steps.
58 """
60 alpha = None
61 beta = None
63 def __init__(self, params, level):
64 """
65 Initialization routine for the base sweeper.
67 Multistep methods work by constructing Euleresque steps with the right hand side composed of a weighted sum of solutions and right hand side evaluations at previous steps.
68 The coefficients in the weighted sum for the solutions are called alpha and the ones for the right hand sides are called beta in this implementation.
69 So for an N step method, there need to be N alpha values and N + 1 beta values, where the last one is on the left hand side for implicit methods.
70 The first element in the coefficients belongs to the value furthest in the past and vice versa. Values from previous time steps are stored in a `Cache` object.
71 Be careful with the sign of the alpha values. You can look at the implementations of the Euler methods for guidance.
73 Class attributes:
74 alpha (list): Coefficients for the solutions of previous steps
75 beta (list): Coefficients for the right hand side evaluations
77 Args:
78 params (dict): parameter object
79 level (pySDC.Level.level): the level that uses this sweeper
80 """
81 import logging
82 from pySDC.core.collocation import CollBase
84 # set up logger
85 self.logger = logging.getLogger('sweeper')
87 self.params = _Pars(params)
89 # check if some parameters are set which only apply to actual sweepers
90 for key in ['initial_guess', 'collocation_class', 'num_nodes', 'quad_type']:
91 if key in params:
92 self.logger.warning(f'"{key}" will be ignored by multistep sweeper')
94 # we need a dummy collocation object to instantiate the levels.
95 self.coll = CollBase(num_nodes=1, quad_type='RADAU-RIGHT')
97 self.__level = level
99 self.parallelizable = False
101 # proprietary variables for the multistep methods
102 self.steps = len(self.alpha)
103 self.cache = Cache(self.steps)
105 @property
106 def level(self):
107 """
108 Returns the current level
110 Returns:
111 pySDC.Level.level: Current level
112 """
113 return self.__level
115 @level.setter
116 def level(self, lvl):
117 """
118 Sets a reference to the current level (done in the initialization of the level)
120 Args:
121 lvl (pySDC.Level.level): Current level
122 """
123 assert isinstance(lvl, Level), f"You tried to set the sweeper's level with an instance of {type(lvl)}!"
125 self.__level = lvl
127 def predict(self):
128 """
129 Add the initial conditions to the cache if needed.
131 Default prediction for the sweepers, only copies the values to all collocation nodes
132 and evaluates the RHS of the ODE there
133 """
134 lvl = self.level
136 if all(me is None for me in self.cache.t):
137 prob = lvl.prob
138 lvl.f[0] = prob.eval_f(lvl.u[0], lvl.time)
139 self.cache.update(lvl.time, lvl.u[0], lvl.f[0])
141 lvl.status.unlocked = True
142 lvl.status.updated = True
144 def compute_residual(self, stage=None):
145 """
146 Do nothing.
148 Args:
149 stage (str): The current stage of the step the level belongs to
150 """
151 lvl = self.level
152 lvl.status.residual = 0.0
153 lvl.status.updated = False
155 return None
157 def compute_end_point(self):
158 """
159 The solution is stored in the single node that we have.
160 """
161 self.level.uend = self.level.u[-1]
163 def update_nodes(self):
164 """
165 Compute the solution to the current step. If the cache is not filled, we compute a provisional solution with a different method.
166 """
168 lvl = self.level
169 prob = lvl.prob
170 time = lvl.time + lvl.dt
172 if None in self.cache.t:
173 self.generate_starting_values()
175 else:
176 # build the right hand side from the previous solutions
177 rhs = prob.dtype_u(prob.init)
178 dts = [self.cache.t[i + 1] - self.cache.t[i] for i in range(self.steps - 1)] + [
179 time - self.cache.t[-1]
180 ] # TODO: Does this work for adaptive step sizes?
181 for i in range(len(self.alpha)):
182 rhs -= self.alpha[i] * self.cache.u[i]
183 rhs += dts[i] * self.beta[i] * self.cache.f[i]
185 # compute the solution to the current step "implicit Euler style"
186 lvl.u[1] = prob.solve_system(rhs, lvl.dt * self.beta[-1], self.cache.u[-1], time)
188 lvl.f[1] = prob.eval_f(lvl.u[1], time)
189 self.cache.update(time, lvl.u[1], lvl.f[1])
191 def generate_starting_values(self):
192 """
193 Compute solutions to the steps when not enough previous values are available for the multistep method.
194 The initial conditions are added in `predict` since this is not bespoke behaviour to any method.
195 """
196 raise NotImplementedError(
197 "No implementation for generating solutions when not enough previous values are available!"
198 )
201class AdamsBashforthExplicit1Step(MultiStep):
202 """
203 One-step Adams-Bashforth method, which is just forward Euler.
204 """
206 alpha = [-1.0]
207 beta = [1.0, 0.0]
210class BackwardEuler(MultiStep):
211 """
212 Backward Euler written as a one-step implicit multistep method.
213 Almost as old, impressive and beloved as Koelner Dom.
214 """
216 alpha = [-1.0]
217 beta = [0.0, 1.0]
220class AdamsMoultonImplicit1Step(MultiStep):
221 """
222 Trapezoidal method dressed up as a multistep method.
223 """
225 alpha = [-1.0]
226 beta = [0.5, 0.5]
229class AdamsMoultonImplicit2Step(MultiStep):
230 """
231 Third-order implicit two-step Adams-Moulton method, started with a trapezoidal rule step.
232 """
234 alpha = [0.0, -1.0]
235 beta = [-1.0 / 12.0, 8.0 / 12.0, 5.0 / 12.0]
237 def generate_starting_values(self):
238 """
239 Generate starting value by trapezoidal rule.
240 """
241 lvl = self.level
242 prob = lvl.prob
243 time = lvl.time + lvl.dt
245 # do trapezoidal rule
246 rhs = lvl.u[0] + lvl.dt / 2 * lvl.f[0]
248 lvl.u[1] = prob.solve_system(rhs, lvl.dt / 2.0, lvl.u[0], time)