diff --git a/pySDC/implementations/problem_classes/Battery.py b/pySDC/implementations/problem_classes/Battery.py index e7b590f1e6..c0b0641804 100644 --- a/pySDC/implementations/problem_classes/Battery.py +++ b/pySDC/implementations/problem_classes/Battery.py @@ -1,7 +1,7 @@ import numpy as np -from pySDC.core.Errors import ProblemError -from pySDC.core.Problem import ptype +from pySDC.core.Errors import ParameterError, ProblemError +from pySDC.core.Problem import ptype, WorkCounter from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh @@ -66,10 +66,37 @@ def __init__(self, ncapacitors=2, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1. nvars = n + 1 if C is None: - C = np.array([1.0, 1.0]) + if ncapacitors == 1: + C = np.array([1.0]) + elif ncapacitors == 2: + C = np.array([1.0, 1.0]) + else: + raise ParameterError(f"No default value for C is set up for ncapacitors={ncapacitors}") + else: + msgC = "ERROR for capacitance C: C has to be an np.ndarray and/or length of array needs to be equal to number of capacitances" + assert all([type(C) == np.ndarray, np.shape(C)[0] == n]), msgC if V_ref is None: - V_ref = np.array([1.0, 1.0]) + if ncapacitors == 1: + V_ref = np.array([1.0]) + elif ncapacitors == 2: + V_ref = np.array([1.0, 1.0]) + else: + raise ParameterError(f"No default value for V_ref is set up for ncapacitors={ncapacitors}") + else: + assertions_V_ref_1 = [ + type(V_ref) == np.ndarray, + np.shape(V_ref)[0] == n, + ] + msg1 = "ERROR for reference value V_ref: V_ref has to be an np.ndarray and/or length of array needs to be equal to number of capacitances" + assert all(assertions_V_ref_1), msg1 + + assertions_V_ref_2 = [ + (alpha > V_ref[k] for k in range(n)), + (V_ref[k] > 0 for k in range(n)), + ] + msg2 = "ERROR for V_ref: At least one of V_ref is less than zero and/or alpha!" + assert all(assertions_V_ref_2), msg2 # invoke super init, passing number of dofs, dtype_u and dtype_f super().__init__(init=(nvars, None, np.dtype('float64'))) @@ -230,6 +257,7 @@ def get_switching_info(self, u, t): switch_detected = False m_guess = -100 break_flag = False + k_detected = 1 for m in range(1, len(u)): for k in range(1, self.nvars): h_prev_node = u[m - 1][k] - self.V_ref[k - 1] @@ -244,9 +272,7 @@ def get_switching_info(self, u, t): if break_flag: break - state_function = ( - [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(len(u))] if switch_detected else [] - ) + state_function = [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(len(u))] return switch_detected, m_guess, state_function @@ -303,9 +329,12 @@ class battery(battery_n_capacitors): ---- This class has the same attributes as the class it inherits from. """ - dtype_f = imex_mesh + def __init__(self, ncapacitors=1, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1.2, V_ref=None): + """Initialization routine""" + super().__init__(ncapacitors, Vs, Rs, C, R, L, alpha, V_ref) + def eval_f(self, u, t): """ Routine to evaluate the right-hand side of the problem. @@ -423,10 +452,8 @@ class battery_implicit(battery): Attributes ---------- - newton_itercount: int - Counts the number of Newton iterations. - newton_ncalls: int - Counts the number of how often Newton is called in the simulation of the problem. + work_counters : WorkCounter + Counts different things, here: Number of Newton iterations is counted. """ dtype_f = mesh @@ -443,17 +470,10 @@ def __init__( newton_maxiter=100, newton_tol=1e-11, ): - if C is None: - C = np.array([1.0]) - - if V_ref is None: - V_ref = np.array([1.0]) - super().__init__(ncapacitors, Vs, Rs, C, R, L, alpha, V_ref) self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', localVars=locals(), readOnly=True) - self.newton_itercount = 0 - self.newton_ncalls = 0 + self.work_counters['newton'] = WorkCounter() def eval_f(self, u, t): """ @@ -544,6 +564,7 @@ def solve_system(self, rhs, factor, u0, t): # increase iteration count n += 1 + self.work_counters['newton']() if np.isnan(res) and self.stop_at_nan: raise ProblemError('Newton got nan after %i iterations, aborting...' % n) @@ -553,9 +574,6 @@ def solve_system(self, rhs, factor, u0, t): if n == self.newton_maxiter: self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res)) - self.newton_ncalls += 1 - self.newton_itercount += n - me = self.dtype_u(self.init) me[:] = u[:] diff --git a/pySDC/playgrounds/EnergyGrids/battery_adaptivity.py b/pySDC/playgrounds/EnergyGrids/battery_adaptivity.py index c8de68f799..60090e7750 100644 --- a/pySDC/playgrounds/EnergyGrids/battery_adaptivity.py +++ b/pySDC/playgrounds/EnergyGrids/battery_adaptivity.py @@ -2,15 +2,19 @@ import dill from pySDC.helpers.stats_helper import get_sorted -from pySDC.implementations.collocations import Collocation -from pySDC.projects.PinTSimE.switch_controller_nonMPI import switch_controller_nonMPI +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.implementations.problem_classes.Battery import battery from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order from pySDC.playgrounds.EnergyGrids.log_data_battery import log_data_battery -from pySDC.projects.PinTSimE.piline_model import setup_mpl import pySDC.helpers.plot_helper as plt_helper +from pySDC.implementations.hooks.log_solution import LogSolution + +from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator +from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity +from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + def main(): """ @@ -18,55 +22,65 @@ def main(): """ # initialize level parameters - level_params = dict() - level_params['restol'] = -1e-10 - level_params['e_tol'] = 1e-5 - level_params['dt'] = 2e-2 + level_params = { + 'restol': -1, + 'dt': 2e-2, + } # initialize sweeper parameters - sweeper_params = dict() - - sweeper_params['quad_type'] = 'LOBATTO' - sweeper_params['num_nodes'] = 5 - sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part - sweeper_params['initial_guess'] = 'spread' - - # initialize problem parameters - problem_params = dict() - problem_params['Vs'] = 5.0 - problem_params['Rs'] = 0.5 - problem_params['C'] = 1 - problem_params['R'] = 1 - problem_params['L'] = 1 - problem_params['alpha'] = 3 # 10 - problem_params['V_ref'] = 1 + sweeper_params = { + 'quad_type': 'LOBATTO', + 'num_nodes': 5, + 'QI': 'LU', + 'initial_guess': 'spread', + } # initialize step parameters - step_params = dict() - step_params['maxiter'] = 5 + step_params = { + 'maxiter': 10, + } # initialize controller parameters - controller_params = dict() - controller_params['use_adaptivity'] = True - controller_params['use_switch_estimator'] = True - controller_params['logger_level'] = 20 - controller_params['hook_class'] = log_data_battery + controller_params = { + 'logger_level': 30, + 'hook_class': [LogSolution], + 'mssdc_jac': False, + } + + # convergence controllers + convergence_controllers = {} + switch_estimator_params = { + 'tol': 1e-10, + 'alpha': 1.0, + } + convergence_controllers.update({SwitchEstimator: switch_estimator_params}) + adaptivity_params = { + 'e_tol': 1e-7, + } + convergence_controllers.update({Adaptivity: adaptivity_params}) + restarting_params = { + 'max_restarts': 50, + 'crash_after_max_restarts': False, + } + convergence_controllers.update({BasicRestartingNonMPI: restarting_params}) # fill description dictionary for easy step instantiation - description = dict() - description['problem_class'] = battery - description['problem_params'] = problem_params - description['sweeper_class'] = imex_1st_order - description['sweeper_params'] = sweeper_params - description['level_params'] = level_params - description['step_params'] = step_params + description = { + 'problem_class': battery, + 'problem_params': {}, # use default problem params + 'sweeper_class': imex_1st_order, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } # set time parameters t0 = 0.0 - Tend = 4 + Tend = 0.5 # instantiate controller - controller = switch_controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + 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 @@ -91,19 +105,17 @@ def plot_voltages(cwd='./'): stats = dill.load(f) f.close() - # convert filtered statistics to list of iterations count, sorted by process - cL = get_sorted(stats, type='current L', sortby='time', recomputed=False) - vC = get_sorted(stats, type='voltage C', sortby='time', recomputed=False) - - times = [v[0] for v in cL] + u = get_sorted(stats, type='u', sortby='time', recomputed=False) + t = np.array([me[0] for me in u]) + iL = np.array([me[1][0] for me in u]) + vC = np.array([me[1][1] for me in u]) dt = np.array(get_sorted(stats, type='dt', recomputed=False)) list_gs = get_sorted(stats, type='restart') - setup_mpl() fig, ax = plt_helper.plt.subplots(1, 1) - ax.plot(times, [v[1] for v in cL], label='$i_L$') - ax.plot(times, [v[1] for v in vC], label='$v_C$') + ax.plot(t, iL, label='$i_L$') + ax.plot(t, vC, label='$v_C$') ax.set_xlabel('Time', fontsize=20) ax.set_ylabel('Energy', fontsize=20) for element in list_gs: diff --git a/pySDC/playgrounds/EnergyGrids/log_data.py b/pySDC/playgrounds/EnergyGrids/log_data.py deleted file mode 100644 index 87ad2b5a48..0000000000 --- a/pySDC/playgrounds/EnergyGrids/log_data.py +++ /dev/null @@ -1,40 +0,0 @@ -from pySDC.core.Hooks import hooks - - -class log_data(hooks): - def post_step(self, step, level_number): - - super(log_data, self).post_step(step, level_number) - - # some abbreviations - L = step.levels[level_number] - - L.sweep.compute_end_point() - - self.add_to_stats( - process=step.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='v1', - value=L.uend[0], - ) - self.add_to_stats( - process=step.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='v2', - value=L.uend[1], - ) - self.add_to_stats( - process=step.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='p3', - value=L.uend[2], - ) diff --git a/pySDC/playgrounds/EnergyGrids/log_data_battery.py b/pySDC/playgrounds/EnergyGrids/log_data_battery.py deleted file mode 100644 index 6579c1fdae..0000000000 --- a/pySDC/playgrounds/EnergyGrids/log_data_battery.py +++ /dev/null @@ -1,50 +0,0 @@ -from pySDC.core.Hooks import hooks - - -class log_data_battery(hooks): - def post_step(self, step, level_number): - - super(log_data_battery, self).post_step(step, level_number) - - # some abbreviations - L = step.levels[level_number] - - L.sweep.compute_end_point() - - self.add_to_stats( - process=step.status.slot, - time=L.time + L.dt, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='current L', - value=L.uend[0], - ) - self.add_to_stats( - process=step.status.slot, - time=L.time + L.dt, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='voltage C', - value=L.uend[1], - ) - self.increment_stats( - process=step.status.slot, - time=L.time + L.dt, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='restart', - value=1, - initialize=0, - ) - self.add_to_stats( - process=step.status.slot, - time=L.time + L.dt, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='dt', - value=L.dt, - ) diff --git a/pySDC/playgrounds/EnergyGrids/playground.py b/pySDC/playgrounds/EnergyGrids/playground.py index 7517d9148b..a8730108a4 100644 --- a/pySDC/playgrounds/EnergyGrids/playground.py +++ b/pySDC/playgrounds/EnergyGrids/playground.py @@ -2,12 +2,11 @@ import dill from pySDC.helpers.stats_helper import get_sorted -from pySDC.implementations.collocations import Collocation from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.implementations.problem_classes.Piline import piline from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -from pySDC.playgrounds.EnergyGrids.log_data import log_data +from pySDC.implementations.hooks.log_solution import LogSolution import pySDC.helpers.plot_helper as plt_helper @@ -45,7 +44,7 @@ def main(): # initialize controller parameters controller_params = dict() controller_params['logger_level'] = 20 - controller_params['hook_class'] = log_data + controller_params['hook_class'] = [LogSolution] # fill description dictionary for easy step instantiation description = dict() @@ -82,16 +81,16 @@ def plot_voltages(cwd='./'): f.close() # convert filtered statistics to list of iterations count, sorted by process - v1 = get_sorted(stats, type='v1', sortby='time') - v2 = get_sorted(stats, type='v2', sortby='time') - p3 = get_sorted(stats, type='p3', sortby='time') - - times = [v[0] for v in v1] + u_val = get_sorted(stats, type='u', sortby='time') + t = np.array([me[0] for me in u_val]) + v1 = np.array([me[1][0] for me in u_val]) + v2 = np.array([me[1][1] for me in u_val]) + p3 = np.array([me[1][2] for me in u_val]) # plt_helper.setup_mpl() - plt_helper.plt.plot(times, [v[1] for v in v1], label='v1') - plt_helper.plt.plot(times, [v[1] for v in v2], label='v2') - plt_helper.plt.plot(times, [v[1] for v in p3], label='p3') + plt_helper.plt.plot(t, v1, label='v1') + plt_helper.plt.plot(t, v2, label='v2') + plt_helper.plt.plot(t, p3, label='p3') plt_helper.plt.legend() plt_helper.plt.show() diff --git a/pySDC/playgrounds/EnergyGrids/playground_battery.py b/pySDC/playgrounds/EnergyGrids/playground_battery.py deleted file mode 100644 index 968d0541b1..0000000000 --- a/pySDC/playgrounds/EnergyGrids/playground_battery.py +++ /dev/null @@ -1,142 +0,0 @@ -import numpy as np -import dill - -from pySDC.helpers.stats_helper import get_sorted - -# from pySDC.helpers.visualization_tools import show_residual_across_simulation -from pySDC.implementations.collocations import Collocation -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.implementations.problem_classes.Battery import battery -from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order - -from pySDC.playgrounds.EnergyGrids.log_data_battery import log_data_battery -import pySDC.helpers.plot_helper as plt_helper - - -def main(): - """ - A simple test program to do SDC/PFASST runs for the battery drain model - """ - - # initialize level parameters - level_params = dict() - level_params['restol'] = -1e-10 - level_params['e_tol'] = 1e-5 - level_params['dt'] = 2e-2 - - # initialize sweeper parameters - sweeper_params = dict() - - sweeper_params['quad_type'] = 'LOBATTO' - sweeper_params['num_nodes'] = 5 - sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part - sweeper_params['initial_guess'] = 'spread' - - # initialize problem parameters - problem_params = dict() - problem_params['Vs'] = 5.0 - problem_params['Rs'] = 0.5 - problem_params['C'] = 1 - problem_params['R'] = 1 - problem_params['L'] = 1 - problem_params['alpha'] = 10 - problem_params['V_ref'] = 1 - - # initialize step parameters - step_params = dict() - step_params['maxiter'] = 5 - - # initialize controller parameters - controller_params = dict() - controller_params['use_adaptivity'] = False - controller_params['use_switch_estimator'] = True - controller_params['logger_level'] = 20 - controller_params['hook_class'] = log_data_battery - - # fill description dictionary for easy step instantiation - description = dict() - description['problem_class'] = battery # pass problem class - description['problem_params'] = problem_params # pass problem parameters - description['sweeper_class'] = imex_1st_order # pass sweeper - description['sweeper_params'] = sweeper_params # pass sweeper parameters - description['level_params'] = level_params # pass level parameters - description['step_params'] = step_params - - # set time parameters - t0 = 0.0 - Tend = 2.4 - - # 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) - - # call main function to get things done... - uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - - # fname = 'data/piline.dat' - fname = 'battery.dat' - f = open(fname, 'wb') - dill.dump(stats, f) - f.close() - - -def plot_voltages(cwd='./'): - """ - Routine to plot the numerical solution of the model - """ - - f = open(cwd + 'battery.dat', 'rb') - stats = dill.load(f) - f.close() - - # convert filtered statistics to list of iterations count, sorted by process - cL = get_sorted(stats, type='current L', sortby='time') - vC = get_sorted(stats, type='voltage C', sortby='time') - - times = [v[0] for v in cL] - - # plt_helper.setup_mpl() - plt_helper.plt.plot(times, [v[1] for v in cL], label='current L') - plt_helper.plt.plot(times, [v[1] for v in vC], label='voltage C') - plt_helper.plt.legend() - - plt_helper.plt.show() - - -def plot_residuals(cwd='./'): - """ - Routine to plot the residuals for one block of each iteration and each process - """ - - f = open(cwd + 'battery.dat', 'rb') - stats = dill.load(f) - f.close() - - # filter statistics by number of iterations - iter_counts = get_sorted(stats, type='niter', sortby='time') - - # compute and print statistics - min_iter = 20 - max_iter = 0 - f = open('battery_out.txt', 'w') - for item in iter_counts: - out = 'Number of iterations for time %4.2f: %1i' % item - f.write(out + '\n') - print(out) - min_iter = min(min_iter, item[1]) - max_iter = max(max_iter, item[1]) - f.close() - - # call helper routine to produce residual plot - - fname = 'battery_residuals.png' - show_residual_across_simulation(stats=stats, fname=fname) - - -if __name__ == "__main__": - main() - plot_voltages() - # plot_residuals() diff --git a/pySDC/playgrounds/EnergyGrids/playground_buck.py b/pySDC/playgrounds/EnergyGrids/playground_buck.py index 5de076edf0..26272e7c78 100644 --- a/pySDC/playgrounds/EnergyGrids/playground_buck.py +++ b/pySDC/playgrounds/EnergyGrids/playground_buck.py @@ -1,14 +1,12 @@ import numpy as np import dill -from scipy.integrate import solve_ivp from pySDC.helpers.stats_helper import get_sorted -from pySDC.implementations.collocations import Collocation from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.implementations.problem_classes.BuckConverter import buck_converter from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -from pySDC.playgrounds.EnergyGrids.log_data import log_data +from pySDC.implementations.hooks.log_solution import LogSolution import pySDC.helpers.plot_helper as plt_helper @@ -48,7 +46,7 @@ def main(): # initialize controller parameters controller_params = dict() controller_params['logger_level'] = 20 - controller_params['hook_class'] = log_data + controller_params['hook_class'] = [LogSolution] # fill description dictionary for easy step instantiation description = dict() @@ -86,16 +84,16 @@ def plot_voltages(cwd='./'): f.close() # convert filtered statistics to list of iterations count, sorted by process - v1 = get_sorted(stats, type='v1', sortby='time') - v2 = get_sorted(stats, type='v2', sortby='time') - p3 = get_sorted(stats, type='p3', sortby='time') - - times = [v[0] for v in v1] + u_val = get_sorted(stats, type='u', sortby='time') + t = np.array([me[0] for me in u_val]) + v1 = np.array([me[1][0] for me in u_val]) + v2 = np.array([me[1][1] for me in u_val]) + p3 = np.array([me[1][2] for me in u_val]) # plt_helper.setup_mpl() - plt_helper.plt.plot(times, [v[1] for v in v1], label='v1') - plt_helper.plt.plot(times, [v[1] for v in v2], label='v2') - plt_helper.plt.plot(times, [v[1] for v in p3], label='p3') + plt_helper.plt.plot(t, v1, label='v1') + plt_helper.plt.plot(t, v2, label='v2') + plt_helper.plt.plot(t, p3, label='p3') plt_helper.plt.legend() plt_helper.plt.show() diff --git a/pySDC/projects/PinTSimE/battery_2capacitors_model.py b/pySDC/projects/PinTSimE/battery_2capacitors_model.py deleted file mode 100644 index c1769e232b..0000000000 --- a/pySDC/projects/PinTSimE/battery_2capacitors_model.py +++ /dev/null @@ -1,306 +0,0 @@ -import numpy as np -import dill -from pathlib import Path - -from pySDC.helpers.stats_helper import get_sorted -from pySDC.implementations.problem_classes.Battery import battery_n_capacitors -from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.projects.PinTSimE.battery_model import ( - controller_run, - generate_description, - get_recomputed, - proof_assertions_description, -) -from pySDC.projects.PinTSimE.piline_model import setup_mpl -import pySDC.helpers.plot_helper as plt_helper - -from pySDC.core.Hooks import hooks -from pySDC.implementations.hooks.log_solution import LogSolution -from pySDC.implementations.hooks.default_hook import DefaultHooks - -from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator - - -class LogEvent(hooks): - """ - Logs the problem dependent state function of the battery drain model. - """ - - def post_step(self, step, level_number): - super(LogEvent, self).post_step(step, level_number) - - L = step.levels[level_number] - P = L.prob - - L.sweep.compute_end_point() - - self.add_to_stats( - process=step.status.slot, - time=L.time + L.dt, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='state_function_1', - value=L.uend[1] - P.V_ref[0], - ) - self.add_to_stats( - process=step.status.slot, - time=L.time + L.dt, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='state_function_2', - value=L.uend[2] - P.V_ref[1], - ) - - -def run(): - """ - Executes the simulation for the battery model using the IMEX sweeper and plot the results - as _model_solution_.png - """ - - dt = 1e-2 - t0 = 0.0 - Tend = 3.5 - - problem_classes = [battery_n_capacitors] - sweeper_classes = [imex_1st_order] - num_nodes = 4 - restol = -1 - maxiter = 8 - - ncapacitors = 2 - alpha = 5.0 - V_ref = np.array([1.0, 1.0]) - C = np.array([1.0, 1.0]) - - problem_params = dict() - problem_params['ncapacitors'] = ncapacitors - problem_params['C'] = C - problem_params['alpha'] = alpha - problem_params['V_ref'] = V_ref - - recomputed = False - use_switch_estimator = [True] - max_restarts = 1 - tol_event = 1e-8 - - hook_class = [DefaultHooks, LogSolution, LogEvent] - - for problem, sweeper in zip(problem_classes, sweeper_classes): - for use_SE in use_switch_estimator: - description, controller_params = generate_description( - dt, - problem, - sweeper, - num_nodes, - hook_class, - False, - use_SE, - problem_params, - restol, - maxiter, - max_restarts, - tol_event, - ) - - # Assertions - proof_assertions_description(description, False, use_SE) - - proof_assertions_time(dt, Tend, V_ref, alpha) - - stats = controller_run(description, controller_params, False, use_SE, t0, Tend) - - check_solution(stats, dt, use_SE) - - plot_voltages(description, problem.__name__, sweeper.__name__, recomputed, use_SE, False) - - -def plot_voltages(description, problem, sweeper, recomputed, use_switch_estimator, use_adaptivity, cwd='./'): - """ - Routine to plot the numerical solution of the model. - - Parameters - ---------- - description : dict - Contains all information for a controller run. - problem : pySDC.core.Problem.ptype - Problem class that wants to be simulated. - sweeper : pySDC.core.Sweeper.sweeper - Sweeper class for solving the problem class numerically. - recomputed : bool - Flag if the values after a restart are used or before. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. - use_adaptivity : bool - Flag if adaptivity wants to be used or not. - cwd : str - Current working directory. - """ - - f = open(cwd + 'data/{}_{}_USE{}_USA{}.dat'.format(problem, sweeper, use_switch_estimator, use_adaptivity), 'rb') - stats = dill.load(f) - f.close() - - # convert filtered statistics to list of iterations count, sorted by process - cL = np.array([me[1][0] for me in get_sorted(stats, type='u', recomputed=recomputed)]) - vC1 = np.array([me[1][1] for me in get_sorted(stats, type='u', recomputed=recomputed)]) - vC2 = np.array([me[1][2] for me in get_sorted(stats, type='u', recomputed=recomputed)]) - - t = np.array([me[0] for me in get_sorted(stats, type='u', recomputed=recomputed)]) - - setup_mpl() - fig, ax = plt_helper.plt.subplots(1, 1, figsize=(4.5, 3)) - ax.plot(t, cL, label='$i_L$') - ax.plot(t, vC1, label='$v_{C_1}$') - ax.plot(t, vC2, label='$v_{C_2}$') - - if use_switch_estimator: - switches = get_recomputed(stats, type='switch', sortby='time') - if recomputed is not None: - assert len(switches) >= 2, f"Expected at least 2 switches, got {len(switches)}!" - t_switches = [v[1] for v in switches] - - for i in range(len(t_switches)): - ax.axvline(x=t_switches[i], linestyle='--', color='k', label='Switch {}'.format(i + 1)) - - ax.legend(frameon=False, fontsize=12, loc='upper right') - - ax.set_xlabel('Time') - ax.set_ylabel('Energy') - - fig.savefig('data/battery_2capacitors_model_solution.png', dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig) - - -def check_solution(stats, dt, use_switch_estimator): - """ - Function that checks the solution based on a hardcoded reference solution. - Based on check_solution function from @brownbaerchen. - - Parameters - ---------- - stats : dict - Raw statistics from a controller run. - dt : float - Initial time step. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. - """ - - data = get_data_dict(stats, use_switch_estimator) - - if use_switch_estimator: - msg = f'Error when using the switch estimator for battery_2capacitors for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 1.1783297877614183, - 'vC1': 0.9999999999967468, - 'vC2': 0.999999999996747, - 'state_function_1': -3.2531755067566337e-12, - 'state_function_2': -3.2529534621517087e-12, - 'restarts': 2.0, - 'sum_niters': 2824.0, - } - elif dt == 4e-1: - expected = { - 'cL': 1.5039617338098907, - 'vC1': 0.9999999968387812, - 'vC2': 0.9999999968387812, - 'state_function_1': -3.161218842251401e-09, - 'state_function_2': -3.161218842251401e-09, - 'restarts': 10.0, - 'sum_niters': 200.0, - } - elif dt == 4e-2: - expected = { - 'cL': 1.2707220273133215, - 'vC1': 1.0000000041344774, - 'vC2': 0.999999999632751, - 'state_function_1': 4.134477427086836e-09, - 'state_function_2': -3.672490089812186e-10, - 'restarts': 6.0, - 'sum_niters': 792.0, - } - - got = { - 'cL': data['cL'][-1], - 'vC1': data['vC1'][-1], - 'vC2': data['vC2'][-1], - 'state_function_1': data['state_function_1'][-1], - 'state_function_2': data['state_function_2'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } - - for key in expected.keys(): - assert np.isclose( - expected[key], got[key], rtol=1e-4 - ), f'{msg} Expected {key}={expected[key]:.4e}, got {key}={got[key]:.4e}' - - -def get_data_dict(stats, use_switch_estimator, recomputed=False): - """ - Converts the statistics in a useful data dictionary so that it can be easily checked in the check_solution function. - Based on @brownbaerchen's get_data function. - - Parameters - ---------- - stats : dict - Raw statistics from a controller run. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. - recomputed : bool - Flag if the values after a restart are used or before. - - Returns - ------- - data : dict - Contains all information as the statistics dict. - """ - - data = dict() - data['cL'] = np.array([me[1][0] for me in get_sorted(stats, type='u', recomputed=False, sortby='time')]) - data['vC1'] = np.array([me[1][1] for me in get_sorted(stats, type='u', recomputed=False, sortby='time')]) - data['vC2'] = np.array([me[1][2] for me in get_sorted(stats, type='u', recomputed=False, sortby='time')]) - data['state_function_1'] = np.array(get_sorted(stats, type='state_function_1', sortby='time', recomputed=False))[ - :, 1 - ] - data['state_function_2'] = np.array(get_sorted(stats, type='state_function_2', sortby='time', recomputed=False))[ - :, 1 - ] - data['restarts'] = np.sum(np.array(get_sorted(stats, type='restart', recomputed=None, sortby='time'))[:, 1]) - data['sum_niters'] = np.sum(np.array(get_sorted(stats, type='niter', recomputed=None, sortby='time'))[:, 1]) - - return data - - -def proof_assertions_time(dt, Tend, V_ref, alpha): - """ - Function to proof the assertions regarding the time domain (in combination with the specific problem). - - Parameters - ---------- - dt : float - Time step for computation. - Tend : float - End time. - V_ref : np.ndarray - Reference values (problem parameter). - alpha : np.float - Multiple used for initial conditions (problem_parameter). - """ - - assert ( - Tend == 3.5 and V_ref[0] == 1.0 and V_ref[1] == 1.0 and alpha == 5.0 - ), "Error! Do not use other parameters for V_ref[:] != 1.0, alpha != 1.2, Tend != 0.3 due to hardcoded reference!" - - assert ( - dt == 1e-2 or dt == 4e-1 or dt == 4e-2 - ), "Error! Do not use other time steps dt != 4e-1 or dt != 4e-2 or dt != 4e-3 due to hardcoded references!" - - -if __name__ == "__main__": - run() diff --git a/pySDC/projects/PinTSimE/battery_model.py b/pySDC/projects/PinTSimE/battery_model.py index 7de772ba4f..20f2e11f1c 100644 --- a/pySDC/projects/PinTSimE/battery_model.py +++ b/pySDC/projects/PinTSimE/battery_model.py @@ -1,14 +1,12 @@ import numpy as np -import dill from pathlib import Path from pySDC.helpers.stats_helper import sort_stats, filter_stats, get_sorted -from pySDC.implementations.problem_classes.Battery import battery, battery_implicit +from pySDC.implementations.problem_classes.Battery import battery, battery_implicit, battery_n_capacitors from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.projects.PinTSimE.piline_model import setup_mpl import pySDC.helpers.plot_helper as plt_helper from pySDC.core.Hooks import hooks @@ -20,14 +18,16 @@ from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI +from pySDC.projects.PinTSimE.hardcoded_solutions import testSolution -class LogEvent(hooks): + +class LogEventBattery(hooks): """ Logs the problem dependent state function of the battery drain model. """ def post_step(self, step, level_number): - super(LogEvent, self).post_step(step, level_number) + super().post_step(step, level_number) L = step.levels[level_number] P = L.prob @@ -41,15 +41,17 @@ def post_step(self, step, level_number): iter=0, sweep=L.status.sweep, type='state_function', - value=L.uend[1] - P.V_ref[0], + value=L.uend[1:] - P.V_ref[:], ) -def generate_description( +def generateDescription( dt, problem, sweeper, num_nodes, + quad_type, + QI, hook_class, use_adaptivity, use_switch_estimator, @@ -58,30 +60,32 @@ def generate_description( maxiter, max_restarts=None, tol_event=1e-10, + alpha=1.0, ): - """ + r""" Generate a description for the battery models for a controller run. Parameters ---------- dt : float Time step for computation. - problem : pySDC.core.Problem.ptype + problem : pySDC.core.Problem Problem class that wants to be simulated. - sweeper : pySDC.core.Sweeper.sweeper + sweeper : pySDC.core.Sweeper Sweeper class for solving the problem class numerically. num_nodes : int Number of collocation nodes. - hook_class : pySDC.core.Hooks - Logged data for a problem. + quad_type : str + Type of quadrature nodes, e.g. ``'LOBATTO'`` or ``'RADAU-RIGHT'``. + QI : str + Type of preconditioner used in SDC, e.g. ``'IE'`` or ``'LU'``. + hook_class : List of pySDC.core.Hooks + Logged data for a problem, e.g., hook_class=[LogSolution] logs the solution ``'u'`` + during the simulation. use_adaptivity : bool Flag if the adaptivity wants to be used or not. use_switch_estimator : bool Flag if the switch estimator wants to be used or not. - ncapacitors : int - Number of capacitors used for the battery_model. - alpha : float - Multiple used for the initial conditions (problem_parameter). problem_params : dict Dictionary containing the problem parameters. restol : float @@ -91,7 +95,9 @@ def generate_description( max_restarts : int, optional Maximum number of restarts per step. tol_event : float, optional - Tolerance for switch estimation to terminate. + Tolerance for event detection to terminate. + alpha : float, optional + Factor that indicates how the new step size in the Switch Estimator is reduced. Returns ------- @@ -102,59 +108,69 @@ def generate_description( """ # initialize level parameters - level_params = dict() - level_params['restol'] = -1 if use_adaptivity else restol - level_params['dt'] = dt + level_params = { + 'restol': -1 if use_adaptivity else restol, + 'dt': dt, + } + if use_adaptivity: + assert restol == -1, "Please set restol to -1 or omit it" # initialize sweeper parameters - sweeper_params = dict() - sweeper_params['quad_type'] = 'LOBATTO' - sweeper_params['num_nodes'] = num_nodes - sweeper_params['QI'] = 'IE' - sweeper_params['initial_guess'] = 'spread' + sweeper_params = { + 'quad_type': quad_type, + 'num_nodes': num_nodes, + 'QI': QI, + 'initial_guess': 'spread', + } # initialize step parameters - step_params = dict() - step_params['maxiter'] = maxiter + step_params = { + 'maxiter': maxiter, + } + assert 'errtol' not in step_params.keys(), 'No exact solution known to compute error' # initialize controller parameters - controller_params = dict() - controller_params['logger_level'] = 30 - controller_params['hook_class'] = hook_class - controller_params['mssdc_jac'] = False + controller_params = { + 'logger_level': 30, + 'hook_class': hook_class, + 'mssdc_jac': False, + } # convergence controllers - convergence_controllers = dict() + convergence_controllers = {} if use_switch_estimator: - switch_estimator_params = {} - switch_estimator_params['tol'] = tol_event + switch_estimator_params = { + 'tol': tol_event, + 'alpha': alpha, + } convergence_controllers.update({SwitchEstimator: switch_estimator_params}) - if use_adaptivity: - adaptivity_params = dict() - adaptivity_params['e_tol'] = 1e-7 + adaptivity_params = { + 'e_tol': 1e-7, + } convergence_controllers.update({Adaptivity: adaptivity_params}) - if max_restarts is not None: - convergence_controllers[BasicRestartingNonMPI] = { + restarting_params = { 'max_restarts': max_restarts, 'crash_after_max_restarts': False, } + convergence_controllers.update({BasicRestartingNonMPI: restarting_params}) # fill description dictionary for easy step instantiation - description = dict() - description['problem_class'] = problem # pass problem class - description['problem_params'] = problem_params # pass problem parameters - description['sweeper_class'] = sweeper # pass sweeper - description['sweeper_params'] = sweeper_params # pass sweeper parameters - description['level_params'] = level_params # pass level parameters - description['step_params'] = step_params - description['convergence_controllers'] = convergence_controllers + description = { + 'problem_class': problem, + 'problem_params': problem_params, + 'sweeper_class': sweeper, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } return description, controller_params -def controller_run(description, controller_params, use_adaptivity, use_switch_estimator, t0, Tend): +def controllerRun(description, controller_params, t0, Tend, exact_event_time_avail=False): """ Executes a controller run for a problem defined in the description. @@ -164,10 +180,12 @@ def controller_run(description, controller_params, use_adaptivity, use_switch_es Contains all information for a controller run. controller_params : dict Parameters needed for a controller run. - use_adaptivity : bool - Flag if the adaptivity wants to be used or not. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. + t0 : float + Starting time of simulation. + Tend : float + End time of simulation. + exact_event_time_avail : bool, optional + Indicates if exact event time of a problem is available. Returns ------- @@ -181,504 +199,402 @@ def controller_run(description, controller_params, use_adaptivity, use_switch_es # 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 if exact_event_time_avail else None # call main function to get things done... uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - problem = description['problem_class'] - sweeper = description['sweeper_class'] + return stats, t_switch_exact - Path("data").mkdir(parents=True, exist_ok=True) - fname = 'data/{}_{}_USE{}_USA{}.dat'.format( - problem.__name__, sweeper.__name__, use_switch_estimator, use_adaptivity - ) - f = open(fname, 'wb') - dill.dump(stats, f) - f.close() - return stats +def main(): + r""" + Executes the simulation. - -def run(): - """ - Executes the simulation for the battery model using two different sweepers and plot the results - as _model_solution_.png + Note + ---- + Hardcoded solutions for battery models in `pySDC.projects.PinTSimE.hardcoded_solutions` are only computed for + ``dt_list=[1e-2, 1e-3]`` and ``M_fix=4``. Hence changing ``dt_list`` and ``M_fix`` to different values could arise + an ``AssertionError``. """ - dt = 1e-2 - t0 = 0.0 - Tend = 0.3 - - problem_classes = [battery, battery_implicit] - sweeper_classes = [imex_1st_order, generic_implicit] - num_nodes = 4 - restol = -1 - maxiter = 8 - - ncapacitors = 1 - alpha = 1.2 - V_ref = np.array([1.0]) - C = np.array([1.0]) - - problem_params = dict() - problem_params['ncapacitors'] = ncapacitors - problem_params['C'] = C - problem_params['alpha'] = alpha - problem_params['V_ref'] = V_ref - - max_restarts = 1 - recomputed = False - use_switch_estimator = [True] + # defines parameters for sweeper + M_fix = 4 + sweeper_params = { + 'num_nodes': M_fix, + 'quad_type': 'LOBATTO', + 'QI': 'IE', + } + + # defines parameters for event detection, restol, and max. number of iterations + handling_params = { + 'restol': -1, + 'maxiter': 8, + 'max_restarts': 50, + 'recomputed': False, + 'tol_event': 1e-10, + 'alpha': 1.0, + 'exact_event_time_avail': None, + } + + all_params = { + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, + } + + hook_class = [LogSolution, LogEventBattery, LogEmbeddedErrorEstimate, LogStepSize] + + use_detection = [True] use_adaptivity = [True] - hook_class = [LogSolution, LogEvent, LogEmbeddedErrorEstimate, LogStepSize] - - for problem, sweeper in zip(problem_classes, sweeper_classes): - for use_SE in use_switch_estimator: - for use_A in use_adaptivity: - tol_event = 1e-10 if sweeper.__name__ == 'generic_implicit' else 1e-17 - - description, controller_params = generate_description( - dt, - problem, - sweeper, - num_nodes, - hook_class, - use_A, - use_SE, - problem_params, - restol, - maxiter, - max_restarts, - tol_event, - ) - - # Assertions - proof_assertions_description(description, use_A, use_SE) - - proof_assertions_time(dt, Tend, V_ref, alpha) - - stats = controller_run(description, controller_params, use_A, use_SE, t0, Tend) - - check_solution(stats, dt, problem.__name__, use_A, use_SE) + for problem, sweeper in zip([battery, battery_implicit], [imex_1st_order, generic_implicit]): + for defaults in [False, True]: + # for hardcoded solutions problem parameter defaults should match with parameters here + if defaults: + params_battery_1capacitor = { + 'ncapacitors': 1, + } + else: + params_battery_1capacitor = { + 'ncapacitors': 1, + 'C': np.array([1.0]), + 'alpha': 1.2, + 'V_ref': np.array([1.0]), + } - plot_voltages(description, problem.__name__, sweeper.__name__, recomputed, use_SE, use_A) + all_params.update({'problem_params': params_battery_1capacitor}) + + _ = runSimulation( + problem=problem, + sweeper=sweeper, + all_params=all_params, + use_adaptivity=use_adaptivity, + use_detection=use_detection, + hook_class=hook_class, + interval=(0.0, 0.3), + dt_list=[1e-2, 1e-3], + nnodes=[M_fix], + ) + + # defines parameters for the problem class + params_battery_2capacitors = { + 'ncapacitors': 2, + 'C': np.array([1.0, 1.0]), + 'alpha': 1.2, + 'V_ref': np.array([1.0, 1.0]), + } + + all_params.update({'problem_params': params_battery_2capacitors}) + + _ = runSimulation( + problem=battery_n_capacitors, + sweeper=imex_1st_order, + all_params=all_params, + use_adaptivity=use_adaptivity, + use_detection=use_detection, + hook_class=hook_class, + interval=(0.0, 0.5), + dt_list=[1e-2, 1e-3], + nnodes=[sweeper_params['num_nodes']], + ) -def plot_voltages(description, problem, sweeper, recomputed, use_switch_estimator, use_adaptivity, cwd='./'): - """ - Routine to plot the numerical solution of the model. +def runSimulation(problem, sweeper, all_params, use_adaptivity, use_detection, hook_class, interval, dt_list, nnodes): + r""" + Script that executes the simulation for a given problem class for given parameters defined by the user. Parameters ---------- - description : dict - Contains all information for a controller run. - problem : pySDC.core.Problem.ptype - Problem class that wants to be simulated. - sweeper : pySDC.core.Sweeper.sweeper - Sweeper class for solving the problem class numerically. - recomputed : bool - Flag if the values after a restart are used or before. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. - use_adaptivity : bool - Flag if adaptivity wants to be used or not. - cwd : str - Current working directory. + problem : pySDC.core.Problem + Problem class to be simulated. + sweeper : pySDC.core.Sweeper + Sweeper that is used to simulate the problem class. + all_params : dict + Dictionary contains the problem parameters for ``problem``, the sweeper parameters for ``sweeper``, + and handling parameters needed for event detection, i.e., ``max_restarts``, ``recomputed``, ``tol_event``, + ``alpha``, and ``exact_event_time_available``. + use_adaptivity : list of bool + Indicates whether adaptivity is used in the simulation or not. Here a list is used to iterate over the + different cases, i.e., ``use_adaptivity=[True, False]``. + use_detection : list of bool + Indicates whether event detection is used in the simulation or not. Here a list is used to iterate over the + different cases, i.e., ``use_detection=[True, False]``. + hook_class : list of pySDC.core.Hooks + List containing the different hook classes to log data during the simulation, i.e., ``hook_class=[LogSolution]`` + logs the solution ``u``. + interval : tuple + Simulation interval. + dt_list : list of float + List containing different step sizes where the solution is computed. + nnodes : list of int + The solution can be computed for different number of collocation nodes. """ - f = open(cwd + 'data/{}_{}_USE{}_USA{}.dat'.format(problem, sweeper, use_switch_estimator, use_adaptivity), 'rb') - stats = dill.load(f) - f.close() + Path("data").mkdir(parents=True, exist_ok=True) - # convert filtered statistics to list of iterations count, sorted by process - cL = np.array([me[1][0] for me in get_sorted(stats, type='u', recomputed=recomputed)]) - vC = np.array([me[1][1] for me in get_sorted(stats, type='u', recomputed=recomputed)]) + prob_cls_name = problem.__name__ - t = np.array([me[0] for me in get_sorted(stats, type='u', recomputed=recomputed)]) + u_num = {} - setup_mpl() - fig, ax = plt_helper.plt.subplots(1, 1, figsize=(3, 3)) - ax.set_title('Simulation of {} using {}'.format(problem, sweeper), fontsize=10) - ax.plot(t, cL, label=r'$i_L$') - ax.plot(t, vC, label=r'$v_C$') + for dt in dt_list: + u_num[dt] = {} - if use_switch_estimator: - switches = get_recomputed(stats, type='switch', sortby='time') + for M in nnodes: + u_num[dt][M] = {} - assert len(switches) >= 1, 'No switches found!' - t_switch = [v[1] for v in switches] + for use_SE in use_detection: + u_num[dt][M][use_SE] = {} - ax.axvline(x=t_switch[-1], linestyle='--', linewidth=0.8, color='r', label='Switch') + for use_A in use_adaptivity: + u_num[dt][M][use_SE][use_A] = {} - if use_adaptivity: - dt = np.array(get_sorted(stats, type='dt', recomputed=False)) + problem_params = all_params['problem_params'] + sweeper_params = all_params['sweeper_params'] + handling_params = all_params['handling_params'] - dt_ax = ax.twinx() - dt_ax.plot(dt[:, 0], dt[:, 1], linestyle='-', linewidth=0.8, color='k', label=r'$\Delta t$') - dt_ax.set_ylabel(r'$\Delta t$', fontsize=8) - dt_ax.legend(frameon=False, fontsize=8, loc='center right') + # plotting results for fixed M requires that M_fix is included in nnodes! + M_fix = sweeper_params['num_nodes'] + assert ( + M_fix in nnodes + ), f"For fixed number of collocation nodes {M_fix} no solution will be computed!" - ax.axhline(y=1.0, linestyle='--', linewidth=0.8, color='g', label='$V_{ref}$') + restol = -1 if use_A else handling_params['restol'] - ax.legend(frameon=False, fontsize=8, loc='upper right') + description, controller_params = generateDescription( + dt=dt, + problem=problem, + sweeper=sweeper, + num_nodes=M, + quad_type=sweeper_params['quad_type'], + QI=sweeper_params['QI'], + hook_class=hook_class, + use_adaptivity=use_A, + use_switch_estimator=use_SE, + problem_params=problem_params, + restol=restol, + maxiter=handling_params['maxiter'], + max_restarts=handling_params['max_restarts'], + tol_event=handling_params['tol_event'], + alpha=handling_params['alpha'], + ) - ax.set_xlabel('Time', fontsize=8) - ax.set_ylabel('Energy', fontsize=8) + stats, t_switch_exact = controllerRun( + description=description, + controller_params=controller_params, + t0=interval[0], + Tend=interval[-1], + exact_event_time_avail=handling_params['exact_event_time_avail'], + ) - fig.savefig('data/{}_model_solution_{}.png'.format(problem, sweeper), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig) + u_num[dt][M][use_SE][use_A] = getDataDict( + stats, prob_cls_name, use_A, use_SE, handling_params['recomputed'], t_switch_exact + ) + plotSolution(u_num[dt][M][use_SE][use_A], prob_cls_name, use_A, use_SE) -def check_solution(stats, dt, problem, use_adaptivity, use_switch_estimator): + testSolution(u_num[dt][M_fix][use_SE][use_A], prob_cls_name, dt, use_A, use_SE) + + return u_num + + +def getUnknownLabels(prob_cls_name): """ - Function that checks the solution based on a hardcoded reference solution. - Based on check_solution function from @brownbaerchen. + Returns the unknown for a problem and corresponding labels for a plot. Parameters ---------- - stats : dict - Raw statistics from a controller run. - dt : float - Initial time step. - problem : problem_class.__name__ - The problem_class that is numerically solved - use_switch_estimator : bool - Indicates if switch detection is used or not. - use_adaptivity : bool - Indicate if adaptivity is used or not. - """ + prob_cls_name : str + Name of the problem class. - data = get_data_dict(stats, use_adaptivity, use_switch_estimator) - - if problem == 'battery': - if use_switch_estimator and use_adaptivity: - msg = f'Error when using switch estimator and adaptivity for battery for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 0.5446532674094873, - 'vC': 0.9999999999883544, - 'dt': 0.01, - 'e_em': 2.220446049250313e-16, - 'state_function': -1.1645573394503117e-11, - 'restarts': 3.0, - 'sum_niters': 136.0, - } - elif dt == 1e-3: - expected = { - 'cL': 0.539386744746365, - 'vC': 0.9999999710472945, - 'dt': 0.005520873635314061, - 'e_em': 2.220446049250313e-16, - 'state_function': -2.8952705455331795e-08, - 'restarts': 11.0, - 'sum_niters': 264.0, - } + Returns + ------- + unknowns : list of str + Contains the names of unknowns. + unknowns_labels : list of str + Contains the labels of unknowns for plotting. + """ - got = { - 'cL': data['cL'][-1], - 'vC': data['vC'][-1], - 'dt': data['dt'][-1], - 'e_em': data['e_em'][-1], - 'state_function': data['state_function'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } - elif use_switch_estimator and not use_adaptivity: - msg = f'Error when using switch estimator for battery for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 0.5456190026495924, - 'vC': 0.999166666941434, - 'state_function': -0.0008333330585660326, - 'restarts': 4.0, - 'sum_niters': 296.0, - } - elif dt == 1e-3: - expected = { - 'cL': 0.5403849766797957, - 'vC': 0.9999166666752302, - 'state_function': -8.33333247698409e-05, - 'restarts': 2.0, - 'sum_niters': 2424.0, - } + unknowns = { + 'battery': ['iL', 'vC'], + 'battery_implicit': ['iL', 'vC'], + 'battery_n_capacitors': ['iL', 'vC1', 'vC2'], + 'DiscontinuousTestODE': ['u'], + 'piline': ['vC1', 'vC2', 'iLp'], + 'buck_converter': ['vC1', 'vC2', 'iLp'], + } - got = { - 'cL': data['cL'][-1], - 'vC': data['vC'][-1], - 'state_function': data['state_function'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } - - elif not use_switch_estimator and use_adaptivity: - msg = f'Error when using adaptivity for battery for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 0.4433805288639916, - 'vC': 0.90262388393713, - 'dt': 0.18137307612335937, - 'e_em': 2.7177844974524135e-09, - 'restarts': 0.0, - 'sum_niters': 24.0, - } - elif dt == 1e-3: - expected = { - 'cL': 0.3994744179584864, - 'vC': 0.9679037468770668, - 'dt': 0.1701392217033212, - 'e_em': 2.0992988458701234e-09, - 'restarts': 0.0, - 'sum_niters': 32.0, - } + unknowns_labels = { + 'battery': [r'$i_L$', r'$v_C$'], + 'battery_implicit': [r'$i_L$', r'$v_C$'], + 'battery_n_capacitors': [r'$i_L$', r'$v_{C_1}$', r'$v_{C_2}$'], + 'DiscontinuousTestODE': [r'$u$'], + 'piline': [r'$v_{C_1}$', r'$v_{C_2}$', r'$i_{L_\pi}$'], + 'buck_converter': [r'$v_{C_1}$', r'$v_{C_2}$', r'$i_{L_\pi}$'], + } - got = { - 'cL': data['cL'][-1], - 'vC': data['vC'][-1], - 'dt': data['dt'][-1], - 'e_em': data['e_em'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } - - elif problem == 'battery_implicit': - if use_switch_estimator and use_adaptivity: - msg = f'Error when using switch estimator and adaptivity for battery_implicit for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 0.5446675396652545, - 'vC': 0.9999999999883541, - 'dt': 0.01, - 'e_em': 2.220446049250313e-16, - 'state_function': -1.1645906461410505e-11, - 'restarts': 3.0, - 'sum_niters': 136.0, - } - elif dt == 1e-3: - expected = { - 'cL': 0.5393867447463223, - 'vC': 0.9999999710472952, - 'dt': 0.005520876908755634, - 'e_em': 2.220446049250313e-16, - 'state_function': -2.895270478919798e-08, - 'restarts': 11.0, - 'sum_niters': 264.0, - } + return unknowns[prob_cls_name], unknowns_labels[prob_cls_name] - got = { - 'cL': data['cL'][-1], - 'vC': data['vC'][-1], - 'dt': data['dt'][-1], - 'e_em': data['e_em'][-1], - 'state_function': data['state_function'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } - elif use_switch_estimator and not use_adaptivity: - msg = f'Error when using switch estimator for battery_implicit for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 0.5456190026495138, - 'vC': 0.9991666669414431, - 'state_function': -0.0008333330585569287, - 'restarts': 4.0, - 'sum_niters': 296.0, - } - elif dt == 1e-3: - expected = { - 'cL': 0.5403849766797896, - 'vC': 0.9999166666752302, - 'state_function': -8.33333247698409e-05, - 'restarts': 2.0, - 'sum_niters': 2424.0, - } - got = { - 'cL': data['cL'][-1], - 'vC': data['vC'][-1], - 'state_function': data['state_function'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } - - elif not use_switch_estimator and use_adaptivity: - msg = f'Error when using adaptivity for battery_implicit for dt={dt:.1e}:' - if dt == 1e-2: - expected = { - 'cL': 0.4694087102919169, - 'vC': 0.9026238839371302, - 'dt': 0.18137307612335937, - 'e_em': 2.3469713394952407e-09, - 'restarts': 0.0, - 'sum_niters': 24.0, - } - elif dt == 1e-3: - expected = { - 'cL': 0.39947441811958956, - 'vC': 0.9679037468770735, - 'dt': 0.1701392217033212, - 'e_em': 1.147640815712947e-09, - 'restarts': 0.0, - 'sum_niters': 32.0, - } +def plotStylingStuff(): # pragma: no cover + """ + Returns plot stuff such as colors, line styles for making plots more pretty. + """ - got = { - 'cL': data['cL'][-1], - 'vC': data['vC'][-1], - 'dt': data['dt'][-1], - 'e_em': data['e_em'][-1], - 'restarts': data['restarts'], - 'sum_niters': data['sum_niters'], - } + colors = { + False: { + False: 'dodgerblue', + True: 'navy', + }, + True: { + False: 'linegreen', + True: 'darkgreen', + }, + } - for key in expected.keys(): - err_msg = f'{msg} Expected {key}={expected[key]:.4e}, got {key}={got[key]:.4e}' - if key == 'cL': - assert abs(expected[key] - got[key]) <= 1e-2, err_msg - else: - assert np.isclose(expected[key], got[key], rtol=1e-3), err_msg + return colors -def get_data_dict(stats, use_adaptivity, use_switch_estimator, recomputed=False): - """ - Converts the statistics in a useful data dictionary so that it can be easily checked in the check_solution function. - Based on @brownbaerchen's get_data function. +def plotSolution(u_num, prob_cls_name, use_adaptivity, use_detection): # pragma: no cover + r""" + Plots the numerical solution for one simulation run. Parameters ---------- - stats : dict - Raw statistics from a controller run. + u_num : dict + Contains numerical solution with corresponding times for different problem_classes, and + labels for different unknowns of the problem. + prob_cls_name : str + Name of the problem class to be plotted. use_adaptivity : bool - Flag if adaptivity wants to be used or not. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. - recomputed : bool - Flag if the values after a restart are used or before. - - Returns - ------- - data : dict - Contains all information as the statistics dict. + Indicates whether adaptivity is used in the simulation or not. """ - data = dict() - data['cL'] = np.array([me[1][0] for me in get_sorted(stats, type='u', sortby='time', recomputed=recomputed)]) - data['vC'] = np.array([me[1][1] for me in get_sorted(stats, type='u', sortby='time', recomputed=recomputed)]) + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5)) + + unknowns = u_num['unknowns'] + unknowns_labels = u_num['unknowns_labels'] + for unknown, unknown_label in zip(unknowns, unknowns_labels): + ax.plot(u_num['t'], u_num[unknown], label=unknown_label) + + if use_detection: + t_switches = u_num['t_switches'] + for i in range(len(t_switches)): + ax.axvline(x=t_switches[i], linestyle='--', linewidth=0.8, color='r', label='Event {}'.format(i + 1)) + if use_adaptivity: - data['dt'] = np.array(get_sorted(stats, type='dt', sortby='time', recomputed=recomputed))[:, 1] - data['e_em'] = np.array( - get_sorted(stats, type='error_embedded_estimate', sortby='time', recomputed=recomputed) - )[:, 1] - if use_switch_estimator: - data['state_function'] = np.array( - get_sorted(stats, type='state_function', sortby='time', recomputed=recomputed) - )[:, 1] - if use_adaptivity or use_switch_estimator: - data['restarts'] = np.sum(np.array(get_sorted(stats, type='restart', recomputed=None, sortby='time'))[:, 1]) - data['sum_niters'] = np.sum(np.array(get_sorted(stats, type='niter', recomputed=None, sortby='time'))[:, 1]) + dt_ax = ax.twinx() + dt = u_num['dt'] + dt_ax.plot(dt[:, 0], dt[:, 1], linestyle='-', linewidth=0.8, color='k', label=r'$\Delta t$') + dt_ax.set_ylabel(r'$\Delta t$', fontsize=16) + dt_ax.legend(frameon=False, fontsize=12, loc='center right') + + ax.legend(frameon=False, fontsize=12, loc='upper right') + ax.set_xlabel(r'$t$', fontsize=16) + ax.set_ylabel(r'$u(t)$', fontsize=16) - return data + fig.savefig('data/{}_model_solution.png'.format(prob_cls_name), dpi=300, bbox_inches='tight') + plt_helper.plt.close(fig) -def get_recomputed(stats, type, sortby): - """ - Function that filters statistics after a recomputation. It stores all value of a type before restart. If there are multiple values - with same time point, it only stores the elements with unique times. +def getDataDict(stats, prob_cls_name, use_adaptivity, use_detection, recomputed, t_switch_exact): + r""" + Extracts statistics and store it in a dictionary. In this routine, from ``stats`` different data are extracted + such as + + - each component of solution ``'u'`` and corresponding time domain ``'t'``, + - the unknowns of the problem ``'unknowns'``, + - the unknowns of the problem as labels for plotting ``'unknowns_labels'``, + - global error ``'e_global'`` after each step, + - events found by event detection ``'t_switches''``, + - exact event time ``'t_switch_exact'``, + - event error ``'e_event'``, + - state function ``'state_function'``, + - embedded error estimate computing when using adaptivity ``'e_em'``, + - (adjusted) step sizes ``'dt'``, + - sum over restarts ``'sum_restarts'``, + - and the sum over all iterations ``'sum_niters'``. + + Note + ---- + In order to use these data, corresponding hook classes has to be defined before the simulation. Otherwise, no values can + be obtained. + + The global error does only make sense when an exact solution for the problem is available. Since ``'e_global'`` is stored + for each problem class, only for ``DiscontinuousTestODE`` the global error is taken into account when testing the solution. + + Also the event error ``'e_event'`` can only be computed if an exact event time is available. Since the function + ``controllerRun`` returns ``t_switch_exact=None`` when no exact event time is available, in order to compute the event error, + it has to be proven whether the list (in case of more than one event) contains ``None`` or not. Parameters ---------- stats : dict - Raw statistics from a controller run. - type : str - The type the be filtered. - sortby : str - String to specify which key to use for sorting. + Raw statistics of one simulation run. + prob_cls_name : str + Name of the problem class. + use_adaptivity : bool + Indicates whether adaptivity is used in the simulation or not. + use_detection : bool + Indicates whether event detection is used or not. + recomputed : bool + Indicates if values after successfully steps are used or not. + t_switch_exact : float + Exact event time of the problem. Returns ------- - sorted_list : list - List of filtered statistics. + res : dict + Dictionary with extracted data separated with reasonable keys. """ - sorted_nested_list = [] - times_unique = np.unique([me[0] for me in get_sorted(stats, type=type)]) - filtered_list = [ - filter_stats( - stats, - time=t_unique, - num_restarts=max([me.num_restarts for me in filter_stats(stats, type=type, time=t_unique).keys()]), - type=type, - ) - for t_unique in times_unique - ] - for item in filtered_list: - sorted_nested_list.append(sort_stats(item, sortby=sortby)) - sorted_list = [item for sub_item in sorted_nested_list for item in sub_item] - return sorted_list + res = {} + unknowns, unknowns_labels = getUnknownLabels(prob_cls_name) + # numerical solution + u_val = get_sorted(stats, type='u', sortby='time', recomputed=recomputed) + res['t'] = np.array([item[0] for item in u_val]) + for i, label in enumerate(unknowns): + res[label] = np.array([item[1][i] for item in u_val]) -def proof_assertions_description(description, use_adaptivity, use_switch_estimator): - """ - Function to proof the assertions in the description. + res['unknowns'] = unknowns + res['unknowns_labels'] = unknowns_labels - Parameters - ---------- - description : dict - Contains all information for a controller run. - use_adaptivity : bool - Flag if adaptivity wants to be used or not. - use_switch_estimator : bool - Flag if the switch estimator wants to be used or not. - """ + # global error + res['e_global'] = np.array(get_sorted(stats, type='e_global_post_step', sortby='time', recomputed=recomputed)) - n = description['problem_params']['ncapacitors'] - assert ( - description['problem_params']['alpha'] > description['problem_params']['V_ref'][k] for k in range(n) - ), 'Please set "alpha" greater than values of "V_ref"' - assert type(description['problem_params']['V_ref']) == np.ndarray, '"V_ref" needs to be an np.ndarray' - assert type(description['problem_params']['C']) == np.ndarray, '"C" needs to be an np.ndarray ' - assert ( - np.shape(description['problem_params']['V_ref'])[0] == n - ), 'Number of reference values needs to be equal to number of condensators' - assert ( - np.shape(description['problem_params']['C'])[0] == n - ), 'Number of capacitance values needs to be equal to number of condensators' - - assert ( - description['problem_params']['V_ref'][k] > 0 for k in range(n) - ), 'Please set values of "V_ref" greater than 0' - - assert 'errtol' not in description['step_params'].keys(), 'No exact solution known to compute error' - assert 'alpha' in description['problem_params'].keys(), 'Please supply "alpha" in the problem parameters' - assert 'V_ref' in description['problem_params'].keys(), 'Please supply "V_ref" in the problem parameters' + # event time(s) found by event detection + if use_detection: + switches = get_sorted(stats, type='switch', sortby='time', recomputed=recomputed) + assert len(switches) >= 1, 'No events found!' + t_switches = [t[1] for t in switches] + res['t_switches'] = t_switches - if use_adaptivity: - assert description['level_params']['restol'] == -1, "Please set restol to -1 or omit it" + t_switch_exact = [t_switch_exact] + res['t_switch_exact'] = t_switch_exact + if not all(t is None for t in t_switch_exact): + event_err = [ + abs(num_item - ex_item) for (num_item, ex_item) in zip(res['t_switches'], res['t_switch_exact']) + ] + res['e_event'] = event_err -def proof_assertions_time(dt, Tend, V_ref, alpha): - """ - Function to proof the assertions regarding the time domain (in combination with the specific problem). + h_val = get_sorted(stats, type='state_function', sortby='time', recomputed=recomputed) + h = np.array([np.abs(val[1]) for val in h_val]) + res['state_function'] = h - Parameters - ---------- - dt : float - Time step for computation. - Tend : float - End time. - V_ref : np.ndarray - Reference values (problem parameter). - alpha : float - Multiple used for initial conditions (problem_parameter). - """ + # embedded error and adapted step sizes + if use_adaptivity: + res['e_em'] = np.array(get_sorted(stats, type='error_embedded_estimate', sortby='time', recomputed=recomputed)) + res['dt'] = np.array(get_sorted(stats, type='dt', recomputed=recomputed)) - assert dt < Tend, "Time step is too large for the time domain!" + # sum over restarts + if use_adaptivity or use_detection: + res['sum_restarts'] = np.sum(np.array(get_sorted(stats, type='restart', recomputed=None, sortby='time'))[:, 1]) - assert ( - Tend == 0.3 and V_ref[0] == 1.0 and alpha == 1.2 - ), "Error! Do not use other parameters for V_ref != 1.0, alpha != 1.2, Tend != 0.3 due to hardcoded reference!" - assert dt == 1e-2, "Error! Do not use another time step dt!= 1e-2!" + # sum over all iterations + res['sum_niters'] = np.sum(np.array(get_sorted(stats, type='niter', recomputed=None, sortby='time'))[:, 1]) + return res if __name__ == "__main__": - run() + main() diff --git a/pySDC/projects/PinTSimE/buck_model.py b/pySDC/projects/PinTSimE/buck_model.py index 4e5da79565..3a2423d339 100644 --- a/pySDC/projects/PinTSimE/buck_model.py +++ b/pySDC/projects/PinTSimE/buck_model.py @@ -1,131 +1,64 @@ -import numpy as np -import dill -from pathlib import Path - -from pySDC.helpers.stats_helper import get_sorted -from pySDC.core.Collocation import CollBase as Collocation -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.implementations.problem_classes.BuckConverter import buck_converter from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -from pySDC.projects.PinTSimE.piline_model import log_data, setup_mpl -import pySDC.helpers.plot_helper as plt_helper - - -def main(): - """ - A simple test program to do SDC/PFASST runs for the buck converter model - """ - - # initialize level parameters - level_params = dict() - level_params['restol'] = 1e-12 - level_params['dt'] = 1e-5 - - # initialize sweeper parameters - sweeper_params = dict() - sweeper_params['collocation_class'] = Collocation - sweeper_params['quad_type'] = 'LOBATTO' - sweeper_params['num_nodes'] = 5 - sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part - - # initialize problem parameters - problem_params = dict() - problem_params['duty'] = 0.5 # duty cycle - problem_params['fsw'] = 1e3 # switching freqency - - # initialize step parameters - step_params = dict() - step_params['maxiter'] = 20 - - # initialize controller parameters - controller_params = dict() - controller_params['logger_level'] = 30 - controller_params['hook_class'] = log_data - - # fill description dictionary for easy step instantiation - description = dict() - description['problem_class'] = buck_converter # pass problem class - description['problem_params'] = problem_params # pass problem parameters - description['sweeper_class'] = imex_1st_order # pass sweeper - description['sweeper_params'] = sweeper_params # pass sweeper parameters - description['level_params'] = level_params # pass level parameters - description['step_params'] = step_params - - assert 'errtol' not in description['step_params'].keys(), 'No exact solution known to compute error' - assert 'duty' in description['problem_params'].keys(), 'Please supply "duty" in the problem parameters' - assert 'fsw' in description['problem_params'].keys(), 'Please supply "fsw" in the problem parameters' - - assert 0 <= problem_params['duty'] <= 1, 'Please set "duty" greater than or equal to 0 and less than or equal to 1' - # set time parameters - t0 = 0.0 - Tend = 2e-2 +from pySDC.projects.PinTSimE.battery_model import runSimulation - # instantiate controller - controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) +from pySDC.implementations.hooks.log_solution import LogSolution - # get initial values on finest level - P = controller.MS[0].levels[0].prob - uinit = P.u_exact(t0) - # call main function to get things done... - uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - - Path("data").mkdir(parents=True, exist_ok=True) - fname = 'data/buck.dat' - f = open(fname, 'wb') - dill.dump(stats, f) - f.close() - - # filter statistics by number of iterations - iter_counts = get_sorted(stats, type='niter', sortby='time') - - # compute and print statistics - min_iter = 20 - max_iter = 0 - - f = open('data/buck_out.txt', 'w') - niters = np.array([item[1] for item in iter_counts]) - out = ' Mean number of iterations: %4.2f' % np.mean(niters) - f.write(out + '\n') - print(out) - for item in iter_counts: - out = 'Number of iterations for time %4.2f: %1i' % item - f.write(out + '\n') - print(out) - min_iter = min(min_iter, item[1]) - max_iter = max(max_iter, item[1]) - - assert np.mean(niters) <= 8, "Mean number of iterations is too high, got %s" % np.mean(niters) - f.close() - - plot_voltages() - - -def plot_voltages(cwd='./'): - f = open(cwd + 'data/buck.dat', 'rb') - stats = dill.load(f) - f.close() - - # convert filtered statistics to list of iterations count, sorted by process - v1 = get_sorted(stats, type='v1', sortby='time') - v2 = get_sorted(stats, type='v2', sortby='time') - p3 = get_sorted(stats, type='p3', sortby='time') - - times = [v[0] for v in v1] - - setup_mpl() - fig, ax = plt_helper.plt.subplots(1, 1, figsize=(4.5, 3)) - ax.plot(times, [v[1] for v in v1], linewidth=1, label=r'$v_{C_1}$') - ax.plot(times, [v[1] for v in v2], linewidth=1, label=r'$v_{C_2}$') - ax.plot(times, [v[1] for v in p3], linewidth=1, label=r'$i_{L_\pi}$') - ax.legend(frameon=False, fontsize=12, loc='upper right') - - ax.set_xlabel('Time') - ax.set_ylabel('Energy') +def main(): + r""" + Executes the simulation. + + Note + ---- + Hardcoded solutions for battery models in `pySDC.projects.PinTSimE.hardcoded_solutions` are only computed for + ``dt_list=[1e-2, 1e-3]`` and ``M_fix=4``. Hence changing ``dt_list`` and ``M_fix`` to different values could arise + an ``AssertionError``. + """ - fig.savefig('data/buck_model_solution.png', dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig) + # --- defines parameters for sweeper ---- + M_fix = 5 + sweeper_params = { + 'num_nodes': M_fix, + 'quad_type': 'LOBATTO', + 'QI': 'IE', + } + + # --- defines parameters for event detection, max. number of iterations and restol ---- + handling_params = { + 'restol': 1e-12, + 'maxiter': 8, + 'max_restarts': None, + 'recomputed': None, + 'tol_event': None, + 'alpha': None, + 'exact_event_time_avail': None, + } + + # ---- all parameters are stored in this dictionary - note: defaults are used for the problem ---- + all_params = { + 'problem_params': {}, + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, + } + + hook_class = [LogSolution] + + use_detection = [False] + use_adaptivity = [False] + + _ = runSimulation( + problem=buck_converter, + sweeper=imex_1st_order, + all_params=all_params, + use_adaptivity=use_adaptivity, + use_detection=use_detection, + hook_class=hook_class, + interval=(0.0, 2e-2), + dt_list=[1e-5, 2e-5], + nnodes=[M_fix], + ) if __name__ == "__main__": diff --git a/pySDC/projects/PinTSimE/discontinuous_test_ODE.py b/pySDC/projects/PinTSimE/discontinuous_test_ODE.py index a8b343bf92..46bee5c4d9 100644 --- a/pySDC/projects/PinTSimE/discontinuous_test_ODE.py +++ b/pySDC/projects/PinTSimE/discontinuous_test_ODE.py @@ -1,27 +1,25 @@ import numpy as np -import dill from pathlib import Path from pySDC.helpers.stats_helper import get_sorted from pySDC.implementations.problem_classes.DiscontinuousTestODE import DiscontinuousTestODE from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI -from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator -from pySDC.projects.PinTSimE.battery_model import get_recomputed, generate_description +from pySDC.projects.PinTSimE.battery_model import runSimulation + import pySDC.helpers.plot_helper as plt_helper + from pySDC.core.Hooks import hooks from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep from pySDC.implementations.hooks.log_solution import LogSolution -class LogEvent(hooks): +class LogEventDiscontinuousTestODE(hooks): """ Logs the problem dependent state function of the discontinuous test ODE. """ def post_step(self, step, level_number): - super(LogEvent, self).post_step(step, level_number) + super().post_step(step, level_number) L = step.levels[level_number] @@ -39,239 +37,64 @@ def post_step(self, step, level_number): def main(): + r""" + Executes the simulation. + + Note + ---- + Hardcoded solutions for battery models in `pySDC.projects.PinTSimE.hardcoded_solutions` are only computed for + ``dt_list=[1e-2, 1e-3]`` and ``M_fix=4``. Hence changing ``dt_list`` and ``M_fix`` to different values could arise + an ``AssertionError``. """ - Executes the main stuff of the file. - """ - - Path("data").mkdir(parents=True, exist_ok=True) - - hookclass = [LogEvent, LogSolution, LogGlobalErrorPostStep] - - problem_class = DiscontinuousTestODE - - sweeper = generic_implicit - nnodes = [2, 3, 4] - maxiter = 8 - newton_tol = 1e-11 - - problem_params = dict() - problem_params['newton_maxiter'] = 50 - problem_params['newton_tol'] = newton_tol - - use_detection = [True, False] - - t0 = 1.4 - Tend = 1.7 - dt_list = [1e-2, 1e-3] - - for dt in dt_list: - for num_nodes in nnodes: - for use_SE in use_detection: - print(f'Controller run -- Simulation for step size: {dt}') - - restol = 1e-14 - recomputed = False if use_SE else None - - description, controller_params = generate_description( - dt, - problem_class, - sweeper, - num_nodes, - hookclass, - False, - use_SE, - problem_params, - restol, - maxiter, - max_restarts=None, - tol_event=1e-10, - ) - - proof_assertions(description, t0, Tend, recomputed, use_SE) - - stats, t_switch_exact = controller_run(t0, Tend, controller_params, description) - - if use_SE: - switches = get_recomputed(stats, type='switch', sortby='time') - assert len(switches) >= 1, 'No events found!' - test_event_error(stats, dt, t_switch_exact, num_nodes) - - test_error(stats, dt, num_nodes, use_SE, recomputed) - - -def controller_run(t0, Tend, controller_params, description): - """ - Executes a controller run for time interval to be specified in the arguments. - - Parameters - ---------- - t0 : float - Initial time of simulation. - Tend : float - End time of simulation. - controller_params : dict - Parameters needed for the controller. - description : dict - Contains all information for a controller run. - - Returns - ------- - stats : dict - Raw statistics from a controller run. - """ - - # instantiate the 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 - -def test_event_error(stats, dt, t_switch_exact, num_nodes): - """ - Tests the error between the exact event time and the event time founded by switch estimation. - - The errors to the exact event time are very small. The higher the number of collocation nodes - is the smaller the error to the exact event time is. - - Parameter - --------- - stats : dict - Raw statistics from a controller run. - dt : float - Current step size. - t_switch_exact : float - Exact event time of the problem. - num_nodes : int - Number of collocation nodes used. - """ - - switches = get_recomputed(stats, type='switch', sortby='time') - assert len(switches) >= 1, 'No switches found!' - t_switches = [v[1] for v in switches] - - # dict with hardcoded solution for event time error - t_event_err = { - 1e-2: { - 2: 1.7020665390443668e-06, - 3: 5.532252433937401e-10, - 4: 6.2776006615195e-11, - }, - 1e-3: { - 2: 1.7060500789867206e-08, - 3: 5.890081755666188e-10, - 4: 0.00057634035536136, - }, + # --- defines parameters for the problem class ---- + problem_params = { + 'newton_maxiter': 50, + 'newton_tol': 1e-11, } - t_event_err_got = abs(t_switch_exact - t_switches[-1]) - t_event_err_expected = t_event_err[dt][num_nodes] - - msg = f'Expected event time error {t_event_err_expected:.5f}, got {t_event_err_got:.5f}' - assert np.isclose(t_event_err_got, t_event_err_expected, atol=5e-3), msg - - -def test_error(stats, dt, num_nodes, use_SE, recomputed): - """ - Tests the error between the exact event solution and the numerical solution founded. - - In the dictionary containing the errors it can be clearly seen that errors are inherently reduced - using the switch estimator to predict the event and adapt the time step to resolve the event in - an more accurate way! - - Parameter - --------- - stats : dict - Raw statistics from a controller run. - dt : float - Current step size. - num_nodes : int - Number of collocation nodes used. - use_SE : bool - Indicates whether switch detection is used or not. - recomputed : bool - Indicates whether the values after a restart will be used. - """ - - err = get_sorted(stats, type='e_global_post_step', sortby='time', recomputed=recomputed) - err_norm = max([me[1] for me in err]) - - u_err = { - True: { - 1e-2: { - 2: 8.513409107457903e-06, - 3: 4.046930790480019e-09, - 4: 3.8459546658486943e-10, - }, - 1e-3: { - 2: 8.9185828500149e-08, - 3: 4.459276503609999e-09, - 4: 0.00015611434317808204, - }, - }, - False: { - 1e-2: { - 2: 0.014137551021780936, - 3: 0.009855041165877765, - 4: 0.006289698596543047, - }, - 1e-3: { - 2: 0.002674332734426521, - 3: 0.00057634035536136, - 4: 0.00015611434317808204, - }, - }, + # --- defines parameters for sweeper ---- + M_fix = 3 + sweeper_params = { + 'num_nodes': M_fix, + 'quad_type': 'LOBATTO', + 'QI': 'IE', } - u_err_expected = u_err[use_SE][dt][num_nodes] - u_err_got = err_norm - - msg = f'Expected event time error {u_err_expected:.7f}, got {u_err_got:.7f}' - assert np.isclose(u_err_got - u_err_expected, 0, atol=1e-11), msg - - -def proof_assertions(description, t0, Tend, recomputed, use_detection): - """ - Tests the parameters if they would not change. - - Parameters - ---------- - description : dict - Contains all information for a controller run. - t0 : float - Starting time. - Tend : float - End time. - recomputed : bool - Indicates whether the values after a restart are considered. - use_detection : bool - Indicates whether switch estimation is used. - """ - - newton_tol = description['problem_params']['newton_tol'] - msg = 'Newton tolerance should be set as small as possible to get best possible resolution of solution' - assert newton_tol <= 1e-8, msg - - assert t0 >= 1.0, 'Problem is only defined for t >= 1!' - assert Tend >= np.log(5), f'To investigate event, please set Tend larger than {np.log(5):.5f}' + # --- defines parameters for event detection ---- + handling_params = { + 'restol': 1e-13, + 'maxiter': 8, + 'max_restarts': 50, + 'recomputed': False, + 'tol_event': 1e-12, + 'alpha': 1.0, + 'exact_event_time_avail': True, + } - num_nodes = description['sweeper_params']['num_nodes'] - for M in [2, 3, 4]: - if num_nodes not in [2, 3, 4]: - assert num_nodes == M, f'Hardcoded solutions are only for M={M}!' + # ---- all parameters are stored in this dictionary ---- + all_params = { + 'problem_params': problem_params, + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, + } - sweeper = description['sweeper_class'].__name__ - assert sweeper == 'generic_implicit', 'Only generic_implicit sweeper is tested!' + hook_class = [LogEventDiscontinuousTestODE, LogSolution, LogGlobalErrorPostStep] - if use_detection: - assert recomputed == False, 'Be aware that recomputed is set to False by using switch detection!' + use_detection = [True, False] + use_adaptivity = [False] + + _ = runSimulation( + problem=DiscontinuousTestODE, + sweeper=generic_implicit, + all_params=all_params, + use_adaptivity=use_adaptivity, + use_detection=use_detection, + hook_class=hook_class, + interval=(1.0, 2.0), + dt_list=[1e-2, 1e-3], + nnodes=[M_fix], + ) if __name__ == "__main__": diff --git a/pySDC/projects/PinTSimE/estimation_check.py b/pySDC/projects/PinTSimE/estimation_check.py index 6fef419cb0..61fb858584 100644 --- a/pySDC/projects/PinTSimE/estimation_check.py +++ b/pySDC/projects/PinTSimE/estimation_check.py @@ -1,685 +1,394 @@ import numpy as np -import dill from pathlib import Path -from pySDC.helpers.stats_helper import get_sorted -from pySDC.implementations.problem_classes.Battery import battery, battery_implicit +from pySDC.helpers.stats_helper import sort_stats, filter_stats, get_sorted +from pySDC.implementations.problem_classes.Battery import battery_n_capacitors, battery, battery_implicit from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.projects.PinTSimE.battery_model import ( - controller_run, - check_solution, - generate_description, - get_recomputed, - LogEvent, - proof_assertions_description, -) - -from pySDC.projects.PinTSimE.piline_model import setup_mpl + +from pySDC.projects.PinTSimE.battery_model import runSimulation, plotStylingStuff + import pySDC.helpers.plot_helper as plt_helper -from pySDC.core.Hooks import hooks +from pySDC.projects.PinTSimE.battery_model import LogEventBattery from pySDC.implementations.hooks.log_solution import LogSolution from pySDC.implementations.hooks.log_step_size import LogStepSize from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate -from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator -from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity - -def run(cwd='./'): - """ - Routine to check the differences between using a switch estimator or not +def run_estimation_check(): + r""" + Generates plots to visualise results applying the Switch Estimator and Adaptivity to the battery models + containing. - Parameters - ---------- - cwd : str - Current working directory. + Note + ---- + Hardcoded solutions for battery models in `pySDC.projects.PinTSimE.hardcoded_solutions` are only computed for + ``dt_list=[1e-2, 1e-3]`` and ``M_fix=4``. Hence changing ``dt_list`` and ``M_fix`` to different values could arise + an ``AssertionError``. """ - dt_list = [1e-2, 1e-3] - t0 = 0.0 - Tend = 0.3 - - problem_classes = [battery, battery_implicit] - sweeper_classes = [imex_1st_order, generic_implicit] - num_nodes = 4 - restol = -1 - maxiter = 8 - - ncapacitors = 1 - alpha = 1.2 - V_ref = np.array([1.0]) - C = np.array([1.0]) - - problem_params = dict() - problem_params['ncapacitors'] = ncapacitors - problem_params['C'] = C - problem_params['alpha'] = alpha - problem_params['V_ref'] = V_ref - - hook_class = [LogSolution, LogEvent, LogStepSize, LogEmbeddedErrorEstimate] - - max_restarts = 1 - use_switch_estimator = [True, False] + Path("data").mkdir(parents=True, exist_ok=True) + + # --- defines parameters for sweeper ---- + M_fix = 4 + sweeper_params = { + 'num_nodes': M_fix, + 'quad_type': 'LOBATTO', + 'QI': 'IE', + } + + # --- defines parameters for event detection and maximum number of iterations ---- + handling_params = { + 'restol': -1, + 'maxiter': 8, + 'max_restarts': 50, + 'recomputed': False, + 'tol_event': 1e-10, + 'alpha': 1.0, + 'exact_event_time_avail': None, + } + + problem_classes = [battery, battery_implicit, battery_n_capacitors] + prob_class_names = [cls.__name__ for cls in problem_classes] + sweeper_classes = [imex_1st_order, generic_implicit, imex_1st_order] + + # --- defines parameters for battery models ---- + params_battery_1capacitor = { + 'ncapacitors': 1, + 'C': np.array([1.0]), + 'alpha': 1.2, + 'V_ref': np.array([1.0]), + } + + params_battery_2capacitors = { + 'ncapacitors': 2, + 'C': np.array([1.0, 1.0]), + 'alpha': 1.2, + 'V_ref': np.array([1.0, 1.0]), + } + + # --- parameters for each problem class are stored in this dictionary ---- + all_params = { + 'battery': { + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, + 'problem_params': params_battery_1capacitor, + }, + 'battery_implicit': { + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, + 'problem_params': params_battery_1capacitor, + }, + 'battery_n_capacitors': { + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, + 'problem_params': params_battery_2capacitors, + }, + } + + # ---- simulation domain for each problem class ---- + interval = { + 'battery': (0.0, 0.3), + 'battery_implicit': (0.0, 0.3), + 'battery_n_capacitors': (0.0, 0.5), + } + + hook_class = [LogSolution, LogEventBattery, LogEmbeddedErrorEstimate, LogStepSize] + + use_detection = [True, False] use_adaptivity = [True, False] - restarts_SE = [] - restarts_adapt = [] - restarts_SE_adapt = [] - - for problem, sweeper in zip(problem_classes, sweeper_classes): - for dt_item in dt_list: - for use_SE in use_switch_estimator: - for use_A in use_adaptivity: - tol_event = 1e-10 if sweeper.__name__ == 'generic_implicit' else 1e-17 - description, controller_params = generate_description( - dt_item, - problem, - sweeper, - num_nodes, - hook_class, - use_A, - use_SE, - problem_params, - restol, - maxiter, - max_restarts, - tol_event, - ) - # Assertions - proof_assertions_description(description, use_A, use_SE) - - stats = controller_run(description, controller_params, use_A, use_SE, t0, Tend) - - if use_A or use_SE: - check_solution(stats, dt_item, problem.__name__, use_A, use_SE) - - fname = 'data/battery_dt{}_USE{}_USA{}_{}.dat'.format(dt_item, use_SE, use_A, sweeper.__name__) - f = open(fname, 'wb') - dill.dump(stats, f) - f.close() - - if use_SE or use_A: - restarts_sorted = np.array(get_sorted(stats, type='restart', recomputed=None))[:, 1] - if use_SE and not use_A: - restarts_SE.append(np.sum(restarts_sorted)) - - elif not use_SE and use_A: - restarts_adapt.append(np.sum(restarts_sorted)) - - elif use_SE and use_A: - restarts_SE_adapt.append(np.sum(restarts_sorted)) - - accuracy_check(dt_list, problem.__name__, sweeper.__name__, V_ref) - - differences_around_switch( - dt_list, - problem.__name__, - restarts_SE, - restarts_adapt, - restarts_SE_adapt, - sweeper.__name__, - V_ref, + for problem, sweeper, prob_cls_name in zip(problem_classes, sweeper_classes, prob_class_names): + u_num = runSimulation( + problem=problem, + sweeper=sweeper, + all_params=all_params[prob_cls_name], + use_adaptivity=use_adaptivity, + use_detection=use_detection, + hook_class=hook_class, + interval=interval[prob_cls_name], + dt_list=[1e-2, 1e-3], + nnodes=[M_fix], ) - differences_over_time(dt_list, problem.__name__, sweeper.__name__, V_ref) + plotAccuracyCheck(u_num, prob_cls_name, M_fix) - iterations_over_time(dt_list, description['step_params']['maxiter'], problem.__name__, sweeper.__name__) + plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix) - restarts_SE = [] - restarts_adapt = [] - restarts_SE_adapt = [] + plotStateFunctionOverTime(u_num, prob_cls_name, M_fix) -def accuracy_check(dt_list, problem, sweeper, V_ref, cwd='./'): - """ +def plotAccuracyCheck(u_num, prob_cls_name, M_fix): # pragma: no cover + r""" Routine to check accuracy for different step sizes in case of using adaptivity. Parameters ---------- - dt_list : list - List of considered (initial) step sizes. - problem : pySDC.core.Problem.ptype - Problem class used to consider (the class name). - sweeper : pySDC.core.Sweeper.sweeper - Sweeper used to solve (the class name). - V_ref : float - Reference value for the switch. - cwd : str - Current working directory. + u_num : dict + Contains the all the data. Dictionary has the structure ``u_num[dt][M][use_SE][use_A]``, + where for each step size ``dt``, for each number of collocation node ``M``, for each + combination of event detection ``use_SE`` and adaptivity ``use_A`` appropriate stuff is stored. + For more details, see ``pySDC.projects.PinTSimE.battery_model.getDataDict``. + prob_cls_name : str + Name of the problem class. + M_fix : int + Fixed number of collocation nodes the plot is generated for. """ - if len(dt_list) > 1: - setup_mpl() - fig_acc, ax_acc = plt_helper.plt.subplots( - 1, len(dt_list), figsize=(3 * len(dt_list), 3), sharex='col', sharey='row' + colors = plotStylingStuff() + dt_list = u_num.keys() + + use_A = True + for dt in dt_list: + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5), sharex='col', sharey='row') + e_ax = ax.twinx() + for use_SE in u_num[dt][M_fix].keys(): + dt_val = u_num[dt][M_fix][use_SE][use_A]['dt'] + e_em_val = u_num[dt][M_fix][use_SE][use_A]['e_em'] + if use_SE: + t_switches = u_num[dt][M_fix][use_SE][use_A]['t_switches'] + + for i in range(len(t_switches)): + ax.axvline(x=t_switches[i], linestyle='--', color='tomato', label='Event {}'.format(i + 1)) + + ax.plot(dt_val[:, 0], dt_val[:, 1], color=colors[use_SE][use_A], label=r'SE={}, A={}'.format(use_SE, use_A)) + + e_ax.plot(e_em_val[:, 0], e_em_val[:, 1], linestyle='dashdot', color=colors[use_SE][use_A]) + + ax.plot(0, 0, color='black', linestyle='solid', label=r'$\Delta t_\mathrm{adapt}$') + ax.plot(0, 0, color='black', linestyle='dashdot', label=r'$e_{em}$') + + e_ax.set_yscale('log', base=10) + e_ax.set_ylabel(r'Embedded error estimate $e_{em}$', fontsize=16) + e_ax.set_ylim(1e-16, 1e-7) + e_ax.tick_params(labelsize=16) + e_ax.minorticks_off() + + ax.tick_params(axis='both', which='major', labelsize=16) + ax.set_ylim(1e-9, 1e0) + ax.set_yscale('log', base=10) + ax.set_xlabel(r'Time $t$', fontsize=16) + ax.set_ylabel(r'Adapted step sizes $\Delta t_\mathrm{adapt}$', fontsize=16) + ax.grid(visible=True) + ax.minorticks_off() + ax.legend(frameon=True, fontsize=12, loc='center left') + + fig.savefig( + 'data/detection_and_adaptivity_{}_dt={}_M={}.png'.format(prob_cls_name, dt, M_fix), + dpi=300, + bbox_inches='tight', ) + plt_helper.plt.close(fig) - else: - setup_mpl() - fig_acc, ax_acc = plt_helper.plt.subplots(1, 1, figsize=(3, 3), sharex='col', sharey='row') - - count_ax = 0 - for dt_item in dt_list: - f3 = open(cwd + 'data/battery_dt{}_USETrue_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE_adapt = dill.load(f3) - f3.close() - - f4 = open(cwd + 'data/battery_dt{}_USEFalse_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_adapt = dill.load(f4) - f4.close() - - switches_SE_adapt = get_recomputed(stats_SE_adapt, type='switch', sortby='time') - t_switch_SE_adapt = [v[1] for v in switches_SE_adapt] - t_switch_SE_adapt = t_switch_SE_adapt[-1] - - dt_SE_adapt_val = get_sorted(stats_SE_adapt, type='dt', recomputed=False) - dt_adapt_val = get_sorted(stats_adapt, type='dt', recomputed=False) - - e_emb_SE_adapt_val = get_sorted(stats_SE_adapt, type='e_embedded', recomputed=False) - e_emb_adapt_val = get_sorted(stats_adapt, type='e_embedded', recomputed=False) - - times_SE_adapt = [v[0] for v in e_emb_SE_adapt_val] - times_adapt = [v[0] for v in e_emb_adapt_val] - - e_emb_SE_adapt = [v[1] for v in e_emb_SE_adapt_val] - e_emb_adapt = [v[1] for v in e_emb_adapt_val] - - if len(dt_list) > 1: - ax_acc[count_ax].set_title(r'$\Delta t_\mathrm{initial}$=%s' % dt_item) - dt1 = ax_acc[count_ax].plot( - [v[0] for v in dt_SE_adapt_val], - [v[1] for v in dt_SE_adapt_val], - 'ko-', - label=r'SE+A - $\Delta t_\mathrm{adapt}$', - ) - dt2 = ax_acc[count_ax].plot( - [v[0] for v in dt_adapt_val], [v[1] for v in dt_adapt_val], 'g-', label=r'A - $\Delta t_\mathrm{adapt}$' - ) - ax_acc[count_ax].axvline(x=t_switch_SE_adapt, linestyle='--', linewidth=0.5, color='r', label='Switch') - ax_acc[count_ax].tick_params(axis='both', which='major', labelsize=6) - ax_acc[count_ax].set_xlabel('Time', fontsize=6) - if count_ax == 0: - ax_acc[count_ax].set_ylabel(r'$\Delta t_\mathrm{adapt}$', fontsize=6) - - e_ax = ax_acc[count_ax].twinx() - e_plt1 = e_ax.plot(times_SE_adapt, e_emb_SE_adapt, 'k--', label=r'SE+A - $\epsilon_{emb}$') - e_plt2 = e_ax.plot(times_adapt, e_emb_adapt, 'g--', label=r'A - $\epsilon_{emb}$') - e_ax.set_yscale('log', base=10) - e_ax.set_ylim(1e-16, 1e-7) - e_ax.tick_params(labelsize=6) - - lines = dt1 + e_plt1 + dt2 + e_plt2 - labels = [l.get_label() for l in lines] - - ax_acc[count_ax].legend(lines, labels, frameon=False, fontsize=6, loc='upper right') - - else: - ax_acc.set_title(r'$\Delta t_\mathrm{initial}$=%s' % dt_item) - dt1 = ax_acc.plot( - [v[0] for v in dt_SE_adapt_val], - [v[1] for v in dt_SE_adapt_val], - 'ko-', - label=r'SE+A - $\Delta t_\mathrm{adapt}$', - ) - dt2 = ax_acc.plot( - [v[0] for v in dt_adapt_val], - [v[1] for v in dt_adapt_val], - 'go-', - label=r'A - $\Delta t_\mathrm{adapt}$', - ) - ax_acc.axvline(x=t_switch_SE_adapt, linestyle='--', linewidth=0.5, color='r', label='Switch') - ax_acc.tick_params(axis='both', which='major', labelsize=6) - ax_acc.set_xlabel('Time', fontsize=6) - ax_acc.set_ylabel(r'$Delta t_\mathrm{adapt}$', fontsize=6) - - e_ax = ax_acc.twinx() - e_plt1 = e_ax.plot(times_SE_adapt, e_emb_SE_adapt, 'k--', label=r'SE+A - $\epsilon_{emb}$') - e_plt2 = e_ax.plot(times_adapt, e_emb_adapt, 'g--', label=r'A - $\epsilon_{emb}$') - e_ax.set_yscale('log', base=10) - e_ax.tick_params(labelsize=6) - - lines = dt1 + e_plt1 + dt2 + e_plt2 - labels = [l.get_label() for l in lines] - - ax_acc.legend(lines, labels, frameon=False, fontsize=6, loc='upper right') - - count_ax += 1 - - fig_acc.savefig('data/embedded_error_adaptivity_{}.png'.format(sweeper), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig_acc) - - -def differences_around_switch( - dt_list, problem, restarts_SE, restarts_adapt, restarts_SE_adapt, sweeper, V_ref, cwd='./' -): - """ - Routine to plot the differences before, at, and after the switch. Produces the - diffs_estimation_.png file + +def plotStateFunctionAroundEvent(u_num, prob_cls_name, M_fix): # pragma: no cover + r""" + 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. Parameters ---------- - dt_list : list - List of considered (initial) step sizes. - problem : pySDC.core.Problem.ptype - Problem class used to consider (the class name). - restarts_SE : list - Restarts for the solve only using the switch estimator. - restarts_adapt : list - Restarts for the solve of only using adaptivity. - restarts_SE_adapt : list - Restarts for the solve of using both, switch estimator and adaptivity. - sweeper : pySDC.core.Sweeper.sweeper - Sweeper used to solve (the class name). - V_ref float - Reference value for the switch. - cwd : str - Current working directory. + u_num : dict + Contains the all the data. Dictionary has the structure ``u_num[dt][M][use_SE][use_A]``, + where for each step size ``dt``, for each number of collocation node ``M``, for each + combination of event detection ``use_SE`` and adaptivity ``use_A`` appropriate stuff is stored. + For more details, see ``pySDC.projects.PinTSimE.battery_model.getDataDict``. + prob_cls_name : str + Name of the problem class. + M_fix : int + Fixed number of collocation nodes the plot is generated for. """ - diffs_true_at = [] - diffs_false_before = [] - diffs_false_after = [] - - diffs_true_at_adapt = [] - diffs_true_before_adapt = [] - diffs_true_after_adapt = [] - - diffs_false_before_adapt = [] - diffs_false_after_adapt = [] - for dt_item in dt_list: - f1 = open(cwd + 'data/battery_dt{}_USETrue_USAFalse_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE = dill.load(f1) - f1.close() - - f2 = open(cwd + 'data/battery_dt{}_USEFalse_USAFalse_{}.dat'.format(dt_item, sweeper), 'rb') - stats = dill.load(f2) - f2.close() - - f3 = open(cwd + 'data/battery_dt{}_USETrue_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE_adapt = dill.load(f3) - f3.close() - - f4 = open(cwd + 'data/battery_dt{}_USEFalse_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_adapt = dill.load(f4) - f4.close() - - switches_SE = get_recomputed(stats_SE, type='switch', sortby='time') - t_switch = [v[1] for v in switches_SE] - t_switch = t_switch[-1] # battery has only one single switch - - switches_SE_adapt = get_recomputed(stats_SE_adapt, type='switch', sortby='time') - t_switch_SE_adapt = [v[1] for v in switches_SE_adapt] - t_switch_SE_adapt = t_switch_SE_adapt[-1] - - vC_SE = [me[1][1] for me in get_sorted(stats_SE, type='u', recomputed=False)] - vC_adapt = [me[1][1] for me in get_sorted(stats_adapt, type='u', recomputed=False)] - vC_SE_adapt = [me[1][1] for me in get_sorted(stats_SE_adapt, type='u', recomputed=False)] - vC = [me[1][1] for me in get_sorted(stats, type='u', recomputed=False)] - - diff_SE, diff = vC_SE - V_ref[0], vC - V_ref[0] - times_SE = [me[0] for me in get_sorted(stats_SE, type='u', recomputed=False)] - times = [me[0] for me in get_sorted(stats, type='u', recomputed=False)] - - diff_adapt, diff_SE_adapt = vC_adapt - V_ref[0], vC_SE_adapt - V_ref[0] - times_adapt = [me[0] for me in get_sorted(stats_adapt, type='u', recomputed=False)] - times_SE_adapt = [me[0] for me in get_sorted(stats_SE_adapt, type='u', recomputed=False)] - - diffs_true_at.append([diff_SE[m] for m in range(len(times_SE)) if abs(times_SE[m] - t_switch) <= 1e-7][0]) - - diffs_false_before.append( - [diff[m - 1] for m in range(1, len(times)) if times[m - 1] <= t_switch <= times[m]][0] + title_cases = { + 0: 'Using detection', + 1: 'Using adaptivity', + 2: 'Using adaptivity and detection', + } + + dt_list = list(u_num.keys()) + use_detection = u_num[list(dt_list)[0]][M_fix].keys() + use_adaptivity = u_num[list(dt_list)[0]][M_fix][list(use_detection)[0]].keys() + h0 = u_num[list(dt_list)[0]][M_fix][list(use_detection)[0]][list(use_adaptivity)[0]]['state_function'] + n = h0[0].shape[0] + + for i in range(n): + fig, ax = plt_helper.plt.subplots(1, 3, figsize=(12, 4), sharex='col', sharey='row', squeeze=False) + dt_list = list(u_num.keys()) + for use_SE in use_detection: + for use_A in use_adaptivity: + # ---- decide whether state function (w/o handling) has two entries or one; choose correct one with reshaping ---- + h_val_no_handling = [u_num[dt][M_fix][False][False]['state_function'] for dt in dt_list] + h_no_handling = [item[:] if n == 1 else item[:, i] for item in h_val_no_handling] + h_no_handling = [item.reshape((item.shape[0],)) for item in h_no_handling] + + t_no_handling = [u_num[dt][M_fix][False][False]['t'] for dt in dt_list] + + if not use_A and not use_SE: + continue + else: + ind = 0 if (not use_A and use_SE) else (1 if (use_A and not use_SE) else 2) + ax[0, ind].set_title(r'{} for $n={}$'.format(title_cases[ind], i + 1)) + + # ---- same is done here for state function of other cases ---- + h_val = [u_num[dt][M_fix][use_SE][use_A]['state_function'] for dt in dt_list] + h = [item[:] if n == 1 else item[:, i] for item in h_val] + h = [item.reshape((item.shape[0],)) for item in h] + + t = [u_num[dt][M_fix][use_SE][use_A]['t'] for dt in dt_list] + + 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] + + 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 + ], + color='limegreen', + marker='s', + linestyle='solid', + alpha=0.4, + label='At event', + ) + + ax[0, ind].plot( + dt_list, + [ + h_item[m - 1] + for (t_item, h_item, t_switch_item) in zip(t_no_handling, h_no_handling, t_switch) + for m in range(1, len(t_item)) + if t_item[m - 1] < t_switch_item < t_item[m] + ], + color='firebrick', + marker='o', + linestyle='solid', + alpha=0.4, + label='Before event', + ) + + ax[0, ind].plot( + dt_list, + [ + h_item[m] + for (t_item, h_item, t_switch_item) in zip(t_no_handling, h_no_handling, t_switch) + for m in range(1, len(t_item)) + if t_item[m - 1] < t_switch_item < t_item[m] + ], + color='deepskyblue', + marker='*', + linestyle='solid', + alpha=0.4, + label='After event', + ) + + else: + ax[0, ind].plot( + dt_list, + [ + h_item[m - 1] + for (t_item, h_item, t_switch_item) in zip(t, h, t_switch) + for m in range(1, len(t_item)) + if t_item[m - 1] < t_switch_item < t_item[m] + ], + color='firebrick', + marker='o', + linestyle='solid', + alpha=0.4, + label='Before event', + ) + + 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(1, len(t_item)) + if t_item[m - 1] < t_switch_item < t_item[m] + ], + color='deepskyblue', + marker='*', + linestyle='solid', + alpha=0.4, + label='After event', + ) + + ax[0, ind].tick_params(axis='both', which='major', labelsize=16) + ax[0, ind].set_xticks(dt_list) + ax[0, ind].set_xticklabels(dt_list) + ax[0, ind].set_ylim(1e-15, 1e1) + ax[0, ind].set_yscale('log', base=10) + ax[0, ind].set_xlabel(r'Step size $\Delta t$', fontsize=16) + ax[0, 0].set_ylabel(r'Absolute values of h $|h(v_{C_n}(t))|$', fontsize=16) + ax[0, ind].grid(visible=True) + ax[0, ind].minorticks_off() + ax[0, ind].legend(frameon=True, fontsize=12, loc='lower left') + + fig.savefig( + 'data/{}_comparison_event{}_M={}.png'.format(prob_cls_name, i + 1, M_fix), dpi=300, bbox_inches='tight' ) - diffs_false_after.append([diff[m] for m in range(1, len(times)) if times[m - 1] <= t_switch <= times[m]][0]) + plt_helper.plt.close(fig) - for m in range(len(times_SE_adapt)): - if abs(times_SE_adapt[m] - t_switch_SE_adapt) <= 1e-10: - diffs_true_at_adapt.append(diff_SE_adapt[m]) - diffs_true_before_adapt.append(diff_SE_adapt[m - 1]) - diffs_true_after_adapt.append(diff_SE_adapt[m + 1]) - diffs_false_before_adapt.append( - [diff_adapt[m - 1] for m in range(len(times_adapt)) if times_adapt[m - 1] <= t_switch <= times_adapt[m]][0] - ) - diffs_false_after_adapt.append( - [diff_adapt[m] for m in range(len(times_adapt)) if times_adapt[m - 1] <= t_switch <= times_adapt[m]][0] - ) - - setup_mpl() - fig_around, ax_around = plt_helper.plt.subplots(1, 3, figsize=(9, 3), sharex='col', sharey='row') - ax_around[0].set_title("Using SE") - pos11 = ax_around[0].plot(dt_list, diffs_false_before, 'rs-', label='before switch') - pos12 = ax_around[0].plot(dt_list, diffs_false_after, 'bd--', label='after switch') - pos13 = ax_around[0].plot(dt_list, diffs_true_at, 'ko--', label='at switch') - ax_around[0].set_xticks(dt_list) - ax_around[0].set_xticklabels(dt_list) - ax_around[0].tick_params(axis='both', which='major', labelsize=6) - ax_around[0].set_xscale('log', base=10) - ax_around[0].set_yscale('symlog', linthresh=1e-8) - ax_around[0].set_ylim(-1, 1) - ax_around[0].set_xlabel(r'$\Delta t_\mathrm{initial}$', fontsize=6) - ax_around[0].set_ylabel(r'$v_{C}-V_{ref}$', fontsize=6) - - restart_ax0 = ax_around[0].twinx() - restarts_plt0 = restart_ax0.plot(dt_list, restarts_SE, 'cs--', label='Restarts') - restart_ax0.tick_params(labelsize=6) - - lines = pos11 + pos12 + pos13 + restarts_plt0 - labels = [l.get_label() for l in lines] - ax_around[0].legend(lines, labels, frameon=False, fontsize=6, loc='lower right') - - ax_around[1].set_title("Using Adaptivity") - pos21 = ax_around[1].plot(dt_list, diffs_false_before_adapt, 'rs-', label='before switch') - pos22 = ax_around[1].plot(dt_list, diffs_false_after_adapt, 'bd--', label='after switch') - ax_around[1].set_xticks(dt_list) - ax_around[1].set_xticklabels(dt_list) - ax_around[1].tick_params(axis='both', which='major', labelsize=6) - ax_around[1].set_xscale('log', base=10) - ax_around[1].set_yscale('symlog', linthresh=1e-8) - ax_around[1].set_ylim(-1, 1) - ax_around[1].set_xlabel(r'$\Delta t_\mathrm{initial}$', fontsize=6) - - restart_ax1 = ax_around[1].twinx() - restarts_plt1 = restart_ax1.plot(dt_list, restarts_adapt, 'cs--', label='Restarts') - restart_ax1.tick_params(labelsize=6) - - lines = pos21 + pos22 + restarts_plt1 - labels = [l.get_label() for l in lines] - ax_around[1].legend(lines, labels, frameon=False, fontsize=6, loc='lower right') - - ax_around[2].set_title("Using SE + Adaptivity") - pos31 = ax_around[2].plot(dt_list, diffs_true_before_adapt, 'rs-', label='before switch') - pos32 = ax_around[2].plot(dt_list, diffs_true_after_adapt, 'bd--', label='after switch') - pos33 = ax_around[2].plot(dt_list, diffs_true_at_adapt, 'ko--', label='at switch') - ax_around[2].set_xticks(dt_list) - ax_around[2].set_xticklabels(dt_list) - ax_around[2].tick_params(axis='both', which='major', labelsize=6) - ax_around[2].set_xscale('log', base=10) - ax_around[2].set_yscale('symlog', linthresh=1e-8) - ax_around[2].set_ylim(-1, 1) - ax_around[2].set_xlabel(r'$\Delta t_\mathrm{initial}$', fontsize=6) - - restart_ax2 = ax_around[2].twinx() - restarts_plt2 = restart_ax2.plot(dt_list, restarts_SE_adapt, 'cs--', label='Restarts') - restart_ax2.tick_params(labelsize=6) - - lines = pos31 + pos32 + pos33 + restarts_plt2 - labels = [l.get_label() for l in lines] - ax_around[2].legend(frameon=False, fontsize=6, loc='lower right') - - fig_around.savefig('data/diffs_around_switch_{}.png'.format(sweeper), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig_around) - - -def differences_over_time(dt_list, problem, sweeper, V_ref, cwd='./'): - """ - Routine to plot the differences in time using the switch estimator or not. Produces the - difference_estimation_.png file. +def plotStateFunctionOverTime(u_num, prob_cls_name, M_fix): # pragma: no cover + r""" + Routine that plots the state function over time. Parameters ---------- - dt_list : list - List of considered (initial) step sizes. - problem : pySDC.core.Problem.ptype - Problem class used to consider (the class name). - sweeper : pySDC.core.Sweeper.sweeper - Sweeper used to solve (the class name). - V_ref : float - Reference value for the switch. - cwd : str - Current working directory. + u_num : dict + Contains the all the data. Dictionary has the structure ``u_num[dt][M][use_SE][use_A]``, + where for each step size ``dt``, for each number of collocation node ``M``, for each + combination of event detection ``use_SE`` and adaptivity ``use_A`` appropriate stuff is stored. + For more details, see ``pySDC.projects.PinTSimE.battery_model.getDataDict``. + prob_cls_name : str + Indicates the name of the problem class to be considered. + M_fix : int + Fixed number of collocation nodes the plot is generated for. """ - if len(dt_list) > 1: - setup_mpl() - fig_diffs, ax_diffs = plt_helper.plt.subplots( - 2, len(dt_list), figsize=(4 * len(dt_list), 6), sharex='col', sharey='row' - ) - - else: - setup_mpl() - fig_diffs, ax_diffs = plt_helper.plt.subplots(2, 1, figsize=(4, 6)) - - count_ax = 0 - for dt_item in dt_list: - f1 = open(cwd + 'data/battery_dt{}_USETrue_USAFalse_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE = dill.load(f1) - f1.close() - - f2 = open(cwd + 'data/battery_dt{}_USEFalse_USAFalse_{}.dat'.format(dt_item, sweeper), 'rb') - stats = dill.load(f2) - f2.close() - - f3 = open(cwd + 'data/battery_dt{}_USETrue_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE_adapt = dill.load(f3) - f3.close() - - f4 = open(cwd + 'data/battery_dt{}_USEFalse_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_adapt = dill.load(f4) - f4.close() - - switches_SE = get_recomputed(stats_SE, type='switch', sortby='time') - t_switch_SE = [v[1] for v in switches_SE] - t_switch_SE = t_switch_SE[-1] # battery has only one single switch - - switches_SE_adapt = get_recomputed(stats_SE_adapt, type='switch', sortby='time') - t_switch_SE_adapt = [v[1] for v in switches_SE_adapt] - t_switch_SE_adapt = t_switch_SE_adapt[-1] - - dt_adapt = np.array(get_sorted(stats_adapt, type='dt', recomputed=False)) - dt_SE_adapt = np.array(get_sorted(stats_SE_adapt, type='dt', recomputed=False)) - - restart_adapt = np.array(get_sorted(stats_adapt, type='restart', recomputed=None)) - restart_SE_adapt = np.array(get_sorted(stats_SE_adapt, type='restart', recomputed=None)) - - vC_SE = [me[1][1] for me in get_sorted(stats_SE, type='u', recomputed=False)] - vC_adapt = [me[1][1] for me in get_sorted(stats_adapt, type='u', recomputed=False)] - vC_SE_adapt = [me[1][1] for me in get_sorted(stats_SE_adapt, type='u', recomputed=False)] - vC = [me[1][1] for me in get_sorted(stats, type='u', recomputed=False)] - - diff_SE, diff = vC_SE - V_ref[0], vC - V_ref[0] - times_SE = [me[0] for me in get_sorted(stats_SE, type='u', recomputed=False)] - times = [me[0] for me in get_sorted(stats, type='u', recomputed=False)] - - diff_adapt, diff_SE_adapt = vC_adapt - V_ref[0], vC_SE_adapt - V_ref[0] - times_adapt = [me[0] for me in get_sorted(stats_adapt, type='u', recomputed=False)] - times_SE_adapt = [me[0] for me in get_sorted(stats_SE_adapt, type='u', recomputed=False)] - - if len(dt_list) > 1: - ax_diffs[0, count_ax].set_title(r'$\Delta t$=%s' % dt_item) - ax_diffs[0, count_ax].plot(times_SE, diff_SE, label='SE=True, A=False', color='#ff7f0e') - ax_diffs[0, count_ax].plot(times, diff, label='SE=False, A=False', color='#1f77b4') - ax_diffs[0, count_ax].plot(times_adapt, diff_adapt, label='SE=False, A=True', color='red', linestyle='--') - ax_diffs[0, count_ax].plot( - times_SE_adapt, diff_SE_adapt, label='SE=True, A=True', color='limegreen', linestyle='-.' - ) - ax_diffs[0, count_ax].axvline(x=t_switch_SE, linestyle='--', linewidth=0.5, color='k', label='Switch') - ax_diffs[0, count_ax].legend(frameon=False, fontsize=6, loc='lower left') - ax_diffs[0, count_ax].set_yscale('symlog', linthresh=1e-5) - ax_diffs[0, count_ax].tick_params(axis='both', which='major', labelsize=6) - if count_ax == 0: - ax_diffs[0, count_ax].set_ylabel('Difference $v_{C}-V_{ref}$', fontsize=6) - - if count_ax == 0 or count_ax == 1: - ax_diffs[0, count_ax].legend(frameon=False, fontsize=6, loc='upper right') - - else: - ax_diffs[0, count_ax].legend(frameon=False, fontsize=6, loc='upper right') - - ax_diffs[1, count_ax].plot( - dt_adapt[:, 0], dt_adapt[:, 1], label=r'$\Delta t$ - SE=F, A=T', color='red', linestyle='--' - ) - ax_diffs[1, count_ax].plot([None], [None], label='Restart - SE=F, A=T', color='grey', linestyle='-.') - - for i in range(len(restart_adapt)): - if restart_adapt[i, 1] > 0: - ax_diffs[1, count_ax].axvline(restart_adapt[i, 0], color='grey', linestyle='-.') - - ax_diffs[1, count_ax].plot( - dt_SE_adapt[:, 0], - dt_SE_adapt[:, 1], - label=r'$ \Delta t$ - SE=T, A=T', - color='limegreen', - linestyle='-.', - ) - ax_diffs[1, count_ax].plot([None], [None], label='Restart - SE=T, A=T', color='black', linestyle='-.') - - for i in range(len(restart_SE_adapt)): - if restart_SE_adapt[i, 1] > 0: - ax_diffs[1, count_ax].axvline(restart_SE_adapt[i, 0], color='black', linestyle='-.') - - ax_diffs[1, count_ax].set_xlabel('Time', fontsize=6) - ax_diffs[1, count_ax].tick_params(axis='both', which='major', labelsize=6) - if count_ax == 0: - ax_diffs[1, count_ax].set_ylabel(r'$\Delta t_\mathrm{adapted}$', fontsize=6) - - ax_diffs[1, count_ax].set_yscale('log', base=10) - ax_diffs[1, count_ax].legend(frameon=True, fontsize=6, loc='lower left') - - else: - ax_diffs[0].set_title(r'$\Delta t$=%s' % dt_item) - ax_diffs[0].plot(times_SE, diff_SE, label='SE=True', color='#ff7f0e') - ax_diffs[0].plot(times, diff, label='SE=False', color='#1f77b4') - ax_diffs[0].plot(times_adapt, diff_adapt, label='SE=False, A=True', color='red', linestyle='--') - ax_diffs[0].plot(times_SE_adapt, diff_SE_adapt, label='SE=True, A=True', color='limegreen', linestyle='-.') - ax_diffs[0].axvline(x=t_switch_SE, linestyle='--', linewidth=0.5, color='k', label='Switch') - ax_diffs[0].tick_params(axis='both', which='major', labelsize=6) - ax_diffs[0].set_yscale('symlog', linthresh=1e-5) - ax_diffs[0].set_ylabel('Difference $v_{C}-V_{ref}$', fontsize=6) - ax_diffs[0].legend(frameon=False, fontsize=6, loc='center right') - - ax_diffs[1].plot(dt_adapt[:, 0], dt_adapt[:, 1], label='SE=False, A=True', color='red', linestyle='--') - ax_diffs[1].plot( - dt_SE_adapt[:, 0], dt_SE_adapt[:, 1], label='SE=True, A=True', color='limegreen', linestyle='-.' - ) - ax_diffs[1].tick_params(axis='both', which='major', labelsize=6) - ax_diffs[1].set_xlabel('Time', fontsize=6) - ax_diffs[1].set_ylabel(r'$\Delta t_\mathrm{adapted}$', fontsize=6) - ax_diffs[1].set_yscale('log', base=10) - - ax_diffs[1].legend(frameon=False, fontsize=6, loc='upper right') - - count_ax += 1 - - plt_helper.plt.tight_layout() - fig_diffs.savefig('data/diffs_over_time_{}.png'.format(sweeper), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig_diffs) - - -def iterations_over_time(dt_list, maxiter, problem, sweeper, cwd='./'): - """ - Routine to plot the number of iterations over time using switch estimator or not. Produces the - iters_.png file. - - Parameters - ---------- - dt_list : list - List of considered (initial) step sizes. - maxiter : int - Maximum number of iterations. - problem : pySDC.core.Problem.ptype - Problem class used to consider (the class name). - sweeper : pySDC.core.Sweeper.sweeper - Sweeper used to solve (the class name). - cwd : str - Current working directory. - """ + dt_list = u_num.keys() + use_detection = u_num[list(dt_list)[0]][M_fix].keys() + use_adaptivity = u_num[list(dt_list)[0]][M_fix][list(use_detection)[0]].keys() + h0 = u_num[list(dt_list)[0]][M_fix][list(use_detection)[0]][list(use_adaptivity)[0]]['state_function'] + n = h0[0].shape[0] + for dt in dt_list: + figsize = (7.5, 5) if n == 1 else (12, 5) + fig, ax = plt_helper.plt.subplots(1, n, figsize=figsize, sharex='col', sharey='row', squeeze=False) + + for use_SE in use_detection: + for use_A in use_adaptivity: + t = u_num[dt][M_fix][use_SE][use_A]['t'] + h_val = u_num[dt][M_fix][use_SE][use_A]['state_function'] + + linestyle = 'dashdot' if use_A else 'dotted' + for i in range(n): + h = h_val[:] if n == 1 else h_val[:, i] + ax[0, i].set_title(r'$n={}$'.format(i + 1)) + ax[0, i].plot( + t, h, linestyle=linestyle, label='Detection: {}, Adaptivity: {}'.format(use_SE, use_A) + ) - iters_time_SE = [] - iters_time = [] - iters_time_SE_adapt = [] - iters_time_adapt = [] - times_SE = [] - times = [] - times_SE_adapt = [] - times_adapt = [] - t_switches_SE = [] - t_switches_SE_adapt = [] - - for dt_item in dt_list: - f1 = open(cwd + 'data/battery_dt{}_USETrue_USAFalse_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE = dill.load(f1) - f1.close() - - f2 = open(cwd + 'data/battery_dt{}_USEFalse_USAFalse_{}.dat'.format(dt_item, sweeper), 'rb') - stats = dill.load(f2) - f2.close() - - f3 = open(cwd + 'data/battery_dt{}_USETrue_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_SE_adapt = dill.load(f3) - f3.close() - - f4 = open(cwd + 'data/battery_dt{}_USEFalse_USATrue_{}.dat'.format(dt_item, sweeper), 'rb') - stats_adapt = dill.load(f4) - f4.close() - - # consider iterations before restarts to see what happens - iter_counts_SE_val = get_sorted(stats_SE, type='niter') - iter_counts_SE_adapt_val = get_sorted(stats_SE_adapt, type='niter') - iter_counts_adapt_val = get_sorted(stats_adapt, type='niter') - iter_counts_val = get_sorted(stats, type='niter') - - iters_time_SE.append([v[1] for v in iter_counts_SE_val]) - iters_time_SE_adapt.append([v[1] for v in iter_counts_SE_adapt_val]) - iters_time_adapt.append([v[1] for v in iter_counts_adapt_val]) - iters_time.append([v[1] for v in iter_counts_val]) - - times_SE.append([v[0] for v in iter_counts_SE_val]) - times_SE_adapt.append([v[0] for v in iter_counts_SE_adapt_val]) - times_adapt.append([v[0] for v in iter_counts_adapt_val]) - times.append([v[0] for v in iter_counts_val]) - - switches_SE = get_recomputed(stats_SE, type='switch', sortby='time') - t_switch_SE = [v[1] for v in switches_SE] - t_switches_SE.append(t_switch_SE[-1]) - - switches_SE_adapt = get_recomputed(stats_SE_adapt, type='switch', sortby='time') - t_switch_SE_adapt = [v[1] for v in switches_SE_adapt] - t_switches_SE_adapt.append(t_switch_SE_adapt[-1]) - - if len(dt_list) > 1: - setup_mpl() - fig_iter_all, ax_iter_all = plt_helper.plt.subplots( - nrows=1, ncols=len(dt_list), figsize=(2 * len(dt_list) - 1, 3), sharex='col', sharey='row' + ax[0, i].tick_params(axis='both', which='major', labelsize=16) + ax[0, i].set_ylim(1e-15, 1e0) + ax[0, i].set_yscale('log', base=10) + ax[0, i].set_xlabel(r'Time $t$', fontsize=16) + ax[0, 0].set_ylabel(r'Absolute values of h $|h(v_{C_n}(t))|$', fontsize=16) + ax[0, i].grid(visible=True) + ax[0, i].minorticks_off() + ax[0, i].legend(frameon=True, fontsize=12, loc='lower left') + + fig.savefig( + 'data/{}_state_function_over_time_dt={}_M={}.png'.format(prob_cls_name, dt, M_fix), + dpi=300, + bbox_inches='tight', ) - for col in range(len(dt_list)): - ax_iter_all[col].plot(times[col], iters_time[col], label='SE=F, A=F') - ax_iter_all[col].plot(times_SE[col], iters_time_SE[col], label='SE=T, A=F') - ax_iter_all[col].plot(times_SE_adapt[col], iters_time_SE_adapt[col], '--', label='SE=T, A=T') - ax_iter_all[col].plot(times_adapt[col], iters_time_adapt[col], '--', label='SE=F, A=T') - ax_iter_all[col].axvline(x=t_switches_SE[col], linestyle='--', linewidth=0.5, color='k', label='Switch') - ax_iter_all[col].set_title(r'$\Delta t_\mathrm{initial}$=%s' % dt_list[col]) - ax_iter_all[col].set_ylim(0, maxiter + 2) - ax_iter_all[col].set_xlabel('Time', fontsize=6) - ax_iter_all[col].tick_params(axis='both', which='major', labelsize=6) - - if col == 0: - ax_iter_all[col].set_ylabel('Number iterations', fontsize=6) - - ax_iter_all[col].legend(frameon=False, fontsize=6, loc='upper right') - else: - setup_mpl() - fig_iter_all, ax_iter_all = plt_helper.plt.subplots(nrows=1, ncols=1, figsize=(3, 3)) - - ax_iter_all.plot(times[0], iters_time[0], label='SE=False') - ax_iter_all.plot(times_SE[0], iters_time_SE[0], label='SE=True') - ax_iter_all.plot(times_SE_adapt[0], iters_time_SE_adapt[0], '--', label='SE=T, A=T') - ax_iter_all.plot(times_adapt[0], iters_time_adapt[0], '--', label='SE=F, A=T') - ax_iter_all.axvline(x=t_switches_SE[0], linestyle='--', linewidth=0.5, color='k', label='Switch') - ax_iter_all.set_title(r'$\Delta t_\mathrm{initial}$=%s' % dt_list[0]) - ax_iter_all.set_ylim(0, maxiter + 2) - ax_iter_all.set_xlabel('Time', fontsize=6) - ax_iter_all.tick_params(axis='both', which='major', labelsize=6) - - ax_iter_all.set_ylabel('Number iterations', fontsize=6) - ax_iter_all.legend(frameon=False, fontsize=6, loc='upper right') - - plt_helper.plt.tight_layout() - fig_iter_all.savefig('data/iters_{}.png'.format(sweeper), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig_iter_all) + plt_helper.plt.close(fig) if __name__ == "__main__": - run() + run_estimation_check() diff --git a/pySDC/projects/PinTSimE/estimation_check_2capacitors.py b/pySDC/projects/PinTSimE/estimation_check_2capacitors.py deleted file mode 100644 index ac269d744c..0000000000 --- a/pySDC/projects/PinTSimE/estimation_check_2capacitors.py +++ /dev/null @@ -1,245 +0,0 @@ -import numpy as np -import dill -from pathlib import Path - -from pySDC.helpers.stats_helper import get_sorted -from pySDC.implementations.problem_classes.Battery import battery_n_capacitors -from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.projects.PinTSimE.battery_model import controller_run, generate_description, get_recomputed - -from pySDC.projects.PinTSimE.battery_2capacitors_model import ( - LogEvent, - check_solution, - proof_assertions_description, - proof_assertions_time, -) - -from pySDC.projects.PinTSimE.piline_model import setup_mpl -import pySDC.helpers.plot_helper as plt_helper - -from pySDC.implementations.hooks.log_solution import LogSolution - -from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator - - -def run(cwd='./'): - """ - Routine to check the differences between using a switch estimator or not. - - Parameters - ---------- - cwd : str - Current working directory. - """ - - dt_list = [4e-1, 4e-2] - t0 = 0.0 - Tend = 3.5 - - problem_classes = [battery_n_capacitors] - sweeper_classes = [imex_1st_order] - num_nodes = 4 - restol = -1 - maxiter = 8 - - ncapacitors = 2 - alpha = 5.0 - V_ref = np.array([1.0, 1.0]) - C = np.array([1.0, 1.0]) - - problem_params = dict() - problem_params['ncapacitors'] = ncapacitors - problem_params['C'] = C - problem_params['alpha'] = alpha - problem_params['V_ref'] = V_ref - - hook_class = [LogSolution, LogEvent] - - use_switch_estimator = [True, False] - restarts_all = [] - restarts_dict = dict() - for problem, sweeper in zip(problem_classes, sweeper_classes): - for dt_item in dt_list: - for use_SE in use_switch_estimator: - description, controller_params = generate_description( - dt_item, - problem, - sweeper, - num_nodes, - hook_class, - False, - use_SE, - problem_params, - restol, - maxiter, - 1, - 1e-8, - ) - - # Assertions - proof_assertions_description(description, False, use_SE) - - proof_assertions_time(dt_item, Tend, V_ref, alpha) - - stats = controller_run(description, controller_params, False, use_SE, t0, Tend) - - if use_SE: - switches = get_recomputed(stats, type='switch', sortby='time') - assert len(switches) >= 2, f"Expected at least 2 switches for dt: {dt_item}, got {len(switches)}!" - - check_solution(stats, dt_item, use_SE) - - fname = 'data/{}_dt{}_USE{}.dat'.format(problem.__name__, dt_item, use_SE) - f = open(fname, 'wb') - dill.dump(stats, f) - f.close() - - if use_SE: - restarts_dict[dt_item] = np.array(get_sorted(stats, type='restart', recomputed=None)) - restarts = restarts_dict[dt_item][:, 1] - restarts_all.append(np.sum(restarts)) - print(f"Restarts for dt: {dt_item:.2e} -- {np.sum(restarts):.0f}") - - V_ref = description['problem_params']['V_ref'] - - val_switch_all = [] - diff_true_all1 = [] - diff_false_all_before1 = [] - diff_false_all_after1 = [] - diff_true_all2 = [] - diff_false_all_before2 = [] - diff_false_all_after2 = [] - restarts_dt_switch1 = [] - restarts_dt_switch2 = [] - for dt_item in dt_list: - f1 = open(cwd + 'data/{}_dt{}_USETrue.dat'.format(problem.__name__, dt_item), 'rb') - stats_true = dill.load(f1) - f1.close() - - f2 = open(cwd + 'data/{}_dt{}_USEFalse.dat'.format(problem.__name__, dt_item), 'rb') - stats_false = dill.load(f2) - f2.close() - - switches = get_recomputed(stats_true, type='switch', sortby='time') - t_switch = [v[1] for v in switches] - val_switch_all.append([t_switch[0], t_switch[1]]) - - vC1_true = [me[1][1] for me in get_sorted(stats_true, type='u', recomputed=False)] - vC2_true = [me[1][2] for me in get_sorted(stats_true, type='u', recomputed=False)] - vC1_false = [me[1][1] for me in get_sorted(stats_false, type='u', recomputed=False)] - vC2_false = [me[1][2] for me in get_sorted(stats_false, type='u', recomputed=False)] - - diff_true1 = vC1_true - V_ref[0] - diff_true2 = vC2_true - V_ref[1] - diff_false1 = vC1_false - V_ref[0] - diff_false2 = vC2_false - V_ref[1] - - t_true = [me[0] for me in get_sorted(stats_true, type='u', recomputed=False)] - t_false = [me[0] for me in get_sorted(stats_false, type='u', recomputed=False)] - - diff_true_all1.append([diff_true1[m] for m in range(len(t_true)) if abs(t_true[m] - t_switch[0]) <= 1e-17]) - diff_true_all2.append([diff_true2[np.argmin([abs(me - t_switch[1]) for me in t_true])]]) - - diff_false_all_before1.append( - [diff_false1[m - 1] for m in range(1, len(t_false)) if t_false[m - 1] < t_switch[0] < t_false[m]] - ) - diff_false_all_after1.append( - [diff_false1[m] for m in range(1, len(t_false)) if t_false[m - 1] < t_switch[0] < t_false[m]] - ) - - diff_false_all_before2.append( - [diff_false2[m - 1] for m in range(1, len(t_false)) if t_false[m - 1] < t_switch[1] < t_false[m]] - ) - diff_false_all_after2.append( - [diff_false2[m] for m in range(1, len(t_false)) if t_false[m - 1] < t_switch[1] < t_false[m]] - ) - - restarts_dt = restarts_dict[dt_item] - for i in range(len(restarts_dt[:, 0])): - if np.isclose(restarts_dt[i, 0], t_switch[0], atol=1e-15): - restarts_dt_switch1.append(np.sum(restarts_dt[0 : i - 1, 1])) - break - - for i in range(len(restarts_dt[:, 0])): - if np.isclose(restarts_dt[i, 0], t_switch[1], atol=1e-15): - restarts_dt_switch2.append(np.sum(restarts_dt[i - 2 :, 1])) - break - - setup_mpl() - fig1, ax1 = plt_helper.plt.subplots(1, 1, figsize=(4.5, 3)) - ax1.set_title('Time evolution of $v_{C_{1}}-V_{ref1}$') - ax1.plot(t_true, diff_true1, label='SE=True', color='#ff7f0e') - ax1.plot(t_false, diff_false1, label='SE=False', color='#1f77b4') - ax1.axvline(x=t_switch[0], linestyle='--', color='k', label='Switch1') - ax1.legend(frameon=False, fontsize=10, loc='lower left') - ax1.set_yscale('symlog', linthresh=1e-5) - ax1.set_xlabel('Time') - - fig1.savefig('data/difference_estimation_vC1_dt{}.png'.format(dt_item), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig1) - - setup_mpl() - fig2, ax2 = plt_helper.plt.subplots(1, 1, figsize=(4.5, 3)) - ax2.set_title('Time evolution of $v_{C_{2}}-V_{ref2}$') - ax2.plot(t_true, diff_true2, label='SE=True', color='#ff7f0e') - ax2.plot(t_false, diff_false2, label='SE=False', color='#1f77b4') - ax2.axvline(x=t_switch[1], linestyle='--', color='k', label='Switch2') - ax2.legend(frameon=False, fontsize=10, loc='lower left') - ax2.set_yscale('symlog', linthresh=1e-5) - ax2.set_xlabel('Time') - - fig2.savefig('data/difference_estimation_vC2_dt{}.png'.format(dt_item), dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig2) - - setup_mpl() - fig1, ax1 = plt_helper.plt.subplots(1, 1, figsize=(3, 3)) - ax1.set_title("Difference $v_{C_{1}}-V_{ref1}$") - pos1 = ax1.plot(dt_list, diff_false_all_before1, 'rs-', label='SE=False - before switch1') - pos2 = ax1.plot(dt_list, diff_false_all_after1, 'bd-', label='SE=False - after switch1') - pos3 = ax1.plot(dt_list, diff_true_all1, 'kd-', label='SE=True') - ax1.set_xticks(dt_list) - ax1.set_xticklabels(dt_list) - ax1.set_xscale('log', base=10) - ax1.set_yscale('symlog', linthresh=1e-10) - ax1.set_ylim(-2, 2) - ax1.set_xlabel(r'$\Delta t$') - - restart_ax = ax1.twinx() - restarts = restart_ax.plot(dt_list, restarts_dt_switch1, 'cs--', label='Restarts') - restart_ax.set_ylabel('Restarts') - - lines = pos1 + pos2 + pos3 + restarts - labels = [l.get_label() for l in lines] - ax1.legend(lines, labels, frameon=False, fontsize=8, loc='center right') - - fig1.savefig('data/diffs_estimation_vC1.png', dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig1) - - setup_mpl() - fig2, ax2 = plt_helper.plt.subplots(1, 1, figsize=(3, 3)) - ax2.set_title("Difference $v_{C_{2}}-V_{ref2}$") - pos1 = ax2.plot(dt_list, diff_false_all_before2, 'rs-', label='SE=False - before switch2') - pos2 = ax2.plot(dt_list, diff_false_all_after2, 'bd-', label='SE=False - after switch2') - pos3 = ax2.plot(dt_list, diff_true_all2, 'kd-', label='SE=True') - ax2.set_xticks(dt_list) - ax2.set_xticklabels(dt_list) - ax2.set_xscale('log', base=10) - ax2.set_yscale('symlog', linthresh=1e-10) - ax2.set_ylim(-2, 2) - ax2.set_xlabel(r'$\Delta t$') - - restart_ax = ax2.twinx() - restarts = restart_ax.plot(dt_list, restarts_dt_switch2, 'cs--', label='Restarts') - restart_ax.set_ylabel('Restarts') - - lines = pos1 + pos2 + pos3 + restarts - labels = [l.get_label() for l in lines] - ax2.legend(lines, labels, frameon=False, fontsize=8, loc='center right') - - fig2.savefig('data/diffs_estimation_vC2.png', dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig2) - - -if __name__ == "__main__": - run() diff --git a/pySDC/projects/PinTSimE/hardcoded_solutions.py b/pySDC/projects/PinTSimE/hardcoded_solutions.py new file mode 100644 index 0000000000..27733a12ec --- /dev/null +++ b/pySDC/projects/PinTSimE/hardcoded_solutions.py @@ -0,0 +1,480 @@ +import numpy as np + +from pySDC.core.Errors import ParameterError + + +def testSolution(u_num, prob_cls_name, dt, use_adaptivity, use_detection): + r""" + Test for numerical solution if values satisfy hardcoded values. + + Note + ---- + Only the items are tested which does make sense for any problem. For instance, `getDataDict` stores the global error + after each step **for each problem class**. Since an exact solution is only available for e.g. ``DiscontinuousTestODE``, + a test is only be done for this problem class. + + Hardcoded solutions are computed for only one collocation node ``M_fix``, which is specified for each problem in the + related files, see ``pySDC.projects.PinTSimE.battery_model``, ``pySDC.projects.PinTSimE.estimation_check`` and + ``pySDC.projects.PinTSimE.discontinuous_test_ODE``. + + Parameters + ---------- + u_num : dict + Contains the numerical solution together with event time found by event detection, step sizes adjusted + via adaptivity and/or switch estimation. + prob_cls_name : str + Indicates which problem class is tested. + dt : float + (Initial) step sizes used for the simulation. + use_adaptivity : bool + Indicates whether adaptivity is used in the simulation or not. + use_detection : bool + Indicates whether discontinuity handling is used in the simulation or not. + """ + + unknowns = u_num['unknowns'] + u_num_tmp = {unknown: u_num[unknown][-1] for unknown in unknowns} + + got = {} + got = u_num_tmp + + if prob_cls_name == 'battery': + if use_adaptivity and 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, + } + 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, + } + got.update( + { + 't_switches': u_num['t_switches'], + 'dt': u_num['dt'][-1, 1], + 'e_em': u_num['e_em'][-1, 1], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + elif not use_adaptivity and 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, + } + elif dt == 1e-3: + expected = { + 'iL': 0.5408462989990014, + 'vC': 1.0, + 't_switches': [0.18232155679395023], + 'sum_restarts': 2.0, + 'sum_niters': 2424.0, + } + got.update( + { + 't_switches': u_num['t_switches'], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + if use_adaptivity and not use_detection: + msg = f"Error when using adaptivity for {prob_cls_name} for dt={dt:.1e}:" + if dt == 1e-2: + expected = { + 'iL': 0.4433805288639916, + 'vC': 0.90262388393713, + 'dt': 0.18137307612335937, + 'e_em': 2.7177844974524135e-09, + 'sum_restarts': 0.0, + 'sum_niters': 24.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.3994744179584864, + 'vC': 0.9679037468770668, + 'dt': 0.1701392217033212, + 'e_em': 2.0992988458701234e-09, + 'sum_restarts': 0.0, + 'sum_niters': 32.0, + } + got.update( + { + 'dt': u_num['dt'][-1, 1], + 'e_em': u_num['e_em'][-1, 1], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + elif not use_adaptivity and not use_detection: + msg = f'Error for {prob_cls_name} for dt={dt:.1e}:' + if dt == 1e-2: + expected = { + 'iL': 0.5456625861551172, + 'vC': 0.9973251377556902, + 'sum_niters': 240.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.538639340748308, + 'vC': 0.9994050706403905, + 'sum_niters': 2400.0, + } + got.update( + { + 'sum_niters': u_num['sum_niters'], + } + ) + + elif prob_cls_name == 'battery_implicit': + if use_adaptivity and 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, + 'e_em': 2.220446049250313e-16, + 'sum_restarts': 3.0, + 'sum_niters': 56.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.5393867577837699, + 'vC': 1.0000000000250129, + 't_switches': [0.1823215568036829], + 'dt': 0.015641237833012522, + 'e_em': 2.220446049250313e-16, + 'sum_restarts': 14.0, + 'sum_niters': 328.0, + } + got.update( + { + 't_switches': u_num['t_switches'], + 'dt': u_num['dt'][-1, 1], + 'e_em': u_num['e_em'][-1, 1], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + elif not use_adaptivity and 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, + } + elif dt == 1e-3: + expected = { + 'iL': 0.5407340515794409, + 'vC': 1.0000000000010945, + 't_switches': [0.182321556796257], + 'sum_restarts': 3.0, + 'sum_niters': 2440.0, + } + got.update( + { + 't_switches': u_num['t_switches'], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + if use_adaptivity and not use_detection: + msg = f"Error when using adaptivity for {prob_cls_name} for dt={dt:.1e}:" + if dt == 1e-2: + expected = { + 'iL': 0.4694087102919169, + 'vC': 0.9026238839418407, + 'dt': 0.18137307612335937, + 'e_em': 2.3469713394952407e-09, + 'sum_restarts': 0.0, + 'sum_niters': 24.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.39947441811958956, + 'vC': 0.9679037468856341, + 'dt': 0.1701392217033212, + 'e_em': 1.147640815712947e-09, + 'sum_restarts': 0.0, + 'sum_niters': 32.0, + } + got.update( + { + 'dt': u_num['dt'][-1, 1], + 'e_em': u_num['e_em'][-1, 1], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + elif not use_adaptivity and not use_detection: + msg = f'Error for {prob_cls_name} for dt={dt:.1e}:' + if dt == 1e-2: + expected = { + 'iL': 0.5456915668459889, + 'vC': 0.9973251377578705, + 'sum_niters': 240.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.5386399890100035, + 'vC': 0.999405070641239, + 'sum_niters': 2400.0, + } + got.update( + { + 'sum_niters': u_num['sum_niters'], + } + ) + + elif prob_cls_name == 'battery_n_capacitors': + if use_adaptivity and 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, + } + 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, + } + got.update( + { + 't_switches': u_num['t_switches'], + 'dt': u_num['dt'][-1, 1], + 'e_em': u_num['e_em'][-1, 1], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + elif not use_adaptivity and 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, + } + 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, + } + got.update( + { + 't_switches': u_num['t_switches'], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + if use_adaptivity and not use_detection: + msg = f"Error when using adaptivity for {prob_cls_name} for dt={dt:.1e}:" + if dt == 1e-2: + expected = { + 'iL': 0.15890544838473294, + 'vC1': 0.8806086293336285, + 'vC2': 0.9915019206727803, + 'dt': 0.38137307612335936, + 'e_em': 4.145817911194172e-09, + 'sum_restarts': 0.0, + 'sum_niters': 24.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.15422467570971707, + 'vC1': 0.8756872272783145, + 'vC2': 0.9971015415168025, + 'dt': 0.3701392217033212, + 'e_em': 3.6970297934146856e-09, + 'sum_restarts': 0.0, + 'sum_niters': 32.0, + } + got.update( + { + 'dt': u_num['dt'][-1, 1], + 'e_em': u_num['e_em'][-1, 1], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + } + ) + + elif not use_adaptivity and not use_detection: + msg = f'Error for {prob_cls_name} for dt={dt:.1e}:' + if dt == 1e-2: + expected = { + 'iL': 0.5939796175551723, + 'vC1': 0.9973251377556902, + 'vC2': 0.9973251377466236, + 'sum_niters': 400.0, + } + elif dt == 1e-3: + expected = { + 'iL': 0.6107051960885036, + 'vC1': 0.9994050706403905, + 'vC2': 0.9997382611893499, + 'sum_niters': 4000.0, + } + got.update( + { + 'sum_niters': u_num['sum_niters'], + } + ) + + elif prob_cls_name == 'DiscontinuousTestODE': + if not use_adaptivity and 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], + } + 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], + } + got.update( + { + 't_switches': u_num['t_switches'], + 'sum_restarts': u_num['sum_restarts'], + 'sum_niters': u_num['sum_niters'], + 'e_global': u_num['e_global'][-1, 1], + 'e_event': u_num['e_event'], + } + ) + + elif not use_adaptivity and not use_detection: + msg = f"Error for {prob_cls_name} for dt={dt:.1e}:" + if dt == 1e-2: + expected = { + 'u': 5.9805345175338225, + 'sum_niters': 527.0, + 'e_global': 0.009855041056925806, + } + elif dt == 1e-3: + expected = { + 'u': 5.9737411566014105, + 'sum_niters': 2226.0, + 'e_global': 0.0005763403865515215, + } + got.update( + { + 'sum_niters': u_num['sum_niters'], + 'e_global': u_num['e_global'][-1, 1], + } + ) + + elif prob_cls_name == 'piline': + if not use_adaptivity and not use_detection: + msg = f"Error for {prob_cls_name} for dt={dt:.1e}:" + if dt == 5e-2: + expected = { + 'vC1': 83.97535096068474, + 'vC2': 80.52142367760014, + 'iLp': 16.096505806313573, + 'sum_niters': 2045.0, + } + elif dt == 1e-2: + expected = { + 'vC1': 83.97462422132232, + 'vC2': 80.52135774600202, + 'iLp': 16.09884649820726, + 'sum_niters': 7206.0, + } + got.update( + { + 'sum_niters': u_num['sum_niters'], + } + ) + + elif prob_cls_name == 'buck_converter': + if not use_adaptivity and not 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, + } + elif dt == 1e-5: + expected = { + 'vC1': 9.891508522329485, + 'vC2': 4.70939963429714, + 'iLp': -0.32177442457657557, + 'sum_niters': 8262.0, + } + got.update( + { + 'sum_niters': u_num['sum_niters'], + } + ) + + else: + raise ParameterError(f"For {prob_cls_name} there is no test implemented here!") + + 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 + 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/piline_model.py b/pySDC/projects/PinTSimE/piline_model.py index d576cfcf5d..fc2f5f56f6 100644 --- a/pySDC/projects/PinTSimE/piline_model.py +++ b/pySDC/projects/PinTSimE/piline_model.py @@ -1,183 +1,64 @@ -import matplotlib as mpl -import numpy as np -import dill -from pathlib import Path - -mpl.use('Agg') - -from pySDC.helpers.stats_helper import get_sorted -from pySDC.core.Collocation import CollBase as Collocation -from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.implementations.problem_classes.Piline import piline from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -import pySDC.helpers.plot_helper as plt_helper - -from pySDC.core.Hooks import hooks - -class log_data(hooks): - def post_step(self, step, level_number): - super(log_data, self).post_step(step, level_number) +from pySDC.projects.PinTSimE.battery_model import runSimulation - # some abbreviations - L = step.levels[level_number] - - L.sweep.compute_end_point() - - self.add_to_stats( - process=step.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='v1', - value=L.uend[0], - ) - self.add_to_stats( - process=step.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='v2', - value=L.uend[1], - ) - self.add_to_stats( - process=step.status.slot, - time=L.time, - level=L.level_index, - iter=0, - sweep=L.status.sweep, - type='p3', - value=L.uend[2], - ) +from pySDC.implementations.hooks.log_solution import LogSolution def main(): + r""" + Executes the simulation. + + Note + ---- + Hardcoded solutions for battery models in `pySDC.projects.PinTSimE.hardcoded_solutions` are only computed for + ``dt_list=[1e-2, 1e-3]`` and ``M_fix=4``. Hence changing ``dt_list`` and ``M_fix`` to different values could arise + an ``AssertionError``. """ - A simple test program to do SDC/PFASST runs for the Piline model - """ - - # initialize level parameters - level_params = dict() - level_params['restol'] = 1e-10 - level_params['dt'] = 0.25 - - # initialize sweeper parameters - sweeper_params = dict() - sweeper_params['collocation_class'] = Collocation - sweeper_params['quad_type'] = 'LOBATTO' - sweeper_params['num_nodes'] = 5 - # sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part - - # initialize step parameters - step_params = dict() - step_params['maxiter'] = 20 - - # initialize controller parameters - controller_params = dict() - controller_params['logger_level'] = 30 - controller_params['hook_class'] = log_data - - # fill description dictionary for easy step instantiation - # keep in mind: default params are used for the problem, but can be changed - description = dict() - description['problem_class'] = piline # pass problem class - description['sweeper_class'] = imex_1st_order # pass sweeper - description['sweeper_params'] = sweeper_params # pass sweeper parameters - description['level_params'] = level_params # pass level parameters - description['step_params'] = step_params # pass step parameters - - assert 'errtol' not in description['step_params'].keys(), "No exact or reference solution known to compute error" - - # set time parameters - t0 = 0.0 - Tend = 15 - num_procs = 1 - - # instantiate controller - controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) - - # get initial values on finest level - P = controller.MS[0].levels[0].prob - uinit = P.u_exact(t0) - - # call main function to get things done... - uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - - Path("data").mkdir(parents=True, exist_ok=True) - fname = 'data/piline.dat' - f = open(fname, 'wb') - dill.dump(stats, f) - f.close() - - # filter statistics by number of iterations - iter_counts = get_sorted(stats, type='niter', sortby='time') - - # compute and print statistics - min_iter = 20 - max_iter = 0 - - f = open('data/piline_out.txt', 'w') - niters = np.array([item[1] for item in iter_counts]) - out = ' Mean number of iterations: %4.2f' % np.mean(niters) - f.write(out + '\n') - print(out) - for item in iter_counts: - out = 'Number of iterations for time %4.2f: %1i' % item - f.write(out + '\n') - print(out) - min_iter = min(min_iter, item[1]) - max_iter = max(max_iter, item[1]) - - assert np.mean(niters) <= 10, "Mean number of iterations is too high, got %s" % np.mean(niters) - f.close() - - plot_voltages() - - -def plot_voltages(cwd='./'): - f = open(cwd + 'data/piline.dat', 'rb') - stats = dill.load(f) - f.close() - - # convert filtered statistics to list of iterations count, sorted by process - v1 = get_sorted(stats, type='v1', sortby='time') - v2 = get_sorted(stats, type='v2', sortby='time') - p3 = get_sorted(stats, type='p3', sortby='time') - - times = [v[0] for v in v1] - - setup_mpl() - fig, ax = plt_helper.plt.subplots(1, 1, figsize=(4.5, 3)) - ax.plot(times, [v[1] for v in v1], linewidth=1, label=r'$v_{C_1}$') - ax.plot(times, [v[1] for v in v2], linewidth=1, label=r'$v_{C_2}$') - ax.plot(times, [v[1] for v in p3], linewidth=1, label=r'$i_{L_\pi}$') - ax.legend(frameon=False, fontsize=12, loc='center right') - - ax.set_xlabel('Time') - ax.set_ylabel('Energy') - - fig.savefig('data/piline_model_solution.png', dpi=300, bbox_inches='tight') - plt_helper.plt.close(fig) + # --- defines parameters for sweeper ---- + M_fix = 3 + sweeper_params = { + 'num_nodes': M_fix, + 'quad_type': 'LOBATTO', + 'QI': 'IE', + } -def setup_mpl(fontsize=8): - plt_helper.setup_mpl(reset=True) + # --- defines parameters for event detection ---- + handling_params = { + 'restol': 1e-12, + 'maxiter': 8, + 'max_restarts': None, + 'recomputed': None, + 'tol_event': None, + 'alpha': None, + 'exact_event_time_avail': None, + } - style_options = { - "font.family": "sans-serif", - "font.serif": "Computer Modern Sans Serif", - "font.sans-serif": "Computer Modern Sans Serif", - "font.monospace": "Computer Modern Sans Serif", - "axes.labelsize": 12, # LaTeX default is 10pt font. - "legend.fontsize": 13, # Make the legend/label fonts a little smaller - "axes.xmargin": 0.03, - "axes.ymargin": 0.03, - "lines.linewidth": 1, # Make the plot lines a little smaller + # ---- all parameters are stored in this dictionary - note: defaults are used for the problem ---- + all_params = { + 'problem_params': {}, + 'sweeper_params': sweeper_params, + 'handling_params': handling_params, } - mpl.rcParams.update(style_options) + hook_class = [LogSolution] + + use_detection = [False] + use_adaptivity = [False] + + _ = runSimulation( + problem=piline, + sweeper=imex_1st_order, + all_params=all_params, + use_adaptivity=use_adaptivity, + use_detection=use_detection, + hook_class=hook_class, + interval=(0.0, 15.0), + dt_list=[5e-2, 1e-2], + nnodes=[M_fix], + ) if __name__ == "__main__": diff --git a/pySDC/projects/PinTSimE/switch_estimator.py b/pySDC/projects/PinTSimE/switch_estimator.py index d9059b91c1..3c438c134a 100644 --- a/pySDC/projects/PinTSimE/switch_estimator.py +++ b/pySDC/projects/PinTSimE/switch_estimator.py @@ -13,8 +13,12 @@ class SwitchEstimator(ConvergenceController): """ def setup(self, controller, params, description): - """ - Function sets default variables to handle with the switch at the beginning. + 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 Parameters ---------- @@ -39,8 +43,9 @@ def setup(self, controller, params, description): ) defaults = { - 'control_order': 100, + 'control_order': 0, 'nodes': coll.nodes, + 'tol_zero': 1e-13, } return {**defaults, **params} @@ -92,11 +97,20 @@ def get_new_step_size(self, controller, S, **kwargs): ) # when the state function is already close to zero the event is already resolved well - if abs(state_function[-1]) <= self.params.tol: - self.log("Is already close enough to one of the end point!", S) + 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] + boundary = 'left' + elif abs(state_function[-1]) <= self.params.tol_zero: + boundary = 'right' + t_switch = 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_event_time( - controller.hooks[0], S.status.slot, L.time, L.level_index, L.status.sweep, t_interp[-1] + controller.hooks[0], S.status.slot, L.time, L.level_index, L.status.sweep, t_switch ) + L.prob.count_switches() self.status.is_zero = True @@ -104,14 +118,32 @@ def get_new_step_size(self, controller, S, **kwargs): 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) + 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 state_function]), + ) if L.time < self.status.t_switch < L.time + L.dt: - dt_switch = self.status.t_switch - L.time + dt_switch = (self.status.t_switch - L.time) * self.params.alpha if ( abs(self.status.t_switch - L.time) <= self.params.tol or abs((L.time + L.dt) - self.status.t_switch) <= self.params.tol ): - self.log(f"Switch located at time {self.status.t_switch:.12f}", S) + self.log(f"Switch located at time {self.status.t_switch:.15f}", S) L.prob.t_switch = self.status.t_switch self.log_event_time( controller.hooks[0], @@ -125,7 +157,7 @@ def get_new_step_size(self, controller, S, **kwargs): L.prob.count_switches() else: - self.log(f"Located Switch at time {self.status.t_switch:.12f} is outside the range", S) + self.log(f"Located Switch at time {self.status.t_switch:.15f} is outside the range", S) # when an event is found, step size matching with this event should be preferred dt_planned = L.status.dt_new if L.status.dt_new is not None else L.params.dt @@ -137,8 +169,7 @@ def get_new_step_size(self, controller, S, **kwargs): else: # event occurs on L.time or L.time + L.dt; no restart necessary boundary = 'left boundary' if self.status.t_switch == L.time else 'right boundary' - self.log(f"Estimated switch {self.status.t_switch:.12f} occurs at {boundary}", S) - + self.log(f"Estimated switch {self.status.t_switch:.15f} occurs at {boundary}", S) self.log_event_time( controller.hooks[0], S.status.slot, @@ -240,10 +271,10 @@ def get_switch(t_interp, state_function, m_guess): Returns ------- t_switch : float - Time point of the founded switch. + Time point of found event. """ - Interpolator = sp.interpolate.BarycentricInterpolator(t_interp, state_function) + LagrangeInterpolator = LagrangeInterpolation(t_interp, state_function) def p(t): """ @@ -259,7 +290,7 @@ def p(t): p(t) : float The value of the interpolated function at time t. """ - return Interpolator.__call__(t) + return LagrangeInterpolator.eval(t) def fprime(t): """ @@ -275,13 +306,12 @@ def fprime(t): dp : float Derivative of interpolation p at time t. """ - dt = 1e-8 - dp = (p(t + dt) - p(t)) / dt + dt_FD = 1e-10 + dp = (p(t + dt_FD) - p(t)) / dt_FD # forward difference return dp - newton_tol, newton_maxiter = 1e-8, 50 + newton_tol, newton_maxiter = 1e-15, 100 t_switch = newton(t_interp[m_guess], p, fprime, newton_tol, newton_maxiter) - return t_switch @staticmethod @@ -353,7 +383,49 @@ def newton(x0, p, fprime, newton_tol, newton_maxiter): n += 1 root = x0 - msg = "Newton's method took {} iterations".format(n) - print(msg) return root + + +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): + """ + Evaluates the Lagrange interpolation at time t. + + Parameters + ---------- + t : float + Time where interpolation is computed. + + Returns + ------- + p : float + Value of interpolant at time t. + """ + 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_battery_2capacitors_model.py b/pySDC/tests/test_projects/test_pintsime/test_battery_2capacitors_model.py deleted file mode 100644 index cb7dccd0cf..0000000000 --- a/pySDC/tests/test_projects/test_pintsime/test_battery_2capacitors_model.py +++ /dev/null @@ -1,8 +0,0 @@ -import pytest - - -@pytest.mark.base -def test_main(): - from pySDC.projects.PinTSimE.battery_2capacitors_model import run - - run() diff --git a/pySDC/tests/test_projects/test_pintsime/test_battery_model.py b/pySDC/tests/test_projects/test_pintsime/test_battery_model.py index 7dc540868a..1a59a44045 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_battery_model.py +++ b/pySDC/tests/test_projects/test_pintsime/test_battery_model.py @@ -3,6 +3,6 @@ @pytest.mark.base def test_main(): - from pySDC.projects.PinTSimE.battery_model import run + from pySDC.projects.PinTSimE.battery_model import main - run() + main() diff --git a/pySDC/tests/test_projects/test_pintsime/test_estimation_check.py b/pySDC/tests/test_projects/test_pintsime/test_estimation_check.py index 01ec3fed73..25a2813a41 100644 --- a/pySDC/tests/test_projects/test_pintsime/test_estimation_check.py +++ b/pySDC/tests/test_projects/test_pintsime/test_estimation_check.py @@ -3,6 +3,6 @@ @pytest.mark.base def test_main(): - from pySDC.projects.PinTSimE.estimation_check import run + from pySDC.projects.PinTSimE.estimation_check import run_estimation_check - run() + run_estimation_check() diff --git a/pySDC/tests/test_projects/test_pintsime/test_estimation_check_2capacitors.py b/pySDC/tests/test_projects/test_pintsime/test_estimation_check_2capacitors.py deleted file mode 100644 index 09751a7897..0000000000 --- a/pySDC/tests/test_projects/test_pintsime/test_estimation_check_2capacitors.py +++ /dev/null @@ -1,8 +0,0 @@ -import pytest - - -@pytest.mark.base -def test_main(): - from pySDC.projects.PinTSimE.estimation_check_2capacitors import run - - run()