TL: first revision for the parallel SDC theory paper

CITATION.cff

@@ -30,9 +30,9 @@ authors:
     affiliation: "Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH, 52425 Jülich, Germany"
 
 
-version: 5.5.0
+version: 5.5.1
 doi: 10.5281/zenodo.594191
-date-released: 2024-06-24
+date-released: 2024-09-03
 keywords: 
   - "parallel-in-time"
   - "spectral deferred corrections"

docs/source/conf.py

@@ -72,7 +72,7 @@
 # The short X.Y version.
 version = '5.5'
 # The full version, including alpha/beta/rc tags.
-release = '5.5.0'
+release = '5.5.1'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

pySDC/core/sweeper.py

@@ -263,3 +263,15 @@ def level(self, L):
     @property
     def rank(self):
         return 0
+
+    def updateVariableCoeffs(self, k):
+        """
+        Eventually update QDelta implicit coefficients if variable ...
+
+        Parameters
+        ----------
+        k : int
+            Index of the sweep (0 for initial sweep, 1 for the first one, ...).
+        """
+        if self.params.QI == 'MIN-SR-FLEX':
+            self.QI = self.get_Qdelta_implicit(qd_type="MIN-SR-FLEX", k=k)

pySDC/implementations/sweeper_classes/generic_implicit.py

@@ -66,8 +66,7 @@ def update_nodes(self):
         M = self.coll.num_nodes
 
         # update the MIN-SR-FLEX preconditioner
-        if self.params.QI == 'MIN-SR-FLEX':
-            self.QI = self.get_Qdelta_implicit(qd_type="MIN-SR-FLEX", k=L.status.sweep)
+        self.updateVariableCoeffs(L.status.sweep)
 
         # gather all terms which are known already (e.g. from the previous iteration)
         # this corresponds to u0 + QF(u^k) - QdF(u^k) + tau

pySDC/implementations/sweeper_classes/generic_implicit_MPI.py

@@ -201,7 +201,8 @@ def update_nodes(self):
         # only if the level has been touched before
         assert L.status.unlocked
 
-        # get number of collocation nodes for easier access
+        # update the MIN-SR-FLEX preconditioner
+        self.updateVariableCoeffs(L.status.sweep)
 
         # gather all terms which are known already (e.g. from the previous iteration)
         # this corresponds to u0 + QF(u^k) - QdF(u^k) + tau

pySDC/projects/parallelSDC_reloaded/environment.yml

@@ -9,6 +9,8 @@ dependencies:
   - matplotlib>=3.0
   - dill>=0.2.6
   - scipy>=0.17.1
+  - mpich
+  - mpi4py>=3.0.0
   - pip
   - pip:
     - qmat>=0.1.8

pySDC/projects/parallelSDC_reloaded/scripts/fig06_allenCahnMPI.py

