diff --git a/pySDC/core/sweeper.py b/pySDC/core/sweeper.py index 09ad1c5da4..ff8fe29c8c 100644 --- a/pySDC/core/sweeper.py +++ b/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) diff --git a/pySDC/implementations/sweeper_classes/generic_implicit.py b/pySDC/implementations/sweeper_classes/generic_implicit.py index 6263a9ad76..a70b1045c4 100644 --- a/pySDC/implementations/sweeper_classes/generic_implicit.py +++ b/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 diff --git a/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py b/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py index 26ace3eab5..802dc13731 100644 --- a/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py +++ b/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 diff --git a/pySDC/playgrounds/Diagonal/dahlquist.py b/pySDC/playgrounds/Diagonal/dahlquist.py index b9065fe5e5..dbedcee89d 100644 --- a/pySDC/playgrounds/Diagonal/dahlquist.py +++ b/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 diff --git a/pySDC/playgrounds/Diagonal/optim.py b/pySDC/playgrounds/Diagonal/optim.py index c60461cd2d..4d2253a203 100644 --- a/pySDC/playgrounds/Diagonal/optim.py +++ b/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]: diff --git a/pySDC/playgrounds/dedalus/advection.py b/pySDC/playgrounds/dedalus/advection.py index 94169cfbcc..e487458ff0 100644 --- a/pySDC/playgrounds/dedalus/advection.py +++ b/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 diff --git a/pySDC/playgrounds/dedalus/problem.py b/pySDC/playgrounds/dedalus/problem.py index fce924a454..a474d776c1 100644 --- a/pySDC/playgrounds/dedalus/problem.py +++ b/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): diff --git a/pySDC/playgrounds/dedalus/sdc.py b/pySDC/playgrounds/dedalus/sdc.py index 15895c5670..b56b39c610 100644 --- a/pySDC/playgrounds/dedalus/sdc.py +++ b/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): diff --git a/pySDC/projects/parallelSDC_reloaded/environment.yml b/pySDC/projects/parallelSDC_reloaded/environment.yml index 5202baed8c..c3ab643ada 100644 --- a/pySDC/projects/parallelSDC_reloaded/environment.yml +++ b/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 diff --git a/pySDC/projects/parallelSDC_reloaded/scripts/fig06_allenCahnMPI.py b/pySDC/projects/parallelSDC_reloaded/scripts/fig06_allenCahnMPI.py new file mode 100644 index 0000000000..6b1557b61c --- /dev/null +++ b/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) diff --git a/pySDC/projects/parallelSDC_reloaded/scripts/fig06_allenCahnMPI_plot.py b/pySDC/projects/parallelSDC_reloaded/scripts/fig06_allenCahnMPI_plot.py new file mode 100644 index 0000000000..2b8631070e --- /dev/null +++ b/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) diff --git a/pySDC/projects/parallelSDC_reloaded/scripts/fig06_compTime.json b/pySDC/projects/parallelSDC_reloaded/scripts/fig06_compTime.json new file mode 100644 index 0000000000..c4971916a0 --- /dev/null +++ b/pySDC/projects/parallelSDC_reloaded/scripts/fig06_compTime.json @@ -0,0 +1,162 @@ +{ + "MIN-SR-S": { + "errors": [ + 0.08446537094737497, + 0.02891511011839226, + 0.009359670441056173, + 0.0017675799093321094, + 0.0005909668233092055, + 0.000306143216765387, + 0.0002337400892894797 + ], + "costs": [ + 0.9752397537231445, + 1.5291516780853271, + 2.4031898975372314, + 4.854462146759033, + 7.893767356872559, + 15.199107885360718, + 35.19144558906555 + ] + }, + "MIN-SR-S_MPI": { + "errors": [ + 0.08446537106419487, + 0.028915109678786025, + 0.00935967050769007, + 0.0017675799045971865, + 0.0005909666985126798, + 0.0003061431833988446, + 0.00023374006937152005 + ], + "costs": [ + 0.327730655670166, + 0.5010201930999756, + 0.7522416114807129, + 1.5414834022521973, + 2.3344194889068604, + 4.074150085449219, + 10.160308122634888 + ] + }, + "MIN-SR-FLEX": { + "errors": [ + 0.03961908028763937, + 0.003963578959233919, + 0.0003691116001822194, + 7.944059091555469e-05, + 0.00015503113030299905, + 0.00019389551968445538, + 0.00021499305226875678 + ], + "costs": [ + 1.2781708240509033, + 1.6502108573913574, + 2.552668333053589, + 5.057770013809204, + 8.727833271026611, + 15.858268976211548, + 38.207245111465454 + ] + }, + "MIN-SR-FLEX_MPI": { + "errors": [ + 0.03961908010579853, + 0.003964059881633282, + 0.0003691115106070605, + 7.94406499361501e-05, + 0.00015503115406990015, + 0.00019389561909587447, + 0.00021499306618857255 + ], + "costs": [ + 0.5572898387908936, + 0.6142678260803223, + 0.8634381294250488, + 1.6381163597106934, + 2.969147205352783, + 4.86967134475708, + 12.136021375656128 + ] + }, + "VDHS": { + "errors": [ + 0.5345708075891601, + 0.20140479812575454, + 0.05620974256627712, + 0.008404588766557962, + 0.0016989846252961774, + 0.00016807940141687374, + 0.0001875379340551794 + ], + "costs": [ + 1.2595515251159668, + 1.8388524055480957, + 2.7718563079833984, + 5.170547246932983, + 9.220052480697632, + 15.198433876037598, + 36.616854667663574 + ] + }, + "VDHS_MPI": { + "errors": [ + 0.5345708077129743, + 0.2014047980511802, + 0.056209742298451895, + 0.0084045887253569, + 0.001698984580450411, + 0.00016807941468539286, + 0.0001875379224585031 + ], + "costs": [ + 0.4351925849914551, + 0.630584716796875, + 0.858457088470459, + 1.5588116645812988, + 2.849358081817627, + 4.095851898193359, + 10.163822889328003 + ] + }, + "ESDIRK43": { + "errors": [ + 0.22484823728586387, + 0.018644034062360124, + 0.0013501854062157983, + 0.00023722708007943585, + 0.00022105426812846953, + 0.00022273332127340696, + 0.0002235453435056003 + ], + "costs": [ + 0.49336743354797363, + 0.7599084377288818, + 1.2564899921417236, + 2.5580623149871826, + 4.78606653213501, + 9.57835078239441, + 16.35263991355896 + ] + }, + "LU": { + "errors": [ + 0.0048080489209774285, + 0.0006654206107703006, + 0.00026623686663334385, + 0.0002278193055409893, + 0.0002251848068728332, + 0.000224433969255424, + 0.0002295136162850787 + ], + "costs": [ + 0.9455454349517822, + 1.4841115474700928, + 2.4991953372955322, + 4.7198779582977295, + 7.599146366119385, + 15.103773593902588, + 35.04678988456726 + ] + } +} \ No newline at end of file diff --git a/pySDC/projects/parallelSDC_reloaded/tests/test_parallelSDC_reloaded.py b/pySDC/projects/parallelSDC_reloaded/tests/test_parallelSDC_reloaded.py index 8614f4716b..ef2d6e6fec 100644 --- a/pySDC/projects/parallelSDC_reloaded/tests/test_parallelSDC_reloaded.py +++ b/pySDC/projects/parallelSDC_reloaded/tests/test_parallelSDC_reloaded.py @@ -113,3 +113,15 @@ def test_script_fig05_allenCahn(): assert fig05_allenCahn.config == [ (*minPrec, "VDHS", "ESDIRK43", "LU"), ] + + +@pytest.mark.base +def test_script_fig06_allenCahn(): + # Test fails for python < 3.8, so avoid it + if sys.version_info.minor < 8: + return + + from pySDC.projects.parallelSDC_reloaded.scripts import fig06_allenCahnMPI, fig06_allenCahnMPI_plot + + assert fig06_allenCahnMPI.nSweeps == 4 + assert fig06_allenCahnMPI_plot.minPrec == ["MIN-SR-NS", "MIN-SR-S", "MIN-SR-FLEX"] diff --git a/pySDC/projects/parallelSDC_reloaded/utils.py b/pySDC/projects/parallelSDC_reloaded/utils.py index 898f0253ad..028d4029ac 100644 --- a/pySDC/projects/parallelSDC_reloaded/utils.py +++ b/pySDC/projects/parallelSDC_reloaded/utils.py @@ -210,9 +210,12 @@ def setupProblem(name, description, dt, **kwargs): raise NotImplementedError(f"problem {name} not implemented") -def solutionSDC(tEnd, nSteps, params, probName, **kwargs): +def solutionSDC(tEnd, nSteps, params, probName, verbose=True, noExcept=False, **kwargs): dt = tEnd / nSteps setupProblem(probName, params, dt, **kwargs) + if noExcept: + params["problem_params"]["stop_at_nan"] = False + params["problem_params"]["stop_at_maxiter"] = False controller = controller_nonMPI(num_procs=1, controller_params={'logger_level': 30}, description=params) @@ -225,7 +228,8 @@ def solutionSDC(tEnd, nSteps, params, probName, **kwargs): uSDC[0] = uInit tVals = np.linspace(0, tEnd, nSteps + 1) tBeg = time() - print(" -- computing numerical solution with pySDC ...") + if verbose: + print(" -- computing numerical solution with pySDC ...") warnings.filterwarnings("ignore") for i in range(nSteps): uTmp[:] = uSDC[i] @@ -243,7 +247,8 @@ def solutionSDC(tEnd, nSteps, params, probName, **kwargs): except KeyError: nNewton = 0 nRHS = prob.work_counters["rhs"].niter - print(f" done, newton:{nNewton}, rhs:{nRHS}, tComp:{tComp}") + if verbose: + print(f" done, newton:{nNewton}, rhs:{nRHS}, tComp:{tComp}") try: parallel = controller.MS[0].levels[0].sweep.parallelizable except AttributeError: # pragma: no cover @@ -270,7 +275,7 @@ def solutionExact(tEnd, nSteps, probName, **kwargs): key = f"{tEnd}_{nSteps}" cacheFile = '_solJacobiEllipticExact.json' - # Eventually load already computed solution from local cache + # Potentially load already computed solution from local cache try: with open(f"{PATH}/{cacheFile}", "r") as f: cache = json.load(f)