Second order sdc (#345)

* citation * Hamiltonia * Cleaned codes * much more clean and more documentation * More Black * Code coverage * Hopefully, it passes Codecov * Totally removed Hamiltonian from coverage report --------- Co-authored-by: Ikrom Akramov <[email protected]>
Parallel-in-Time · Aug 15, 2023 · efe2999 · efe2999
1 parent 196ae36
commit efe2999
Show file tree

Hide file tree

Showing 10 changed files with 727 additions and 609 deletions.
diff --git a/pySDC/implementations/problem_classes/PenningTrap_3D.py b/pySDC/implementations/problem_classes/PenningTrap_3D.py
@@ -8,8 +8,74 @@
 
 # noinspection PyUnusedLocal
 class penningtrap(ptype):
-    """
-    Example implementing particles in a penning trap
+    r"""
+    This class implements a standard Penning trap problem on the time interval :math:`[0, t_{end}]`
+    fully investigated in [1]_. The equations are given by the following equation of motion
+
+    .. math::
+        \frac{dv}{dt}=f(x,v)=\alpha[E(x,t)+v\times B(x,t)],
+
+    .. math::
+        \frac{dx}{dt}=v.
+    with the particles :math:`x, v\in \mathbb{R}^{3}`. For the penning trap problem, the other parameters are given by
+    The contant magnetic field :math:`B=\frac{\omega_{B}}{\alpha}\cdot \hat{e_{z}}\in \mathbb{R}^{3}`
+    along the :math:`z-` axis with the particle's charge-to-mass ratio :math:`\alpha=\frac{q}{m}` so that
+
+    .. math::
+        v\times B=\frac{\omega_{B}}{\alpha}\left(
+            \begin{matrix}
+            0 & 1 & 0\\
+                -1 & 0 & 0\\
+                    0 & 0 & 0
+            \end{matrix}
+            \right)v.
+
+    The electric field :math:`E(x_{i})=E_{ext}(x_{i})+E_{int}(x_{i})\in \mathbb{R}^{3}` where
+
+    .. math::
+        E_{ext}(x_{i})=-\epsilon\frac{\omega_{E}^{2}}{\alpha}\left(
+            \begin{matrix}
+            1 & 0 & 0\\
+                0 & 1 & 0\\
+                    0 & 0 & -2
+            \end{matrix}
+            \right)x_{i}.
+
+    and the inter-particle Coulomb interaction
+
+    .. math::
+        E_{int}(x_{i})=\sum_{k=1, k\neq i}^{N_{particles}}Q_{k}\frac{x_{i}-x_{k}}{(|x_{i}-x_{k}|^{2}+\lambda^{2})^{3/2}}
+
+    with the smoothing parameter :math:`\lambda>0`.
+    The exact solution also given for the single particle penning trap more detailed [1]_ [2]_.
+    For to solve nonlinear eqaution of system, Boris trick is used (see. [2]_)
+
+    Parameters
+
+    ----------
+    omega_B : float
+        Amplitute of magnetic field.
+    omega_E : float
+        Amplitute of electric field.
+    u0 : array
+        initial condition for postion
+
+        initial condition for velocity
+
+        q : value
+
+        m : mass
+    nparts : int
+        The number of particles.
+    sig : float
+        The smoothing parameter :math:`\lambda>0`.
+
+    References
+    ----------
+    .. [1] F. Penning. Die Glimmentladung bei niedrigem Druck zwischen koaxialen Zylindern in einem axialen Magnetfeld.
+        Physica. Vol. 3 (1936).
+    .. [2] Mathias Winkel, Robert Speck and Daniel Ruprecht. A high-order Boris integrator.
+        Journal of Computational Physics (2015).
     """
 
     dtype_u = particles
@@ -76,8 +142,12 @@ def eval_f(self, part, t):
         N = self.nparts
 
         self.work_counters['rhs']()
+        try:
+            penningtrap.Harmonic_oscillator
+            Emat = np.diag([0, 0, -1])
+        except AttributeError:
+            Emat = np.diag([1, 1, -2])
 
-        Emat = np.diag([1, 1, -2])
         f = self.dtype_f(self.init)
 
         f.elec[:] = self.get_interactions(part)

diff --git a/pySDC/projects/Second_orderSDC/README.md b/pySDC/projects/Second_orderSDC/README.md
@@ -1,28 +1,44 @@
-# Spectral Deferred Correction methods for Second-order problems
+# Spectral Deferred Correction Methods for Second-Order Problems
 
-Python code to implement Second-order SDC paper plots.
+Python code for implementing the paper's plots on Second-order SDC methods.
 
 ## Attribution
-You can freely use and reuse this code in line with the BSD license. 
-If you use it (or parts of it) for a publication, please cite:
+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:
 
-- [Publication citation information will be added soon]
+**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)
 
-## How can I reproduce Figures from the publication?
-
-- Fig. 1: Run `dampedharmonic_oscillator_run_stability.py`.
-- Fig. 2: Run `dampedharmonic_oscillator_run_stability.py` with the following settings:
-   - Set `maxiter` to 1, 2, 3, or 4 and run individually.
-- Fig. 3: Run `penningtrap_run_error.py` (Run local convergence) with `dt=0.015625/4` and `axes=[0]`.
-- Fig. 4: Run `penningtrap_run_error.py` (Run local convergence) with `dt=0.015625*4` and `axes=[2]`.
-- Fig. 5: Run `penningtrap_run_error.py` (Run global convergence) with `dt=0.015625*2` and `axes=[0]` and `axes=[2]`.
-   - Note: Each run needs to be performed individually.
-   - The y-axis limits should be manually set in `penningtrap_run_error.py`, specifically in lines 103-104.
-- Fig. 6: Run `penningtrap_run_work_precision.py` with the following settings:
-   - Set `dt=0.015625*2` and `K_iter=[1, 2, 3]` in `axis=[2]`.
-   - Set `K_iter=[2, 4, 6]` in `axis=[0]`.
+## 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:
+   - Set `kappa_max=30` and `mu_max=30`.
+   - Adjust `maxiter` to 1, 2, 3, or 4 and execute each individually.
+- **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:
+   - 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`:
+      - Set `dt=0.015625 / 2` and `num_nodes=2` (Both Horizontal and Vertical axes)
+   - For `M=3`:
+      - Set `dt=0.015625` and `num_nodes=3` (Both Horizontal and Vertical axes)
+   - For `M=4`:
+      - Set `dt=0.015625 * 4` and `num_nodes=4` (Both Horizontal and Vertical axes)
+- **Fig. 6:** Run `penningtrap_run_Hamiltonian_error.py`:
+   - Note: Execute each computation individually in case of crashes.
+   - Data is saved in the data folder.
+- **Fig. 7:** Execute `penningtrap_run_work_precision.py` with the following configurations:
+   - Set `dt=0.015625*2` and `K_iter=(1, 2, 3)` for `axis=(2,)`.
+   - Set `dt=0.015625*4`, and `K_iter=(2, 4, 6)`, and `dt_cont=2` for `axis=(0,)`.
+
 
 ## Who do I talk to?
 

diff --git a/pySDC/projects/Second_orderSDC/dampedharmonic_oscillator_run_stability.py b/pySDC/projects/Second_orderSDC/dampedharmonic_oscillator_run_stability.py
@@ -9,19 +9,8 @@ def dampedharmonic_oscillator_params():
     """
     Runtine to compute modulues of the stability function
 
-    Args:
-        None
-
-
     Returns:
-        numpy.narray: values for the spring pendulum
-        numpy.narray: values for the Friction
-        numpy.narray: number of num_nodes
-        numpy.narray: number of iterations
-        numpy.narray: moduli for the SDC
-        numpy.narray: moduli for the K_{sdc} marrix
-        numpy.narray: moduli for the Pircard iteration
-        numpy.narray: moduli for the K_{sdc} Picard iteration
+        description
     """
 
     # initialize level parameters
@@ -44,7 +33,7 @@ def dampedharmonic_oscillator_params():
 
     # initialize step parameters
     step_params = dict()
-    step_params['maxiter'] = 100
+    step_params['maxiter'] = 50
 
     # fill description dictionary for easy step instantiation
     description = dict()
@@ -71,8 +60,8 @@ def dampedharmonic_oscillator_params():
     # exec(open("check_data_folder.py").read())
     description = dampedharmonic_oscillator_params()
     Stability = Stability_implementation(description, kappa_max=18, mu_max=18, Num_iter=(200, 200))
-    Stability.run_SDC_stability
-    Stability.run_Picard_stability
-    Stability.run_RKN_stability
-    Stability.run_Ksdc
+    Stability.run_SDC_stability()
+    Stability.run_Picard_stability()
+    Stability.run_RKN_stability()
+    Stability.run_Ksdc()
     # Stability.run_Kpicard
diff --git a/pySDC/projects/Second_orderSDC/penningtrap_HookClass.py b/pySDC/projects/Second_orderSDC/penningtrap_HookClass.py
@@ -1,6 +1,5 @@
-import numpy as np
-
 from pySDC.core.Hooks import hooks
+import numpy as np
 
 
 class particles_output(hooks):
@@ -19,72 +18,15 @@ def pre_run(self, step, level_number):
         """
         super(particles_output, self).pre_run(step, level_number)
 
-        # some abbreviations
-        L = step.levels[level_number]
-
-        part = L.u[0]
-        N = L.prob.nparts
-        w = np.array([1, 1, -2])
-
-        # compute (slowly..) the potential at u0
-        fpot = np.zeros(N)
-        for i in range(N):
-            # inner loop, omit ith particle
-            for j in range(0, i):
-                dist2 = np.linalg.norm(part.pos[:, i] - part.pos[:, j], 2) ** 2 + L.prob.sig**2
-                fpot[i] += part.q[j] / np.sqrt(dist2)
-            for j in range(i + 1, N):
-                dist2 = np.linalg.norm(part.pos[:, i] - part.pos[:, j], 2) ** 2 + L.prob.sig**2
-                fpot[i] += part.q[j] / np.sqrt(dist2)
-            fpot[i] -= L.prob.omega_E**2 * part.m[i] / part.q[i] / 2.0 * np.dot(w, part.pos[:, i] * part.pos[:, i])
-
-        # add up kinetic and potntial contributions to total energy
-        epot = 0
-        ekin = 0
-        for n in range(N):
-            epot += part.q[n] * fpot[n]
-            ekin += part.m[n] / 2.0 * np.dot(part.vel[:, n], part.vel[:, n])
-
-        self.add_to_stats(
-            process=step.status.slot,
-            time=L.time,
-            level=L.level_index,
-            iter=0,
-            sweep=L.status.sweep,
-            type="etot",
-            value=epot + ekin,
-        )
-
-
-class convergence_data(hooks):
-    def __init__(self):
-        super(convergence_data, self).__init__()
-
-        self.storage = dict()
-
-        self.values = [
-            "position",
-            "velocity",
-            "position_exact",
-            "velocity_exact",
-            "pos_nodes",
-            "vel_nodes",
-            "pos_nodes_ex",
-            "vel_nodes_ex",
-        ]
-
-        for _, jj in enumerate(self.values):
-            self.storage[jj] = dict()
-
     def post_step(self, step, level_number):
         """
-        Default runtine called after each iteration
+        Default routine called after each iteration
         Args:
             step: the current step
             level_number: the current level number
         """
 
-        super(convergence_data, self).pre_run(step, level_number)
+        super(particles_output, self).post_step(step, level_number)
 
         # some abbreviations
         L = step.levels[level_number]
@@ -93,21 +35,77 @@ def post_step(self, step, level_number):
 
         L.sweep.compute_end_point()
         part = L.uend
+        N = L.prob.nparts
 
-        self.storage["position"][L.time] = part.pos
-        self.storage["velocity"][L.time] = part.vel
-        self.storage["position_exact"][L.time] = L.prob.u_exact(L.time + L.dt).pos
-        self.storage["velocity_exact"][L.time] = L.prob.u_exact(L.time + L.dt).vel
-
-        if L.time + L.dt >= self.Tend:
-            self.add_to_stats(
-                process=step.status.slot,
-                time=L.dt,
-                level=L.level_index,
-                iter=step.status.iter,
-                sweep=L.status.sweep,
-                type="error",
-                value=self.storage,
-            )
-
-        return None
+        # =============================================================================
+        #
+
+        try:  # pragma: no cover
+            L.prob.Harmonic_oscillator
+
+            # add up kinetic and potntial contributions to total energy
+            epot = 0
+            ekin = 0
+            name = str(L.sweep)
+            for n in range(N):
+                epot += 1 / 2.0 * np.dot(part.pos[:, n], part.pos[:, n])
+                ekin += 1 / 2.0 * np.dot(part.vel[:, n], part.vel[:, n])
+                Energy = epot + ekin
+            uinit = L.u[0]
+            H0 = 1 / 2 * (np.dot(uinit.vel[:].T, uinit.vel[:]) + np.dot(uinit.pos[:].T, uinit.pos[:]))
+            Ham = abs(Energy - H0) / abs(H0)
+            if 'RKN' in name:
+                filename = 'data/Ham_RKN' + '{}.txt'.format(step.params.maxiter)
+            else:
+                filename = 'data/Ham_SDC' + '{}{}.txt'.format(L.sweep.coll.num_nodes, step.params.maxiter)
+
+            if L.time == 0.0:
+                file = open(filename, 'w')
+                file.write(str('time') + " | " + str('Hamiltonian error') + '\n')
+            else:
+                file = open(filename, 'a')
+
+            file.write(str(L.time) + " | " + str(Ham) + '\n')
+
+            file.close()
+        except AttributeError:
+            pass
+        # =============================================================================
+
+        part_exact = L.prob.u_exact(L.time + L.dt)
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='position',
+            value=part.pos,
+        )
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='velocity',
+            value=part.vel,
+        )
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='position_exact',
+            value=part_exact.pos,
+        )
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='velocity_exact',
+            value=part_exact.vel,
+        )