diff --git a/pySDC/projects/Second_orderSDC/harmonic_oscillator_params.py b/pySDC/projects/Second_orderSDC/harmonic_oscillator_params.py index 504af7ec95..06382fe2a9 100644 --- a/pySDC/projects/Second_orderSDC/harmonic_oscillator_params.py +++ b/pySDC/projects/Second_orderSDC/harmonic_oscillator_params.py @@ -3,7 +3,6 @@ from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order - def harmonic_oscillator_params(): """ Routine to compute modules of the stability function diff --git a/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py b/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py index add9cdab07..5e8246367a 100644 --- a/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py +++ b/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py @@ -4,13 +4,15 @@ if __name__ == '__main__': # Define lists for the number of nodes and maximum iterations - description=harmonic_oscillator_params() - helper_params={'quad_type_list':['GAUSS', 'LOBATTO'], - 'Num_iter':(2, 2), - 'num_nodes_list':np.arange(3, 6, 1), - 'max_iter_list':np.arange(2, 10, 1)} - description['helper_params']=helper_params - points=((1, 100), (3, 100), (10, 100)) + description = harmonic_oscillator_params() + helper_params = { + 'quad_type_list': ['GAUSS', 'LOBATTO'], + 'Num_iter': (2, 2), + 'num_nodes_list': np.arange(3, 6, 1), + 'max_iter_list': np.arange(2, 10, 1), + } + description['helper_params'] = helper_params + points = ((1, 100), (3, 100), (10, 100)) # Iterate through points and perform stability check for ii in points: compute_and_check_stability(description, ii, check_stability=True) diff --git a/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stab_interval.py b/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stab_interval.py index de5e0778e1..6effbff377 100644 --- a/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stab_interval.py +++ b/pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stab_interval.py @@ -4,13 +4,15 @@ if __name__ == '__main__': # Define lists for the number of nodes and maximum iterations - description=harmonic_oscillator_params() - helper_params={'quad_type_list':['GAUSS'], - 'Num_iter':(2000, 1), - 'num_nodes_list':np.arange(2, 7, 1), - 'max_iter_list':np.arange(1, 11, 1)} - description['helper_params']=helper_params - points=((100, 1e-10),) + description = harmonic_oscillator_params() + helper_params = { + 'quad_type_list': ['GAUSS'], + 'Num_iter': (2000, 1), + 'num_nodes_list': np.arange(2, 7, 1), + 'max_iter_list': np.arange(1, 11, 1), + } + description['helper_params'] = helper_params + points = ((100, 1e-10),) # Iterate through points and perform stability check for ii in points: compute_and_check_stability(description, ii, compute_interval=True) diff --git a/pySDC/projects/Second_orderSDC/penningtrap_Simulation.py b/pySDC/projects/Second_orderSDC/penningtrap_Simulation.py index 27a3cf392c..480cfe064b 100644 --- a/pySDC/projects/Second_orderSDC/penningtrap_Simulation.py +++ b/pySDC/projects/Second_orderSDC/penningtrap_Simulation.py @@ -6,13 +6,16 @@ from pySDC.implementations.sweeper_classes.Runge_Kutta_Nystrom import RKN, Velocity_Verlet from pySDC.projects.Second_orderSDC.plot_helper import PlotManager + class ComputeError(PlotManager): """ This class generates data for plots and computations for Second-order SDC """ def __init__(self, controller_params, description, time_iter=3, K_iter=(1, 2, 3), Tend=2, axes=(1,), cwd=''): - super().__init__(controller_params, description, time_iter=time_iter, K_iter=K_iter, Tend=Tend, axes=axes, cwd='') + super().__init__( + controller_params, description, time_iter=time_iter, K_iter=K_iter, Tend=Tend, axes=axes, cwd='' + ) def run_local_error(self): """ @@ -138,7 +141,7 @@ def compute_global_error_data(self, Picard=False, RK=False, VV=False, work_count "Time_steps/Work_counter | Order_pos | Abs_error_position | Order_vel | Abs_error_velocity\n" ) - cont = 2 if self.time_iter == 3 else 2**abs(3 - self.time_iter) + cont = 2 if self.time_iter == 3 else 2 ** abs(3 - self.time_iter) cont = cont if not Picard else dt_cont for ii in range(0, self.time_iter): @@ -153,7 +156,7 @@ def compute_global_error_data(self, Picard=False, RK=False, VV=False, work_count t0, Tend = 0.0, self.Tend P = controller.MS[0].levels[0].prob - uinit =P.u_init() + uinit = P.u_init() uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) func_eval = P.work_counters['Boris_solver'].niter + P.work_counters['rhs'].niter diff --git a/pySDC/projects/Second_orderSDC/penningtrap_run_Hamiltonian_error.py b/pySDC/projects/Second_orderSDC/penningtrap_run_Hamiltonian_error.py index 254a55626b..cab9299836 100644 --- a/pySDC/projects/Second_orderSDC/penningtrap_run_Hamiltonian_error.py +++ b/pySDC/projects/Second_orderSDC/penningtrap_run_Hamiltonian_error.py @@ -99,7 +99,7 @@ def plot_Hamiltonian_error(K, M, dt): # pragma: no cover """ set_fixed_plot_params() # Define final time - time = 1e+6 + time = 1e6 tn = dt # Find time nodes t = np.arange(0, time + tn, tn) diff --git a/pySDC/projects/Second_orderSDC/stability_simulation.py b/pySDC/projects/Second_orderSDC/stability_simulation.py index 804f1f334c..3f35191405 100644 --- a/pySDC/projects/Second_orderSDC/stability_simulation.py +++ b/pySDC/projects/Second_orderSDC/stability_simulation.py @@ -5,6 +5,8 @@ from pySDC.projects.Second_orderSDC.plot_helper import set_fixed_plot_params from tabulate import tabulate + + class StabilityImplementation: """ Routine to compute the stability domains of different configurations of SDC @@ -209,11 +211,17 @@ def run_RKN_stability(self): # pragma: no cover self.plot_stability(region_RKN.T, title='RKN-4 stability region') -def compute_and_check_stability(description, - point, compute_interval=False, save_interval_file=False, interval_filename='./data/stab_interval.txt', - check_stability=False, store_results=False, results_filename='./data/stab_results.txt', quadrature_list=['GAUSS', 'LOBATTO'] +def compute_and_check_stability( + description, + point, + compute_interval=False, + save_interval_file=False, + interval_filename='./data/stab_interval.txt', + check_stability=False, + store_results=False, + results_filename='./data/stab_results.txt', + quadrature_list=['GAUSS', 'LOBATTO'], ): - # Storage for stability interval and stability check interval_data = [] results_data = [] @@ -229,14 +237,13 @@ def compute_and_check_stability(description, # Create Stability_implementation instance for stability check - stab_model = StabilityImplementation(description, kappa_max=point[0], mu_max=point[1], Num_iter=description['helper_params']['Num_iter']) + stab_model = StabilityImplementation( + description, kappa_max=point[0], mu_max=point[1], Num_iter=description['helper_params']['Num_iter'] + ) if compute_interval: - - # Extract the values where SDC is less than or equal to 1 mask = stab_model.picard <= 1 + 1e-14 for ii in range(len(mask)): - if mask[ii]: kappa_max_interval = stab_model.lambda_kappa[ii] else: @@ -246,8 +253,6 @@ def compute_and_check_stability(description, interval_data.append([quad_type, num_nodes, max_iter, kappa_max_interval]) if check_stability: - - # Check stability and print results if stab_model.SDC[-1, -1] <= 1: stability_result = "Stable" @@ -255,9 +260,9 @@ def compute_and_check_stability(description, stability_result = "Unstable. Increase M or K" # Add row to the results data - results_data.append([ - quad_type, num_nodes, max_iter, point, stability_result, stab_model.SDC[-1, -1] - ]) + results_data.append( + [quad_type, num_nodes, max_iter, point, stability_result, stab_model.SDC[-1, -1]] + ) # Define column names for interval data interval_headers = ["Quad Type", "Num Nodes", "Max Iter", 'kappa_max']