Skip to content

Commit

Permalink
forgot black
Browse files Browse the repository at this point in the history
  • Loading branch information
ikrom96git committed Jan 11, 2024
1 parent 4bf0008 commit 3d04ad1
Show file tree
Hide file tree
Showing 6 changed files with 43 additions and 32 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
16 changes: 9 additions & 7 deletions pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Original file line number Diff line number Diff line change
Expand Up @@ -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)
9 changes: 6 additions & 3 deletions pySDC/projects/Second_orderSDC/penningtrap_Simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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):
Expand All @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
31 changes: 18 additions & 13 deletions pySDC/projects/Second_orderSDC/stability_simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 = []
Expand All @@ -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:
Expand All @@ -246,18 +253,16 @@ 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"
else:
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']
Expand Down

0 comments on commit 3d04ad1

Please sign in to comment.