Bug fix (#396)

* Clear documentation and some edits in the code

* forgot black

* some changes

pySDC/projects/Second_orderSDC/README.rst

@@ -24,11 +24,12 @@ Reproducing Figures from the Publication
 - **Fig. 2:** Run `harmonic_oscillator_run_stability.py` with the following configurations:
    - Set `kappa_max=30` and `mu_max=30`.
    - Adjust `maxiter` to 1, 2, 3, or 4 and execute each individually.
-- **Table 1:** Execute `harmonic_socillator_run_stab_interval.py`:
-   - To save the results set: `save_file=True`
+- **Table 1:** Execute `harmonic_oscillator_run_stab_interval.py`:
+   - For the Picard iteration set: `Picard=True`
+   - To save the results set: `save_interval_file=True`
 
 - Use the script `harmonic_oscillator_run_points.py` to create a table based on given $(\kappa, \mu)$ points. This table assists in determining suitable values for `M`, `K`, and `quadrature nodes` to ensure stability in the SDC method.
-   - To save the results set: `save_file=True`
+   - To save the results set: `save_points_file=True`
 
 - **Fig. 3:** Run `penningtrap_run_error.py` (Run local convergence: `conv.run_local_error()`) with `dt=0.015625/4` and `axes=(0,)`.
 - **Fig. 4:** Run `penningtrap_run_error.py` (Run local convergence: `conv.run_local_error()`) using `dt=0.015625*4` and `axes=(2,)`.

pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py

@@ -2,26 +2,30 @@
 from pySDC.projects.Second_orderSDC.harmonic_oscillator_params import get_default_harmonic_oscillator_description
 from pySDC.projects.Second_orderSDC.stability_simulation import compute_and_generate_table
 
-
 if __name__ == '__main__':
     '''
-    This script generates table for the given points to compare in what quadrature type,
-    number of nodes and number of iterations the SDC iteration is stable or not.
-    Additional parameter:
-        To save the table set: save_points_table=True
-        Change filename: points_table_filename='FILENAME' by default: './data/point_table.txt'
+    This script generates a table to compare stability of the SDC iteration for given points,
+    exploring different quadrature types, number of nodes, and number of iterations.
+
+    Additional parameters in the function `compute_and_generate_table`:
+        - To save the table, set: `save_points_table=True`.
+        - To change the filename, set `points_table_filename='FILENAME'`. Default is './data/point_table.txt'.
     '''
-    # Get default parameters for the harmonic osicillator problem
+    # This code checks if the "data" folder exists or not.
+    exec(open("check_data_folder.py").read())
+    # Get default parameters for the harmonic oscillator problem
     description = get_default_harmonic_oscillator_description()
-    # Additonal params to compute stability points
+
+    # Additional parameters to compute stability points
     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),
+        'quad_type_list': ('GAUSS', 'LOBATTO'),  # List of quadrature types
+        'Num_iter': (2, 2),  # Number of iterations
+        'num_nodes_list': np.arange(3, 6, 1),  # List of number of nodes
+        'max_iter_list': np.arange(2, 10, 1),  # List of maximum iterations
     }
 
-    points = ((1, 100), (3, 100), (10, 100))
+    points = ((1, 100), (3, 100), (10, 100))  # Stability parameters: (kappa, mu)
+
     # Iterate through points and perform stability check
     for ii in points:
-        compute_and_generate_table(description, helper_params, ii, check_stability_point=True)
+        compute_and_generate_table(description, helper_params, ii, check_stability_point=True, save_points_table=False)

pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stab_interval.py

@@ -2,28 +2,40 @@
 from pySDC.projects.Second_orderSDC.harmonic_oscillator_params import get_default_harmonic_oscillator_description
 from pySDC.projects.Second_orderSDC.stability_simulation import compute_and_generate_table
 
-
 if __name__ == '__main__':
     '''
-    To compute maximum stable values for SDC and Picard iterations for the purely oscillatory case with
-    no damping (mu=0)
-    Additional parameter:
-        To save the table set: save_interval_file=True
-        Change filename: interval_filename='FILENAME' by default: './data/stab_interval.txt'
+    Main script to compute maximum stable values for SDC and Picard iterations for the purely oscillatory case with no damping (mu=0).
+
+    Additional parameters in the `compute_and_generate_table` function:
+        - To save the stability table, set `save_interval_file=True`.
+        - To change the filename, set `interval_filename='FILENAME'`. Default is './data/stab_interval.txt'.
+
     Output:
-        it generates to compare with different values of M (number of nodes) and K (maximal number of iterations)
+        The script generates data to compare different values of M (number of nodes) and K (maximal number of iterations).
     '''
-    # Ger default parameters
+    # This code checks if the "data" folder exists or not.
+    exec(open("check_data_folder.py").read())
+    # Get default parameters for the harmonic oscillator
     description = get_default_harmonic_oscillator_description()
