+ Source code for implementations.convergence_controller_classes.estimate_polynomial_error
import numpy as np
from pySDC.core.Lagrange import LagrangeApproximation
@@ -44,9 +44,9 @@ Source code for implementations.convergence_controller_classes.estimate_inte
from pySDC.core.Collocation import CollBase
-
-
[docs]
-
class EstimateInterpolationError(ConvergenceController):
+
+
[docs]
+
class EstimatePolynomialError(ConvergenceController):
"""
Estimate the local error by using all but one collocation node in a polynomial interpolation to that node.
While the converged collocation problem with M nodes gives a order M approximation to this point, the interpolation
@@ -54,11 +54,13 @@
Source code for implementations.convergence_controller_classes.estimate_inte
That is to say this gives an error estimate that is order M. Keep in mind that the collocation problem should be
converged for this and has order up to 2M. Still, the lower order method can be used for time step selection, for
instance.
+ If the last node is not the end point, we can interpolate to that node, which is an order M approximation and compare
+ to the order 2M approximation we get from the extrapolation step.
By default, we interpolate to the second to last node.
"""
-
-
[docs]
+
+
[docs]
def setup(self, controller, params, description, **kwargs):
"""
Args:
@@ -72,9 +74,14 @@
Source code for implementations.convergence_controller_classes.estimate_inte
from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence
+ sweeper_params = description['sweeper_params']
+ num_nodes = sweeper_params['num_nodes']
+ quad_type = sweeper_params['quad_type']
+
defaults = {
'control_order': -75,
- 'estimate_on_node': description['sweeper_params'].get('num_nodes', 2) - 1,
+ 'estimate_on_node': num_nodes + 1 if quad_type == 'GAUSS' else num_nodes - 1,
+ **super().setup(controller, params, description, **kwargs),
}
self.comm = description['sweeper_params'].get('comm', None)
@@ -87,11 +94,18 @@ Source code for implementations.convergence_controller_classes.estimate_inte
controller.add_hook(LogEmbeddedErrorEstimate)
self.check_convergence = CheckConvergence.check_convergence
- return {**defaults, **super().setup(controller, params, description, **kwargs)}
+
if quad_type != 'GAUSS' and defaults['estimate_on_node'] > num_nodes:
+
from pySDC.core.Errors import ParameterError
+
+
raise ParameterError(
+
'You cannot interpolate with lower accuracy to the end point if the end point is a node!'
+
)
+
+
return defaults
-
-
[docs]
+
+
[docs]
def reset_status_variables(self, controller, **kwargs):
"""
Add variable for embedded error
@@ -112,11 +126,12 @@
Source code for implementations.convergence_controller_classes.estimate_inte
where = ["levels", "status"]
for S in steps:
- self.add_variable(S, name='error_embedded_estimate', where=where, init=None)
+
self.add_variable(S, name='error_embedded_estimate', where=where, init=None)
+
self.add_variable(S, name='order_embedded_estimate', where=where, init=None)
-
-
[docs]
+
+
[docs]
def matmul(self, A, b):
"""
Matrix vector multiplication, possibly MPI parallel.
@@ -148,8 +163,8 @@
Source code for implementations.convergence_controller_classes.estimate_inte
return A @ np.asarray(b)
-
-
[docs]
+
+
[docs]
def post_iteration_processing(self, controller, S, **kwargs):
"""
Estimate the error
@@ -165,13 +180,13 @@
Source code for implementations.convergence_controller_classes.estimate_inte
if self.check_convergence(S):
L = S.levels[0]
coll = L.sweep.coll
- nodes = coll.nodes
+ nodes = np.append(np.append(0, coll.nodes), 1.0)
estimate_on_node = self.params.estimate_on_node
interpolator = LagrangeApproximation(
- points=[0] + [nodes[i - 1] for i in range(1, coll.num_nodes + 1) if i != estimate_on_node]
+ points=[nodes[i] for i in range(coll.num_nodes + 1) if i != estimate_on_node]
)
- interpolation_matrix = interpolator.getInterpolationMatrix([nodes[estimate_on_node - 1]])
+ interpolation_matrix = interpolator.getInterpolationMatrix([nodes[estimate_on_node]])
u = [
L.u[i].flatten() if L.u[i] is not None else L.u[i]
for i in range(coll.num_nodes + 1)
@@ -179,16 +194,28 @@ Source code for implementations.convergence_controller_classes.estimate_inte
]
u_inter = self.matmul(interpolation_matrix, u)[0].reshape(L.prob.init[0])
+ # compute end point if needed
+ if estimate_on_node == len(nodes) - 1:
+ if L.uend is None:
+ L.sweep.compute_end_point()
+ high_order_sol = L.uend
+ rank = 0
+ L.status.order_embedded_estimate = coll.num_nodes + 1
+ else:
+ high_order_sol = L.u[estimate_on_node]
+ rank = estimate_on_node - 1
+ L.status.order_embedded_estimate = coll.num_nodes * 1
+
if self.comm:
- buf = np.array(abs(u_inter - L.u[estimate_on_node]) if self.comm.rank == estimate_on_node - 1 else 0.0)
- self.comm.Bcast(buf, root=estimate_on_node - 1)
+ buf = np.array(abs(u_inter - high_order_sol) if self.comm.rank == rank else 0.0)
+ self.comm.Bcast(buf, root=rank)
L.status.error_embedded_estimate = buf
else:
- L.status.error_embedded_estimate = abs(u_inter - L.u[estimate_on_node])
+
L.status.error_embedded_estimate = abs(u_inter - high_order_sol)
-
-
[docs]
+
diff --git a/coverage/d_020efe120a771d8a_hamiltonian_and_energy_output_py.html b/coverage/d_020efe120a771d8a_hamiltonian_and_energy_output_py.html
index daa86f2e47..e44c1b78fe 100644
--- a/coverage/d_020efe120a771d8a_hamiltonian_and_energy_output_py.html
+++ b/coverage/d_020efe120a771d8a_hamiltonian_and_energy_output_py.html
@@ -65,7 +65,7 @@
» next
coverage.py v7.3.1,
- created at 2023-09-22 06:33 +0000
+ created at 2023-09-22 06:42 +0000
diff --git a/coverage/d_020efe120a771d8a_hamiltonian_output_py.html b/coverage/d_020efe120a771d8a_hamiltonian_output_py.html
index 837662567c..02fcbd74db 100644
--- a/coverage/d_020efe120a771d8a_hamiltonian_output_py.html
+++ b/coverage/d_020efe120a771d8a_hamiltonian_output_py.html
@@ -65,7 +65,7 @@
» next
coverage.py v7.3.1,
- created at 2023-09-22 06:33 +0000
+ created at 2023-09-22 06:42 +0000
diff --git a/coverage/d_020efe120a771d8a_harmonic_oscillator_py.html b/coverage/d_020efe120a771d8a_harmonic_oscillator_py.html
index a8aa390d5d..aaa9ce494e 100644
--- a/coverage/d_020efe120a771d8a_harmonic_oscillator_py.html
+++ b/coverage/d_020efe120a771d8a_harmonic_oscillator_py.html
@@ -65,7 +65,7 @@
» next
coverage.py v7.3.1,
- created at 2023-09-22 06:33 +0000
+ created at 2023-09-22 06:42 +0000
diff --git a/coverage/d_020efe120a771d8a_simple_problems_py.html b/coverage/d_020efe120a771d8a_simple_problems_py.html
index b4b7fc58a0..e06899bf37 100644
--- a/coverage/d_020efe120a771d8a_simple_problems_py.html
+++ b/coverage/d_020efe120a771d8a_simple_problems_py.html
@@ -65,7 +65,7 @@
» next
coverage.py v7.3.1,
- created at 2023-09-22 06:33 +0000
+ created at 2023-09-22 06:42 +0000
diff --git a/coverage/d_020efe120a771d8a_solar_system_py.html b/coverage/d_020efe120a771d8a_solar_system_py.html
index b28ccc32a4..fec5dde827 100644
--- a/coverage/d_020efe120a771d8a_solar_system_py.html
+++ b/coverage/d_020efe120a771d8a_solar_system_py.html
@@ -65,7 +65,7 @@
» next
coverage.py v7.3.1,
- created at 2023-09-22 06:33 +0000
+ created at 2023-09-22 06:42 +0000
diff --git a/coverage/d_020efe120a771d8a_stop_at_error_hook_py.html b/coverage/d_020efe120a771d8a_stop_at_error_hook_py.html
index 2db5138a96..011cb93aa0 100644
--- a/coverage/d_020efe120a771d8a_stop_at_error_hook_py.html
+++ b/coverage/d_020efe120a771d8a_stop_at_error_hook_py.html
@@ -65,7 +65,7 @@
» next
coverage.py v7.3.1,
- created at 2023-09-22 06:33 +0000
+ created at 2023-09-22 06:42 +0000
diff --git a/coverage/d_064a9f2a35945611_FaultHooks_py.html b/coverage/d_064a9f2a35945611_FaultHooks_py.html
index c5a3745954..c45e119a35 100644
--- a/coverage/d_064a9f2a35945611_FaultHooks_py.html
+++ b/coverage/d_064a9f2a35945611_FaultHooks_py.html
@@ -65,7 +65,7 @@