Coverage for pySDC/implementations/sweeper_classes/Multistep.py: 95%
83 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1from pySDC.core.sweeper import Sweeper, _Pars
4class Cache(object):
5 """
6 Class for managing solutions and right hand side evaluations of previous steps for the MultiStep "sweeper".
8 Attributes:
9 - u (list): Contains solution from previous steps
10 - f (list): Contains right hand side evaluations from previous steps
11 - t (list): Contains time of previous steps
12 """
14 def __init__(self, num_steps):
15 """
16 Initialization routing
18 Args:
19 num_steps (int): Number of entries for the cache
20 """
21 self.num_steps = num_steps
22 self.u = [None] * num_steps
23 self.f = [None] * num_steps
24 self.t = [None] * num_steps
26 def update(self, t, u, f):
27 """
28 Add a new value to the cache and remove the oldest one.
30 Args:
31 t (float): Time of new step
32 u (dtype_u): Solution of new step
33 f (dtype_f): Right hand side evaluation at new step
34 """
35 self.u[:-1] = self.u[1:]
36 self.f[:-1] = self.f[1:]
37 self.t[:-1] = self.t[1:]
39 self.u[-1] = u
40 self.f[-1] = f
41 self.t[-1] = t
43 def __str__(self):
44 """
45 Print the contents of the cache for debugging purposes.
46 """
47 string = ''
48 for i in range(self.num_steps):
49 string = f'{string} t={self.t[i]}: u={self.u[i]}, f={self.f[i]}'
51 return string
54class MultiStep(Sweeper):
55 alpha = None
56 beta = None
58 def __init__(self, params):
59 """
60 Initialization routine for the base sweeper.
62 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.
63 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.
64 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.
65 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.
66 Be careful with the sign of the alpha values. You can look at the implementations of the Euler methods for guidance.
68 Class attributes:
69 alpha (list): Coefficients for the solutions of previous steps
70 beta (list): Coefficients for the right hand side evaluations
72 Args:
73 params (dict): parameter object
74 """
75 import logging
76 from pySDC.core.collocation import CollBase
78 # set up logger
79 self.logger = logging.getLogger('sweeper')
81 self.params = _Pars(params)
83 # check if some parameters are set which only apply to actual sweepers
84 for key in ['initial_guess', 'collocation_class', 'num_nodes', 'quad_type']:
85 if key in params:
86 self.logger.warning(f'"{key}" will be ignored by multistep sweeper')
88 # we need a dummy collocation object to instantiate the levels.
89 self.coll = CollBase(num_nodes=1, quad_type='RADAU-RIGHT')
91 # This will be set as soon as the sweeper is instantiated at the level
92 self.__level = None
94 self.parallelizable = False
96 # proprietary variables for the multistep methods
97 self.steps = len(self.alpha)
98 self.cache = Cache(self.steps)
100 def predict(self):
101 """
102 Add the initial conditions to the cache if needed.
104 Default prediction for the sweepers, only copies the values to all collocation nodes
105 and evaluates the RHS of the ODE there
106 """
107 lvl = self.level
109 if all(me is None for me in self.cache.t):
110 prob = lvl.prob
111 lvl.f[0] = prob.eval_f(lvl.u[0], lvl.time)
112 self.cache.update(lvl.time, lvl.u[0], lvl.f[0])
114 lvl.status.unlocked = True
115 lvl.status.updated = True
117 def compute_residual(self, stage=None):
118 """
119 Do nothing.
121 Args:
122 stage (str): The current stage of the step the level belongs to
123 """
124 lvl = self.level
125 lvl.status.residual = 0.0
126 lvl.status.updated = False
128 return None
130 def compute_end_point(self):
131 """
132 The solution is stored in the single node that we have.
133 """
134 self.level.uend = self.level.u[-1]
136 def update_nodes(self):
137 """
138 Compute the solution to the current step. If the cache is not filled, we compute a provisional solution with a different method.
139 """
141 lvl = self.level
142 prob = lvl.prob
143 time = lvl.time + lvl.dt
145 if None in self.cache.t:
146 self.generate_starting_values()
148 else:
149 # build the right hand side from the previous solutions
150 rhs = prob.dtype_u(prob.init)
151 dts = [self.cache.t[i + 1] - self.cache.t[i] for i in range(self.steps - 1)] + [
152 time - self.cache.t[-1]
153 ] # TODO: Does this work for adaptive step sizes?
154 for i in range(len(self.alpha)):
155 rhs -= self.alpha[i] * self.cache.u[i]
156 rhs += dts[i] * self.beta[i] * self.cache.f[i]
158 # compute the solution to the current step "implicit Euler style"
159 lvl.u[1] = prob.solve_system(rhs, lvl.dt * self.beta[-1], self.cache.u[-1], time)
161 lvl.f[1] = prob.eval_f(lvl.u[1], time)
162 self.cache.update(time, lvl.u[1], lvl.f[1])
164 def generate_starting_values(self):
165 """
166 Compute solutions to the steps when not enough previous values are available for the multistep method.
167 The initial conditions are added in `predict` since this is not bespoke behaviour to any method.
168 """
169 raise NotImplementedError(
170 "No implementation for generating solutions when not enough previous values are available!"
171 )
174class AdamsBashforthExplicit1Step(MultiStep):
175 """
176 This is just forward Euler.
177 """
179 alpha = [-1.0]
180 beta = [1.0, 0.0]
183class BackwardEuler(MultiStep):
184 """
185 Almost as old, impressive and beloved as Koelner Dom.
186 """
188 alpha = [-1.0]
189 beta = [0.0, 1.0]
192class AdamsMoultonImplicit1Step(MultiStep):
193 """
194 Trapezoidal method dressed up as a multistep method.
195 """
197 alpha = [-1.0]
198 beta = [0.5, 0.5]
201class AdamsMoultonImplicit2Step(MultiStep):
202 """
203 Third order implicit scheme
204 """
206 alpha = [0.0, -1.0]
207 beta = [-1.0 / 12.0, 8.0 / 12.0, 5.0 / 12.0]
209 def generate_starting_values(self):
210 """
211 Generate starting value by trapezoidal rule.
212 """
213 lvl = self.level
214 prob = lvl.prob
215 time = lvl.time + lvl.dt
217 # do trapezoidal rule
218 rhs = lvl.u[0] + lvl.dt / 2 * lvl.f[0]
220 lvl.u[1] = prob.solve_system(rhs, lvl.dt / 2.0, lvl.u[0], time)