+
     # Additional parameters to compute stability interval on the kappa
+    # =============================================================================
+    #     To get exactly the same as table in the paper set:
+    #       'Num_iter': (500, 1) for the SDC iteration
+    #       'Num_iter': (2000, 1) for the Picard iteration
+    # =============================================================================
     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),
+        'quad_type_list': ('GAUSS',),  # Type of quadrature
+        'Num_iter': (500, 1),  # Number of iterations
+        'num_nodes_list': np.arange(2, 7, 1),  # List of number of nodes
+        'max_iter_list': np.arange(1, 11, 1),  # List of maximum iterations
     }
 
-    points = ((100, 1e-10),)
+    points = ((100, 1e-10),)  # Stability parameters: (Num_iter, mu)
+
     # Iterate through points and perform stability check
     for ii in points:
-        compute_and_generate_table(description, helper_params, ii, compute_interval=True)
+        # If you want to get the table for the Picard iteration set Picard=True
+        compute_and_generate_table(
+            description, helper_params, ii, compute_interval=True, Picard=False, save_interval_file=True
+        )

pySDC/projects/Second_orderSDC/penningtrap_run_error.py

@@ -21,11 +21,11 @@
     controller_params, description = penningtrap_params()
     ## =============================================================================
     ##     dt-timestep and num_nodes can be changed here manually
-    description['level_params']['dt'] = 0.015625 * 4
+    description['level_params']['dt'] = 0.015625 / 4
     description['sweeper_params']['num_nodes'] = 4
     ## =============================================================================
     # Give the parameters to the class
-    conv = ComputeError(controller_params, description, time_iter=3, K_iter=(1, 2, 3, 10), axes=(2,))
+    conv = ComputeError(controller_params, description, time_iter=3, K_iter=(1, 2, 3), axes=(0,))
     # Run local convergence order
     # conv.run_local_error()
     # Run global convergence order

pySDC/projects/Second_orderSDC/penningtrap_run_work_precision.py

@@ -21,8 +21,8 @@
     ## =============================================================================
     ##     dt-timestep and Tend can be changed here manually
     Tend = 128 * 0.015625
-    description['level_params']['dt'] = 0.015625 * 4
-    description['sweeper_params']['initial_guess'] = 'spread'  # 'zero', 'spread'
+    description['level_params']['dt'] = 0.015625 * 2
+    description['sweeper_params']['initial_guess'] = 'spread'  #'random' 'zero', 'spread'
     ## =============================================================================
     work_pre = ComputeError(controller_params, description, time_iter=3, Tend=Tend, K_iter=(1, 2, 3), axes=(2,))
     work_pre.run_work_precision(RK=True)

pySDC/projects/Second_orderSDC/stability_simulation.py

@@ -228,7 +228,9 @@ def run_RKN_stability(self):  # pragma: no cover
         self.plot_stability(region_RKN.T, title='RKN-4 stability region')
 
 
-def check_points_and_interval(description, helper_params, point, compute_interval=False, check_stability_point=False):
+def check_points_and_interval(
+    description, helper_params, point, compute_interval=False, check_stability_point=False, Picard=False
+):
     # Storage for stability interval and stability check
     interval_data = []
     points_data = []
@@ -249,7 +251,10 @@ def check_points_and_interval(description, helper_params, point, compute_interva
                 )
                 if compute_interval:
                     # Extract the values where SDC is less than or equal to 1
-                    mask = stab_model.picard <= 1 + 1e-14
+                    if Picard:
+                        mask = stab_model.picard <= 1 + 1e-14
+                    else:
+                        mask = stab_model.SDC <= 1.0
                     for ii in range(len(mask)):
                         if mask[ii]:
                             kappa_max_interval = stab_model.lambda_kappa[ii]
@@ -287,11 +292,14 @@ def compute_and_generate_table(
     save_points_table=False,
     points_table_filename='./data/point_table.txt',
     quadrature_list=('GAUSS', 'LOBATTO'),
+    Picard=False,
 ):  # pragma: no cover
     from tabulate import tabulate
 
     if compute_interval:
-        interval_data = check_points_and_interval(description, helper_params, point, compute_interval=compute_interval)
+        interval_data = check_points_and_interval(
+            description, helper_params, point, compute_interval=compute_interval, Picard=Picard
+        )
     else:
         points_data = check_points_and_interval(
             description, helper_params, point, check_stability_point=check_stability_point
@@ -316,7 +324,10 @@ def compute_and_generate_table(
         print(f"Stability Results Table saved to {points_table_filename}")
 
     if compute_interval:
-        print("Stability Interval Table:")
+        if Picard:
+            print("Picard stability Interval Table:")
+        else:
+            print("SDC stability Interval Table:")
         print(tabulate(interval_data, headers=interval_headers, tablefmt="grid"))
 
     if check_stability_point: