Coverage for pySDC/implementations/convergence_controller_classes/adaptive_collocation.py: 99%
77 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-04 15:08 +0000
1import numpy as np
3from qmat.lagrange import LagrangeApproximation
4from pySDC.core.convergence_controller import ConvergenceController, Status
7class AdaptiveCollocation(ConvergenceController):
8 """
9 This convergence controller allows to change the underlying quadrature between iterations.
11 Supplying multiple quadrature rules will result in a change to a new quadrature whenever the previous one is
12 converged until all methods given are converged, at which point the step is ended.
13 Whenever the quadrature is changed, the solution is interpolated from the old nodes to the new nodes to ensure
14 accelerated convergence compared to starting from the initial conditions.
16 Use this convergence controller by supplying parameters that the sweeper accepts as a list to the `params`.
17 For instance, supplying
19 .. code-block:: python
21 params = {
22 'num_nodes': [2, 3],
23 }
25 will use collocation methods like you passed to the `sweeper_params` in the `description` object, but will change
26 the number of nodes to 2 before the first iteration and to 3 as soon as the 2-node collocation problem is converged.
27 This will override whatever you set for the number of nodes in the `sweeper_params`, but you still need to set
28 something to allow instantiation of the levels before this convergence controller becomes active.
29 Make sure all lists you supply here have the same length.
31 Feel free to set `logger_level = 15` in the controller parameters to get comprehensive text output on what exactly
32 is happening.
34 This convergence controller has various applications.
35 - You could try to obtain speedup by doing some inexactness. It is currently possible to set various residual
36 tolerances, which will be passed to the levels, corresponding to the accuracy with which each collocation problem
37 is solved.
38 - You can compute multiple solutions to the same initial value problem with different order. This allows, for
39 instance, to do adaptive time stepping.
41 When trying to obtain speedup with this, be ware that the interpolation is not for free. In particular, it is
42 necessary to reevaluate the right hand side at all nodes afterwards.
43 """
45 def setup(self, controller, params, description, **kwargs):
46 """
47 Record what variables we want to vary.
49 Args:
50 controller (pySDC.Controller.controller): The controller
51 params (dict): The params passed for this specific convergence controller
52 description (dict): The description object used to instantiate the controller
54 Returns:
55 (dict): The updated params dictionary
56 """
58 defaults = {
59 'control_order': 300,
60 'num_colls': 0,
61 'sweeper_params': description['sweeper_params'],
62 'vary_keys_sweeper': [],
63 'vary_keys_level': [],
64 }
66 # only these keys can be changed by this convergence controller
67 self.allowed_sweeper_keys = ['quad_type', 'num_nodes', 'node_type', 'do_coll_update']
68 self.allowed_level_keys = ['restol']
70 # add the keys to lists so we know what we need to change later
71 for key in params.keys():
72 if type(params[key]) == list:
73 if key in self.allowed_sweeper_keys:
74 defaults['vary_keys_sweeper'] += [key]
75 elif key in self.allowed_level_keys:
76 defaults['vary_keys_level'] += [key]
77 else:
78 raise NotImplementedError(f'Don\'t know what to do with key {key} here!')
80 defaults['num_colls'] = max([defaults['num_colls'], len(params[key])])
82 self.comm = description['sweeper_params'].get('comm', None)
83 if self.comm:
84 from mpi4py import MPI
86 self.MPI_SUM = MPI.SUM
88 return {**defaults, **super().setup(controller, params, description, **kwargs)}
90 def matmul(self, A, b):
91 """
92 Matrix vector multiplication, possibly MPI parallel.
93 The parallel implementation performs a reduce operation in every row of the matrix. While communicating the
94 entire vector once could reduce the number of communications, this way we never need to store the entire vector
95 on any specific rank.
97 Args:
98 A (2d np.ndarray): Matrix
99 b (list): Vector
101 Returns:
102 List: Axb
103 """
104 if self.comm:
105 res = [A[i, 0] * b[0] if b[i] is not None else None for i in range(A.shape[0])]
106 buf = b[0] * 0.0
107 for i in range(1, A.shape[1]):
108 self.comm.Reduce(A[i, self.comm.rank + 1] * b[self.comm.rank + 1], buf, op=self.MPI_SUM, root=i - 1)
109 if i == self.comm.rank + 1:
110 res[i] += buf
111 return res
112 else:
113 return A @ b
115 def switch_sweeper(self, S):
116 """
117 Update to the next sweeper in line.
119 Args:
120 S (pySDC.Step.step): The current step
122 Returns:
123 None
124 """
126 # generate dictionaries with the new parameters
127 new_params_sweeper = {
128 key: self.params.get(key)[self.status.active_coll] for key in self.params.vary_keys_sweeper
129 }
130 sweeper_params = self.params.sweeper_params.copy()
131 update_params_sweeper = {**sweeper_params, **new_params_sweeper}
133 new_params_level = {key: self.params.get(key)[self.status.active_coll] for key in self.params.vary_keys_level}
135 # update sweeper for all levels
136 for L in S.levels:
137 P = L.prob
139 # store solution of current level which will be interpolated to new level
140 u_old = [me.flatten() if me is not None else me for me in L.u]
141 nodes_old = L.sweep.coll.nodes.copy()
143 # change sweeper
144 L.sweep.__init__(update_params_sweeper, L)
146 # reset level to tell it the new structure of the solution
147 L.params.__dict__.update(new_params_level)
148 L.reset_level(reset_status=False)
150 # interpolate solution of old collocation problem to new one
151 nodes_new = L.sweep.coll.nodes.copy()
152 interpolator = LagrangeApproximation(points=np.append(0, nodes_old))
154 u_inter = self.matmul(interpolator.getInterpolationMatrix(np.append(0, nodes_new)), u_old)
156 # assign the interpolated values to the nodes in the level
157 for i in range(0, len(u_inter)):
158 if u_inter[i] is not None:
159 me = P.dtype_u(P.init)
160 me[:] = np.reshape(u_inter[i], P.init[0])
161 L.u[i] = me
163 # reevaluate rhs
164 for i in range(L.sweep.coll.num_nodes + 1):
165 if L.u[i] is not None:
166 L.f[i] = L.prob.eval_f(L.u[i], L.time)
168 # log the new parameters
169 self.log(f'Switching to collocation {self.status.active_coll + 1} of {self.params.num_colls}', S, level=20)
170 msg = 'New quadrature:'
171 for key in list(sweeper_params.keys()) + list(new_params_level.keys()):
172 if key in self.params.vary_keys_sweeper:
173 msg += f'\n--> {key}: {update_params_sweeper[key]}'
174 elif key in self.params.vary_keys_level:
175 msg += f'\n--> {key}: {new_params_level[key]}'
176 else:
177 msg += f'\n {key}: {update_params_sweeper[key]}'
178 self.log(msg, S)
180 def setup_status_variables(self, controller, **kwargs):
181 """
182 Add an index for which collocation method to use.
184 Args:
185 controller (pySDC.Controller.controller): The controller
187 Returns:
188 None
189 """
190 self.status = Status(['active_coll'])
192 def reset_status_variables(self, controller, **kwargs):
193 """
194 Reset the status variables between time steps.
196 Args:
197 controller (pySDC.Controller.controller): The controller
199 Returns:
200 None
201 """
202 self.status.active_coll = 0
204 def post_iteration_processing(self, controller, S, **kwargs):
205 """
206 Switch to the next collocation method if the current one is converged.
208 Args:
209 controller (pySDC.Controller.controller): The controller
210 S (pySDC.Step.step): The current step
212 Returns:
213 None
214 """
215 if (self.status.active_coll < self.params.num_colls - 1) and S.status.done:
216 self.status.active_coll += 1
217 S.status.done = False
218 self.switch_sweeper(S)
220 def post_spread_processing(self, controller, S, **kwargs):
221 """
222 Overwrite the sweeper parameters with the first collocation parameters set up here.
224 Args:
225 controller (pySDC.Controller.controller): The controller
226 S (pySDC.Step.step): The current step
228 Returns:
229 None
230 """
231 self.switch_sweeper(S)
233 def check_parameters(self, controller, params, description, **kwargs):
234 """
235 Check if we allow the scheme to solve the collocation problems to convergence.
237 Args:
238 controller (pySDC.Controller): The controller
239 params (dict): The params passed for this specific convergence controller
240 description (dict): The description object used to instantiate the controller
242 Returns:
243 bool: Whether the parameters are compatible
244 str: The error message
245 """
246 if description["level_params"].get("restol", -1.0) <= 1e-16:
247 return (
248 False,
249 "Switching the collocation problems requires solving them to some tolerance that can be reached. Please set attainable `restol` in the level params",
250 )
252 return True, ""