Coverage for pySDC/implementations/convergence_controller_classes/estimate_polynomial_error.py: 90%
97 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
3from qmat.lagrange import LagrangeApproximation
4from pySDC.core.convergence_controller import ConvergenceController
7class EstimatePolynomialError(ConvergenceController):
8 """
9 Estimate the local error by using all but one collocation node in a polynomial interpolation to that node.
10 While the converged collocation problem with M nodes gives a order M approximation to this point, the interpolation
11 gives only an order M-1 approximation. Hence, we have two solutions with different order, and we know their order.
12 That is to say this gives an error estimate that is order M. Keep in mind that the collocation problem should be
13 converged for this and has order up to 2M. Still, the lower order method can be used for time step selection, for
14 instance.
15 If the last node is not the end point, we can interpolate to that node, which is an order M approximation and compare
16 to the order 2M approximation we get from the extrapolation step.
17 By default, we interpolate to the second to last node, or to the end point for Gauss nodes, where it is not a node.
18 """
20 def setup(self, controller, params, description, **kwargs):
21 """
22 Args:
23 controller (pySDC.Controller.controller): The controller
24 params (dict): The params passed for this specific convergence controller
25 description (dict): The description object used to instantiate the controller
27 Returns:
28 (dict): The updated params dictionary
29 """
30 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
31 from pySDC.core.check_convergence import CheckConvergence
33 sweeper_params = description['sweeper_params']
34 num_nodes = sweeper_params['num_nodes']
35 quad_type = sweeper_params['quad_type']
37 defaults = {
38 'control_order': -75,
39 'estimate_on_node': num_nodes + 1 if quad_type == 'GAUSS' else num_nodes - 1,
40 'rel_error': False,
41 **super().setup(controller, params, description, **kwargs),
42 }
43 self.comm = description['sweeper_params'].get('comm', None)
45 if self.comm:
46 from mpi4py import MPI
48 self.prepare_MPI_datatypes()
49 self.MPI_SUM = MPI.SUM
51 controller.add_hook(LogEmbeddedErrorEstimate)
52 self.check_convergence = CheckConvergence.check_convergence
54 if quad_type != 'GAUSS' and defaults['estimate_on_node'] > num_nodes:
55 from pySDC.core.errors import ParameterError
57 raise ParameterError(
58 'You cannot interpolate with lower accuracy to the end point if the end point is a node!'
59 )
61 self.interpolation_matrix = None
63 return defaults
65 def reset_status_variables(self, *args, **kwargs):
66 """
67 Add variable for embedded error
69 Returns:
70 None
71 """
72 self.add_status_variable_to_level('error_embedded_estimate')
73 self.add_status_variable_to_level('order_embedded_estimate')
75 def matmul(self, A, b, xp=np):
76 """
77 Matrix vector multiplication, possibly MPI parallel.
78 The parallel implementation performs a reduce operation in every row of the matrix. While communicating the
79 entire vector once could reduce the number of communications, this way we never need to store the entire vector
80 on any specific rank.
82 Args:
83 A (2d np.ndarray): Matrix
84 b (list): Vector
85 xp: Either numpy or cupy
87 Returns:
88 List: Axb
89 """
91 if self.comm:
92 res = [A[i, 0] * b[0] if b[i] is not None else None for i in range(A.shape[0])]
93 buf = b[0] * 0.0
94 for i in range(0, A.shape[0]):
95 index = self.comm.rank + (1 if self.comm.rank < self.params.estimate_on_node - 1 else 0)
96 # `zeros_like` gives a bare array on CuPy, where NumPy would keep the datatype:
97 # CuPy has no `subok`. That is fine here, since this is only a send buffer and
98 # needs a device pointer rather than a datatype. Do not "tidy" it into
99 # `res[0] * 0.0`, which would keep the type but propagate `NaN` -- and `res` is a
100 # residual, so it is `NaN` exactly when a solve has already gone wrong.
101 send_buf = (
102 (A[i, index] * b[index])
103 if self.comm.rank != self.params.estimate_on_node - 1
104 else xp.zeros_like(res[0])
105 )
106 self.comm.Allreduce(send_buf, buf, op=self.MPI_SUM)
107 res[i] += buf
108 return res
109 else:
110 # `asarray` stacks the vector into a single array, and CuPy builds that through
111 # `copy`, which asserts -- with no message -- when it is handed back an `ndarray`
112 # subclass instead of the base type. `cupy_mesh.copy` returns `type(self)` on
113 # purpose, so that a copied mesh stays a mesh, which makes every datatype here such a
114 # subclass. NumPy has no such restriction, but dropping the wrapper is right either
115 # way: what this wants is the numbers, not the datatype.
116 return A @ xp.asarray([entry.view(xp.ndarray) for entry in b])
118 def get_interpolated_solution(self, L, xp):
119 """
120 Get the interpolated solution for numpy or cupy data types
122 Args:
123 L (pySDC.level): The level holding the solutions at the collocation nodes
124 xp: The array module of the data, numpy or cupy
125 """
126 coll = L.sweep.coll
128 u = [
129 L.u[i].flatten() if L.u[i] is not None else L.u[i]
130 for i in range(coll.num_nodes + 1)
131 if i != self.params.estimate_on_node
132 ]
133 return self.matmul(self.interpolation_matrix, u, xp=xp)[0].reshape(L.prob.init[0])
135 def post_iteration_processing(self, controller, S, **kwargs):
136 """
137 Estimate the error
139 Args:
140 controller (pySDC.Controller.controller): The controller
141 S (pySDC.Step.step): The current step
143 Returns:
144 None
145 """
147 if self.check_convergence(S):
148 L = S.levels[0]
149 coll = L.sweep.coll
150 nodes = np.append(np.append(0, coll.nodes), 1.0)
151 estimate_on_node = self.params.estimate_on_node
153 if hasattr(L.u[0], 'xp'):
154 xp = L.u[0].xp
155 else:
156 xp = np
158 if self.interpolation_matrix is None:
159 interpolator = LagrangeApproximation(
160 points=[nodes[i] for i in range(coll.num_nodes + 1) if i != estimate_on_node]
161 )
162 self.interpolation_matrix = xp.array(interpolator.getInterpolationMatrix([nodes[estimate_on_node]]))
164 u_inter = self.get_interpolated_solution(L, xp)
166 # compute end point if needed
167 if estimate_on_node == len(nodes) - 1:
168 if L.uend is None:
169 L.sweep.compute_end_point()
170 high_order_sol = L.uend
171 rank = 0
172 L.status.order_embedded_estimate = coll.num_nodes + 1
173 else:
174 high_order_sol = L.u[estimate_on_node]
175 rank = estimate_on_node - 1
176 L.status.order_embedded_estimate = coll.num_nodes * 1
178 if not self.comm or self.comm.rank == rank:
179 error = abs(u_inter - high_order_sol)
180 if self.params.rel_error:
181 # the norm of the datatype, as for the error: `u_inter` may be a bare array
182 error /= abs(high_order_sol)
184 if self.comm:
185 buf = np.array(error if self.comm.rank == rank else 0.0)
186 self.comm.Bcast(buf, root=rank)
187 L.status.error_embedded_estimate = float(buf)
188 else:
189 L.status.error_embedded_estimate = error
191 self.debug(
192 f'Obtained error estimate: {L.status.error_embedded_estimate:.2e} of order {L.status.order_embedded_estimate}',
193 S,
194 )
196 def check_parameters(self, controller, params, description, **kwargs):
197 """
198 Check if we allow the scheme to solve the collocation problems to convergence.
200 Args:
201 controller (pySDC.Controller): The controller
202 params (dict): The params passed for this specific convergence controller
203 description (dict): The description object used to instantiate the controller
205 Returns:
206 bool: Whether the parameters are compatible
207 str: The error message
208 """
209 if description['sweeper_params'].get('num_nodes', 0) < 2:
210 return False, 'Need at least two collocation nodes to interpolate to one!'
212 return True, ""
215class EstimatePolynomialErrorFiredrake(EstimatePolynomialError):
216 """
217 Polynomial interpolation error estimate for Firedrake functions, combining the node solutions term by term.
218 """
220 def matmul(self, A, b):
221 """
222 Matrix vector multiplication, possibly MPI parallel.
223 The parallel implementation performs a reduce operation in every row of the matrix. While communicating the
224 entire vector once could reduce the number of communications, this way we never need to store the entire vector
225 on any specific rank.
227 Args:
228 A (2d np.ndarray): Matrix
229 b (list): Vector
231 Returns:
232 List: Axb
233 """
235 if self.comm:
236 res = [A[i, 0] * b[0] if b[i] is not None else None for i in range(A.shape[0])]
237 buf = 0 * b[0]
238 for i in range(0, A.shape[0]):
239 index = self.comm.rank + (1 if self.comm.rank < self.params.estimate_on_node - 1 else 0)
240 send_buf = (
241 (A[i, index] * b[index]) if self.comm.rank != self.params.estimate_on_node - 1 else 0 * res[0]
242 )
243 self.comm.Allreduce(send_buf, buf, op=self.MPI_SUM)
244 res[i] += buf
245 return res
246 else:
247 res = []
248 for i in range(A.shape[0]):
249 res.append(A[i, 0] * b[0])
250 for j in range(1, A.shape[1]):
251 res[-1] += A[i, j] * b[j]
253 return res
255 def get_interpolated_solution(self, L, xp=None):
256 """
257 Get the interpolated solution for Firedrake data types. The nodes hold solutions, not mass matrix applied to
258 them, so interpolating them needs no mass matrix.
260 Args:
261 L (pySDC.level): The level holding the solutions at the collocation nodes
262 xp: Not used, Firedrake data is not an array; accepted because the base class passes it
263 """
264 coll = L.sweep.coll
266 u = [
267 L.u[i] if L.u[i] is not None else L.u[i]
268 for i in range(coll.num_nodes + 1)
269 if i != self.params.estimate_on_node
270 ]
271 return L.prob.dtype_u(self.matmul(self.interpolation_matrix, u)[0])