Review

pySDC/implementations/problem_classes/HarmonicOscillator.py

@@ -28,6 +28,7 @@ class harmonic_oscillator(ptype):
         Phase of the oscillation.
     amp : float, optional
         Amplitude of the oscillation.
+    Source: https://beltoforion.de/en/harmonic_oscillator/
     """
 
     dtype_u = particles

pySDC/projects/Second_orderSDC/README.md

@@ -1,29 +1,52 @@
 # Spectral Deferred Correction Methods for Second-Order Problems
 
+**Title:** Spectral Deferred Correction Methods for Second-order Problems
+
+**Authors:** Ikrom Akramov, Sebastian Götschel, Michael Minion, Daniel Ruprecht, and Robert Speck.
+
+
 Python code for implementing the paper's plots on Second-order SDC methods.
 
 ## Attribution
 You are welcome to use and adapt this code under the terms of the BSD license.
 If you utilize it, either in whole or in part, for a publication, please provide proper citation:
 
-**Title:** Spectral Deferred Correction Methods for Second-order Problems
-
-**Authors:** Ikrom Akramov, Sebastian Götschel, Michael Minion, Daniel Ruprecht, and Robert Speck.
 
-[![DOI](http://example.com)](http://example.com)
+
+   **BibTeX formatted citation:**
+
+          @misc{akramov2023spectral,
+                title={Spectral deferred correction methods for second-order problems}, 
+               author={Ikrom Akramov and Sebastian Götschel and Michael Minion and Daniel Ruprecht and Robert Speck},
+               year={2023},
+               eprint={2310.08352},
+               archivePrefix={arXiv},
+               primaryClass={math.NA}
+
+
+
+
+
+[![arXiv](https://img.shields.io/badge/arXiv-2310.08352-b31b1b.svg)](https://arxiv.org/abs/2310.08352)
 
 ## Reproducing Figures from the Publication
 
-- **Fig. 1:** Execute `dampedharmonic_oscillator_run_stability.py` while setting `kappa_max=18` and `mu_max=18`.
-- **Fig. 2:** Run `dampedharmonic_oscillator_run_stability.py` with the following configurations:
+- **Fig. 1:** Execute `harmonic_oscillator_run_stability.py` while setting `kappa_max=18` and `mu_max=18`.
+- **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`
+
+- 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`
+
 - **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,)`.
 - **Fig. 5:** Run `penningtrap_run_error.py` (Run global convergence: `conv.run_global_error()`) with `dt=0.015625*2`:
    - Note: Perform each run individually: first with `axes=(0,)`, then with `axes=(2,)`.
    - Manually set y-axis limits in `penningtrap_run_error.py`, specifically in lines 147-148.
-- **Table 1:** Execute `penningtrap_run_error.py` (Run global convergence: `conv.run_global_error()`) with the following settings:
+- **Table 2:** Execute `penningtrap_run_error.py` (Run global convergence: `conv.run_global_error()`) with the following settings:
    - Expected order and approximate order are saved in the file: `data/global_order_vs_approx_order.csv`
    - Set: `K_iter=(2, 3, 4, 10)`
    - For `M=2`:

pySDC/projects/Second_orderSDC/check_data_folder.py

@@ -3,7 +3,7 @@
 folder_name = "./data"
 
 # Check if the folder already exists
-if not os.path.isdir(folder_name):  # pragma: no cover
+if not os.path.isdir(folder_name):
     # Create the folder
     os.makedirs(folder_name)
 else:

pySDC/projects/Second_orderSDC/harmonic_oscillator_params.py

@@ -0,0 +1,36 @@
+import numpy as np
+from pySDC.implementations.problem_classes.HarmonicOscillator import harmonic_oscillator
+from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order
+
+
+def get_default_harmonic_oscillator_description():
+    """
+    Routine to compute modules of the stability function
+
+    Returns:
+        description (dict): A dictionary containing parameters for the damped harmonic oscillator problem
+    """
+
+    # Initialize level parameters
+    level_params = {'restol': 1e-16, 'dt': 1.0}
+
+    # Initialize problem parameters for the Damped harmonic oscillator problem
+    problem_params = {'k': 0, 'mu': 0, 'u0': np.array([1, 1])}
+
+    # Initialize sweeper parameters
+    sweeper_params = {'quad_type': 'GAUSS', 'num_nodes': 3, 'do_coll_update': True, 'picard_mats_sweep': True}
+
+    # Initialize step parameters
+    step_params = {'maxiter': 5}
+
+    # Fill description dictionary for easy step instantiation
+    description = {
+        'problem_class': harmonic_oscillator,
+        'problem_params': problem_params,
+        'sweeper_class': boris_2nd_order,
+        'sweeper_params': sweeper_params,
+        'level_params': level_params,
+        'step_params': step_params,
+    }
+
+    return description

pySDC/projects/Second_orderSDC/harmonic_oscillator_run_points.py

@@ -0,0 +1,27 @@
+import numpy as np
+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'
+    '''
+    # Get default parameters for the harmonic osicillator problem
+    description = get_default_harmonic_oscillator_description()
+    # Additonal params 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),
+    }
+
+    points = ((1, 100), (3, 100), (10, 100))
+    # Iterate through points and perform stability check
+    for ii in points:
+        compute_and_generate_table(description, helper_params, ii, check_stability_point=True)

pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stab_interval.py

@@ -0,0 +1,29 @@
+import numpy as np
+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'
+    Output:
+        it generates to compare with different values of M (number of nodes) and K (maximal number of iterations)
+    '''
+    # Ger default parameters
+    description = get_default_harmonic_oscillator_description()
+    # Additional parameters to compute stability interval on the kappa
+    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),
+    }
+
+    points = ((100, 1e-10),)
+    # Iterate through points and perform stability check
+    for ii in points:
+        compute_and_generate_table(description, helper_params, ii, compute_interval=True)

pySDC/projects/Second_orderSDC/harmonic_oscillator_run_stability.py

@@ -0,0 +1,26 @@
+from pySDC.projects.Second_orderSDC.harmonic_oscillator_params import get_default_harmonic_oscillator_description
+from pySDC.projects.Second_orderSDC.stability_simulation import StabilityImplementation
+
+
+if __name__ == '__main__':
+    """
+    To implement Stability region for the Harmonic Oscillator problem
+    Run for
+        SDC stability region: model_stab.run_SDC_stability()
+        Picard stability region: model_stab.run_Picard_stability()
+        Runge-Kutta-Nzström stability region: model_run_RKN_stability()
+    To implement spectral radius of iteration matrix
+    Run:
+        Iteration matrix of SDC method: model_stab.run_Ksdc()
+        Iteration matrix of Picard method: model_stab.run_Kpicard()
+
+    """
+    # Execute the stability analysis for the damped harmonic oscillator
+    description = get_default_harmonic_oscillator_description()
+    model_stab = StabilityImplementation(description, kappa_max=30, mu_max=30, Num_iter=(200, 200))
+
+    model_stab.run_SDC_stability()
+    model_stab.run_Picard_stability()
+    model_stab.run_RKN_stability()
+    model_stab.run_Ksdc()
+    # model_stab.run_Kpicard()