TL: first revision for the parallel SDC theory paper (#481)

* TL: dump version

* TL: added MPI script for pSDC theory paper

* TL: adapted script for parallel computation

* TL: minor detail for output

* TL: mini-fixes and black

* TL: update for plots

* TL: upload reference data from jusuf runs

* TL: minor modifications

* TL: prepared PR

* TL: fix linting error

* TL: added missing dependency for project tests

* TL: another detail

* TL: eventually <=> finally apparently

pySDC/core/sweeper.py

@@ -263,3 +263,15 @@ def level(self, L):
     @property
     def rank(self):
         return 0
+
+    def updateVariableCoeffs(self, k):
+        """
+        Potentially 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/playgrounds/Diagonal/dahlquist.py

@@ -71,7 +71,7 @@ def setParameters(cls, M=None, nodeDistr=None, quadType=None,
             cls.QDeltaE, cls.dtauE = genQDelta(cls.nodes, explSweep, cls.Q)
             cls.explSweep = explSweep
 
-        # Eventually update nSweep, initSweep and forceProlongation
+        # Potentially update nSweep, initSweep and forceProlongation
         cls.initSweep = cls.initSweep if initSweep is None else initSweep
         cls.nSweep = cls.nSweep if nSweep is None else nSweep
         cls.forceProl = cls.forceProl if forceProl is None else forceProl

pySDC/playgrounds/Diagonal/optim.py

@@ -50,7 +50,7 @@ def findRes(x):
 
         funcEval = func(x0)
 
-        # Eventually add local minimum to results
+        # Potentially add local minimum to results
         xOrig = findRes(x0)
         if xOrig:
             if funcEval < res[xOrig]:

pySDC/playgrounds/dedalus/advection.py

@@ -81,7 +81,7 @@ def getErr(nStep):
     sym = '^-' if useSDC else 'o-'
     plt.loglog(dt, err, sym, label=lbl)
 
-    # Eventually plot order curve
+    # Potentially plot order curve
     if name in orderPlot:
         order = orderPlot[name]
         c = err[-1]/dt[-1]**order * 2

pySDC/playgrounds/dedalus/problem.py

@@ -117,7 +117,7 @@ def updateLHS(self, dt, qI, init=False):
         # Update LHS and LHS solvers for each subproblems
         for sp in solver.subproblems:
             if self.init:
-                # Eventually instanciate list of solvers (ony first time step)
+                # Instanciate list of solvers (ony first time step)
                 sp.LHS_solvers = [None] * self.M
                 self.init = False
             for i in range(self.M):

pySDC/playgrounds/dedalus/sdc.py

@@ -173,7 +173,7 @@ def setParameters(cls, nNodes=None, nodeType=None, quadType=None,
                 Q=cls.Q, nodes=cls.nodes, nNodes=nNodes,
                 nodeType=nodeType, quadType=quadType)
 
-        # Eventually update additional parameters
+        # Potentially update additional parameters
         if forceProl is not None: cls.forceProl = forceProl
         if initSweep is not None: cls.initSweep = initSweep
         if not keepNSweeps:
@@ -302,7 +302,7 @@ def _updateLHS(self, dt, init=False):
         # Update LHS and LHS solvers for each subproblems
         for sp in solver.subproblems:
             if init:
-                # Eventually instantiate list of solver (ony first time step)
+                # Potentially instantiate list of solver (ony first time step)
                 sp.LHS_solvers = [[None for _ in range(self.M)] for _ in range(self.nSweeps)]
             for k in range(self.nSweeps):
                 for m in range(self.M):

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,107 @@
+#!/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]
+    if qDelta.startswith("--"):
+        qDelta = "MIN-SR-FLEX"
+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)