@@ -0,0 +1,105 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 11 11:14:01 2024
+
+Figures with experiments on the Allen-Cahn problem (MPI runs)
+"""
+import os
+import sys
+import json
+import numpy as np
+from mpi4py import MPI
+
+from pySDC.projects.parallelSDC_reloaded.utils import solutionExact, getParamsSDC, solutionSDC, getParamsRK
+from pySDC.implementations.sweeper_classes.generic_implicit_MPI import generic_implicit_MPI
+
+PATH = '/' + os.path.join(*__file__.split('/')[:-1])
+SCRIPT = __file__.split('/')[-1].split('.')[0]
+
+COMM_WORLD = MPI.COMM_WORLD
+
+# SDC parameters
+nNodes = 4
+quadType = 'RADAU-RIGHT'
+nodeType = 'LEGENDRE'
+parEfficiency = 0.8  # 1/nNodes
+nSweeps = 4
+
+# Problem parameters
+pName = "ALLEN-CAHN"
+tEnd = 50
+pParams = {
+    "periodic": False,
+    "nvars": 2**11 - 1,
+    "epsilon": 0.04,
+}
+
+# -----------------------------------------------------------------------------
+# %% Convergence and error VS cost plots
+# -----------------------------------------------------------------------------
+nStepsList = np.array([5, 10, 20, 50, 100, 200, 500])
+dtVals = tEnd / nStepsList
+
+
+def getError(uNum, uRef):
+    if uNum is None:
+        return np.inf
+    return np.linalg.norm(uRef[-1, :] - uNum[-1, :], ord=2)
+
+
+def getCost(counters):
+    _, _, tComp = counters
+    return tComp
+
+
+try:
+    qDelta = sys.argv[1]
+except IndexError:
+    qDelta = "MIN-SR-FLEX"
+
+try:
+    params = getParamsRK(qDelta)
+except KeyError:
+    params = getParamsSDC(quadType=quadType, numNodes=nNodes, nodeType=nodeType, qDeltaI=qDelta, nSweeps=nSweeps)
+
+useMPI = False
+if COMM_WORLD.Get_size() == 4 and qDelta in ["MIN-SR-NS", "MIN-SR-S", "MIN-SR-FLEX", "VDHS"]:  # pragma: no cover
+    params['sweeper_class'] = generic_implicit_MPI
+    useMPI = True
+
+errors = []
+costs = []
+
+root = COMM_WORLD.Get_rank() == 0
+if root:
+    print(f"Running simulation with {qDelta}")
+
+for nSteps in nStepsList:
+    if root:
+        uRef = solutionExact(tEnd, nSteps, pName, **pParams)
+
+    uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, pName, verbose=root, noExcept=True, **pParams)
+
+    if root:
+        err = getError(uSDC, uRef)
+        errors.append(err)
+
+        cost = getCost(counters)
+        costs.append(cost)
+
+if COMM_WORLD.Get_rank() == 0:
+    errors = [float(e) for e in errors]
+
+    print("errors : ", errors)
+    print("tComps : ", costs)
+    fileName = f"{PATH}/fig06_compTime.json"
+    timings = {}
+    if os.path.isfile(fileName):
+        with open(fileName, "r") as f:
+            timings = json.load(f)
+
+    timings[qDelta + "_MPI" * useMPI] = {"errors": errors, "costs": costs}
+
+    with open(fileName, 'w') as f:
+        json.dump(timings, f, indent=4)

pySDC/projects/parallelSDC_reloaded/scripts/fig06_allenCahnMPI_plot.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 11 11:14:01 2024
+
+Figures with experiments on the Allen-Cahn problem (MPI runs)
+"""
+import os
+import json
+import numpy as np
+
+from pySDC.projects.parallelSDC_reloaded.utils import plt
+
+PATH = '/' + os.path.join(*__file__.split('/')[:-1])
+SCRIPT = __file__.split('/')[-1].split('.')[0]
+
+fileName = f"{PATH}/fig06_compTime.json"
+timings = {}
+with open(fileName, "r") as f:
+    timings = json.load(f)
+
+minPrec = ["MIN-SR-NS", "MIN-SR-S", "MIN-SR-FLEX"]
+
+symList = ['^', '>', '<', 'o', 's', '*', 'p']
+config = [
+    (*minPrec, "VDHS", "ESDIRK43", "LU"),
+]
+nStepsList = np.array([5, 10, 20, 50, 100, 200, 500])
+tEnd = 50
+dtVals = tEnd / nStepsList
+
+# -----------------------------------------------------------------------------
+# %% Error VS cost plots
+# -----------------------------------------------------------------------------
+for qDeltaList in config:
+    figNameConv = f"{SCRIPT}_conv_1"
+    figNameCost = f"{SCRIPT}_cost_1"
+
+    for qDelta, sym in zip(qDeltaList, symList):
+        if qDelta == "MIN-SR-NS":
+            res = timings["MIN-SR-S_MPI"].copy()
+            res["errors"] = [np.nan for _ in res["errors"]]
+        elif qDelta in ["MIN-SR-S", "MIN-SR-FLEX", "VDHS"]:
+            res = timings[f"{qDelta}_MPI"]
+        else:
+            res = timings[qDelta]
+
+        errors = res["errors"]
+        costs = res["costs"]
+
+        ls = '-' if qDelta.startswith("MIN-SR-") else "--"
+
+        plt.figure(figNameConv)
+        plt.loglog(dtVals, errors, sym + ls, label=qDelta)
+
+        plt.figure(figNameCost)
+        plt.loglog(costs, errors, sym + ls, label=qDelta)
+
+    for figName in [figNameConv, figNameCost]:
+        plt.figure(figName)
+        plt.gca().set(
+            xlabel="Computation Time" if "cost" in figName else r"$\Delta {t}$",
+            ylabel=r"$L_2$ error at $T$",
+        )
+        plt.legend()
+        plt.grid(True)
+        plt.tight_layout()
+        plt.savefig(f"{PATH}/{figName}.pdf")
+
+# -----------------------------------------------------------------------------
+# %% Speedup tables
+# -----------------------------------------------------------------------------
+header = f"{'dt size':12} |"
+header += '|'.join(f"  {dt:1.1e}  " for dt in dtVals)
+print(header)
+print("-" * len(header))
+meanEff = 0
+for qDelta in ["MIN-SR-S", "MIN-SR-FLEX", "VDHS"]:
+    seq = timings[f"{qDelta}"]
+    par = timings[f"{qDelta}_MPI"]
+    assert np.allclose(
+        seq["errors"], par["errors"], atol=1e-6
+    ), f"parallel and sequential errors are not close for {qDelta}"
+    tComp = seq["costs"]
+    tCompMPI = par["costs"]
+    meanEff += np.mean([tS / tP for tP, tS in zip(tCompMPI, tComp)])
+    print(f"{qDelta:12} |" + '|'.join(f" {tS/tP:1.1f} ({tS/tP/4*100:1.0f}%) " for tP, tS in zip(tCompMPI, tComp)))
+meanEff /= 3
+print("Mean parallel efficiency :", meanEff / 4)