Coverage for pySDC/implementations/sweeper_classes/Multistep.py: 93%
91 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-18 07:07 +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 alpha = None
57 beta = None
59 def __init__(self, params, level):
60 """
61 Initialization routine for the base sweeper.
63 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.
64 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.
65 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.
66 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.
67 Be careful with the sign of the alpha values. You can look at the implementations of the Euler methods for guidance.
69 Class attributes:
70 alpha (list): Coefficients for the solutions of previous steps
71 beta (list): Coefficients for the right hand side evaluations
73 Args:
74 params (dict): parameter object
75 level (pySDC.Level.level): the level that uses this sweeper
76 """
77 import logging
78 from pySDC.core.collocation import CollBase
80 # set up logger
81 self.logger = logging.getLogger('sweeper')
83 self.params = _Pars(params)
85 # check if some parameters are set which only apply to actual sweepers
86 for key in ['initial_guess', 'collocation_class', 'num_nodes', 'quad_type']:
87 if key in params:
88 self.logger.warning(f'"{key}" will be ignored by multistep sweeper')
90 # we need a dummy collocation object to instantiate the levels.
91 self.coll = CollBase(num_nodes=1, quad_type='RADAU-RIGHT')
93 self.__level = level
95 self.parallelizable = False
97 # proprietary variables for the multistep methods
98 self.steps = len(self.alpha)
99 self.cache = Cache(self.steps)
101 @property
102 def level(self):
103 """
104 Returns the current level
106 Returns:
107 pySDC.Level.level: Current level
108 """
109 return self.__level
111 @level.setter
112 def level(self, lvl):
113 """
114 Sets a reference to the current level (done in the initialization of the level)
116 Args:
117 lvl (pySDC.Level.level): Current level
118 """
119 assert isinstance(lvl, Level), f"You tried to set the sweeper's level with an instance of {type(lvl)}!"
121 self.__level = lvl
123 def predict(self):
124 """
125 Add the initial conditions to the cache if needed.
127 Default prediction for the sweepers, only copies the values to all collocation nodes
128 and evaluates the RHS of the ODE there
129 """
130 lvl = self.level
132 if all(me is None for me in self.cache.t):
133 prob = lvl.prob
134 lvl.f[0] = prob.eval_f(lvl.u[0], lvl.time)
135 self.cache.update(lvl.time, lvl.u[0], lvl.f[0])
137 lvl.status.unlocked = True
138 lvl.status.updated = True
140 def compute_residual(self, stage=None):
141 """
142 Do nothing.
144 Args:
145 stage (str): The current stage of the step the level belongs to
146 """
147 lvl = self.level
148 lvl.status.residual = 0.0
149 lvl.status.updated = False
151 return None
153 def compute_end_point(self):
154 """
155 The solution is stored in the single node that we have.
156 """
157 self.level.uend = self.level.u[-1]
159 def update_nodes(self):
160 """
161 Compute the solution to the current step. If the cache is not filled, we compute a provisional solution with a different method.
162 """
164 lvl = self.level
165 prob = lvl.prob
166 time = lvl.time + lvl.dt
168 if None in self.cache.t:
169 self.generate_starting_values()
171 else:
172 # build the right hand side from the previous solutions
173 rhs = prob.dtype_u(prob.init)
174 dts = [self.cache.t[i + 1] - self.cache.t[i] for i in range(self.steps - 1)] + [
175 time - self.cache.t[-1]
176 ] # TODO: Does this work for adaptive step sizes?
177 for i in range(len(self.alpha)):
178 rhs -= self.alpha[i] * self.cache.u[i]
179 rhs += dts[i] * self.beta[i] * self.cache.f[i]
181 # compute the solution to the current step "implicit Euler style"
182 lvl.u[1] = prob.solve_system(rhs, lvl.dt * self.beta[-1], self.cache.u[-1], time)
184 lvl.f[1] = prob.eval_f(lvl.u[1], time)
185 self.cache.update(time, lvl.u[1], lvl.f[1])
187 def generate_starting_values(self):
188 """
189 Compute solutions to the steps when not enough previous values are available for the multistep method.
190 The initial conditions are added in `predict` since this is not bespoke behaviour to any method.
191 """
192 raise NotImplementedError(
193 "No implementation for generating solutions when not enough previous values are available!"
194 )
197class AdamsBashforthExplicit1Step(MultiStep):
198 """
199 This is just forward Euler.
200 """
202 alpha = [-1.0]
203 beta = [1.0, 0.0]
206class BackwardEuler(MultiStep):
207 """
208 Almost as old, impressive and beloved as Koelner Dom.
209 """
211 alpha = [-1.0]
212 beta = [0.0, 1.0]
215class AdamsMoultonImplicit1Step(MultiStep):
216 """
217 Trapezoidal method dressed up as a multistep method.
218 """
220 alpha = [-1.0]
221 beta = [0.5, 0.5]
224class AdamsMoultonImplicit2Step(MultiStep):
225 """
226 Third order implicit scheme
227 """
229 alpha = [0.0, -1.0]
230 beta = [-1.0 / 12.0, 8.0 / 12.0, 5.0 / 12.0]
232 def generate_starting_values(self):
233 """
234 Generate starting value by trapezoidal rule.
235 """
236 lvl = self.level
237 prob = lvl.prob
238 time = lvl.time + lvl.dt
240 # do trapezoidal rule
241 rhs = lvl.u[0] + lvl.dt / 2 * lvl.f[0]
243 lvl.u[1] = prob.solve_system(rhs, lvl.dt / 2.0, lvl.u[0], time)