From 92cf3855ff85e67b0b9a6b5636d796f1a41431a4 Mon Sep 17 00:00:00 2001 From: lisawim Date: Wed, 25 Oct 2023 12:55:25 +0200 Subject: [PATCH 01/29] Added piecewise linear interpolation to SwitchEstimator --- pySDC/projects/PinTSimE/switch_estimator.py | 41 +++++++-------------- 1 file changed, 14 insertions(+), 27 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 3c438c134a..a075189eb9 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -274,7 +274,7 @@ def get_switch(t_interp, state_function, m_guess): Time point of found event. """ - LagrangeInterpolator = LagrangeInterpolation(t_interp, state_function) + LinearInterpolator = LinearInterpolation(t_interp, state_function) def p(t): """ @@ -290,7 +290,7 @@ def p(t): p(t) : float The value of the interpolated function at time t. """ - return LagrangeInterpolator.eval(t) + return LinearInterpolator.eval(t) def fprime(t): """ @@ -387,45 +387,32 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): return root -class LagrangeInterpolation(object): +class LinearInterpolation(object): def __init__(self, ti, yi): """Initialization routine""" self.ti = np.asarray(ti) self.yi = np.asarray(yi) self.n = len(ti) - def get_Lagrange_polynomial(self, t, i): - """ - Computes the basis of the i-th Lagrange polynomial. - - Parameters - ---------- - t : float - Time where the polynomial is computed at. - i : int - Index of the Lagrange polynomial - - Returns - ------- - product : float - The product of the bases. - """ - product = np.prod([(t - self.ti[k]) / (self.ti[i] - self.ti[k]) for k in range(self.n) if k != i]) - return product - def eval(self, t): - """ - Evaluates the Lagrange interpolation at time t. + r""" + Evaluates the linear interpolant at time :math:`t`. Parameters ---------- t : float - Time where interpolation is computed. + Time where the interpolant is evaluated. Returns ------- p : float - Value of interpolant at time t. + Value of the interpolant at time :math:`t`. """ - p = np.sum([self.yi[i] * self.get_Lagrange_polynomial(t, i) for i in range(self.n)]) + + for i in range(1, self.n): + if t >= self.ti[i - 1] or t <= self.ti[i]: + p = self.yi[i - 1] + (self.yi[i] - self.yi[i - 1]) / (self.ti[i] - self.ti[i - 1]) * (t - self.ti[i - 1]) + else: + p = 0 + raise ParameterError('Interpolant is only evaluated for time occuring between time nodes') return p From 4040463f027a86d880e47d256b26e511e0a4bd89 Mon Sep 17 00:00:00 2001 From: lisawim Date: Wed, 25 Oct 2023 12:56:36 +0200 Subject: [PATCH 02/29] Started with test for SwitchEstimator [WIP] --- .../test_pintsime/test_SwitchEstimator.py | 191 ++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100644 pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py new file mode 100644 index 0000000000..c68c22a936 --- /dev/null +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -0,0 +1,191 @@ +import numpy as np +import pytest + +from pySDC.implementations.problem_classes.DiscontinuousTestODE import DiscontinuousTestODE +from pySDC.core.Problem import ptype +from pySDC.implementations.datatype_classes.mesh import mesh + + +class ExactDiscontinuousTestODE(DiscontinuousTestODE): + r""" + Dummy problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestODE``. + """ + + def eval_f(self, u, t): + """ + Derivative. + + Parameters + ---------- + u : dtype_u + Exact value of u. + t : _type_ + Time :math:`t`. + + Returns + ------- + f : dtype_f + Derivative. + """ + + f = self.dtype_f(self.init) + + t_switch = np.inf if self.t_switch is None else self.t_switch + h = u[0] - 5 + if h >= 0 or t >= t_switch: + f[:] = 1 + else: + f[:] = np.exp(t) + return f + + def solve_system(self, rhs, factor, u0, t): + """ + Just return the exact solution... + + Parameters + ---------- + rhs : dtype_f + Right-hand side for the linear system. + factor : float + Abbrev. for the local stepsize (or any other factor required). + u0 : dtype_u + Initial guess for the iterative solver. + t : float + Current time (e.g. for time-dependent BCs). + + Returns + ------- + me : dtype_u + The solution as mesh. + """ + + return self.u_exact(t) + + +def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Tend, dt, tol): + """ + Simulates one run of a discontinuous test problem class for one specific tolerance specified in + the parameters. For testing, only one time step should be considered. + + Parameters + ---------- + problem : pySDC.core.Problem + Problem class to be simulated. + sweeper: pySDC.core.Sweeper + Sweeper used for the simulation. + quad_type : str + Type of quadrature used. + num_nodes : int + Number of collocation nodes. + t0 : float + Starting time. + Tend : float + End time. + dt : float + Time step size. + tol : float + Tolerance for the switch estimator. + """ + + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + + # initialize level parameters + level_params = { + 'restol': 1e-13, + 'dt': dt, + } + + # initialize sweeper parameters + sweeper_params = { + 'quad_type': quad_type, + 'num_nodes': num_nodes, + 'QI': 'IE', + 'initial_guess': 'spread', + } + + # initialize step parameters + step_params = { + 'maxiter': 10, + } + + # initialize controller parameters + controller_params = { + 'logger_level': 30, + 'hook_class': LogSolution, + } + + # convergence controllers + convergence_controllers = dict() + switch_estimator_params = { + 'tol': tol, + 'alpha': 1.0, + } + convergence_controllers.update({SwitchEstimator: switch_estimator_params}) + + convergence_controllers[BasicRestartingNonMPI] = { + 'max_restarts': 3, + 'crash_after_max_restarts': False, + } + + # fill description dictionary for easy step instantiation + description = { + 'problem_class': problem, # pass problem class + 'problem_params': dict(), # problem_params, # pass problem parameters + 'sweeper_class': sweeper, # pass sweeper + 'sweeper_params': sweeper_params, # pass sweeper parameters + 'level_params': level_params, # pass level parameters + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } + + # set time parameters + t0 = t0 + Tend = Tend + + # instantiate controller + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + # get initial values on finest level + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + t_switch_exact = P.t_switch_exact + + # call main function to get things done... + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + + return stats, t_switch_exact + + +@pytest.mark.base +@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) +@pytest.mark.parametrize('num_nodes', [3, 4, 5]) +def test_all_tolerances_ODE(tol, num_nodes): + """Tests a problem for different tolerances for the switch estimator.""" + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.helpers.stats_helper import get_sorted + + t0 = 1.6 + Tend = 1.62 + problem = ExactDiscontinuousTestODE + + stats, t_switch_exact = discontinuousTestProblem_run( + problem=problem, + sweeper=generic_implicit, + quad_type='LOBATTO', + num_nodes=num_nodes, + t0=t0, + Tend=Tend, + dt=2e-2, + tol=tol, + ) + + # in this specific example only one event has to be found + switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] + assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' + + t_switch = switches[-1] + event_err = abs(t_switch - t_switch_exact) + assert event_err < 1e-14, f'Event time error for is not small enough!' From 2a797177b7582b29f089f9b116403771c083046a Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 26 Oct 2023 10:26:21 +0200 Subject: [PATCH 03/29] Test to proof sum_restarts when event occuring at boundary --- .../problem_classes/DiscontinuousTestODE.py | 2 +- pySDC/projects/PinTSimE/switch_estimator.py | 23 +- .../test_pintsime/test_SwitchEstimator.py | 211 +++++++++++++++--- 3 files changed, 197 insertions(+), 39 deletions(-) diff --git a/pySDC/implementations/problem_classes/DiscontinuousTestODE.py b/pySDC/implementations/problem_classes/DiscontinuousTestODE.py index 1179cb07c6..a0edd0db6f 100644 --- a/pySDC/implementations/problem_classes/DiscontinuousTestODE.py +++ b/pySDC/implementations/problem_classes/DiscontinuousTestODE.py @@ -216,7 +216,7 @@ def get_switching_info(self, u, t): m_guess = m - 1 break - state_function = [u[m][0] - 5 for m in range(len(u))] if switch_detected else [] + state_function = [u[m][0] - 5 for m in range(len(u))] return switch_detected, m_guess, state_function def count_switches(self): diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index a075189eb9..33c9749c06 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -1,6 +1,7 @@ import numpy as np import scipy as sp +from pySDC.core.Errors import ParameterError from pySDC.core.Collocation import CollBase from pySDC.core.ConvergenceController import ConvergenceController, Status from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence @@ -46,6 +47,7 @@ def setup(self, controller, params, description): 'control_order': 0, 'nodes': coll.nodes, 'tol_zero': 1e-13, + 't_interp': [], } return {**defaults, **params} @@ -86,27 +88,27 @@ def get_new_step_size(self, controller, S, **kwargs): """ L = S.levels[0] - + print('conv', CheckConvergence.check_convergence(S)) if CheckConvergence.check_convergence(S): self.status.switch_detected, m_guess, state_function = L.prob.get_switching_info(L.u, L.time) - + print('detected:', self.status.switch_detected, state_function) if self.status.switch_detected: - t_interp = [L.time + L.dt * self.params.nodes[m] for m in range(len(self.params.nodes))] - t_interp, state_function = self.adapt_interpolation_info( - L.time, L.sweep.coll.left_is_node, t_interp, state_function + self.params.t_interp = [L.time + L.dt * self.params.nodes[m] for m in range(len(self.params.nodes))] + self.params.t_interp, state_function = self.adapt_interpolation_info( + L.time, L.sweep.coll.left_is_node, self.params.t_interp, state_function ) - + # print(t_interp, state_function) # when the state function is already close to zero the event is already resolved well if abs(state_function[-1]) <= self.params.tol_zero or abs(state_function[0]) <= self.params.tol_zero: if abs(state_function[0]) <= self.params.tol_zero: - t_switch = t_interp[0] + t_switch = self.params.t_interp[0] boundary = 'left' elif abs(state_function[-1]) <= self.params.tol_zero: boundary = 'right' - t_switch = t_interp[-1] + t_switch = self.params.t_interp[-1] msg = f"The value of state function is close to zero, thus event time is already close enough to the {boundary} end point!" - self.log(msg, S) + # self.log(msg, S) self.log_event_time( controller.hooks[0], S.status.slot, L.time, L.level_index, L.status.sweep, t_switch ) @@ -116,7 +118,7 @@ def get_new_step_size(self, controller, S, **kwargs): # intermediate value theorem states that a root is contained in current step if state_function[0] * state_function[-1] < 0 and self.status.is_zero is None: - self.status.t_switch = self.get_switch(t_interp, state_function, m_guess) + self.status.t_switch = self.get_switch(self.params.t_interp, state_function, m_guess) controller.hooks[0].add_to_stats( process=S.status.slot, @@ -178,6 +180,7 @@ def get_new_step_size(self, controller, S, **kwargs): L.status.sweep, self.status.t_switch, ) + print('Occurs at boundary') L.prob.count_switches() self.status.switch_detected = False diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index c68c22a936..02620ebff6 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -11,6 +11,14 @@ class ExactDiscontinuousTestODE(DiscontinuousTestODE): Dummy problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestODE``. """ + def __init__(self, newton_maxiter=100, newton_tol=1e-8): + """Initialization routine""" + super().__init__(newton_maxiter, newton_tol) + + self.t_switch_exact = np.log(5) + self.t_switch = None + self.nswitches = 0 + def eval_f(self, u, t): """ Derivative. @@ -88,10 +96,12 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten """ from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI - from pySDC.implementations.hooks.log_solution import LogSolution from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.implementations.hooks.log_restarts import LogRestarts + # initialize level parameters level_params = { 'restol': 1e-13, @@ -114,7 +124,7 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten # initialize controller parameters controller_params = { 'logger_level': 30, - 'hook_class': LogSolution, + 'hook_class': [LogSolution, LogRestarts], } # convergence controllers @@ -160,32 +170,177 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten @pytest.mark.base -@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) -@pytest.mark.parametrize('num_nodes', [3, 4, 5]) -def test_all_tolerances_ODE(tol, num_nodes): - """Tests a problem for different tolerances for the switch estimator.""" +def test_adapt_interpolation_info(): + """ + Test if the ``adapt_interpolation_info`` method of ``SwitchEstimator`` does what it is supposed to do. + """ + from pySDC.core.Step import step + from pySDC.core.Controller import controller from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit - from pySDC.helpers.stats_helper import get_sorted + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI t0 = 1.6 - Tend = 1.62 - problem = ExactDiscontinuousTestODE - - stats, t_switch_exact = discontinuousTestProblem_run( - problem=problem, - sweeper=generic_implicit, - quad_type='LOBATTO', - num_nodes=num_nodes, - t0=t0, - Tend=Tend, - dt=2e-2, - tol=tol, - ) - - # in this specific example only one event has to be found - switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] - assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' - - t_switch = switches[-1] - event_err = abs(t_switch - t_switch_exact) - assert event_err < 1e-14, f'Event time error for is not small enough!' + Tend = ExactDiscontinuousTestODE().t_switch_exact + eps = 1e-14 # choose this eps to enforce a sign chance in state function + dt = (Tend - t0) + eps + + level_params = { + 'dt': dt, + 'restol': 1e-13, + } + + sweeper_params = { + 'quad_type': 'LOBATTO', + 'num_nodes': 3, + 'QI': 'IE', + 'initial_guess': 'spread', + } + + step_params = { + 'maxiter': 10, + } + + # convergence controllers + convergence_controllers = dict() + switch_estimator_params = { + 'tol': 1e-10, + 'alpha': 1.0, + } + convergence_controllers.update({SwitchEstimator: switch_estimator_params}) + + convergence_controllers[BasicRestartingNonMPI] = { + 'max_restarts': 3, + 'crash_after_max_restarts': False, + } + + description = { + 'problem_class': ExactDiscontinuousTestODE, + 'problem_params': dict(), + 'sweeper_class': generic_implicit, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } + + # set up step + S = step(description=description) + L = S.levels[0] + P = L.prob + + # initialise switch estimator + switch_estimator_params = { + 'tol': 1e-10, + 'alpha': 1.0, + } + SE = SwitchEstimator(controller=controller, params=switch_estimator_params, description=description) + + SE.setup_status_variables(controller) + + # set initial time in the status of the level + L.status.time = t0 + + # compute initial value (using the exact function here) + L.u[0] = P.u_exact(L.time) + + # call prediction function to initialise nodes + L.sweep.predict() + + # compute the residual (we may be done already!) + L.sweep.compute_residual() + + # reset iteration counter + S.status.iter = 0 + + # run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough + while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol: + # this is where the nodes are actually updated according to the SDC formulas + L.sweep.update_nodes() + + # compute/update the residual + L.sweep.compute_residual() + + # increment the iteration counter + S.status.iter += 1 + + # # get state function + # _, _, state_function = P.get_switching_info(L.u, L.status.time) + # t_interp = SE.params.t_interp + # left_is_node = L.sweep.coll.left_is_node + # print('in test:', state_function) + # # t_interp_update, state_function_update = SE.adapt_interpolation_info( + # # t=L.status.time, left_is_node=left_is_node, t_interp=t_interp, state_function=state_function + # # ) + + # SE.get_new_step_size(controller, S) + # t_interp2 = SE.params.t_interp + # print(t_interp) + + +# @pytest.mark.base +# @pytest.mark.parametrize('num_nodes', [3, 4, 5]) +# def test_detection_at_boundary(num_nodes): +# """ +# This test checks whether a restart is executed or not when the event exactly occurs at the boundary. In this case, +# no restart should be done because occuring the event at the boundary means that the event is already resolved well, +# i.e., the state function there should have a value close to zero. +# """ +# from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit +# from pySDC.helpers.stats_helper import get_sorted + +# problem = ExactDiscontinuousTestODE +# t0 = 1.6 +# Tend = ExactDiscontinuousTestODE().t_switch_exact +# dt = Tend - t0 + +# stats, _ = discontinuousTestProblem_run( +# problem=problem, +# sweeper=generic_implicit, +# quad_type='LOBATTO', +# num_nodes=num_nodes, +# t0=t0, +# Tend=Tend, +# dt=dt, +# tol=1e-10, +# ) + +# sum_restarts = np.sum(np.array(get_sorted(stats, type='restart', sortby='time', recomputed=None))[:, 1]) +# assert sum_restarts == 0, 'Event occurs at boundary, but restart(s) are executed anyway!' + + +# @pytest.mark.base +# @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) +# @pytest.mark.parametrize('num_nodes', [3, 4, 5]) +# def test_all_tolerances_ODE(tol, num_nodes): +# r""" +# Tests a problem for different tolerances for the switch estimator. Here, a dummy problem of +# ``DiscontinuousTestODE`` is used, where the dynamics of the differential equation is replaced by its +# exact dynamics to see if the switch estimator predicts the event correctly. +# """ +# from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit +# from pySDC.helpers.stats_helper import get_sorted + +# t0 = 1.6 +# Tend = 1.62 +# problem = ExactDiscontinuousTestODE + +# stats, t_switch_exact = discontinuousTestProblem_run( +# problem=problem, +# sweeper=generic_implicit, +# quad_type='LOBATTO', +# num_nodes=num_nodes, +# t0=t0, +# Tend=Tend, +# dt=2e-2, +# tol=tol, +# ) + +# # in this specific example only one event has to be found +# switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] +# assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' + +# t_switch = switches[-1] +# event_err = abs(t_switch - t_switch_exact) +# assert event_err < 1e-14, f'Event time error for is not small enough!' From 4d98155e1212b176960f086539233a56f5ee57b4 Mon Sep 17 00:00:00 2001 From: lisawim Date: Fri, 27 Oct 2023 07:35:15 +0200 Subject: [PATCH 04/29] Started with test to check adapt_interpolation_info [WIP] --- .../test_pintsime/test_SwitchEstimator.py | 306 +++++++++--------- 1 file changed, 153 insertions(+), 153 deletions(-) diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 02620ebff6..6f8752b9f9 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -169,101 +169,101 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten return stats, t_switch_exact -@pytest.mark.base -def test_adapt_interpolation_info(): - """ - Test if the ``adapt_interpolation_info`` method of ``SwitchEstimator`` does what it is supposed to do. - """ - from pySDC.core.Step import step - from pySDC.core.Controller import controller - from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit - from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI - from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator - from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI - - t0 = 1.6 - Tend = ExactDiscontinuousTestODE().t_switch_exact - eps = 1e-14 # choose this eps to enforce a sign chance in state function - dt = (Tend - t0) + eps - - level_params = { - 'dt': dt, - 'restol': 1e-13, - } - - sweeper_params = { - 'quad_type': 'LOBATTO', - 'num_nodes': 3, - 'QI': 'IE', - 'initial_guess': 'spread', - } - - step_params = { - 'maxiter': 10, - } - - # convergence controllers - convergence_controllers = dict() - switch_estimator_params = { - 'tol': 1e-10, - 'alpha': 1.0, - } - convergence_controllers.update({SwitchEstimator: switch_estimator_params}) - - convergence_controllers[BasicRestartingNonMPI] = { - 'max_restarts': 3, - 'crash_after_max_restarts': False, - } - - description = { - 'problem_class': ExactDiscontinuousTestODE, - 'problem_params': dict(), - 'sweeper_class': generic_implicit, - 'sweeper_params': sweeper_params, - 'level_params': level_params, - 'step_params': step_params, - 'convergence_controllers': convergence_controllers, - } - - # set up step - S = step(description=description) - L = S.levels[0] - P = L.prob - - # initialise switch estimator - switch_estimator_params = { - 'tol': 1e-10, - 'alpha': 1.0, - } - SE = SwitchEstimator(controller=controller, params=switch_estimator_params, description=description) - - SE.setup_status_variables(controller) - - # set initial time in the status of the level - L.status.time = t0 - - # compute initial value (using the exact function here) - L.u[0] = P.u_exact(L.time) - - # call prediction function to initialise nodes - L.sweep.predict() - - # compute the residual (we may be done already!) - L.sweep.compute_residual() - - # reset iteration counter - S.status.iter = 0 - - # run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough - while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol: - # this is where the nodes are actually updated according to the SDC formulas - L.sweep.update_nodes() - - # compute/update the residual - L.sweep.compute_residual() +# @pytest.mark.base +# def test_adapt_interpolation_info(): +# """ +# Test if the ``adapt_interpolation_info`` method of ``SwitchEstimator`` does what it is supposed to do. +# """ +# from pySDC.core.Step import step +# from pySDC.core.Controller import controller +# from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit +# from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI +# from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator +# from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI - # increment the iteration counter - S.status.iter += 1 +# t0 = 1.6 +# Tend = ExactDiscontinuousTestODE().t_switch_exact +# eps = 1e-14 # choose this eps to enforce a sign chance in state function +# dt = (Tend - t0) + eps + +# level_params = { +# 'dt': dt, +# 'restol': 1e-13, +# } + +# sweeper_params = { +# 'quad_type': 'LOBATTO', +# 'num_nodes': 3, +# 'QI': 'IE', +# 'initial_guess': 'spread', +# } + +# step_params = { +# 'maxiter': 10, +# } + +# # convergence controllers +# convergence_controllers = dict() +# switch_estimator_params = { +# 'tol': 1e-10, +# 'alpha': 1.0, +# } +# convergence_controllers.update({SwitchEstimator: switch_estimator_params}) + +# convergence_controllers[BasicRestartingNonMPI] = { +# 'max_restarts': 3, +# 'crash_after_max_restarts': False, +# } + +# description = { +# 'problem_class': ExactDiscontinuousTestODE, +# 'problem_params': dict(), +# 'sweeper_class': generic_implicit, +# 'sweeper_params': sweeper_params, +# 'level_params': level_params, +# 'step_params': step_params, +# 'convergence_controllers': convergence_controllers, +# } + +# # set up step +# S = step(description=description) +# L = S.levels[0] +# P = L.prob + +# # initialise switch estimator +# switch_estimator_params = { +# 'tol': 1e-10, +# 'alpha': 1.0, +# } +# SE = SwitchEstimator(controller=controller, params=switch_estimator_params, description=description) + +# SE.setup_status_variables(controller) + +# # set initial time in the status of the level +# L.status.time = t0 + +# # compute initial value (using the exact function here) +# L.u[0] = P.u_exact(L.time) + +# # call prediction function to initialise nodes +# L.sweep.predict() + +# # compute the residual (we may be done already!) +# L.sweep.compute_residual() + +# # reset iteration counter +# S.status.iter = 0 + +# # run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough +# while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol: +# # this is where the nodes are actually updated according to the SDC formulas +# L.sweep.update_nodes() + +# # compute/update the residual +# L.sweep.compute_residual() + +# # increment the iteration counter +# S.status.iter += 1 # # get state function # _, _, state_function = P.get_switching_info(L.u, L.status.time) @@ -279,68 +279,68 @@ def test_adapt_interpolation_info(): # print(t_interp) -# @pytest.mark.base -# @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -# def test_detection_at_boundary(num_nodes): -# """ -# This test checks whether a restart is executed or not when the event exactly occurs at the boundary. In this case, -# no restart should be done because occuring the event at the boundary means that the event is already resolved well, -# i.e., the state function there should have a value close to zero. -# """ -# from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -# from pySDC.helpers.stats_helper import get_sorted +@pytest.mark.base +@pytest.mark.parametrize('num_nodes', [3, 4, 5]) +def test_detection_at_boundary(num_nodes): + """ + This test checks whether a restart is executed or not when the event exactly occurs at the boundary. In this case, + no restart should be done because occuring the event at the boundary means that the event is already resolved well, + i.e., the state function there should have a value close to zero. + """ + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.helpers.stats_helper import get_sorted -# problem = ExactDiscontinuousTestODE -# t0 = 1.6 -# Tend = ExactDiscontinuousTestODE().t_switch_exact -# dt = Tend - t0 + problem = ExactDiscontinuousTestODE + t0 = 1.6 + Tend = ExactDiscontinuousTestODE().t_switch_exact + dt = Tend - t0 -# stats, _ = discontinuousTestProblem_run( -# problem=problem, -# sweeper=generic_implicit, -# quad_type='LOBATTO', -# num_nodes=num_nodes, -# t0=t0, -# Tend=Tend, -# dt=dt, -# tol=1e-10, -# ) + stats, _ = discontinuousTestProblem_run( + problem=problem, + sweeper=generic_implicit, + quad_type='LOBATTO', + num_nodes=num_nodes, + t0=t0, + Tend=Tend, + dt=dt, + tol=1e-10, + ) -# sum_restarts = np.sum(np.array(get_sorted(stats, type='restart', sortby='time', recomputed=None))[:, 1]) -# assert sum_restarts == 0, 'Event occurs at boundary, but restart(s) are executed anyway!' + sum_restarts = np.sum(np.array(get_sorted(stats, type='restart', sortby='time', recomputed=None))[:, 1]) + assert sum_restarts == 0, 'Event occurs at boundary, but restart(s) are executed anyway!' -# @pytest.mark.base -# @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) -# @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -# def test_all_tolerances_ODE(tol, num_nodes): -# r""" -# Tests a problem for different tolerances for the switch estimator. Here, a dummy problem of -# ``DiscontinuousTestODE`` is used, where the dynamics of the differential equation is replaced by its -# exact dynamics to see if the switch estimator predicts the event correctly. -# """ -# from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -# from pySDC.helpers.stats_helper import get_sorted +@pytest.mark.base +@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) +@pytest.mark.parametrize('num_nodes', [3, 4, 5]) +def test_all_tolerances_ODE(tol, num_nodes): + r""" + Tests a problem for different tolerances for the switch estimator. Here, a dummy problem of + ``DiscontinuousTestODE`` is used, where the dynamics of the differential equation is replaced by its + exact dynamics to see if the switch estimator predicts the event correctly. + """ + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.helpers.stats_helper import get_sorted -# t0 = 1.6 -# Tend = 1.62 -# problem = ExactDiscontinuousTestODE - -# stats, t_switch_exact = discontinuousTestProblem_run( -# problem=problem, -# sweeper=generic_implicit, -# quad_type='LOBATTO', -# num_nodes=num_nodes, -# t0=t0, -# Tend=Tend, -# dt=2e-2, -# tol=tol, -# ) - -# # in this specific example only one event has to be found -# switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] -# assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' - -# t_switch = switches[-1] -# event_err = abs(t_switch - t_switch_exact) -# assert event_err < 1e-14, f'Event time error for is not small enough!' + t0 = 1.6 + Tend = 1.62 + problem = ExactDiscontinuousTestODE + + stats, t_switch_exact = discontinuousTestProblem_run( + problem=problem, + sweeper=generic_implicit, + quad_type='LOBATTO', + num_nodes=num_nodes, + t0=t0, + Tend=Tend, + dt=2e-2, + tol=tol, + ) + + # in this specific example only one event has to be found + switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] + assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' + + t_switch = switches[-1] + event_err = abs(t_switch - t_switch_exact) + assert event_err < 1e-14, f'Event time error for is not small enough!' From 42aaaf598a30d4befdc30e57a64f40a0f88f126a Mon Sep 17 00:00:00 2001 From: lisawim Date: Fri, 27 Oct 2023 18:32:17 +0200 Subject: [PATCH 05/29] Added test for SE.adapt_interpolation_info() --- pySDC/projects/PinTSimE/switch_estimator.py | 34 ++- .../test_pintsime/test_SwitchEstimator.py | 279 ++++++++++-------- 2 files changed, 174 insertions(+), 139 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 33c9749c06..b232a53697 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -48,6 +48,7 @@ def setup(self, controller, params, description): 'nodes': coll.nodes, 'tol_zero': 1e-13, 't_interp': [], + 'state_function': [], } return {**defaults, **params} @@ -88,27 +89,30 @@ def get_new_step_size(self, controller, S, **kwargs): """ L = S.levels[0] - print('conv', CheckConvergence.check_convergence(S)) + if CheckConvergence.check_convergence(S): - self.status.switch_detected, m_guess, state_function = L.prob.get_switching_info(L.u, L.time) - print('detected:', self.status.switch_detected, state_function) + self.status.switch_detected, m_guess, self.params.state_function = L.prob.get_switching_info(L.u, L.time) + if self.status.switch_detected: self.params.t_interp = [L.time + L.dt * self.params.nodes[m] for m in range(len(self.params.nodes))] - self.params.t_interp, state_function = self.adapt_interpolation_info( - L.time, L.sweep.coll.left_is_node, self.params.t_interp, state_function + self.params.t_interp, self.params.state_function = self.adapt_interpolation_info( + L.time, L.sweep.coll.left_is_node, self.params.t_interp, self.params.state_function ) - # print(t_interp, state_function) + # when the state function is already close to zero the event is already resolved well - if abs(state_function[-1]) <= self.params.tol_zero or abs(state_function[0]) <= self.params.tol_zero: - if abs(state_function[0]) <= self.params.tol_zero: + if ( + abs(self.params.state_function[-1]) <= self.params.tol_zero + or abs(self.params.state_function[0]) <= self.params.tol_zero + ): + if abs(self.params.state_function[0]) <= self.params.tol_zero: t_switch = self.params.t_interp[0] boundary = 'left' - elif abs(state_function[-1]) <= self.params.tol_zero: + elif abs(self.params.state_function[-1]) <= self.params.tol_zero: boundary = 'right' t_switch = self.params.t_interp[-1] msg = f"The value of state function is close to zero, thus event time is already close enough to the {boundary} end point!" - # self.log(msg, S) + self.log(msg, S) self.log_event_time( controller.hooks[0], S.status.slot, L.time, L.level_index, L.status.sweep, t_switch ) @@ -117,8 +121,8 @@ def get_new_step_size(self, controller, S, **kwargs): self.status.is_zero = True # intermediate value theorem states that a root is contained in current step - if state_function[0] * state_function[-1] < 0 and self.status.is_zero is None: - self.status.t_switch = self.get_switch(self.params.t_interp, state_function, m_guess) + if self.params.state_function[0] * self.params.state_function[-1] < 0 and self.status.is_zero is None: + self.status.t_switch = self.get_switch(self.params.t_interp, self.params.state_function, m_guess) controller.hooks[0].add_to_stats( process=S.status.slot, @@ -136,7 +140,7 @@ def get_new_step_size(self, controller, S, **kwargs): iter=0, sweep=L.status.sweep, type='h_all', - value=max([abs(item) for item in state_function]), + value=max([abs(item) for item in self.params.state_function]), ) if L.time < self.status.t_switch < L.time + L.dt: dt_switch = (self.status.t_switch - L.time) * self.params.alpha @@ -414,7 +418,9 @@ def eval(self, t): for i in range(1, self.n): if t >= self.ti[i - 1] or t <= self.ti[i]: - p = self.yi[i - 1] + (self.yi[i] - self.yi[i - 1]) / (self.ti[i] - self.ti[i - 1]) * (t - self.ti[i - 1]) + p = self.yi[i - 1] + (self.yi[i] - self.yi[i - 1]) / (self.ti[i] - self.ti[i - 1]) * ( + t - self.ti[i - 1] + ) else: p = 0 raise ParameterError('Interpolant is only evaluated for time occuring between time nodes') diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 6f8752b9f9..c66d8f9542 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -45,7 +45,7 @@ def eval_f(self, u, t): else: f[:] = np.exp(t) return f - + def solve_system(self, rhs, factor, u0, t): """ Just return the exact solution... @@ -67,16 +67,18 @@ def solve_system(self, rhs, factor, u0, t): The solution as mesh. """ - return self.u_exact(t) + return self.u_exact(t) -def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Tend, dt, tol): - """ - Simulates one run of a discontinuous test problem class for one specific tolerance specified in - the parameters. For testing, only one time step should be considered. +def get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, tol): + r""" + This function prepares the controller for one simulation run. Based on function in + ``pySDC.tests.test_convergence_controllers.test_adaptivity.py`` Parameters ---------- + switch_estimator : pySDC.core.ConvergenceController + Switch estimator. problem : pySDC.core.Problem Problem class to be simulated. sweeper: pySDC.core.Sweeper @@ -94,9 +96,7 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten tol : float Tolerance for the switch estimator. """ - from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI - from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI from pySDC.implementations.hooks.log_solution import LogSolution @@ -133,7 +133,7 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten 'tol': tol, 'alpha': 1.0, } - convergence_controllers.update({SwitchEstimator: switch_estimator_params}) + convergence_controllers.update({switch_estimator: switch_estimator_params}) convergence_controllers[BasicRestartingNonMPI] = { 'max_restarts': 3, @@ -151,13 +151,45 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten 'convergence_controllers': convergence_controllers, } + # instantiate controller + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + return controller + + +def discontinuousTestProblem_run(switch_estimator, problem, sweeper, quad_type, num_nodes, t0, Tend, dt, tol): + """ + Simulates one run of a discontinuous test problem class for one specific tolerance specified in + the parameters. For testing, only one time step should be considered. + + Parameters + ---------- + switch_estimator : pySDC.core.ConvergenceController + Switch Estimator. + problem : pySDC.core.Problem + Problem class to be simulated. + sweeper: pySDC.core.Sweeper + Sweeper used for the simulation. + quad_type : str + Type of quadrature used. + num_nodes : int + Number of collocation nodes. + t0 : float + Starting time. + Tend : float + End time. + dt : float + Time step size. + tol : float + Tolerance for the switch estimator. + """ + + controller = get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, tol) + # set time parameters t0 = t0 Tend = Tend - # instantiate controller - controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) - # get initial values on finest level P = controller.MS[0].levels[0].prob uinit = P.u_exact(t0) @@ -169,114 +201,86 @@ def discontinuousTestProblem_run(problem, sweeper, quad_type, num_nodes, t0, Ten return stats, t_switch_exact -# @pytest.mark.base -# def test_adapt_interpolation_info(): -# """ -# Test if the ``adapt_interpolation_info`` method of ``SwitchEstimator`` does what it is supposed to do. -# """ -# from pySDC.core.Step import step -# from pySDC.core.Controller import controller -# from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -# from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI -# from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator -# from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI - -# t0 = 1.6 -# Tend = ExactDiscontinuousTestODE().t_switch_exact -# eps = 1e-14 # choose this eps to enforce a sign chance in state function -# dt = (Tend - t0) + eps - -# level_params = { -# 'dt': dt, -# 'restol': 1e-13, -# } - -# sweeper_params = { -# 'quad_type': 'LOBATTO', -# 'num_nodes': 3, -# 'QI': 'IE', -# 'initial_guess': 'spread', -# } - -# step_params = { -# 'maxiter': 10, -# } - -# # convergence controllers -# convergence_controllers = dict() -# switch_estimator_params = { -# 'tol': 1e-10, -# 'alpha': 1.0, -# } -# convergence_controllers.update({SwitchEstimator: switch_estimator_params}) - -# convergence_controllers[BasicRestartingNonMPI] = { -# 'max_restarts': 3, -# 'crash_after_max_restarts': False, -# } - -# description = { -# 'problem_class': ExactDiscontinuousTestODE, -# 'problem_params': dict(), -# 'sweeper_class': generic_implicit, -# 'sweeper_params': sweeper_params, -# 'level_params': level_params, -# 'step_params': step_params, -# 'convergence_controllers': convergence_controllers, -# } - -# # set up step -# S = step(description=description) -# L = S.levels[0] -# P = L.prob - -# # initialise switch estimator -# switch_estimator_params = { -# 'tol': 1e-10, -# 'alpha': 1.0, -# } -# SE = SwitchEstimator(controller=controller, params=switch_estimator_params, description=description) - -# SE.setup_status_variables(controller) - -# # set initial time in the status of the level -# L.status.time = t0 - -# # compute initial value (using the exact function here) -# L.u[0] = P.u_exact(L.time) - -# # call prediction function to initialise nodes -# L.sweep.predict() - -# # compute the residual (we may be done already!) -# L.sweep.compute_residual() - -# # reset iteration counter -# S.status.iter = 0 - -# # run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough -# while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol: -# # this is where the nodes are actually updated according to the SDC formulas -# L.sweep.update_nodes() - -# # compute/update the residual -# L.sweep.compute_residual() - -# # increment the iteration counter -# S.status.iter += 1 - - # # get state function - # _, _, state_function = P.get_switching_info(L.u, L.status.time) - # t_interp = SE.params.t_interp - # left_is_node = L.sweep.coll.left_is_node - # print('in test:', state_function) - # # t_interp_update, state_function_update = SE.adapt_interpolation_info( - # # t=L.status.time, left_is_node=left_is_node, t_interp=t_interp, state_function=state_function - # # ) - - # SE.get_new_step_size(controller, S) - # t_interp2 = SE.params.t_interp - # print(t_interp) +@pytest.mark.base +@pytest.mark.parametrize('quad_type', ['LOBATTO', 'RADAU-RIGHT']) +def test_adapt_interpolation_info(quad_type): + r""" + Tests if the method ``adapt_interpolation_info`` does what it is supposed to do. + + - For ``quad_type='RADAU-RIGHT'``, the value at ``t0`` has to be added to the list of + ``state_function``, since it is no collocation node but can be used for the interpolation + anyway. + + - For ``quad_type='LOBATTO'``, the first value in ``state_function`` has to be removed + since ``t0`` is also a collocation node here and the state function would include double values. + + Parameters + ---------- + quad_type : str + Type of quadrature used. + """ + + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + + t0 = 1.6 + Tend = ExactDiscontinuousTestODE().t_switch_exact + eps = 1e-13 # choose this eps to enforce a sign chance in state function + dt = (Tend - t0) + eps + + tol = 1e-10 + num_nodes = 3 + + controller = get_controller( + switch_estimator=SwitchEstimator, + problem=ExactDiscontinuousTestODE, + sweeper=generic_implicit, + quad_type=quad_type, + num_nodes=num_nodes, + dt=dt, + tol=tol, + ) + + S = controller.MS[0] + L = S.levels[0] + P = L.prob + + # instance of switch estimator + SE = controller.convergence_controllers[ + np.arange(len(controller.convergence_controllers))[ + [type(me).__name__ == SwitchEstimator.__name__ for me in controller.convergence_controllers] + ][0] + ] + + S.status.slot = 0 + L.status.time = t0 + S.status.iter = 10 + L.status.residual = 0.0 + L.u[0] = P.u_exact(L.status.time) + + L.sweep.predict() + + # perform one update to get state function with different signs + L.sweep.update_nodes() + + SE.get_new_step_size(controller, S) + + t_interp, state_function = SE.params.t_interp, SE.params.state_function + + assert len(t_interp) == len( + state_function + ), 'Length of interpolation values does not match with length of list containing the state function' + + if quad_type == 'LOBATTO': + assert t_interp[0] != t_interp[1], 'Starting time from interpolation axis is not removed!' + assert ( + len(t_interp) == num_nodes + ), f'Number of values on interpolation axis does not match. Expected {num_nodes}, got {len(t_interp)}' + + elif quad_type == 'RADAU-RIGHT': + assert ( + len(t_interp) == num_nodes + 1 + ), f'Number of values on interpolation axis does not match. Expected {num_nodes + 1}, got {len(t_interp)}' @pytest.mark.base @@ -286,8 +290,15 @@ def test_detection_at_boundary(num_nodes): This test checks whether a restart is executed or not when the event exactly occurs at the boundary. In this case, no restart should be done because occuring the event at the boundary means that the event is already resolved well, i.e., the state function there should have a value close to zero. + + Parameters + ---------- + num_nodes : int + Number of collocation nodes. """ + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.helpers.stats_helper import get_sorted problem = ExactDiscontinuousTestODE @@ -296,6 +307,7 @@ def test_detection_at_boundary(num_nodes): dt = Tend - t0 stats, _ = discontinuousTestProblem_run( + switch_estimator=SwitchEstimator, problem=problem, sweeper=generic_implicit, quad_type='LOBATTO', @@ -313,13 +325,29 @@ def test_detection_at_boundary(num_nodes): @pytest.mark.base @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -def test_all_tolerances_ODE(tol, num_nodes): +@pytest.mark.parametrize('quad_type', ['LOBATTO', 'RADAU-RIGHT']) +def test_all_tolerances_ODE(tol, num_nodes, quad_type): r""" - Tests a problem for different tolerances for the switch estimator. Here, a dummy problem of - ``DiscontinuousTestODE`` is used, where the dynamics of the differential equation is replaced by its - exact dynamics to see if the switch estimator predicts the event correctly. + Here, the switch estimator is applied to a dummy problem of ``DiscontinuousTestODE``, + where the dynamics of the differential equation is replaced by its exact dynamics to see if + the switch estimator predicts the event correctly. The problem is tested for a combination + of different tolerances ``tol`` and different number of collocation nodes ``num_nodes``. + + Since the problem only uses the exact dynamics, the event should be predicted very accurately + by the switch estimator. + + Parameters + ---------- + tol : float + Tolerance for switch estimator. + num_nodes : int + Number of collocation nodes. + quad_type : str + Type of quadrature. """ + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.helpers.stats_helper import get_sorted t0 = 1.6 @@ -327,9 +355,10 @@ def test_all_tolerances_ODE(tol, num_nodes): problem = ExactDiscontinuousTestODE stats, t_switch_exact = discontinuousTestProblem_run( + switch_estimator=SwitchEstimator, problem=problem, sweeper=generic_implicit, - quad_type='LOBATTO', + quad_type=quad_type, num_nodes=num_nodes, t0=t0, Tend=Tend, From 2cad220571c30292cd5fd001b0b485cfef44ab7c Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 22 Feb 2024 18:13:48 +0100 Subject: [PATCH 06/29] Update linear interpolation + logging + changing tolerances --- pySDC/projects/PinTSimE/switch_estimator.py | 102 ++++++++++---------- 1 file changed, 53 insertions(+), 49 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index b232a53697..0b0dc70403 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -124,24 +124,16 @@ def get_new_step_size(self, controller, S, **kwargs): if self.params.state_function[0] * self.params.state_function[-1] < 0 and self.status.is_zero is None: self.status.t_switch = self.get_switch(self.params.t_interp, self.params.state_function, m_guess) - controller.hooks[0].add_to_stats( - process=S.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='switch_all', - value=self.status.t_switch, - ) - controller.hooks[0].add_to_stats( - process=S.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='h_all', - value=max([abs(item) for item in self.params.state_function]), + self.logging_during_estimation( + controller.hooks[0], + S.status.slot, + L.time, + L.level_index, + L.status.sweep, + self.status.t_switch, + self.params.state_function ) + if L.time < self.status.t_switch < L.time + L.dt: dt_switch = (self.status.t_switch - L.time) * self.params.alpha @@ -184,7 +176,6 @@ def get_new_step_size(self, controller, S, **kwargs): L.status.sweep, self.status.t_switch, ) - print('Occurs at boundary') L.prob.count_switches() self.status.switch_detected = False @@ -261,6 +252,27 @@ def log_event_time(controller_hooks, process, time, level, sweep, t_switch): value=t_switch, ) + @staticmethod + def logging_during_estimation(controller_hooks, process, time, level, sweep, t_switch, state_function): + controller_hooks.add_to_stats( + process=process, + time=time, + level=level, + iter=0, + sweep=sweep, + type='switch_all', + value=t_switch, + ) + controller_hooks.add_to_stats( + process=process, + time=time, + level=level, + iter=0, + sweep=sweep, + type='h_all', + value=max([abs(item) for item in state_function]), + ) + @staticmethod def get_switch(t_interp, state_function, m_guess): """ @@ -280,28 +292,12 @@ def get_switch(t_interp, state_function, m_guess): t_switch : float Time point of found event. """ - LinearInterpolator = LinearInterpolation(t_interp, state_function) - - def p(t): - """ - Simplifies the call of the interpolant. - - Parameters - ---------- - t : float - Time t at which the interpolant is called. - - Returns - ------- - p(t) : float - The value of the interpolated function at time t. - """ - return LinearInterpolator.eval(t) + p = lambda t: LinearInterpolator.eval(t) def fprime(t): """ - Computes the derivative of the scalar interpolant using finite differences. + Computes the derivative of the scalar interpolant using centered finite difference. Parameters ---------- @@ -313,11 +309,11 @@ def fprime(t): dp : float Derivative of interpolation p at time t. """ - dt_FD = 1e-10 - dp = (p(t + dt_FD) - p(t)) / dt_FD # forward difference + dt_FD = 1e-11 + dp = (p(t + dt_FD) - p(t - dt_FD)) / dt_FD return dp - newton_tol, newton_maxiter = 1e-15, 100 + newton_tol, newton_maxiter = 1e-12, 100 t_switch = newton(t_interp[m_guess], p, fprime, newton_tol, newton_maxiter) return t_switch @@ -382,15 +378,21 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): n = 0 while n < newton_maxiter: - if abs(p(x0)) < newton_tol or np.isnan(p(x0)) and np.isnan(fprime(x0)): + res = abs(p(x0)) + if res < newton_tol or np.isnan(p(x0)) and np.isnan(fprime(x0)) or np.isclose(fprime(x0), 0.0): break x0 -= 1.0 / fprime(x0) * p(x0) n += 1 - root = x0 + if n == newton_maxiter: + msg = 'Newton did not converge after %i iterations, error is %s' % (n, res) + else: + msg = f'Newton did converge after {n} iterations, error is {res}' + print(msg) + root = x0 return root @@ -404,6 +406,7 @@ def __init__(self, ti, yi): def eval(self, t): r""" Evaluates the linear interpolant at time :math:`t`. + Binary search is used to find the index in which piece :math:`t` lies. Parameters ---------- @@ -416,12 +419,13 @@ def eval(self, t): Value of the interpolant at time :math:`t`. """ - for i in range(1, self.n): - if t >= self.ti[i - 1] or t <= self.ti[i]: - p = self.yi[i - 1] + (self.yi[i] - self.yi[i - 1]) / (self.ti[i] - self.ti[i - 1]) * ( - t - self.ti[i - 1] - ) - else: - p = 0 - raise ParameterError('Interpolant is only evaluated for time occuring between time nodes') + ind = np.searchsorted(self.ti, t) + if ind == 0: + p = self.yi[0] + elif ind == len(self.ti): + p = self.yi[-1] + else: + p = self.yi[ind - 1] + (self.yi[ind] - self.yi[ind - 1]) / (self.ti[ind] - self.ti[ind - 1]) * ( + t - self.ti[ind - 1] + ) return p From d3c3217e1a1e65668a52d92682e254178f1ce7ad Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 22 Feb 2024 18:14:40 +0100 Subject: [PATCH 07/29] Test for linear interpolation + update of other test --- .../test_pintsime/test_SwitchEstimator.py | 49 +++++++++++++++++-- 1 file changed, 45 insertions(+), 4 deletions(-) diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index c66d8f9542..13789256b5 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -2,8 +2,6 @@ import pytest from pySDC.implementations.problem_classes.DiscontinuousTestODE import DiscontinuousTestODE -from pySDC.core.Problem import ptype -from pySDC.implementations.datatype_classes.mesh import mesh class ExactDiscontinuousTestODE(DiscontinuousTestODE): @@ -131,7 +129,7 @@ def get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, convergence_controllers = dict() switch_estimator_params = { 'tol': tol, - 'alpha': 1.0, + 'alpha': 0.95, } convergence_controllers.update({switch_estimator: switch_estimator_params}) @@ -201,6 +199,49 @@ def discontinuousTestProblem_run(switch_estimator, problem, sweeper, quad_type, return stats, t_switch_exact +def getTestDict(): + testDict = { + 0 : { + 't': [0, 2], + 'p': [1, 5], + 'tMid': [1], + 'pMid': [3], + }, + 1 : { + 't': [1, 2, 4, 5], + 'p': [6, 4, 3, 7], + 'tMid': [1.5, 3], + 'pMid': [5, 3.5], + }, + } + return testDict + + +@pytest.mark.base +@pytest.mark.parametrize("key", [0, 1]) +def testInterpolationValues(key): + """ + Test for linear interpolation class in switch estimator. Linear interpolation is tested against + values in testDict that contains values for interpolation computed by hand. + Further, NumPy's routine interp is used as reference to have a second test on hand. + """ + from pySDC.projects.PinTSimE.switch_estimator import LinearInterpolation + + testDict = getTestDict() + testSet = testDict[key] + t, p = testSet['t'], testSet['p'] + tMid, pMid = testSet['tMid'], testSet['pMid'] + LinearInterpolator = LinearInterpolation(t, p) + + for m in range(len(t)): + assert LinearInterpolator.eval(t[m]) == p[m] + assert np.interp(t[m], t, p) == LinearInterpolator.eval(t[m]) + + for m in range(len(tMid)): + assert LinearInterpolator.eval(tMid[m]) == pMid[m] + assert np.interp(tMid[m], t, p) == LinearInterpolator.eval(tMid[m]) + + @pytest.mark.base @pytest.mark.parametrize('quad_type', ['LOBATTO', 'RADAU-RIGHT']) def test_adapt_interpolation_info(quad_type): @@ -372,4 +413,4 @@ def test_all_tolerances_ODE(tol, num_nodes, quad_type): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert event_err < 1e-14, f'Event time error for is not small enough!' + assert np.isclose(event_err, 0, atol=3e-12), f'Event time error is not small enough!' From 19ab5bab612cad500a4b6514d3c8a223d46676b6 Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 22 Feb 2024 18:28:43 +0100 Subject: [PATCH 08/29] Correction for finite difference + adaption tolerance --- pySDC/projects/PinTSimE/switch_estimator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 0b0dc70403..7cfcbac669 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -309,8 +309,8 @@ def fprime(t): dp : float Derivative of interpolation p at time t. """ - dt_FD = 1e-11 - dp = (p(t + dt_FD) - p(t - dt_FD)) / dt_FD + dt_FD = 1e-12 + dp = (p(t + dt_FD) - p(t - dt_FD)) / (2 * dt_FD) return dp newton_tol, newton_maxiter = 1e-12, 100 @@ -364,7 +364,7 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): p : callable Interpolated function where Newton's method is applied at. fprime : callable - Approximated erivative of p using finite differences. + Approximated derivative of p using finite differences. newton_tol : float Tolerance for termination. newton_maxiter : int From 2268cfa807d7e17e82fd76d1a21073ff53689844 Mon Sep 17 00:00:00 2001 From: lisawim Date: Fri, 23 Feb 2024 17:18:31 +0100 Subject: [PATCH 09/29] Added test for DAE case for SE --- pySDC/projects/PinTSimE/switch_estimator.py | 2 +- .../test_pintsime/test_SwitchEstimator.py | 143 +++++++++++++++++- 2 files changed, 139 insertions(+), 6 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 7cfcbac669..b2cba197ca 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -131,7 +131,7 @@ def get_new_step_size(self, controller, S, **kwargs): L.level_index, L.status.sweep, self.status.t_switch, - self.params.state_function + self.params.state_function, ) if L.time < self.status.t_switch < L.time + L.dt: diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 13789256b5..2833f26eff 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -2,11 +2,12 @@ import pytest from pySDC.implementations.problem_classes.DiscontinuousTestODE import DiscontinuousTestODE +from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE class ExactDiscontinuousTestODE(DiscontinuousTestODE): r""" - Dummy problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestODE``. + Dummy ODE problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestODE``. """ def __init__(self, newton_maxiter=100, newton_tol=1e-8): @@ -25,7 +26,7 @@ def eval_f(self, u, t): ---------- u : dtype_u Exact value of u. - t : _type_ + t : float Time :math:`t`. Returns @@ -68,6 +69,86 @@ def solve_system(self, rhs, factor, u0, t): return self.u_exact(t) +class ExactDiscontinuousTestDAE(DiscontinuousTestDAE): + r""" + Dummy DAE problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestDAE``. + """ + + def __init__(self, newton_tol=1e-8): + """Initialization routine""" + super().__init__(newton_tol) + + self.t_switch_exact = np.arccosh(50) + self.t_switch = None + self.nswitches = 0 + + def eval_f(self, u, du, t): + r""" + Returns the exact right-hand side of the problem. The first components in both + cases of ``f`` are set to 1 do avoid getting a zero residual (enforced the sweeper to stop + since "convergence is reached"). + + Parameters + ---------- + u : dtype_u + Exact value of u. + du : dtype_u + Current values of the derivative of the numerical solution at time t. + t : float + Time :math:`t`. + + Returns + ------- + f : dtype_f + Right-hand side. + """ + + f = self.dtype_f(self.init) + + t_switch = np.inf if self.t_switch is None else self.t_switch + + h = 2 * u[0] - 100 + f = self.dtype_f(self.init) + + if h >= 0 or t >= t_switch: + f[:] = ( + 1, + 0, + ) + else: + f[:] = ( + 1, + 0, + ) + return f + + def solve_system(self, impl_sys, u0, t): + r""" + Just returns the derivative of the exact solution. + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + if t <= self.t_switch_exact: + me[:] = (np.sinh(t), np.cosh(t)) + else: + me[:] = (np.sinh(self.t_switch_exact), np.cosh(self.t_switch_exact)) + return me + + def get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, tol): r""" This function prepares the controller for one simulation run. Based on function in @@ -116,7 +197,7 @@ def get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, # initialize step parameters step_params = { - 'maxiter': 10, + 'maxiter': 1, } # initialize controller parameters @@ -201,13 +282,13 @@ def discontinuousTestProblem_run(switch_estimator, problem, sweeper, quad_type, def getTestDict(): testDict = { - 0 : { + 0: { 't': [0, 2], 'p': [1, 5], 'tMid': [1], 'pMid': [3], }, - 1 : { + 1: { 't': [1, 2, 4, 5], 'p': [6, 4, 3, 7], 'tMid': [1.5, 3], @@ -414,3 +495,55 @@ def test_all_tolerances_ODE(tol, num_nodes, quad_type): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) assert np.isclose(event_err, 0, atol=3e-12), f'Event time error is not small enough!' + + +@pytest.mark.base +@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) +@pytest.mark.parametrize('num_nodes', [3, 4, 5]) +def test_all_tolerances_DAE(tol, num_nodes): + r""" + In this test, the switch estimator is applied to a DAE dummy problem of ``DiscontinuousTestDAE``, + where the dynamics of the differential equation is replaced by its exact dynamics to see if + the switch estimator predicts the event correctly. The problem is tested for a combination + of different tolerances ``tol`` and different number of collocation nodes ``num_nodes``. + + Since the problem only uses the exact dynamics, the event should be predicted very accurately + by the switch estimator. + + Parameters + ---------- + tol : float + Tolerance for switch estimator. + num_nodes : int + Number of collocation nodes. + quad_type : str + Type of quadrature. + """ + + from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.helpers.stats_helper import get_sorted + + t0 = 4.6 + Tend = 4.62 + problem = ExactDiscontinuousTestDAE + + stats, t_switch_exact = discontinuousTestProblem_run( + switch_estimator=SwitchEstimator, + problem=problem, + sweeper=fully_implicit_DAE, + quad_type='RADAU-RIGHT', + num_nodes=num_nodes, + t0=t0, + Tend=Tend, + dt=2e-2, + tol=tol, + ) + + # in this specific example only one event has to be found + switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] + assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' + + t_switch = switches[-1] + event_err = abs(t_switch - t_switch_exact) + assert np.isclose(event_err, 0, atol=1e-14), f'Event time error {event_err} is not small enough!' From 9a6350c1915154b6bf38d632ffceeb447abf6ca8 Mon Sep 17 00:00:00 2001 From: lisawim Date: Sat, 24 Feb 2024 07:24:39 +0100 Subject: [PATCH 10/29] Choice of FD seems to be important for performance of SE --- pySDC/projects/PinTSimE/switch_estimator.py | 19 +++++++++++++++---- .../test_pintsime/test_SwitchEstimator.py | 4 ++-- 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index b2cba197ca..6765699114 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -294,6 +294,10 @@ def get_switch(t_interp, state_function, m_guess): """ LinearInterpolator = LinearInterpolation(t_interp, state_function) p = lambda t: LinearInterpolator.eval(t) + if state_function[0] < 0 and state_function[-1] >= 0: + chooseFD = 'centered' + else: + chooseFD = 'backward' def fprime(t): """ @@ -309,8 +313,15 @@ def fprime(t): dp : float Derivative of interpolation p at time t. """ - dt_FD = 1e-12 - dp = (p(t + dt_FD) - p(t - dt_FD)) / (2 * dt_FD) + + if chooseFD == 'centered': + dt_FD = 1e-12 + dp = (p(t + dt_FD) - p(t - dt_FD)) / (2 * dt_FD) + elif chooseFD == 'backward': + dt_FD = 1e-10 + dp = (p(t) - p(t - dt_FD)) / (dt_FD) + else: + raise NotImplementedError return dp newton_tol, newton_maxiter = 1e-12, 100 @@ -387,9 +398,9 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): n += 1 if n == newton_maxiter: - msg = 'Newton did not converge after %i iterations, error is %s' % (n, res) + msg = f'Newton did not converge after {n} iterations, error is {res}' else: - msg = f'Newton did converge after {n} iterations, error is {res}' + msg = f'Newton did converge after {n} iterations, error for root {x0} is {res}' print(msg) root = x0 diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 2833f26eff..00d9a6204d 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -494,7 +494,7 @@ def test_all_tolerances_ODE(tol, num_nodes, quad_type): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=3e-12), f'Event time error is not small enough!' + assert np.isclose(event_err, 0, atol=1.2e-11), f'Event time error {event_err} is not small enough!' @pytest.mark.base @@ -546,4 +546,4 @@ def test_all_tolerances_DAE(tol, num_nodes): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=1e-14), f'Event time error {event_err} is not small enough!' + assert np.isclose(event_err, 0, atol=8e-14), f'Event time error {event_err} is not small enough!' From e1295a326260fff3a58e5926c4bddbe579d3dcde Mon Sep 17 00:00:00 2001 From: lisawim Date: Sat, 24 Feb 2024 09:39:53 +0100 Subject: [PATCH 11/29] Removed attributes from dummy probs (since the parent classes have it) --- .../test_projects/test_pintsime/test_SwitchEstimator.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 00d9a6204d..ac99344c63 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -14,10 +14,6 @@ def __init__(self, newton_maxiter=100, newton_tol=1e-8): """Initialization routine""" super().__init__(newton_maxiter, newton_tol) - self.t_switch_exact = np.log(5) - self.t_switch = None - self.nswitches = 0 - def eval_f(self, u, t): """ Derivative. @@ -78,10 +74,6 @@ def __init__(self, newton_tol=1e-8): """Initialization routine""" super().__init__(newton_tol) - self.t_switch_exact = np.arccosh(50) - self.t_switch = None - self.nswitches = 0 - def eval_f(self, u, du, t): r""" Returns the exact right-hand side of the problem. The first components in both From 421428ef37793bcef15ac19923ff9da883ea1334 Mon Sep 17 00:00:00 2001 From: lisawim Date: Sun, 25 Feb 2024 09:20:14 +0100 Subject: [PATCH 12/29] Test for dummy problems + using functions from battery_model.py --- pySDC/projects/PinTSimE/battery_model.py | 17 +- .../test_pintsime/test_SwitchEstimator.py | 389 ++++++++++-------- 2 files changed, 230 insertions(+), 176 deletions(-) diff --git a/pySDC/projects/PinTSimE/battery_model.py b/pySDC/projects/PinTSimE/battery_model.py index 20f2e11f1c..7301aceaf4 100644 --- a/pySDC/projects/PinTSimE/battery_model.py +++ b/pySDC/projects/PinTSimE/battery_model.py @@ -106,7 +106,7 @@ def generateDescription( controller_params : dict Parameters needed for a controller run. """ - + print(use_adaptivity) # initialize level parameters level_params = { 'restol': -1 if use_adaptivity else restol, @@ -167,10 +167,13 @@ def generateDescription( 'convergence_controllers': convergence_controllers, } - return description, controller_params + # instantiate controller + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + return description, controller_params, controller -def controllerRun(description, controller_params, t0, Tend, exact_event_time_avail=False): +def controllerRun(description, controller_params, controller, t0, Tend, exact_event_time_avail=False): """ Executes a controller run for a problem defined in the description. @@ -180,6 +183,8 @@ def controllerRun(description, controller_params, t0, Tend, exact_event_time_ava Contains all information for a controller run. controller_params : dict Parameters needed for a controller run. + controller : pySDC.core.Controller + Controller to do the stuff. t0 : float Starting time of simulation. Tend : float @@ -193,9 +198,6 @@ def controllerRun(description, controller_params, t0, Tend, exact_event_time_ava Raw statistics from a controller run. """ - # instantiate controller - controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) - # get initial values on finest level P = controller.MS[0].levels[0].prob uinit = P.u_exact(t0) @@ -360,7 +362,7 @@ def runSimulation(problem, sweeper, all_params, use_adaptivity, use_detection, h restol = -1 if use_A else handling_params['restol'] - description, controller_params = generateDescription( + description, controller_params, controller = generateDescription( dt=dt, problem=problem, sweeper=sweeper, @@ -381,6 +383,7 @@ def runSimulation(problem, sweeper, all_params, use_adaptivity, use_detection, h stats, t_switch_exact = controllerRun( description=description, controller_params=controller_params, + controller=controller, t0=interval[0], Tend=interval[-1], exact_event_time_avail=handling_params['exact_event_time_avail'], diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index ac99344c63..954198dd01 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -114,7 +114,7 @@ def eval_f(self, u, du, t): ) return f - def solve_system(self, impl_sys, u0, t): + def solve_system(self, impl_sys, u0, t, **kwargs): r""" Just returns the derivative of the exact solution. @@ -141,137 +141,6 @@ def solve_system(self, impl_sys, u0, t): return me -def get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, tol): - r""" - This function prepares the controller for one simulation run. Based on function in - ``pySDC.tests.test_convergence_controllers.test_adaptivity.py`` - - Parameters - ---------- - switch_estimator : pySDC.core.ConvergenceController - Switch estimator. - problem : pySDC.core.Problem - Problem class to be simulated. - sweeper: pySDC.core.Sweeper - Sweeper used for the simulation. - quad_type : str - Type of quadrature used. - num_nodes : int - Number of collocation nodes. - t0 : float - Starting time. - Tend : float - End time. - dt : float - Time step size. - tol : float - Tolerance for the switch estimator. - """ - from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI - from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI - - from pySDC.implementations.hooks.log_solution import LogSolution - from pySDC.implementations.hooks.log_restarts import LogRestarts - - # initialize level parameters - level_params = { - 'restol': 1e-13, - 'dt': dt, - } - - # initialize sweeper parameters - sweeper_params = { - 'quad_type': quad_type, - 'num_nodes': num_nodes, - 'QI': 'IE', - 'initial_guess': 'spread', - } - - # initialize step parameters - step_params = { - 'maxiter': 1, - } - - # initialize controller parameters - controller_params = { - 'logger_level': 30, - 'hook_class': [LogSolution, LogRestarts], - } - - # convergence controllers - convergence_controllers = dict() - switch_estimator_params = { - 'tol': tol, - 'alpha': 0.95, - } - convergence_controllers.update({switch_estimator: switch_estimator_params}) - - convergence_controllers[BasicRestartingNonMPI] = { - 'max_restarts': 3, - 'crash_after_max_restarts': False, - } - - # fill description dictionary for easy step instantiation - description = { - 'problem_class': problem, # pass problem class - 'problem_params': dict(), # problem_params, # pass problem parameters - 'sweeper_class': sweeper, # pass sweeper - 'sweeper_params': sweeper_params, # pass sweeper parameters - 'level_params': level_params, # pass level parameters - 'step_params': step_params, - 'convergence_controllers': convergence_controllers, - } - - # instantiate controller - controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) - - return controller - - -def discontinuousTestProblem_run(switch_estimator, problem, sweeper, quad_type, num_nodes, t0, Tend, dt, tol): - """ - Simulates one run of a discontinuous test problem class for one specific tolerance specified in - the parameters. For testing, only one time step should be considered. - - Parameters - ---------- - switch_estimator : pySDC.core.ConvergenceController - Switch Estimator. - problem : pySDC.core.Problem - Problem class to be simulated. - sweeper: pySDC.core.Sweeper - Sweeper used for the simulation. - quad_type : str - Type of quadrature used. - num_nodes : int - Number of collocation nodes. - t0 : float - Starting time. - Tend : float - End time. - dt : float - Time step size. - tol : float - Tolerance for the switch estimator. - """ - - controller = get_controller(switch_estimator, problem, sweeper, quad_type, num_nodes, dt, tol) - - # set time parameters - t0 = t0 - Tend = Tend - - # get initial values on finest level - P = controller.MS[0].levels[0].prob - uinit = P.u_exact(t0) - t_switch_exact = P.t_switch_exact - - # call main function to get things done... - uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - - return stats, t_switch_exact - - def getTestDict(): testDict = { 0: { @@ -290,6 +159,68 @@ def getTestDict(): return testDict +@pytest.mark.base +def testExactDummyProblems(): + r""" + Test for dummy problems. The test verifies that the dummy problems exactly returns the dynamics of + the parent class. ``eval_f`` of ``ExactDiscontinuousTestDAE`` is not tested here, since it only returns + a random right-hand side to enforce the sweeper to do not stop to compute. + """ + + from pySDC.implementations.datatype_classes.mesh import mesh + + childODE = ExactDiscontinuousTestODE(**{}) + parentODE = DiscontinuousTestODE(**{}) + assert childODE.t_switch_exact == parentODE.t_switch_exact, f"Exact event times between classes does not match!" + + t0 = 1.0 + dt = 0.1 + u0 = parentODE.u_exact(t0) + rhs = u0.copy() + + uSolve = childODE.solve_system(rhs, dt, u0, t0) + uExact = parentODE.u_exact(t0) + assert np.allclose(uSolve, uExact) + + # same test for event time + tExactEventODE = parentODE.t_switch_exact + dt = 0.1 + u0Event = parentODE.u_exact(tExactEventODE) + rhsEvent = u0.copy() + + uSolveEvent = childODE.solve_system(rhsEvent, dt, u0Event, tExactEventODE) + uExactEvent = parentODE.u_exact(tExactEventODE) + assert np.allclose(uSolveEvent, uExactEvent) + + fExactOde = childODE.eval_f(u0, t0) + fOde = parentODE.eval_f(u0, t0) + assert np.allclose(fExactOde, fOde), f"Right-hand sides do not match!" + + fExactOdeEvent = childODE.eval_f(u0Event, tExactEventODE) + fOdeEvent = parentODE.eval_f(u0Event, tExactEventODE) + assert np.allclose(fExactOdeEvent, fOdeEvent), f"Right-hand sides at event do not match!" + + childDAE = ExactDiscontinuousTestDAE(**{}) + parentDAE = DiscontinuousTestDAE(**{}) + assert childDAE.t_switch_exact == parentDAE.t_switch_exact, f"Exact event times between classes does not match!" + + t0 = 0.0 + u0 = mesh(childDAE.init) + duExactDAE = childDAE.solve_system(impl_sys=None, u0=u0, t=t0) + duDAE = mesh(parentDAE.init) + duDAE[:] = (np.sinh(t0), np.cosh(t0)) + assert np.allclose(duExactDAE, duDAE), f"Method solve_system of dummy problem do not return exact derivative!" + + tExactEventDAE = parentDAE.t_switch_exact + u0Event = mesh(childDAE.init) + duExactDAEEvent = childDAE.solve_system(impl_sys=None, u0=u0Event, t=tExactEventDAE) + duDAEEvent = mesh(parentDAE.init) + duDAEEvent[:] = (np.sinh(tExactEventDAE), np.cosh(tExactEventDAE)) + assert np.allclose( + duExactDAEEvent, duDAEEvent + ), f"Method solve_system of dummy problem do not return exact derivative at event!" + + @pytest.mark.base @pytest.mark.parametrize("key", [0, 1]) def testInterpolationValues(key): @@ -317,7 +248,7 @@ def testInterpolationValues(key): @pytest.mark.base @pytest.mark.parametrize('quad_type', ['LOBATTO', 'RADAU-RIGHT']) -def test_adapt_interpolation_info(quad_type): +def testAdaptInterpolationInfo(quad_type): r""" Tests if the method ``adapt_interpolation_info`` does what it is supposed to do. @@ -334,25 +265,47 @@ def test_adapt_interpolation_info(quad_type): Type of quadrature used. """ + from pySDC.projects.PinTSimE.battery_model import generateDescription from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + problem = ExactDiscontinuousTestODE + problem_params = dict() t0 = 1.6 - Tend = ExactDiscontinuousTestODE().t_switch_exact + Tend = problem(**problem_params).t_switch_exact eps = 1e-13 # choose this eps to enforce a sign chance in state function dt = (Tend - t0) + eps - tol = 1e-10 + sweeper = generic_implicit num_nodes = 3 + QI = 'IE' - controller = get_controller( - switch_estimator=SwitchEstimator, - problem=ExactDiscontinuousTestODE, - sweeper=generic_implicit, - quad_type=quad_type, - num_nodes=num_nodes, + restol = 1e-13 + tol = 1e-10 + alpha = 0.95 + maxiter = 1 + max_restarts = 3 + useA = False + useSE = True + + hook_class = [] + + _, _, controller = generateDescription( dt=dt, - tol=tol, + problem=problem, + sweeper=sweeper, + num_nodes=num_nodes, + quad_type=quad_type, + QI=QI, + hook_class=hook_class, + use_adaptivity=useA, + use_switch_estimator=useSE, + problem_params=problem_params, + restol=restol, + maxiter=maxiter, + max_restarts=max_restarts, + tol_event=tol, + alpha=alpha, ) S = controller.MS[0] @@ -399,7 +352,7 @@ def test_adapt_interpolation_info(quad_type): @pytest.mark.base @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -def test_detection_at_boundary(num_nodes): +def testDetectionBoundary(num_nodes): """ This test checks whether a restart is executed or not when the event exactly occurs at the boundary. In this case, no restart should be done because occuring the event at the boundary means that the event is already resolved well, @@ -411,25 +364,58 @@ def test_detection_at_boundary(num_nodes): Number of collocation nodes. """ + from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit - from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.implementations.hooks.log_restarts import LogRestarts from pySDC.helpers.stats_helper import get_sorted problem = ExactDiscontinuousTestODE + problem_params = dict() t0 = 1.6 - Tend = ExactDiscontinuousTestODE().t_switch_exact + Tend = problem(**problem_params).t_switch_exact dt = Tend - t0 - stats, _ = discontinuousTestProblem_run( - switch_estimator=SwitchEstimator, + sweeper = generic_implicit + QI = 'IE' + quad_type = 'LOBATTO' + + restol = 1e-13 + tol = 1e-10 + alpha = 0.95 + maxiter = 1 + max_restarts = 3 + useA = False + useSE = True + exact_event_time_avail = True + + hook_class = [LogSolution, LogRestarts] + + description, controller_params, controller = generateDescription( + dt=dt, problem=problem, - sweeper=generic_implicit, - quad_type='LOBATTO', + sweeper=sweeper, num_nodes=num_nodes, + quad_type=quad_type, + QI=QI, + hook_class=hook_class, + use_adaptivity=useA, + use_switch_estimator=useSE, + problem_params=problem_params, + restol=restol, + maxiter=maxiter, + max_restarts=max_restarts, + tol_event=tol, + alpha=alpha, + ) + + stats, _ = controllerRun( + description=description, + controller_params=controller_params, + controller=controller, t0=t0, Tend=Tend, - dt=dt, - tol=1e-10, + exact_event_time_avail=exact_event_time_avail, ) sum_restarts = np.sum(np.array(get_sorted(stats, type='restart', sortby='time', recomputed=None))[:, 1]) @@ -440,7 +426,7 @@ def test_detection_at_boundary(num_nodes): @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) @pytest.mark.parametrize('num_nodes', [3, 4, 5]) @pytest.mark.parametrize('quad_type', ['LOBATTO', 'RADAU-RIGHT']) -def test_all_tolerances_ODE(tol, num_nodes, quad_type): +def testDetectionODE(tol, num_nodes, quad_type): r""" Here, the switch estimator is applied to a dummy problem of ``DiscontinuousTestODE``, where the dynamics of the differential equation is replaced by its exact dynamics to see if @@ -461,23 +447,55 @@ def test_all_tolerances_ODE(tol, num_nodes, quad_type): """ from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit - from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.helpers.stats_helper import get_sorted + from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.implementations.hooks.log_restarts import LogRestarts + problem = ExactDiscontinuousTestODE + problem_params = dict() t0 = 1.6 Tend = 1.62 - problem = ExactDiscontinuousTestODE + dt = Tend - t0 - stats, t_switch_exact = discontinuousTestProblem_run( - switch_estimator=SwitchEstimator, + sweeper = generic_implicit + QI = 'IE' + + restol = 1e-13 + alpha = 0.95 + maxiter = 1 + max_restarts = 3 + useA = False + useSE = True + exact_event_time_avail = True + + hook_class = [LogSolution, LogRestarts] + + description, controller_params, controller = generateDescription( + dt=dt, problem=problem, - sweeper=generic_implicit, - quad_type=quad_type, + sweeper=sweeper, num_nodes=num_nodes, + quad_type=quad_type, + QI=QI, + hook_class=hook_class, + use_adaptivity=useA, + use_switch_estimator=useSE, + problem_params=problem_params, + restol=restol, + maxiter=maxiter, + max_restarts=max_restarts, + tol_event=tol, + alpha=alpha, + ) + + stats, t_switch_exact = controllerRun( + description=description, + controller_params=controller_params, + controller=controller, t0=t0, Tend=Tend, - dt=2e-2, - tol=tol, + exact_event_time_avail=exact_event_time_avail, ) # in this specific example only one event has to be found @@ -492,7 +510,7 @@ def test_all_tolerances_ODE(tol, num_nodes, quad_type): @pytest.mark.base @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -def test_all_tolerances_DAE(tol, num_nodes): +def testDetectionDAE(tol, num_nodes): r""" In this test, the switch estimator is applied to a DAE dummy problem of ``DiscontinuousTestDAE``, where the dynamics of the differential equation is replaced by its exact dynamics to see if @@ -513,23 +531,56 @@ def test_all_tolerances_DAE(tol, num_nodes): """ from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE - from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.helpers.stats_helper import get_sorted + from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.implementations.hooks.log_restarts import LogRestarts + problem = ExactDiscontinuousTestDAE + problem_params = dict() t0 = 4.6 Tend = 4.62 - problem = ExactDiscontinuousTestDAE + dt = Tend - t0 - stats, t_switch_exact = discontinuousTestProblem_run( - switch_estimator=SwitchEstimator, + sweeper = fully_implicit_DAE + QI = 'IE' + quad_type = 'RADAU-RIGHT' + + restol = 1e-13 + alpha = 0.95 + maxiter = 1 + max_restarts = 3 + useA = False + useSE = True + exact_event_time_avail = True + + hook_class = [LogSolution, LogRestarts] + + description, controller_params, controller = generateDescription( + dt=dt, problem=problem, - sweeper=fully_implicit_DAE, - quad_type='RADAU-RIGHT', + sweeper=sweeper, num_nodes=num_nodes, + quad_type=quad_type, + QI=QI, + hook_class=hook_class, + use_adaptivity=useA, + use_switch_estimator=useSE, + problem_params=problem_params, + restol=restol, + maxiter=maxiter, + max_restarts=max_restarts, + tol_event=tol, + alpha=alpha, + ) + + stats, t_switch_exact = controllerRun( + description=description, + controller_params=controller_params, + controller=controller, t0=t0, Tend=Tend, - dt=2e-2, - tol=tol, + exact_event_time_avail=exact_event_time_avail, ) # in this specific example only one event has to be found From f4a5294df38d53f58234c6e0516190fd4168f610 Mon Sep 17 00:00:00 2001 From: lisawim Date: Sun, 25 Feb 2024 09:27:29 +0100 Subject: [PATCH 13/29] Moved standard params for test to function --- .../test_pintsime/test_SwitchEstimator.py | 47 ++++++++----------- 1 file changed, 20 insertions(+), 27 deletions(-) diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 954198dd01..b3f6ecde50 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -159,6 +159,20 @@ def getTestDict(): return testDict +def getParamsRun(): + r""" + Returns parameters for conroller run that are used in each test. + """ + restol = 1e-13 + alpha = 0.95 + maxiter = 1 + max_restarts = 3 + useA = False + useSE = True + exact_event_time_avail = True + return restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail + + @pytest.mark.base def testExactDummyProblems(): r""" @@ -280,13 +294,9 @@ def testAdaptInterpolationInfo(quad_type): num_nodes = 3 QI = 'IE' - restol = 1e-13 tol = 1e-10 - alpha = 0.95 - maxiter = 1 - max_restarts = 3 - useA = False - useSE = True + + restol, alpha, maxiter, max_restarts, useA, useSE, _ = getParamsRun() hook_class = [] @@ -380,14 +390,9 @@ def testDetectionBoundary(num_nodes): QI = 'IE' quad_type = 'LOBATTO' - restol = 1e-13 tol = 1e-10 - alpha = 0.95 - maxiter = 1 - max_restarts = 3 - useA = False - useSE = True - exact_event_time_avail = True + + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -461,13 +466,7 @@ def testDetectionODE(tol, num_nodes, quad_type): sweeper = generic_implicit QI = 'IE' - restol = 1e-13 - alpha = 0.95 - maxiter = 1 - max_restarts = 3 - useA = False - useSE = True - exact_event_time_avail = True + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -546,13 +545,7 @@ def testDetectionDAE(tol, num_nodes): QI = 'IE' quad_type = 'RADAU-RIGHT' - restol = 1e-13 - alpha = 0.95 - maxiter = 1 - max_restarts = 3 - useA = False - useSE = True - exact_event_time_avail = True + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() hook_class = [LogSolution, LogRestarts] From 019e78b94db2b88da0e10f39e453156606aad83a Mon Sep 17 00:00:00 2001 From: lisawim Date: Sun, 25 Feb 2024 12:38:19 +0100 Subject: [PATCH 14/29] Updated hardcoded solutions for battery models --- pySDC/projects/PinTSimE/battery_model.py | 8 +- .../projects/PinTSimE/hardcoded_solutions.py | 148 +++++++++--------- pySDC/projects/PinTSimE/switch_estimator.py | 2 +- .../test_pintsime/test_SwitchEstimator.py | 2 +- 4 files changed, 80 insertions(+), 80 deletions(-) diff --git a/pySDC/projects/PinTSimE/battery_model.py b/pySDC/projects/PinTSimE/battery_model.py index 7301aceaf4..7b6b860677 100644 --- a/pySDC/projects/PinTSimE/battery_model.py +++ b/pySDC/projects/PinTSimE/battery_model.py @@ -106,7 +106,7 @@ def generateDescription( controller_params : dict Parameters needed for a controller run. """ - print(use_adaptivity) + # initialize level parameters level_params = { 'restol': -1 if use_adaptivity else restol, @@ -235,7 +235,7 @@ def main(): 'max_restarts': 50, 'recomputed': False, 'tol_event': 1e-10, - 'alpha': 1.0, + 'alpha': 0.95, 'exact_event_time_avail': None, } @@ -246,8 +246,8 @@ def main(): hook_class = [LogSolution, LogEventBattery, LogEmbeddedErrorEstimate, LogStepSize] - use_detection = [True] - use_adaptivity = [True] + use_detection = [True, False] + use_adaptivity = [True, False] for problem, sweeper in zip([battery, battery_implicit], [imex_1st_order, generic_implicit]): for defaults in [False, True]: diff --git a/pySDC/projects/PinTSimE/hardcoded_solutions.py b/pySDC/projects/PinTSimE/hardcoded_solutions.py index 27733a12ec..f23ed8fb68 100644 --- a/pySDC/projects/PinTSimE/hardcoded_solutions.py +++ b/pySDC/projects/PinTSimE/hardcoded_solutions.py @@ -43,23 +43,23 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator and adaptivity for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.5614559718189012, - 'vC': 1.0053361988800296, - 't_switches': [0.18232155679214296], - 'dt': 0.11767844320785703, - 'e_em': 7.811640223565064e-12, - 'sum_restarts': 3.0, - 'sum_niters': 56.0, + 'iL': 0.5393867578149492, + 'vC': 0.9999999999980306, + 't_switches': [0.18232155679314516, 0.18232155674931952], + 'dt': 0.017755006862796352, + 'e_em': 2.220446049250313e-16, + 'sum_restarts': 17.0, + 'sum_niters': 312.0, } elif dt == 1e-3: expected = { - 'iL': 0.5393867578949986, - 'vC': 1.0000000000165197, - 't_switches': [0.18232155677793654], - 'dt': 0.015641173481932502, - 'e_em': 2.220446049250313e-16, - 'sum_restarts': 14.0, - 'sum_niters': 328.0, + 'iL': 0.5393867577908672, + 'vC': 0.9999999999869332, + 't_switches': [0.182321556754745], + 'dt': 0.056647727943027276, + 'e_em': 4.285460875053104e-14, + 'sum_restarts': 15.0, + 'sum_niters': 296.0, } got.update( { @@ -75,19 +75,19 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.549122133626298, - 'vC': 0.9999999999999998, - 't_switches': [0.1823215567939562], - 'sum_restarts': 4.0, - 'sum_niters': 296.0, + 'iL': 0.5490992955276419, + 'vC': 0.9999999999816187, + 't_switches': [0.182321556761159], + 'sum_restarts': 13.0, + 'sum_niters': 416.0, } elif dt == 1e-3: expected = { - 'iL': 0.5408462989990014, - 'vC': 1.0, - 't_switches': [0.18232155679395023], - 'sum_restarts': 2.0, - 'sum_niters': 2424.0, + 'iL': 0.5403849766797957, + 'vC': 0.9999166666718436, + 't_switches': [0.18232155679295683, 0.1823215567608114], + 'sum_restarts': 11.0, + 'sum_niters': 2536.0, } got.update( { @@ -151,23 +151,23 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator and adaptivity for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.5614559717904407, - 'vC': 1.0053361988803866, - 't_switches': [0.18232155679736195], - 'dt': 0.11767844320263804, + 'iL': 0.5393867577386883, + 'vC': 1.0000000000016669, + 't_switches': [0.18232155680505413, 0.1823215568065723], + 'dt': 0.017701859620598615, 'e_em': 2.220446049250313e-16, - 'sum_restarts': 3.0, - 'sum_niters': 56.0, + 'sum_restarts': 17.0, + 'sum_niters': 312.0, } elif dt == 1e-3: expected = { - 'iL': 0.5393867577837699, - 'vC': 1.0000000000250129, - 't_switches': [0.1823215568036829], - 'dt': 0.015641237833012522, + 'iL': 0.5393867577174891, + 'vC': 0.9999999999869332, + 't_switches': [0.18232155677090756], + 'dt': 0.05664772792686473, 'e_em': 2.220446049250313e-16, - 'sum_restarts': 14.0, - 'sum_niters': 328.0, + 'sum_restarts': 15.0, + 'sum_niters': 296.0, } got.update( { @@ -183,19 +183,19 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.5456190026227917, - 'vC': 0.999166666666676, - 't_switches': [0.18232155679663525], - 'sum_restarts': 4.0, - 'sum_niters': 296.0, + 'iL': 0.549099295497225, + 'vC': 0.9999999999816187, + 't_switches': [0.18232155676455783], + 'sum_restarts': 13.0, + 'sum_niters': 416.0, } elif dt == 1e-3: expected = { - 'iL': 0.5407340515794409, - 'vC': 1.0000000000010945, - 't_switches': [0.182321556796257], - 'sum_restarts': 3.0, - 'sum_niters': 2440.0, + 'iL': 0.5403849766692057, + 'vC': 0.9999166666718439, + 't_switches': [0.18232155679484557, 0.1823215567627017], + 'sum_restarts': 11.0, + 'sum_niters': 2536.0, } got.update( { @@ -259,25 +259,25 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator and adaptivity for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.6244130166029733, - 'vC1': 0.999647921822499, - 'vC2': 1.0000000000714673, - 't_switches': [0.18232155679216916, 0.3649951297770592], - 'dt': 0.01, - 'e_em': 2.220446049250313e-16, - 'sum_restarts': 19.0, - 'sum_niters': 400.0, + 'iL': 0.6125019898490534, + 'vC1': 1.000000000007182, + 'vC2': 1.0000000000479325, + 't_switches': [0.18232155680297893, 0.36464311356978835], + 'dt': 0.02973949262636494, + 'e_em': 3.3306690738754696e-16, + 'sum_restarts': 24.0, + 'sum_niters': 496.0, } elif dt == 1e-3: expected = { - 'iL': 0.6112496171462107, - 'vC1': 0.9996894956748836, - 'vC2': 1.0, - 't_switches': [0.1823215567907929, 0.3649535697059346], - 'dt': 0.07298158272977251, - 'e_em': 2.703393064962256e-13, - 'sum_restarts': 11.0, - 'sum_niters': 216.0, + 'iL': 0.6125019898915002, + 'vC1': 1.0000000000057128, + 'vC2': 1.000000000050465, + 't_switches': [0.1823215567993942, 0.36464311356557344], + 'dt': 0.0241594380776875, + 'e_em': 2.220446049250313e-16, + 'sum_restarts': 24.0, + 'sum_niters': 504.0, } got.update( { @@ -293,21 +293,21 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.6314080101219072, - 'vC1': 0.9999999999999998, - 'vC2': 0.9999999999999996, - 't_switches': [0.1823215567939562, 0.3646431135879125], - 'sum_restarts': 8.0, - 'sum_niters': 512.0, + 'iL': 0.6313858468881018, + 'vC1': 1.0000000000155491, + 'vC2': 1.0000000000316576, + 't_switches': [0.18232155681737047, 0.36464311357240825], + 'sum_restarts': 20.0, + 'sum_niters': 680.0, } elif dt == 1e-3: expected = { - 'iL': 0.6152346866530549, - 'vC1': 1.0, - 'vC2': 1.0, - 't_switches': [0.18232155679395023, 0.3646431135879003], - 'sum_restarts': 4.0, - 'sum_niters': 4048.0, + 'iL': 0.6151254297607596, + 'vC1': 1.0000000000846776, + 'vC2': 1.0000000000993976, + 't_switches': [0.18232155679295683, 0.3646431135032376], + 'sum_restarts': 17.0, + 'sum_niters': 4232.0, } got.update( { diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 6765699114..1e080a50b2 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -401,7 +401,7 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): msg = f'Newton did not converge after {n} iterations, error is {res}' else: msg = f'Newton did converge after {n} iterations, error for root {x0} is {res}' - print(msg) + # print(msg) root = x0 return root diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index b3f6ecde50..9c5937bade 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -77,7 +77,7 @@ def __init__(self, newton_tol=1e-8): def eval_f(self, u, du, t): r""" Returns the exact right-hand side of the problem. The first components in both - cases of ``f`` are set to 1 do avoid getting a zero residual (enforced the sweeper to stop + cases of ``f`` are set to 1 to avoid getting a zero residual (enforced the sweeper to stop since "convergence is reached"). Parameters From 4bd86c43b25ad4deeab1ea55c8784fe48374c466 Mon Sep 17 00:00:00 2001 From: lisawim Date: Sun, 25 Feb 2024 18:37:05 +0100 Subject: [PATCH 15/29] Updated hardcoded solutions for DiscontinuousTestODE --- pySDC/projects/PinTSimE/battery_model.py | 8 ++++++ pySDC/projects/PinTSimE/buck_model.py | 1 + .../PinTSimE/discontinuous_test_ODE.py | 3 +- pySDC/projects/PinTSimE/estimation_check.py | 16 ++++++++--- .../projects/PinTSimE/hardcoded_solutions.py | 28 +++++++++---------- pySDC/projects/PinTSimE/piline_model.py | 1 + pySDC/projects/PinTSimE/switch_estimator.py | 23 +++++++++------ .../test_pintsime/test_SwitchEstimator.py | 15 ++++++---- 8 files changed, 62 insertions(+), 33 deletions(-) diff --git a/pySDC/projects/PinTSimE/battery_model.py b/pySDC/projects/PinTSimE/battery_model.py index 7b6b860677..0cfac58156 100644 --- a/pySDC/projects/PinTSimE/battery_model.py +++ b/pySDC/projects/PinTSimE/battery_model.py @@ -61,6 +61,7 @@ def generateDescription( max_restarts=None, tol_event=1e-10, alpha=1.0, + typeFD='backward', ): r""" Generate a description for the battery models for a controller run. @@ -98,6 +99,9 @@ def generateDescription( Tolerance for event detection to terminate. alpha : float, optional Factor that indicates how the new step size in the Switch Estimator is reduced. + typeFD : str, optional + Type of finite difference in switch estimator. For detailed documentation, see + switch_estimator.py. Returns ------- @@ -139,9 +143,11 @@ def generateDescription( # convergence controllers convergence_controllers = {} if use_switch_estimator: + typeFD = typeFD if typeFD is not None else 'backward' switch_estimator_params = { 'tol': tol_event, 'alpha': alpha, + 'typeFD': typeFD, } convergence_controllers.update({SwitchEstimator: switch_estimator_params}) if use_adaptivity: @@ -237,6 +243,7 @@ def main(): 'tol_event': 1e-10, 'alpha': 0.95, 'exact_event_time_avail': None, + 'typeFD': 'backward', } all_params = { @@ -378,6 +385,7 @@ def runSimulation(problem, sweeper, all_params, use_adaptivity, use_detection, h max_restarts=handling_params['max_restarts'], tol_event=handling_params['tol_event'], alpha=handling_params['alpha'], + typeFD=handling_params['typeFD'], ) stats, t_switch_exact = controllerRun( diff --git a/pySDC/projects/PinTSimE/buck_model.py b/pySDC/projects/PinTSimE/buck_model.py index 3a2423d339..4b1e3a27dc 100644 --- a/pySDC/projects/PinTSimE/buck_model.py +++ b/pySDC/projects/PinTSimE/buck_model.py @@ -34,6 +34,7 @@ def main(): 'tol_event': None, 'alpha': None, 'exact_event_time_avail': None, + 'typeFD': None, } # ---- all parameters are stored in this dictionary - note: defaults are used for the problem ---- diff --git a/pySDC/projects/PinTSimE/discontinuous_test_ODE.py b/pySDC/projects/PinTSimE/discontinuous_test_ODE.py index 46bee5c4d9..a264536444 100644 --- a/pySDC/projects/PinTSimE/discontinuous_test_ODE.py +++ b/pySDC/projects/PinTSimE/discontinuous_test_ODE.py @@ -68,8 +68,9 @@ def main(): 'max_restarts': 50, 'recomputed': False, 'tol_event': 1e-12, - 'alpha': 1.0, + 'alpha': 0.95, 'exact_event_time_avail': True, + 'typeFD': 'backward', } # ---- all parameters are stored in this dictionary ---- diff --git a/pySDC/projects/PinTSimE/estimation_check.py b/pySDC/projects/PinTSimE/estimation_check.py index 61fb858584..d0fdb02655 100644 --- a/pySDC/projects/PinTSimE/estimation_check.py +++ b/pySDC/projects/PinTSimE/estimation_check.py @@ -45,8 +45,9 @@ def run_estimation_check(): 'max_restarts': 50, 'recomputed': False, 'tol_event': 1e-10, - 'alpha': 1.0, + 'alpha': 0.95, 'exact_event_time_avail': None, + 'typeFD': 'backward', } problem_classes = [battery, battery_implicit, battery_n_capacitors] @@ -114,7 +115,7 @@ def run_estimation_check(): plotAccuracyCheck(u_num, prob_cls_name, M_fix) - plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix) + # plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix) plotStateFunctionOverTime(u_num, prob_cls_name, M_fix) @@ -187,6 +188,9 @@ def plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix): # pragma: no cov Routine that plots the state function at time before the event, exactly at the event, and after the event. Note that this routine does make sense only for a state function that remains constant after the event. + TODO: Function still does not work as expected. Every time when the switch estimator is adapted, the tolerances + does not suit anymore! + Parameters ---------- u_num : dict @@ -239,15 +243,19 @@ def plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix): # pragma: no cov if use_SE: t_switches = [u_num[dt][M_fix][use_SE][use_A]['t_switches'] for dt in dt_list] - t_switch = [t_event[i] for t_event in t_switches] + for t_switch_item in t_switches: + for m in range(1, len(t_switch_item)): + if np.isclose(t_switch_item[m - 1], t_switch_item[m], atol=1e-10): + t_switch_item.pop(m - 1) + t_switch = [t_event[i] for t_event in t_switches] ax[0, ind].plot( dt_list, [ h_item[m] for (t_item, h_item, t_switch_item) in zip(t, h, t_switch) for m in range(len(t_item)) - if abs(t_item[m] - t_switch_item) <= 1e-14 + if abs(t_item[m] - t_switch_item) <= 2.7961188919789493e-11 ], color='limegreen', marker='s', diff --git a/pySDC/projects/PinTSimE/hardcoded_solutions.py b/pySDC/projects/PinTSimE/hardcoded_solutions.py index f23ed8fb68..76359d2bad 100644 --- a/pySDC/projects/PinTSimE/hardcoded_solutions.py +++ b/pySDC/projects/PinTSimE/hardcoded_solutions.py @@ -375,21 +375,21 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'u': 5.9941358952954955, - 't_switches': [1.6094379124671208], - 'sum_restarts': 25.0, - 'sum_niters': 710.0, - 'e_global': 8.195133460731086e-11, - 'e_event': [3.302047524300633e-11], + 'u': 5.994135895337547, + 'e_global': 7.765876830490015e-11, + 't_switches': [1.6094379124455047], + 'e_event': [1.1404432953554533e-11], + 'sum_restarts': 22.0, + 'sum_niters': 686.0, } elif dt == 1e-3: expected = { - 'u': 5.971767837651004, - 't_switches': [1.6094379124247695], - 'sum_restarts': 23.0, - 'sum_niters': 2388.0, - 'e_global': 2.3067769916451653e-11, - 'e_event': [9.330758388159666e-12], + 'u': 5.972186947618947, + 'e_global': 0.00041911003080929987, + 't_switches': [1.6094379124055875], + 'e_event': [2.851274771842327e-11], + 'sum_restarts': 17.0, + 'sum_niters': 2338.0, } got.update( { @@ -406,14 +406,14 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): if dt == 1e-2: expected = { 'u': 5.9805345175338225, - 'sum_niters': 527.0, 'e_global': 0.009855041056925806, + 'sum_niters': 527.0, } elif dt == 1e-3: expected = { 'u': 5.9737411566014105, - 'sum_niters': 2226.0, 'e_global': 0.0005763403865515215, + 'sum_niters': 2226.0, } got.update( { diff --git a/pySDC/projects/PinTSimE/piline_model.py b/pySDC/projects/PinTSimE/piline_model.py index fc2f5f56f6..be135eb0ed 100644 --- a/pySDC/projects/PinTSimE/piline_model.py +++ b/pySDC/projects/PinTSimE/piline_model.py @@ -34,6 +34,7 @@ def main(): 'tol_event': None, 'alpha': None, 'exact_event_time_avail': None, + 'typeFD': None, } # ---- all parameters are stored in this dictionary - note: defaults are used for the problem ---- diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 1e080a50b2..ab33711e6b 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -122,7 +122,9 @@ def get_new_step_size(self, controller, S, **kwargs): # intermediate value theorem states that a root is contained in current step if self.params.state_function[0] * self.params.state_function[-1] < 0 and self.status.is_zero is None: - self.status.t_switch = self.get_switch(self.params.t_interp, self.params.state_function, m_guess) + self.status.t_switch = self.get_switch( + self.params.t_interp, self.params.state_function, m_guess, self.params.typeFD + ) self.logging_during_estimation( controller.hooks[0], @@ -274,8 +276,8 @@ def logging_during_estimation(controller_hooks, process, time, level, sweep, t_s ) @staticmethod - def get_switch(t_interp, state_function, m_guess): - """ + def get_switch(t_interp, state_function, m_guess, choiceFD='backward'): + r""" Routine to do the interpolation and root finding stuff. Parameters @@ -286,6 +288,10 @@ def get_switch(t_interp, state_function, m_guess): Contains values of state function at these collocation nodes. m_guess : float Index at which the difference drops below zero. + choiceFD : str + Finite difference to be used to approximate derivative of + state function. Can be ``'forward'``, ``'centered'``, or + ``'backward'``. Default is ``'backward'`` Returns ------- @@ -294,10 +300,6 @@ def get_switch(t_interp, state_function, m_guess): """ LinearInterpolator = LinearInterpolation(t_interp, state_function) p = lambda t: LinearInterpolator.eval(t) - if state_function[0] < 0 and state_function[-1] >= 0: - chooseFD = 'centered' - else: - chooseFD = 'backward' def fprime(t): """ @@ -314,10 +316,13 @@ def fprime(t): Derivative of interpolation p at time t. """ - if chooseFD == 'centered': + if choiceFD == 'forward': + dt_FD = 1e-10 + dp = (p(t + dt_FD) - p(t)) / (dt_FD) + elif choiceFD == 'centered': dt_FD = 1e-12 dp = (p(t + dt_FD) - p(t - dt_FD)) / (2 * dt_FD) - elif chooseFD == 'backward': + elif choiceFD == 'backward': dt_FD = 1e-10 dp = (p(t) - p(t - dt_FD)) / (dt_FD) else: diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 9c5937bade..79ba9b5dad 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -170,7 +170,8 @@ def getParamsRun(): useA = False useSE = True exact_event_time_avail = True - return restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail + typeFD = 'centered' + return restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD @pytest.mark.base @@ -296,7 +297,7 @@ def testAdaptInterpolationInfo(quad_type): tol = 1e-10 - restol, alpha, maxiter, max_restarts, useA, useSE, _ = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, _, typeFD = getParamsRun() hook_class = [] @@ -316,6 +317,7 @@ def testAdaptInterpolationInfo(quad_type): max_restarts=max_restarts, tol_event=tol, alpha=alpha, + typeFD=typeFD, ) S = controller.MS[0] @@ -392,7 +394,7 @@ def testDetectionBoundary(num_nodes): tol = 1e-10 - restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -412,6 +414,7 @@ def testDetectionBoundary(num_nodes): max_restarts=max_restarts, tol_event=tol, alpha=alpha, + typeFD=typeFD, ) stats, _ = controllerRun( @@ -466,7 +469,7 @@ def testDetectionODE(tol, num_nodes, quad_type): sweeper = generic_implicit QI = 'IE' - restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -486,6 +489,7 @@ def testDetectionODE(tol, num_nodes, quad_type): max_restarts=max_restarts, tol_event=tol, alpha=alpha, + typeFD=typeFD, ) stats, t_switch_exact = controllerRun( @@ -545,7 +549,7 @@ def testDetectionDAE(tol, num_nodes): QI = 'IE' quad_type = 'RADAU-RIGHT' - restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -565,6 +569,7 @@ def testDetectionDAE(tol, num_nodes): max_restarts=max_restarts, tol_event=tol, alpha=alpha, + typeFD=typeFD, ) stats, t_switch_exact = controllerRun( From 487573ab6433600f1d93f5abf3e8189d8a3804d6 Mon Sep 17 00:00:00 2001 From: lisawim Date: Sun, 25 Feb 2024 18:43:51 +0100 Subject: [PATCH 16/29] Updated docu in SE for FDs --- pySDC/projects/PinTSimE/switch_estimator.py | 29 +++++++++++++++++---- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index ab33711e6b..3365d3e93e 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -17,9 +17,11 @@ def setup(self, controller, params, description): r""" Function sets default variables to handle with the event at the beginning. The default params are: - - control_order : controls the order of the SE's call of convergence controllers - - coll.nodes : defines the collocation nodes for interpolation - - tol_zero : inner tolerance for SE; state function has to satisfy it to terminate + - control_order : controls the order of the SE's call of convergence controllers. + - coll.nodes : defines the collocation nodes for interpolation. + - tol_zero : inner tolerance for SE; state function has to satisfy it to terminate. + - t_interp : interpolation axis with time points. + - state_function : List of values from state function. Parameters ---------- @@ -302,8 +304,25 @@ def get_switch(t_interp, state_function, m_guess, choiceFD='backward'): p = lambda t: LinearInterpolator.eval(t) def fprime(t): - """ - Computes the derivative of the scalar interpolant using centered finite difference. + r""" + Computes the derivative of the scalar interpolant using finite difference. + Here, different finite differences can be used. The type of FD can be set by + setting ``typeFD`` in switch estimator parameters. There are three choices possible: + + - ``typeFD='backward'`` for :math:`h=10^{-10}`: + + .. math:: + \frac{dp}{dt} \approx \frac{p(t) - p(t - h)}{h} + + - ``typeFD='centered'`` for :math:`h=10^{-12}`: + + .. math:: + \frac{dp}{dt} \approx \frac{p(t + h) - p(t - h)}{2h} + + - ``typeFD='forward'`` for :math:`h=10^{-10}`: + + .. math:: + \frac{dp}{dt} \approx \frac{p(t + h) - p(t)}{h} Parameters ---------- From 28425dae5313b4b9bae713bee142d8e877def2ac Mon Sep 17 00:00:00 2001 From: lisawim Date: Mon, 26 Feb 2024 13:15:46 +0100 Subject: [PATCH 17/29] Lagrange Interpolation works better with baclward FD and alpha=0.9 --- .../PinTSimE/paper_PSCC2024/log_event.py | 4 +- pySDC/projects/PinTSimE/switch_estimator.py | 51 +++--- .../test_pintsime/test_SwitchEstimator.py | 170 +++--------------- 3 files changed, 54 insertions(+), 171 deletions(-) diff --git a/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py b/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py index 6db702465a..8a372babc4 100644 --- a/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py +++ b/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py @@ -8,7 +8,7 @@ class LogEventDiscontinuousTestDAE(hooks): """ def post_step(self, step, level_number): - super(LogEventDiscontinuousTestDAE, self).post_step(step, level_number) + super().post_step(step, level_number) L = step.levels[level_number] @@ -32,7 +32,7 @@ class LogEventWSCC9(hooks): """ def post_step(self, step, level_number): - super(LogEventWSCC9, self).post_step(step, level_number) + super().post_step(step, level_number) L = step.levels[level_number] P = L.prob diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 3365d3e93e..51b96613b9 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -300,7 +300,7 @@ def get_switch(t_interp, state_function, m_guess, choiceFD='backward'): t_switch : float Time point of found event. """ - LinearInterpolator = LinearInterpolation(t_interp, state_function) + LinearInterpolator = LagrangeInterpolation(t_interp, state_function) p = lambda t: LinearInterpolator.eval(t) def fprime(t): @@ -342,13 +342,15 @@ def fprime(t): dt_FD = 1e-12 dp = (p(t + dt_FD) - p(t - dt_FD)) / (2 * dt_FD) elif choiceFD == 'backward': - dt_FD = 1e-10 - dp = (p(t) - p(t - dt_FD)) / (dt_FD) + dt_FD = 1e-12 + # dp = (p(t) - p(t - dt_FD)) / (dt_FD) + # dp = (11 * p(t) - 18 * p(t - dt_FD) + 9 * p(t - 2 * dt_FD) - 2 * p(t - 3 * dt_FD)) / (6 * dt_FD) + dp = (25 * p(t) - 48 * p(t - dt_FD) + 36 * p(t - 2 * dt_FD) - 16 * p(t - 3 * dt_FD) + 3 * p(t - 4 * dt_FD)) / (12 * dt_FD) else: raise NotImplementedError return dp - newton_tol, newton_maxiter = 1e-12, 100 + newton_tol, newton_maxiter = 1e-15, 100 t_switch = newton(t_interp[m_guess], p, fprime, newton_tol, newton_maxiter) return t_switch @@ -431,36 +433,45 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): return root -class LinearInterpolation(object): +class LagrangeInterpolation(object): def __init__(self, ti, yi): """Initialization routine""" self.ti = np.asarray(ti) self.yi = np.asarray(yi) self.n = len(ti) + def get_Lagrange_polynomial(self, t, i): + """ + Computes the basis of the i-th Lagrange polynomial. + + Parameters + ---------- + t : float + Time where the polynomial is computed at. + i : int + Index of the Lagrange polynomial + + Returns + ------- + product : float + The product of the bases. + """ + product = np.prod([(t - self.ti[k]) / (self.ti[i] - self.ti[k]) for k in range(self.n) if k != i]) + return product + def eval(self, t): - r""" - Evaluates the linear interpolant at time :math:`t`. - Binary search is used to find the index in which piece :math:`t` lies. + """ + Evaluates the Lagrange interpolation at time t. Parameters ---------- t : float - Time where the interpolant is evaluated. + Time where interpolation is computed. Returns ------- p : float - Value of the interpolant at time :math:`t`. + Value of interpolant at time t. """ - - ind = np.searchsorted(self.ti, t) - if ind == 0: - p = self.yi[0] - elif ind == len(self.ti): - p = self.yi[-1] - else: - p = self.yi[ind - 1] + (self.yi[ind] - self.yi[ind - 1]) / (self.ti[ind] - self.ti[ind - 1]) * ( - t - self.ti[ind - 1] - ) + p = np.sum([self.yi[i] * self.get_Lagrange_polynomial(t, i) for i in range(self.n)]) return p diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 79ba9b5dad..2d01cfc701 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -65,105 +65,11 @@ def solve_system(self, rhs, factor, u0, t): return self.u_exact(t) -class ExactDiscontinuousTestDAE(DiscontinuousTestDAE): - r""" - Dummy DAE problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestDAE``. - """ - - def __init__(self, newton_tol=1e-8): - """Initialization routine""" - super().__init__(newton_tol) - - def eval_f(self, u, du, t): - r""" - Returns the exact right-hand side of the problem. The first components in both - cases of ``f`` are set to 1 to avoid getting a zero residual (enforced the sweeper to stop - since "convergence is reached"). - - Parameters - ---------- - u : dtype_u - Exact value of u. - du : dtype_u - Current values of the derivative of the numerical solution at time t. - t : float - Time :math:`t`. - - Returns - ------- - f : dtype_f - Right-hand side. - """ - - f = self.dtype_f(self.init) - - t_switch = np.inf if self.t_switch is None else self.t_switch - - h = 2 * u[0] - 100 - f = self.dtype_f(self.init) - - if h >= 0 or t >= t_switch: - f[:] = ( - 1, - 0, - ) - else: - f[:] = ( - 1, - 0, - ) - return f - - def solve_system(self, impl_sys, u0, t, **kwargs): - r""" - Just returns the derivative of the exact solution. - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - if t <= self.t_switch_exact: - me[:] = (np.sinh(t), np.cosh(t)) - else: - me[:] = (np.sinh(self.t_switch_exact), np.cosh(self.t_switch_exact)) - return me - - -def getTestDict(): - testDict = { - 0: { - 't': [0, 2], - 'p': [1, 5], - 'tMid': [1], - 'pMid': [3], - }, - 1: { - 't': [1, 2, 4, 5], - 'p': [6, 4, 3, 7], - 'tMid': [1.5, 3], - 'pMid': [5, 3.5], - }, - } - return testDict - - def getParamsRun(): r""" Returns parameters for conroller run that are used in each test. """ - restol = 1e-13 + restol = -1 alpha = 0.95 maxiter = 1 max_restarts = 3 @@ -215,51 +121,6 @@ def testExactDummyProblems(): fOdeEvent = parentODE.eval_f(u0Event, tExactEventODE) assert np.allclose(fExactOdeEvent, fOdeEvent), f"Right-hand sides at event do not match!" - childDAE = ExactDiscontinuousTestDAE(**{}) - parentDAE = DiscontinuousTestDAE(**{}) - assert childDAE.t_switch_exact == parentDAE.t_switch_exact, f"Exact event times between classes does not match!" - - t0 = 0.0 - u0 = mesh(childDAE.init) - duExactDAE = childDAE.solve_system(impl_sys=None, u0=u0, t=t0) - duDAE = mesh(parentDAE.init) - duDAE[:] = (np.sinh(t0), np.cosh(t0)) - assert np.allclose(duExactDAE, duDAE), f"Method solve_system of dummy problem do not return exact derivative!" - - tExactEventDAE = parentDAE.t_switch_exact - u0Event = mesh(childDAE.init) - duExactDAEEvent = childDAE.solve_system(impl_sys=None, u0=u0Event, t=tExactEventDAE) - duDAEEvent = mesh(parentDAE.init) - duDAEEvent[:] = (np.sinh(tExactEventDAE), np.cosh(tExactEventDAE)) - assert np.allclose( - duExactDAEEvent, duDAEEvent - ), f"Method solve_system of dummy problem do not return exact derivative at event!" - - -@pytest.mark.base -@pytest.mark.parametrize("key", [0, 1]) -def testInterpolationValues(key): - """ - Test for linear interpolation class in switch estimator. Linear interpolation is tested against - values in testDict that contains values for interpolation computed by hand. - Further, NumPy's routine interp is used as reference to have a second test on hand. - """ - from pySDC.projects.PinTSimE.switch_estimator import LinearInterpolation - - testDict = getTestDict() - testSet = testDict[key] - t, p = testSet['t'], testSet['p'] - tMid, pMid = testSet['tMid'], testSet['pMid'] - LinearInterpolator = LinearInterpolation(t, p) - - for m in range(len(t)): - assert LinearInterpolator.eval(t[m]) == p[m] - assert np.interp(t[m], t, p) == LinearInterpolator.eval(t[m]) - - for m in range(len(tMid)): - assert LinearInterpolator.eval(tMid[m]) == pMid[m] - assert np.interp(tMid[m], t, p) == LinearInterpolator.eval(tMid[m]) - @pytest.mark.base @pytest.mark.parametrize('quad_type', ['LOBATTO', 'RADAU-RIGHT']) @@ -511,9 +372,10 @@ def testDetectionODE(tol, num_nodes, quad_type): @pytest.mark.base -@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(8, 13)]) +@pytest.mark.parametrize('dt', np.logspace(-2.5, -1.5, num=4)) +@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(9, 13)]) @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -def testDetectionDAE(tol, num_nodes): +def testDetectionDAE(dt, tol, num_nodes): r""" In this test, the switch estimator is applied to a DAE dummy problem of ``DiscontinuousTestDAE``, where the dynamics of the differential equation is replaced by its exact dynamics to see if @@ -534,24 +396,31 @@ def testDetectionDAE(tol, num_nodes): """ from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE + from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE from pySDC.helpers.stats_helper import get_sorted from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun from pySDC.implementations.hooks.log_solution import LogSolution from pySDC.implementations.hooks.log_restarts import LogRestarts + from pySDC.projects.PinTSimE.paper_PSCC2024.log_event import LogEventDiscontinuousTestDAE - problem = ExactDiscontinuousTestDAE + problem = DiscontinuousTestDAE problem_params = dict() t0 = 4.6 Tend = 4.62 - dt = Tend - t0 sweeper = fully_implicit_DAE - QI = 'IE' + QI = 'LU' quad_type = 'RADAU-RIGHT' - restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD = getParamsRun() + _, _, _, _, useA, useSE, exact_event_time_avail, _ = getParamsRun() - hook_class = [LogSolution, LogRestarts] + restol=1e-13 + maxiter = 60 + typeFD = 'backward' + max_restarts = 20 + alpha = 0.9 + + hook_class = [LogSolution, LogRestarts, LogEventDiscontinuousTestDAE] description, controller_params, controller = generateDescription( dt=dt, @@ -583,8 +452,11 @@ def testDetectionDAE(tol, num_nodes): # in this specific example only one event has to be found switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] - assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol}!' + assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol} and M={num_nodes}!' t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=8e-14), f'Event time error {event_err} is not small enough!' + assert np.isclose(event_err, 0, atol=3e-7), f'Event time error {event_err} is not small enough!' + + h = np.array([np.abs(val[1]) for val in get_sorted(stats, type='state_function', sortby='time', recomputed=False)]) + assert np.isclose(h[-1], 0.0, atol=3e-6), f'State function is not close to zero; value is {h[-1]}' \ No newline at end of file From b5fadcc3f7ecc857a6c3f69c49d1e44177e67156 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 27 Feb 2024 14:15:36 +0100 Subject: [PATCH 18/29] Added test for state function + global error --- .../test_pintsime/test_SwitchEstimator.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 2d01cfc701..dff481df49 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -401,6 +401,7 @@ def testDetectionDAE(dt, tol, num_nodes): from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun from pySDC.implementations.hooks.log_solution import LogSolution from pySDC.implementations.hooks.log_restarts import LogRestarts + from pySDC.projects.DAE.misc.HookClass_DAE import error_hook from pySDC.projects.PinTSimE.paper_PSCC2024.log_event import LogEventDiscontinuousTestDAE problem = DiscontinuousTestDAE @@ -418,9 +419,9 @@ def testDetectionDAE(dt, tol, num_nodes): maxiter = 60 typeFD = 'backward' max_restarts = 20 - alpha = 0.9 + alpha = 0.94 - hook_class = [LogSolution, LogRestarts, LogEventDiscontinuousTestDAE] + hook_class = [LogSolution, LogRestarts, LogEventDiscontinuousTestDAE, error_hook] description, controller_params, controller = generateDescription( dt=dt, @@ -456,7 +457,12 @@ def testDetectionDAE(dt, tol, num_nodes): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=3e-7), f'Event time error {event_err} is not small enough!' + assert np.isclose(event_err, 0, atol=1e-10), f'Event time error {event_err} is not small enough!' - h = np.array([np.abs(val[1]) for val in get_sorted(stats, type='state_function', sortby='time', recomputed=False)]) - assert np.isclose(h[-1], 0.0, atol=3e-6), f'State function is not close to zero; value is {h[-1]}' \ No newline at end of file + h = np.array([val[1] for val in get_sorted(stats, type='state_function', sortby='time', recomputed=False)]) + if h[-1] < 0: + assert h[-1] > -1e-10, f"State function has large negative value -> SE does switch too early!" + assert np.isclose(abs(h[-1]), 0.0, atol=1e-11), f'State function is not close to zero; value is {h[-1]}' + + e_global = np.array(get_sorted(stats, type='error_post_step', sortby='time', recomputed=False)) + assert np.isclose(e_global[-1, 1], 0.0, atol=1e-11), f"Error at end time is too large! Expected {1e-11}, got {e_global[-1, 1]}" From 563bcfd2920fab7199b979eac89ebb40f7b7541a Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 5 Mar 2024 10:14:59 +0100 Subject: [PATCH 19/29] Updated LogEvent hooks --- pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py b/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py index 8a372babc4..02e724d209 100644 --- a/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py +++ b/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py @@ -21,7 +21,7 @@ def post_step(self, step, level_number): iter=0, sweep=L.status.sweep, type='state_function', - value=2 * L.uend[0] - 100, + value=2 * L.uend.diff[0] - 100, ) @@ -46,5 +46,5 @@ def post_step(self, step, level_number): iter=0, sweep=L.status.sweep, type='state_function', - value=L.uend[10 * P.m] - P.psv_max, + value=L.uend.diff[10 * P.m] - P.psv_max, ) From 22f2bcc8a6cf512a002cdaf75f668860d4be1a5b Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 5 Mar 2024 10:15:02 +0100 Subject: [PATCH 20/29] Updated hardcoded solutions again --- pySDC/projects/PinTSimE/battery_model.py | 10 +- .../PinTSimE/discontinuous_test_ODE.py | 2 +- pySDC/projects/PinTSimE/estimation_check.py | 2 +- .../projects/PinTSimE/hardcoded_solutions.py | 161 +++++++++--------- pySDC/projects/PinTSimE/switch_estimator.py | 47 ++--- .../test_pintsime/test_SwitchEstimator.py | 42 +++-- 6 files changed, 115 insertions(+), 149 deletions(-) diff --git a/pySDC/projects/PinTSimE/battery_model.py b/pySDC/projects/PinTSimE/battery_model.py index 0cfac58156..0949e6ffdf 100644 --- a/pySDC/projects/PinTSimE/battery_model.py +++ b/pySDC/projects/PinTSimE/battery_model.py @@ -61,7 +61,6 @@ def generateDescription( max_restarts=None, tol_event=1e-10, alpha=1.0, - typeFD='backward', ): r""" Generate a description for the battery models for a controller run. @@ -99,9 +98,6 @@ def generateDescription( Tolerance for event detection to terminate. alpha : float, optional Factor that indicates how the new step size in the Switch Estimator is reduced. - typeFD : str, optional - Type of finite difference in switch estimator. For detailed documentation, see - switch_estimator.py. Returns ------- @@ -143,11 +139,9 @@ def generateDescription( # convergence controllers convergence_controllers = {} if use_switch_estimator: - typeFD = typeFD if typeFD is not None else 'backward' switch_estimator_params = { 'tol': tol_event, 'alpha': alpha, - 'typeFD': typeFD, } convergence_controllers.update({SwitchEstimator: switch_estimator_params}) if use_adaptivity: @@ -241,9 +235,8 @@ def main(): 'max_restarts': 50, 'recomputed': False, 'tol_event': 1e-10, - 'alpha': 0.95, + 'alpha': 0.96, 'exact_event_time_avail': None, - 'typeFD': 'backward', } all_params = { @@ -385,7 +378,6 @@ def runSimulation(problem, sweeper, all_params, use_adaptivity, use_detection, h max_restarts=handling_params['max_restarts'], tol_event=handling_params['tol_event'], alpha=handling_params['alpha'], - typeFD=handling_params['typeFD'], ) stats, t_switch_exact = controllerRun( diff --git a/pySDC/projects/PinTSimE/discontinuous_test_ODE.py b/pySDC/projects/PinTSimE/discontinuous_test_ODE.py index a264536444..22f36fe426 100644 --- a/pySDC/projects/PinTSimE/discontinuous_test_ODE.py +++ b/pySDC/projects/PinTSimE/discontinuous_test_ODE.py @@ -68,7 +68,7 @@ def main(): 'max_restarts': 50, 'recomputed': False, 'tol_event': 1e-12, - 'alpha': 0.95, + 'alpha': 0.96, 'exact_event_time_avail': True, 'typeFD': 'backward', } diff --git a/pySDC/projects/PinTSimE/estimation_check.py b/pySDC/projects/PinTSimE/estimation_check.py index d0fdb02655..51a332477f 100644 --- a/pySDC/projects/PinTSimE/estimation_check.py +++ b/pySDC/projects/PinTSimE/estimation_check.py @@ -45,7 +45,7 @@ def run_estimation_check(): 'max_restarts': 50, 'recomputed': False, 'tol_event': 1e-10, - 'alpha': 0.95, + 'alpha': 0.96, 'exact_event_time_avail': None, 'typeFD': 'backward', } diff --git a/pySDC/projects/PinTSimE/hardcoded_solutions.py b/pySDC/projects/PinTSimE/hardcoded_solutions.py index 76359d2bad..10c64dc222 100644 --- a/pySDC/projects/PinTSimE/hardcoded_solutions.py +++ b/pySDC/projects/PinTSimE/hardcoded_solutions.py @@ -43,23 +43,23 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator and adaptivity for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.5393867578149492, - 'vC': 0.9999999999980306, - 't_switches': [0.18232155679314516, 0.18232155674931952], - 'dt': 0.017755006862796352, - 'e_em': 2.220446049250313e-16, - 'sum_restarts': 17.0, - 'sum_niters': 312.0, + 'iL': 0.5393867577881641, + 'vC': 0.9999999999913842, + 't_switches': [0.1823215567921536, 0.18232155679215356, 0.18232155679173784], + 'dt': 0.09453745651144455, + 'e_em': 1.7587042933087105e-12, + 'sum_restarts': 15.0, + 'sum_niters': 280.0, } elif dt == 1e-3: expected = { - 'iL': 0.5393867577908672, - 'vC': 0.9999999999869332, - 't_switches': [0.182321556754745], - 'dt': 0.056647727943027276, - 'e_em': 4.285460875053104e-14, - 'sum_restarts': 15.0, - 'sum_niters': 296.0, + 'iL': 0.5393867577223005, + 'vC': 0.9999999999813279, + 't_switches': [0.18232155676894835, 0.1823215567897308, 0.18232155678877865], + 'dt': 0.06467602356229402, + 'e_em': 1.1468603844377867e-13, + 'sum_restarts': 17.0, + 'sum_niters': 368.0, } got.update( { @@ -75,17 +75,17 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.5490992955276419, - 'vC': 0.9999999999816187, - 't_switches': [0.182321556761159], - 'sum_restarts': 13.0, + 'iL': 0.5456190026495924, + 'vC': 0.9991666666670524, + 't_switches': [0.18232155679395579, 0.18232155679395592, 0.18232155679356965], + 'sum_restarts': 14.0, 'sum_niters': 416.0, } elif dt == 1e-3: expected = { 'iL': 0.5403849766797957, - 'vC': 0.9999166666718436, - 't_switches': [0.18232155679295683, 0.1823215567608114], + 'vC': 0.9999166666675774, + 't_switches': [0.18232155679395004, 0.18232155679303919], 'sum_restarts': 11.0, 'sum_niters': 2536.0, } @@ -151,23 +151,23 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator and adaptivity for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.5393867577386883, - 'vC': 1.0000000000016669, - 't_switches': [0.18232155680505413, 0.1823215568065723], - 'dt': 0.017701859620598615, + 'iL': 0.5393867577468375, + 'vC': 0.9999999999980123, + 't_switches': [0.18232155680038617, 0.1823215568023739], + 'dt': 0.08896232033732146, 'e_em': 2.220446049250313e-16, - 'sum_restarts': 17.0, - 'sum_niters': 312.0, + 'sum_restarts': 15.0, + 'sum_niters': 280.0, } elif dt == 1e-3: expected = { - 'iL': 0.5393867577174891, - 'vC': 0.9999999999869332, - 't_switches': [0.18232155677090756], - 'dt': 0.05664772792686473, + 'iL': 0.5393867576415584, + 'vC': 0.9999999999802239, + 't_switches': [0.18232155678530526, 0.1823215568066914, 0.1823215568057151], + 'dt': 0.06333183541149384, 'e_em': 2.220446049250313e-16, - 'sum_restarts': 15.0, - 'sum_niters': 296.0, + 'sum_restarts': 17.0, + 'sum_niters': 368.0, } got.update( { @@ -183,18 +183,18 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.549099295497225, - 'vC': 0.9999999999816187, - 't_switches': [0.18232155676455783], - 'sum_restarts': 13.0, + 'iL': 0.5490992952561473, + 'vC': 0.9999999999982524, + 't_switches': [0.1823215567992934, 0.18232155680104123], + 'sum_restarts': 14.0, 'sum_niters': 416.0, } elif dt == 1e-3: expected = { - 'iL': 0.5403849766692057, - 'vC': 0.9999166666718439, - 't_switches': [0.18232155679484557, 0.1823215567627017], - 'sum_restarts': 11.0, + 'iL': 0.5407340516779595, + 'vC': 0.9999999999936255, + 't_switches': [0.18232155676519302], + 'sum_restarts': 10.0, 'sum_niters': 2536.0, } got.update( @@ -259,25 +259,25 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator and adaptivity for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.6125019898490534, - 'vC1': 1.000000000007182, - 'vC2': 1.0000000000479325, - 't_switches': [0.18232155680297893, 0.36464311356978835], - 'dt': 0.02973949262636494, - 'e_em': 3.3306690738754696e-16, + 'iL': 0.6125019898578352, + 'vC1': 1.0000000000471956, + 'vC2': 1.0000000000165106, + 't_switches': [0.18232155678158268, 0.36464311353802376], + 'dt': 0.0985931246953285, + 'e_em': 2.295386103412511e-12, 'sum_restarts': 24.0, - 'sum_niters': 496.0, + 'sum_niters': 552.0, } elif dt == 1e-3: expected = { - 'iL': 0.6125019898915002, - 'vC1': 1.0000000000057128, - 'vC2': 1.000000000050465, - 't_switches': [0.1823215567993942, 0.36464311356557344], - 'dt': 0.0241594380776875, - 'e_em': 2.220446049250313e-16, - 'sum_restarts': 24.0, - 'sum_niters': 504.0, + 'iL': 0.6125019901065321, + 'vC1': 1.0000000000787372, + 'vC2': 1.000000000028657, + 't_switches': [0.1823215567907939, 0.3646431134803315], + 'dt': 0.07154669986159717, + 'e_em': 2.3381296898605797e-13, + 'sum_restarts': 22.0, + 'sum_niters': 472.0, } got.update( { @@ -293,21 +293,21 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'iL': 0.6313858468881018, - 'vC1': 1.0000000000155491, - 'vC2': 1.0000000000316576, - 't_switches': [0.18232155681737047, 0.36464311357240825], - 'sum_restarts': 20.0, - 'sum_niters': 680.0, + 'iL': 0.6313858468030417, + 'vC1': 1.0000000002414198, + 'vC2': 1.0000000000095093, + 't_switches': [0.18232155679395579, 0.3646431133464922], + 'sum_restarts': 19.0, + 'sum_niters': 664.0, } elif dt == 1e-3: expected = { - 'iL': 0.6151254297607596, - 'vC1': 1.0000000000846776, - 'vC2': 1.0000000000993976, - 't_switches': [0.18232155679295683, 0.3646431135032376], - 'sum_restarts': 17.0, - 'sum_niters': 4232.0, + 'iL': 0.6151254295045797, + 'vC1': 1.0000000000227713, + 'vC2': 1.0000000000329365, + 't_switches': [0.18232155680153855, 0.3646431135651182], + 'sum_restarts': 16.0, + 'sum_niters': 4224.0, } got.update( { @@ -375,21 +375,21 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error when using switch estimator for {prob_cls_name} for dt={dt:.1e}:" if dt == 1e-2: expected = { - 'u': 5.994135895337547, - 'e_global': 7.765876830490015e-11, - 't_switches': [1.6094379124455047], - 'e_event': [1.1404432953554533e-11], + 'u': 5.998326995729771, + 'e_global': 0.0041911003550794135, + 't_switches': [1.6094379124660123], + 'e_event': [3.1912028575220575e-11], 'sum_restarts': 22.0, - 'sum_niters': 686.0, + 'sum_niters': 675.0, } elif dt == 1e-3: expected = { - 'u': 5.972186947618947, - 'e_global': 0.00041911003080929987, - 't_switches': [1.6094379124055875], - 'e_event': [2.851274771842327e-11], - 'sum_restarts': 17.0, - 'sum_niters': 2338.0, + 'u': 5.9721869476933005, + 'e_global': 0.0004191100622819022, + 't_switches': [1.6094379124262566, 1.6094379124260099], + 'e_event': [7.843725668976731e-12], + 'sum_restarts': 20.0, + 'sum_niters': 2352.0, } got.update( { @@ -474,7 +474,10 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): for key in expected.keys(): if key == 't_switches' or key == 'e_event': err_msg = f'{msg} Expected {key}={expected[key]}, got {key}={got[key]}' - assert all(np.isclose(expected[key], got[key], atol=1e-4)) == True, err_msg + if len(expected[key]) == got[key]: + assert np.allclose(expected[key], got[key], atol=1e-4) == True, err_msg + else: + assert np.isclose(expected[key][-1], got[key][-1], atol=1e-4) == True, err_msg else: err_msg = f'{msg} Expected {key}={expected[key]:.4e}, got {key}={got[key]:.4e}' assert np.isclose(expected[key], got[key], atol=1e-4), err_msg diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 51b96613b9..7ea35fe071 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -48,7 +48,7 @@ def setup(self, controller, params, description): defaults = { 'control_order': 0, 'nodes': coll.nodes, - 'tol_zero': 1e-13, + 'tol_zero': 2.5e-12, 't_interp': [], 'state_function': [], } @@ -125,7 +125,7 @@ def get_new_step_size(self, controller, S, **kwargs): # intermediate value theorem states that a root is contained in current step if self.params.state_function[0] * self.params.state_function[-1] < 0 and self.status.is_zero is None: self.status.t_switch = self.get_switch( - self.params.t_interp, self.params.state_function, m_guess, self.params.typeFD + self.params.t_interp, self.params.state_function, m_guess ) self.logging_during_estimation( @@ -278,7 +278,7 @@ def logging_during_estimation(controller_hooks, process, time, level, sweep, t_s ) @staticmethod - def get_switch(t_interp, state_function, m_guess, choiceFD='backward'): + def get_switch(t_interp, state_function, m_guess): r""" Routine to do the interpolation and root finding stuff. @@ -290,10 +290,6 @@ def get_switch(t_interp, state_function, m_guess, choiceFD='backward'): Contains values of state function at these collocation nodes. m_guess : float Index at which the difference drops below zero. - choiceFD : str - Finite difference to be used to approximate derivative of - state function. Can be ``'forward'``, ``'centered'``, or - ``'backward'``. Default is ``'backward'`` Returns ------- @@ -305,24 +301,12 @@ def get_switch(t_interp, state_function, m_guess, choiceFD='backward'): def fprime(t): r""" - Computes the derivative of the scalar interpolant using finite difference. - Here, different finite differences can be used. The type of FD can be set by - setting ``typeFD`` in switch estimator parameters. There are three choices possible: - - - ``typeFD='backward'`` for :math:`h=10^{-10}`: - - .. math:: - \frac{dp}{dt} \approx \frac{p(t) - p(t - h)}{h} - - - ``typeFD='centered'`` for :math:`h=10^{-12}`: + Computes the derivative of the scalar interpolant using finite difference. Here, + the derivative is approximated by the backward difference: .. math:: - \frac{dp}{dt} \approx \frac{p(t + h) - p(t - h)}{2h} - - - ``typeFD='forward'`` for :math:`h=10^{-10}`: + \frac{dp}{dt} \approx \frac{25 p(t) - 48 p(t - h) + 36 p(t - 2 h) - 16 p(t - 3h) + 3 p(t - 4 h)}{12 h} - .. math:: - \frac{dp}{dt} \approx \frac{p(t + h) - p(t)}{h} Parameters ---------- @@ -335,22 +319,13 @@ def fprime(t): Derivative of interpolation p at time t. """ - if choiceFD == 'forward': - dt_FD = 1e-10 - dp = (p(t + dt_FD) - p(t)) / (dt_FD) - elif choiceFD == 'centered': - dt_FD = 1e-12 - dp = (p(t + dt_FD) - p(t - dt_FD)) / (2 * dt_FD) - elif choiceFD == 'backward': - dt_FD = 1e-12 - # dp = (p(t) - p(t - dt_FD)) / (dt_FD) - # dp = (11 * p(t) - 18 * p(t - dt_FD) + 9 * p(t - 2 * dt_FD) - 2 * p(t - 3 * dt_FD)) / (6 * dt_FD) - dp = (25 * p(t) - 48 * p(t - dt_FD) + 36 * p(t - 2 * dt_FD) - 16 * p(t - 3 * dt_FD) + 3 * p(t - 4 * dt_FD)) / (12 * dt_FD) - else: - raise NotImplementedError + dt_FD = 1e-10 + # dp = (p(t) - p(t - dt_FD)) / dt_FD + # dp = (11 * p(t) - 18 * p(t - dt_FD) + 9 * p(t - 2 * dt_FD) - 2 * p(t - 3 * dt_FD)) / (6 * dt_FD) + dp = (25 * p(t) - 48 * p(t - dt_FD) + 36 * p(t - 2 * dt_FD) - 16 * p(t - 3 * dt_FD) + 3 * p(t - 4 * dt_FD)) / (12 * dt_FD) return dp - newton_tol, newton_maxiter = 1e-15, 100 + newton_tol, newton_maxiter = 1e-14, 100 t_switch = newton(t_interp[m_guess], p, fprime, newton_tol, newton_maxiter) return t_switch diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index dff481df49..27ce6d5ee2 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -76,8 +76,7 @@ def getParamsRun(): useA = False useSE = True exact_event_time_avail = True - typeFD = 'centered' - return restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD + return restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail @pytest.mark.base @@ -158,7 +157,7 @@ def testAdaptInterpolationInfo(quad_type): tol = 1e-10 - restol, alpha, maxiter, max_restarts, useA, useSE, _, typeFD = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, _ = getParamsRun() hook_class = [] @@ -178,7 +177,6 @@ def testAdaptInterpolationInfo(quad_type): max_restarts=max_restarts, tol_event=tol, alpha=alpha, - typeFD=typeFD, ) S = controller.MS[0] @@ -255,7 +253,7 @@ def testDetectionBoundary(num_nodes): tol = 1e-10 - restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -275,7 +273,6 @@ def testDetectionBoundary(num_nodes): max_restarts=max_restarts, tol_event=tol, alpha=alpha, - typeFD=typeFD, ) stats, _ = controllerRun( @@ -330,7 +327,7 @@ def testDetectionODE(tol, num_nodes, quad_type): sweeper = generic_implicit QI = 'IE' - restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail, typeFD = getParamsRun() + restol, alpha, maxiter, max_restarts, useA, useSE, exact_event_time_avail = getParamsRun() hook_class = [LogSolution, LogRestarts] @@ -350,7 +347,6 @@ def testDetectionODE(tol, num_nodes, quad_type): max_restarts=max_restarts, tol_event=tol, alpha=alpha, - typeFD=typeFD, ) stats, t_switch_exact = controllerRun( @@ -372,10 +368,10 @@ def testDetectionODE(tol, num_nodes, quad_type): @pytest.mark.base -@pytest.mark.parametrize('dt', np.logspace(-2.5, -1.5, num=4)) -@pytest.mark.parametrize('tol', [10 ** (-m) for m in range(9, 13)]) +# @pytest.mark.parametrize('dt', np.logspace(-2.5, -1.5, num=4)) +# @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(9, 13)]) @pytest.mark.parametrize('num_nodes', [3, 4, 5]) -def testDetectionDAE(dt, tol, num_nodes): +def testDetectionDAE(num_nodes): r""" In this test, the switch estimator is applied to a DAE dummy problem of ``DiscontinuousTestDAE``, where the dynamics of the differential equation is replaced by its exact dynamics to see if @@ -401,27 +397,28 @@ def testDetectionDAE(dt, tol, num_nodes): from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun from pySDC.implementations.hooks.log_solution import LogSolution from pySDC.implementations.hooks.log_restarts import LogRestarts - from pySDC.projects.DAE.misc.HookClass_DAE import error_hook + from pySDC.projects.DAE.misc.HookClass_DAE import LogGlobalErrorPostStepDifferentialVariable from pySDC.projects.PinTSimE.paper_PSCC2024.log_event import LogEventDiscontinuousTestDAE problem = DiscontinuousTestDAE - problem_params = dict() + problem_params = dict({'newton_tol': 1e-6}) t0 = 4.6 Tend = 4.62 + dt = Tend - t0 + tol = 1e-10 sweeper = fully_implicit_DAE QI = 'LU' quad_type = 'RADAU-RIGHT' - _, _, _, _, useA, useSE, exact_event_time_avail, _ = getParamsRun() + _, _, _, _, useA, useSE, exact_event_time_avail = getParamsRun() restol=1e-13 maxiter = 60 - typeFD = 'backward' max_restarts = 20 - alpha = 0.94 + alpha = 0.97 - hook_class = [LogSolution, LogRestarts, LogEventDiscontinuousTestDAE, error_hook] + hook_class = [LogSolution, LogRestarts, LogEventDiscontinuousTestDAE, LogGlobalErrorPostStepDifferentialVariable] description, controller_params, controller = generateDescription( dt=dt, @@ -439,7 +436,6 @@ def testDetectionDAE(dt, tol, num_nodes): max_restarts=max_restarts, tol_event=tol, alpha=alpha, - typeFD=typeFD, ) stats, t_switch_exact = controllerRun( @@ -457,12 +453,12 @@ def testDetectionDAE(dt, tol, num_nodes): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=1e-10), f'Event time error {event_err} is not small enough!' + assert np.isclose(event_err, 0, atol=2.2e-6), f'Event time error {event_err} is not small enough!' h = np.array([val[1] for val in get_sorted(stats, type='state_function', sortby='time', recomputed=False)]) if h[-1] < 0: - assert h[-1] > -1e-10, f"State function has large negative value -> SE does switch too early!" - assert np.isclose(abs(h[-1]), 0.0, atol=1e-11), f'State function is not close to zero; value is {h[-1]}' + assert abs(h[-1]) < 4.7e-10, f"State function has large negative value -> SE does switch too early! Got {h[-1]}" + assert np.isclose(abs(h[-1]), 0.0, atol=5e-10), f'State function is not close to zero; value is {h[-1]}' - e_global = np.array(get_sorted(stats, type='error_post_step', sortby='time', recomputed=False)) - assert np.isclose(e_global[-1, 1], 0.0, atol=1e-11), f"Error at end time is too large! Expected {1e-11}, got {e_global[-1, 1]}" + e_global = np.array(get_sorted(stats, type='e_global_differential_post_step', sortby='time', recomputed=False)) + assert np.isclose(e_global[-1, 1], 0.0, atol=2.4e-10), f"Error at end time is too large! Expected {1e-11}, got {e_global[-1, 1]}" From 8ccbebf29258f81624e6fb6a2827873750abf800 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 5 Mar 2024 11:03:24 +0100 Subject: [PATCH 21/29] Adapted test_problems.py --- .../test_projects/test_DAE/test_problems.py | 25 +++++++++---------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/pySDC/tests/test_projects/test_DAE/test_problems.py b/pySDC/tests/test_projects/test_DAE/test_problems.py index ce03041562..2b3c211f80 100644 --- a/pySDC/tests/test_projects/test_DAE/test_problems.py +++ b/pySDC/tests/test_projects/test_DAE/test_problems.py @@ -428,7 +428,7 @@ def test_DiscontinuousTestDAE_SDC(M): @pytest.mark.base -@pytest.mark.parametrize('M', [2, 3, 4, 5]) +@pytest.mark.parametrize('M', [3, 4, 5]) def test_DiscontinuousTestDAE_SDC_detection(M): """ Test for one SDC run with event detection if the found event is close to the exact value and if the global error @@ -443,17 +443,15 @@ def test_DiscontinuousTestDAE_SDC_detection(M): from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI err_tol = { - 2: 5.3952e-9, - 3: 2.6741e-9, - 4: 1.9163e-8, - 5: 2.4791e-8, + 3: 5.97e-13, + 4: 1.43e-10, + 5: 3.18e-10, } event_err_tol = { - 2: 0.0011, - 3: 0.01, - 4: 0.02, - 5: 0.0101, + 3: 0.02, + 4: 8.22e-13, + 5: 1.64e-12, } level_params = { @@ -481,7 +479,7 @@ def test_DiscontinuousTestDAE_SDC_detection(M): switch_estimator_params = { 'tol': 1e-10, - 'alpha': 0.95, + 'alpha': 0.96, } restarting_params = { @@ -514,8 +512,8 @@ def test_DiscontinuousTestDAE_SDC_detection(M): uex = P.u_exact(Tend) uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - err = abs(uex.diff[0] - uend.diff[0]) + print(err) assert err < err_tol[M], f"ERROR for M={M}: Error is too large! Expected {err_tol[M]}, got {err}" switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) @@ -525,6 +523,7 @@ def test_DiscontinuousTestDAE_SDC_detection(M): t_switch_exact = P.t_switch_exact event_err = abs(t_switch_exact - t_switch) + print(event_err) assert ( event_err < event_err_tol[M] ), f"ERROR for M={M}: Event error is too large! Expected {event_err_tol[M]}, got {event_err}" @@ -660,7 +659,7 @@ def test_WSCC9_SDC_detection(): switch_estimator_params = { 'tol': 1e-10, - 'alpha': 0.95, + 'alpha': 0.97, } restarting_params = { @@ -696,7 +695,7 @@ def test_WSCC9_SDC_detection(): switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) assert len(switches) >= 1, 'ERROR: No events found!' t_switch = [me[1] for me in switches][0] - assert np.isclose(t_switch, 0.5335289411812812, atol=1e-3), 'Found event does not match a threshold!' + assert np.isclose(t_switch, 0.528458886745887, atol=1e-3), f'Found event does not match a threshold! Got {t_switch}' # @pytest.mark.base From 0bfcae7a531c62e75d67670df6db712de67c9988 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 5 Mar 2024 11:03:29 +0100 Subject: [PATCH 22/29] Minor changes --- pySDC/projects/PinTSimE/buck_model.py | 1 - pySDC/projects/PinTSimE/estimation_check.py | 1 - pySDC/projects/PinTSimE/piline_model.py | 1 - pySDC/projects/PinTSimE/switch_estimator.py | 14 ++++++-------- .../test_pintsime/test_SwitchEstimator.py | 10 +++++----- 5 files changed, 11 insertions(+), 16 deletions(-) diff --git a/pySDC/projects/PinTSimE/buck_model.py b/pySDC/projects/PinTSimE/buck_model.py index 4b1e3a27dc..3a2423d339 100644 --- a/pySDC/projects/PinTSimE/buck_model.py +++ b/pySDC/projects/PinTSimE/buck_model.py @@ -34,7 +34,6 @@ def main(): 'tol_event': None, 'alpha': None, 'exact_event_time_avail': None, - 'typeFD': None, } # ---- all parameters are stored in this dictionary - note: defaults are used for the problem ---- diff --git a/pySDC/projects/PinTSimE/estimation_check.py b/pySDC/projects/PinTSimE/estimation_check.py index 51a332477f..7d95358172 100644 --- a/pySDC/projects/PinTSimE/estimation_check.py +++ b/pySDC/projects/PinTSimE/estimation_check.py @@ -47,7 +47,6 @@ def run_estimation_check(): 'tol_event': 1e-10, 'alpha': 0.96, 'exact_event_time_avail': None, - 'typeFD': 'backward', } problem_classes = [battery, battery_implicit, battery_n_capacitors] diff --git a/pySDC/projects/PinTSimE/piline_model.py b/pySDC/projects/PinTSimE/piline_model.py index be135eb0ed..fc2f5f56f6 100644 --- a/pySDC/projects/PinTSimE/piline_model.py +++ b/pySDC/projects/PinTSimE/piline_model.py @@ -34,7 +34,6 @@ def main(): 'tol_event': None, 'alpha': None, 'exact_event_time_avail': None, - 'typeFD': None, } # ---- all parameters are stored in this dictionary - note: defaults are used for the problem ---- diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index 7ea35fe071..ba7c79662b 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -124,9 +124,7 @@ def get_new_step_size(self, controller, S, **kwargs): # intermediate value theorem states that a root is contained in current step if self.params.state_function[0] * self.params.state_function[-1] < 0 and self.status.is_zero is None: - self.status.t_switch = self.get_switch( - self.params.t_interp, self.params.state_function, m_guess - ) + self.status.t_switch = self.get_switch(self.params.t_interp, self.params.state_function, m_guess) self.logging_during_estimation( controller.hooks[0], @@ -296,8 +294,8 @@ def get_switch(t_interp, state_function, m_guess): t_switch : float Time point of found event. """ - LinearInterpolator = LagrangeInterpolation(t_interp, state_function) - p = lambda t: LinearInterpolator.eval(t) + LagrangeInterpolator = LagrangeInterpolation(t_interp, state_function) + p = lambda t: LagrangeInterpolator.eval(t) def fprime(t): r""" @@ -320,9 +318,9 @@ def fprime(t): """ dt_FD = 1e-10 - # dp = (p(t) - p(t - dt_FD)) / dt_FD - # dp = (11 * p(t) - 18 * p(t - dt_FD) + 9 * p(t - 2 * dt_FD) - 2 * p(t - 3 * dt_FD)) / (6 * dt_FD) - dp = (25 * p(t) - 48 * p(t - dt_FD) + 36 * p(t - 2 * dt_FD) - 16 * p(t - 3 * dt_FD) + 3 * p(t - 4 * dt_FD)) / (12 * dt_FD) + dp = ( + 25 * p(t) - 48 * p(t - dt_FD) + 36 * p(t - 2 * dt_FD) - 16 * p(t - 3 * dt_FD) + 3 * p(t - 4 * dt_FD) + ) / (12 * dt_FD) return dp newton_tol, newton_maxiter = 1e-14, 100 diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index 27ce6d5ee2..e97fc0a7dd 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -80,7 +80,7 @@ def getParamsRun(): @pytest.mark.base -def testExactDummyProblems(): +def testExactDummyProblem(): r""" Test for dummy problems. The test verifies that the dummy problems exactly returns the dynamics of the parent class. ``eval_f`` of ``ExactDiscontinuousTestDAE`` is not tested here, since it only returns @@ -368,8 +368,6 @@ def testDetectionODE(tol, num_nodes, quad_type): @pytest.mark.base -# @pytest.mark.parametrize('dt', np.logspace(-2.5, -1.5, num=4)) -# @pytest.mark.parametrize('tol', [10 ** (-m) for m in range(9, 13)]) @pytest.mark.parametrize('num_nodes', [3, 4, 5]) def testDetectionDAE(num_nodes): r""" @@ -413,7 +411,7 @@ def testDetectionDAE(num_nodes): _, _, _, _, useA, useSE, exact_event_time_avail = getParamsRun() - restol=1e-13 + restol = 1e-13 maxiter = 60 max_restarts = 20 alpha = 0.97 @@ -461,4 +459,6 @@ def testDetectionDAE(num_nodes): assert np.isclose(abs(h[-1]), 0.0, atol=5e-10), f'State function is not close to zero; value is {h[-1]}' e_global = np.array(get_sorted(stats, type='e_global_differential_post_step', sortby='time', recomputed=False)) - assert np.isclose(e_global[-1, 1], 0.0, atol=2.4e-10), f"Error at end time is too large! Expected {1e-11}, got {e_global[-1, 1]}" + assert np.isclose( + e_global[-1, 1], 0.0, atol=2.4e-10 + ), f"Error at end time is too large! Expected {1e-11}, got {e_global[-1, 1]}" From 827a8edff40e5b721d33d02f40de9d343ad2bb4c Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 7 Mar 2024 09:33:27 +0100 Subject: [PATCH 23/29] Updated tests --- .../test_projects/test_DAE/test_problems.py | 24 +++++++------------ .../test_pintsime/test_SwitchEstimator.py | 18 +++++++------- 2 files changed, 16 insertions(+), 26 deletions(-) diff --git a/pySDC/tests/test_projects/test_DAE/test_problems.py b/pySDC/tests/test_projects/test_DAE/test_problems.py index 2b3c211f80..e0cf70ccbb 100644 --- a/pySDC/tests/test_projects/test_DAE/test_problems.py +++ b/pySDC/tests/test_projects/test_DAE/test_problems.py @@ -335,7 +335,7 @@ def test_DiscontinuousTestDAE_singularity(): assert np.isclose( abs(f_before_event), 0.0 - ), f"ERROR: Right-hand side after event does not match! Expected {(0.0, 0.0)}, got {f_before_event}" + ), f"ERROR: Right-hand side after event does not match! Expected {(0.0, 0.0)}, got {f_before_event=}" # test for t <= t^* u_event = disc_test_DAE.u_exact(t_event) @@ -424,7 +424,7 @@ def test_DiscontinuousTestDAE_SDC(M): uend, _ = controller.run(u0=uinit, t0=t0, Tend=Tend) err = abs(uex.diff[0] - uend.diff[0]) - assert err < err_tol[M], f"ERROR: Error is too large! Expected {err_tol[M]}, got {err}" + assert err < err_tol[M], f"ERROR: Error is too large! Expected {err_tol[M]=}, got {err=}" @pytest.mark.base @@ -442,20 +442,14 @@ def test_DiscontinuousTestDAE_SDC_detection(M): from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI - err_tol = { - 3: 5.97e-13, - 4: 1.43e-10, - 5: 3.18e-10, - } - event_err_tol = { 3: 0.02, - 4: 8.22e-13, - 5: 1.64e-12, + 4: 5e-10, + 5: 1e-10, } level_params = { - 'restol': 1e-13, + 'restol': 1e-10, 'dt': 1e-2, } @@ -513,8 +507,7 @@ def test_DiscontinuousTestDAE_SDC_detection(M): uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) err = abs(uex.diff[0] - uend.diff[0]) - print(err) - assert err < err_tol[M], f"ERROR for M={M}: Error is too large! Expected {err_tol[M]}, got {err}" + assert err < 2e-9, f"ERROR for M={M}: Error is too large! Expected something lower than {2e-9}, got {err=}" switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) assert len(switches) >= 1, 'ERROR for M={M}: No events found!' @@ -523,10 +516,9 @@ def test_DiscontinuousTestDAE_SDC_detection(M): t_switch_exact = P.t_switch_exact event_err = abs(t_switch_exact - t_switch) - print(event_err) assert ( event_err < event_err_tol[M] - ), f"ERROR for M={M}: Event error is too large! Expected {event_err_tol[M]}, got {event_err}" + ), f"ERROR for M={M}: Event error is too large! Expected {event_err_tol[M]=}, got {event_err=}" @pytest.mark.base @@ -695,7 +687,7 @@ def test_WSCC9_SDC_detection(): switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) assert len(switches) >= 1, 'ERROR: No events found!' t_switch = [me[1] for me in switches][0] - assert np.isclose(t_switch, 0.528458886745887, atol=1e-3), f'Found event does not match a threshold! Got {t_switch}' + assert np.isclose(t_switch, 0.528458886745887, atol=1e-3), f'Found event does not match a threshold! Got {t_switch=}' # @pytest.mark.base diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index e97fc0a7dd..f4958e0b58 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -213,12 +213,12 @@ def testAdaptInterpolationInfo(quad_type): assert t_interp[0] != t_interp[1], 'Starting time from interpolation axis is not removed!' assert ( len(t_interp) == num_nodes - ), f'Number of values on interpolation axis does not match. Expected {num_nodes}, got {len(t_interp)}' + ), f'Number of values on interpolation axis does not match. Expected {num_nodes=}, got {len(t_interp)}' elif quad_type == 'RADAU-RIGHT': assert ( len(t_interp) == num_nodes + 1 - ), f'Number of values on interpolation axis does not match. Expected {num_nodes + 1}, got {len(t_interp)}' + ), f'Number of values on interpolation axis does not match. Expected {num_nodes + 1=}, got {len(t_interp)}' @pytest.mark.base @@ -364,7 +364,7 @@ def testDetectionODE(tol, num_nodes, quad_type): t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=1.2e-11), f'Event time error {event_err} is not small enough!' + assert np.isclose(event_err, 0, atol=1.2e-11), f'Event time error {event_err=} is not small enough!' @pytest.mark.base @@ -411,7 +411,7 @@ def testDetectionDAE(num_nodes): _, _, _, _, useA, useSE, exact_event_time_avail = getParamsRun() - restol = 1e-13 + restol = 1e-11 maxiter = 60 max_restarts = 20 alpha = 0.97 @@ -447,18 +447,16 @@ def testDetectionDAE(num_nodes): # in this specific example only one event has to be found switches = [me[1] for me in get_sorted(stats, type='switch', sortby='time', recomputed=False)] - assert len(switches) >= 1, f'{problem.__name__}: No events found for tol={tol} and M={num_nodes}!' + assert len(switches) >= 1, f'{problem.__name__}: No events found for {tol=} and {num_nodes=}!' t_switch = switches[-1] event_err = abs(t_switch - t_switch_exact) - assert np.isclose(event_err, 0, atol=2.2e-6), f'Event time error {event_err} is not small enough!' + assert np.isclose(event_err, 0, atol=2.2e-6), f'Event time error {event_err=} is not small enough!' h = np.array([val[1] for val in get_sorted(stats, type='state_function', sortby='time', recomputed=False)]) - if h[-1] < 0: - assert abs(h[-1]) < 4.7e-10, f"State function has large negative value -> SE does switch too early! Got {h[-1]}" - assert np.isclose(abs(h[-1]), 0.0, atol=5e-10), f'State function is not close to zero; value is {h[-1]}' + assert np.isclose(abs(h[-1]), 0.0, atol=2e-9), f'State function is not close to zero; value is {h[-1]}' e_global = np.array(get_sorted(stats, type='e_global_differential_post_step', sortby='time', recomputed=False)) assert np.isclose( - e_global[-1, 1], 0.0, atol=2.4e-10 + e_global[-1, 1], 0.0, atol=9.93e-10 ), f"Error at end time is too large! Expected {1e-11}, got {e_global[-1, 1]}" From c3107c568cccd4034137311c5bfa8d3e4fbb938e Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 7 Mar 2024 09:33:36 +0100 Subject: [PATCH 24/29] Speed-up test for buck converter --- pySDC/projects/PinTSimE/buck_model.py | 2 +- pySDC/projects/PinTSimE/hardcoded_solutions.py | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/pySDC/projects/PinTSimE/buck_model.py b/pySDC/projects/PinTSimE/buck_model.py index 3a2423d339..1394d8f397 100644 --- a/pySDC/projects/PinTSimE/buck_model.py +++ b/pySDC/projects/PinTSimE/buck_model.py @@ -55,7 +55,7 @@ def main(): use_adaptivity=use_adaptivity, use_detection=use_detection, hook_class=hook_class, - interval=(0.0, 2e-2), + interval=(0.0, 1e-2), dt_list=[1e-5, 2e-5], nnodes=[M_fix], ) diff --git a/pySDC/projects/PinTSimE/hardcoded_solutions.py b/pySDC/projects/PinTSimE/hardcoded_solutions.py index 10c64dc222..3a8b7b93f6 100644 --- a/pySDC/projects/PinTSimE/hardcoded_solutions.py +++ b/pySDC/projects/PinTSimE/hardcoded_solutions.py @@ -450,17 +450,17 @@ def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): msg = f"Error for {prob_cls_name} for dt={dt:.1e}:" if dt == 2e-5: expected = { - 'vC1': 9.890997780767632, - 'vC2': 4.710415385551326, - 'iLp': -0.315406990615236, - 'sum_niters': 5036.0, + 'vC1': 9.781955920747619, + 'vC2': 6.396971204930281, + 'iLp': -1.1056614708409171, + 'sum_niters': 2519.0, } elif dt == 1e-5: expected = { - 'vC1': 9.891508522329485, - 'vC2': 4.70939963429714, - 'iLp': -0.32177442457657557, - 'sum_niters': 8262.0, + 'vC1': 9.782142840662102, + 'vC2': 6.388775533709242, + 'iLp': -1.0994027552202539, + 'sum_niters': 4242.0, } got.update( { From 74ee1ffe0031b62ef57fa5180920f175bc1028e6 Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 7 Mar 2024 09:37:57 +0100 Subject: [PATCH 25/29] Black.. --- pySDC/tests/test_projects/test_DAE/test_problems.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pySDC/tests/test_projects/test_DAE/test_problems.py b/pySDC/tests/test_projects/test_DAE/test_problems.py index e0cf70ccbb..dbc901cfbb 100644 --- a/pySDC/tests/test_projects/test_DAE/test_problems.py +++ b/pySDC/tests/test_projects/test_DAE/test_problems.py @@ -687,7 +687,9 @@ def test_WSCC9_SDC_detection(): switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) assert len(switches) >= 1, 'ERROR: No events found!' t_switch = [me[1] for me in switches][0] - assert np.isclose(t_switch, 0.528458886745887, atol=1e-3), f'Found event does not match a threshold! Got {t_switch=}' + assert np.isclose( + t_switch, 0.528458886745887, atol=1e-3 + ), f'Found event does not match a threshold! Got {t_switch=}' # @pytest.mark.base From 872e81e5536fbfec9ea6bc96fd2cdd5206788776 Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 7 Mar 2024 09:39:54 +0100 Subject: [PATCH 26/29] Use msg about convergence info in Newton in SE --- pySDC/projects/PinTSimE/switch_estimator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index ba7c79662b..03835712ff 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -400,7 +400,7 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): msg = f'Newton did not converge after {n} iterations, error is {res}' else: msg = f'Newton did converge after {n} iterations, error for root {x0} is {res}' - # print(msg) + print(msg) root = x0 return root From bc751fdca3058971252973d2ab19b4677727782c Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 11 Apr 2024 08:57:59 +0200 Subject: [PATCH 27/29] Moved dummy problem to file --- .../problem_classes/DiscontinuousTestODE.py | 61 ++++++++++++++ .../test_pintsime/test_SwitchEstimator.py | 79 +++---------------- 2 files changed, 72 insertions(+), 68 deletions(-) diff --git a/pySDC/implementations/problem_classes/DiscontinuousTestODE.py b/pySDC/implementations/problem_classes/DiscontinuousTestODE.py index af797577bd..111a8267ea 100644 --- a/pySDC/implementations/problem_classes/DiscontinuousTestODE.py +++ b/pySDC/implementations/problem_classes/DiscontinuousTestODE.py @@ -222,3 +222,64 @@ def count_switches(self): Setter to update the number of switches if one is found. """ self.nswitches += 1 + + +class ExactDiscontinuousTestODE(DiscontinuousTestODE): + r""" + Dummy ODE problem for testing the ``SwitchEstimator`` class. The problem contains the exact dynamics + of the problem class ``DiscontinuousTestODE``. + """ + + def __init__(self, newton_maxiter=100, newton_tol=1e-8): + """Initialization routine""" + super().__init__(newton_maxiter, newton_tol) + + def eval_f(self, u, t): + """ + Derivative. + + Parameters + ---------- + u : dtype_u + Exact value of u. + t : float + Time :math:`t`. + + Returns + ------- + f : dtype_f + Derivative. + """ + + f = self.dtype_f(self.init) + + t_switch = np.inf if self.t_switch is None else self.t_switch + h = u[0] - 5 + if h >= 0 or t >= t_switch: + f[:] = 1 + else: + f[:] = np.exp(t) + return f + + def solve_system(self, rhs, factor, u0, t): + """ + Just return the exact solution... + + Parameters + ---------- + rhs : dtype_f + Right-hand side for the linear system. + factor : float + Abbrev. for the local stepsize (or any other factor required). + u0 : dtype_u + Initial guess for the iterative solver. + t : float + Current time (e.g. for time-dependent BCs). + + Returns + ------- + me : dtype_u + The solution as mesh. + """ + + return self.u_exact(t) diff --git a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py index f4958e0b58..4de3df9db5 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py +++ b/pySDC/tests/test_projects/test_pintsime/test_SwitchEstimator.py @@ -1,69 +1,6 @@ import numpy as np import pytest -from pySDC.implementations.problem_classes.DiscontinuousTestODE import DiscontinuousTestODE -from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE - - -class ExactDiscontinuousTestODE(DiscontinuousTestODE): - r""" - Dummy ODE problem for testing. The problem contains the exact dynamics of the problem class ``DiscontinuousTestODE``. - """ - - def __init__(self, newton_maxiter=100, newton_tol=1e-8): - """Initialization routine""" - super().__init__(newton_maxiter, newton_tol) - - def eval_f(self, u, t): - """ - Derivative. - - Parameters - ---------- - u : dtype_u - Exact value of u. - t : float - Time :math:`t`. - - Returns - ------- - f : dtype_f - Derivative. - """ - - f = self.dtype_f(self.init) - - t_switch = np.inf if self.t_switch is None else self.t_switch - h = u[0] - 5 - if h >= 0 or t >= t_switch: - f[:] = 1 - else: - f[:] = np.exp(t) - return f - - def solve_system(self, rhs, factor, u0, t): - """ - Just return the exact solution... - - Parameters - ---------- - rhs : dtype_f - Right-hand side for the linear system. - factor : float - Abbrev. for the local stepsize (or any other factor required). - u0 : dtype_u - Initial guess for the iterative solver. - t : float - Current time (e.g. for time-dependent BCs). - - Returns - ------- - me : dtype_u - The solution as mesh. - """ - - return self.u_exact(t) - def getParamsRun(): r""" @@ -87,11 +24,14 @@ def testExactDummyProblem(): a random right-hand side to enforce the sweeper to do not stop to compute. """ - from pySDC.implementations.datatype_classes.mesh import mesh + from pySDC.implementations.problem_classes.DiscontinuousTestODE import ( + DiscontinuousTestODE, + ExactDiscontinuousTestODE, + ) childODE = ExactDiscontinuousTestODE(**{}) parentODE = DiscontinuousTestODE(**{}) - assert childODE.t_switch_exact == parentODE.t_switch_exact, f"Exact event times between classes does not match!" + assert childODE.t_switch_exact == parentODE.t_switch_exact, "Exact event times between classes does not match!" t0 = 1.0 dt = 0.1 @@ -114,11 +54,11 @@ def testExactDummyProblem(): fExactOde = childODE.eval_f(u0, t0) fOde = parentODE.eval_f(u0, t0) - assert np.allclose(fExactOde, fOde), f"Right-hand sides do not match!" + assert np.allclose(fExactOde, fOde), "Right-hand sides do not match!" fExactOdeEvent = childODE.eval_f(u0Event, tExactEventODE) fOdeEvent = parentODE.eval_f(u0Event, tExactEventODE) - assert np.allclose(fExactOdeEvent, fOdeEvent), f"Right-hand sides at event do not match!" + assert np.allclose(fExactOdeEvent, fOdeEvent), "Right-hand sides at event do not match!" @pytest.mark.base @@ -143,6 +83,7 @@ def testAdaptInterpolationInfo(quad_type): from pySDC.projects.PinTSimE.battery_model import generateDescription from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.implementations.problem_classes.DiscontinuousTestODE import ExactDiscontinuousTestODE problem = ExactDiscontinuousTestODE problem_params = dict() @@ -237,6 +178,7 @@ def testDetectionBoundary(num_nodes): from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.implementations.problem_classes.DiscontinuousTestODE import ExactDiscontinuousTestODE from pySDC.implementations.hooks.log_solution import LogSolution from pySDC.implementations.hooks.log_restarts import LogRestarts from pySDC.helpers.stats_helper import get_sorted @@ -313,6 +255,7 @@ def testDetectionODE(tol, num_nodes, quad_type): """ from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.implementations.problem_classes.DiscontinuousTestODE import ExactDiscontinuousTestODE from pySDC.helpers.stats_helper import get_sorted from pySDC.projects.PinTSimE.battery_model import generateDescription, controllerRun from pySDC.implementations.hooks.log_solution import LogSolution @@ -399,7 +342,7 @@ def testDetectionDAE(num_nodes): from pySDC.projects.PinTSimE.paper_PSCC2024.log_event import LogEventDiscontinuousTestDAE problem = DiscontinuousTestDAE - problem_params = dict({'newton_tol': 1e-6}) + problem_params = {'newton_tol': 1e-6} t0 = 4.6 Tend = 4.62 dt = Tend - t0 From e8bf17e0db988541e1c00433ad2dfa138d6b73e0 Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 18 Apr 2024 06:52:42 +0200 Subject: [PATCH 28/29] Speed up loop using mask --- pySDC/projects/PinTSimE/estimation_check.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pySDC/projects/PinTSimE/estimation_check.py b/pySDC/projects/PinTSimE/estimation_check.py index 7d95358172..843dbba655 100644 --- a/pySDC/projects/PinTSimE/estimation_check.py +++ b/pySDC/projects/PinTSimE/estimation_check.py @@ -244,8 +244,8 @@ def plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix): # pragma: no cov t_switches = [u_num[dt][M_fix][use_SE][use_A]['t_switches'] for dt in dt_list] for t_switch_item in t_switches: for m in range(1, len(t_switch_item)): - if np.isclose(t_switch_item[m - 1], t_switch_item[m], atol=1e-10): - t_switch_item.pop(m - 1) + mask = np.append([True], np.abs(t_switch_item[1:] - t_switch_item[:-1]) > 1e-10) + t_switch_item = t_switch_item[mask] t_switch = [t_event[i] for t_event in t_switches] ax[0, ind].plot( From 2affa7560cd2d5b98a830197395fc8e73a0ca261 Mon Sep 17 00:00:00 2001 From: lisawim Date: Thu, 18 Apr 2024 06:55:56 +0200 Subject: [PATCH 29/29] Removed loop --- pySDC/projects/PinTSimE/estimation_check.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pySDC/projects/PinTSimE/estimation_check.py b/pySDC/projects/PinTSimE/estimation_check.py index 843dbba655..de05b8d6b3 100644 --- a/pySDC/projects/PinTSimE/estimation_check.py +++ b/pySDC/projects/PinTSimE/estimation_check.py @@ -243,9 +243,8 @@ def plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix): # pragma: no cov if use_SE: t_switches = [u_num[dt][M_fix][use_SE][use_A]['t_switches'] for dt in dt_list] for t_switch_item in t_switches: - for m in range(1, len(t_switch_item)): - mask = np.append([True], np.abs(t_switch_item[1:] - t_switch_item[:-1]) > 1e-10) - t_switch_item = t_switch_item[mask] + mask = np.append([True], np.abs(t_switch_item[1:] - t_switch_item[:-1]) > 1e-10) + t_switch_item = t_switch_item[mask] t_switch = [t_event[i] for t_event in t_switches] ax[0, ind].plot(