diff --git a/.github/coverage/floors.json b/.github/coverage/floors.json index c16b06851f..8e220e39ba 100644 --- a/.github/coverage/floors.json +++ b/.github/coverage/floors.json @@ -116,7 +116,7 @@ "proteus/mprans/MoveMesh.py": 86.4, "proteus/mprans/MoveMeshMonitor.py": 65.9, "proteus/mprans/NCLS.py": 75.4, - "proteus/mprans/NCLS3P.py": 76.2, + "proteus/mprans/NCLS3P.py": 75.6, "proteus/mprans/Pres.py": 76.1, "proteus/mprans/PresInc.py": 76.2, "proteus/mprans/PresInit.py": 23.2, diff --git a/proteus/NumericalFlux.py b/proteus/NumericalFlux.py index 2ed93b556d..e6dd63d0e2 100644 --- a/proteus/NumericalFlux.py +++ b/proteus/NumericalFlux.py @@ -48,6 +48,23 @@ "Stress_SIPG_exterior", ]) +#: Boundary-only (``hasInterior=False``) diffusion fluxes that refuse periodic conditions. NF_base marks the face +#: points of a periodic pair as DOF boundary with no value, and these fluxes then apply their weak-Dirichlet diffusive +#: flux there with the trace equal to the interior value: the residual keeps the one-sided consistent flux +#: -a grad(u).n on each side of the seam (which does not cancel: grad(u) jumps across it) and the Jacobian a penalty +#: the residual does not have. Measured on a C0 P1 periodic Poisson problem (-Laplace u = f on the unit square, +#: periodic in x, u = 0 at y = 0, 1; 8x8 cells), L2 error and Newton iterations: Diffusion_IIPG_exterior 3.2e-2 in 1 +#: without the periodic pair (Dirichlet at x = 0, 1 instead), 3.5e-1 and no convergence in 30 with it; +#: Diffusion_SIPG_exterior 4.3e-2 / 4.4e-1. The Stokes pair applies the same diffusive flux at the same points. Use +#: StrongDirichlet (StrongDirichletFactory) for a periodic C0 problem: 4.7e-2 in 1, as ADRProblem does. Matched by +#: exact class name: subclasses with their own exterior flux are not covered. +PERIODIC_SEAM_REFUSED = frozenset([ + "Diffusion_IIPG_exterior", + "Diffusion_SIPG_exterior", + "Stokes_Advection_DiagonalUpwind_Diffusion_IIPG_exterior", + "Stokes_Advection_DiagonalUpwind_Diffusion_SIPG_exterior", +]) + class NF_base(object): useWeakDirichletConditions=True @@ -92,6 +109,13 @@ def __init__(self,vt, parallelPeriodic=False): import copy name = type(self).__name__ + if getPeriodicBoundaryConditions is not None and name in PERIODIC_SEAM_REFUSED: + raise ValueError( + "%s cannot be used with periodic conditions: it applies its weak Dirichlet diffusive flux at the " + "periodic face points, which then carry a one-sided flux in the residual and a penalty in the " + "Jacobian (a C0 periodic Poisson problem: 10x the error, Newton does not converge). Use " + "NumericalFlux.StrongDirichlet (or StrongDirichletFactory) for a periodic problem; see " + "NumericalFlux.PERIODIC_SEAM_REFUSED" % name) if name in UNUSED_FLUX_CLASSES_REFUSED_IN_PARALLEL: from . import Comm comm = Comm.get() @@ -3026,7 +3050,7 @@ class RusanovNumericalFlux_Diagonal(Advection_DiagonalUpwind): def __init__(self,vt,getPointwiseBoundaryConditions, getAdvectiveFluxBoundaryConditions, getDiffusiveFluxBoundaryConditions, - getPeriodicConditions): + getPeriodicConditions=None): Advection_DiagonalUpwind.__init__(self,vt, getPointwiseBoundaryConditions, getAdvectiveFluxBoundaryConditions, @@ -3270,7 +3294,7 @@ class ConvexOneSonicPointNumericalFlux(Advection_DiagonalUpwind): def __init__(self,vt,getPointwiseBoundaryConditions, getAdvectiveFluxBoundaryConditions, getDiffusiveFluxBoundaryConditions, - getPeriodicBoundaryConditions, + getPeriodicBoundaryConditions=None, sonicPoint=0.0,sonicFlux=0.0): Advection_DiagonalUpwind.__init__(self,vt,getPointwiseBoundaryConditions, getAdvectiveFluxBoundaryConditions, @@ -4722,10 +4746,12 @@ class RusanovNumericalFlux(RusanovNumericalFlux_Diagonal): def __init__(self,vt,getPointwiseBoundaryConditions, getAdvectiveFluxBoundaryConditions, - getDiffusiveFluxBoundaryConditions): + getDiffusiveFluxBoundaryConditions, + getPeriodicConditions=None): RusanovNumericalFlux_Diagonal.__init__(self,vt,getPointwiseBoundaryConditions, getAdvectiveFluxBoundaryConditions, - getDiffusiveFluxBoundaryConditions) + getDiffusiveFluxBoundaryConditions, + getPeriodicConditions) self.ebqe['eigen_bound']=numpy.copy(self.ebqe[('u',0)]) if 'eigen_bound' not in self.vt.q: diff --git a/proteus/NumericalSolution.py b/proteus/NumericalSolution.py index b0db13657b..14de7afd36 100644 --- a/proteus/NumericalSolution.py +++ b/proteus/NumericalSolution.py @@ -1372,6 +1372,9 @@ def saveSystemStepState(self): levels.append(([lm.u[ci].dof.copy() for ci in range(lm.coefficients.nc)], owners, {k: v.copy() for k, v in self._coupledFields(lm)}, m.uList[level].copy())) self._systemStepState.append(levels) + # where each model's last step ended: SO_base.modelStepTaken (Sequential_MinModelStep, the default system + # controller) moves it to the end of the system step as soon as the model's step is taken + self._systemStepModelTimes = [getattr(m.stepController, "t_model_last", None) for m in self.modelList] #: the quadrature-point and face fields that one model reads from another (RANS2P reads VOF's and NCLS's #: ('u', 0) and ('grad(u)', 0), the transport models read the flow model's ('velocity', 0), ...); they are @@ -1409,6 +1412,11 @@ def restoreSystemStepState(self): for owner, st in owners.items(): if st is not None: getattr(lm, owner).set_state(copy.deepcopy(st)) + # a model that had finished in the failed attempt kept its end as the end of its last step, and the retry, + # stepping to an earlier time, skipped it (and every later step of it lagged one step behind) + for m, t in zip(self.modelList, getattr(self, "_systemStepModelTimes", ())): + if t is not None: + m.stepController.t_model_last = t self.endSystemStep(accepted=False) logEvent("restored the models to the start of the system step for the retry", level=3) diff --git a/proteus/ctimeIntegration.pyx b/proteus/ctimeIntegration.pyx index 260a99edd6..6c40542490 100644 --- a/proteus/ctimeIntegration.pyx +++ b/proteus/ctimeIntegration.pyx @@ -141,79 +141,75 @@ cdef extern from "timeIntegration.h": double* Uin, double* Uout) # ------------------------------------------- # -def psiTCtteDT(int nPoints, - double tau, - double dtn, - double dtnm1, - np.ndarray yn, - np.ndarray ypn, - np.ndarray ypnm1, - np.ndarray dtnp1): - cpsiTCtteDT(nPoints, - tau, - dtn, - dtnm1, - yn.data, - ypn.data, - ypnm1.data, - dtnp1.data) +# The Python signatures below are the ones TimeIntegration.py calls, the +# contract of the C module these wrappers replaced in 2019 (#879): the sizes +# come from the array shapes, the limiter flag is the last argument and +# psiTCtteDT returns the new step. The 2019 wrappers asked the caller for the +# sizes instead, so every call in TimeIntegration.py raised TypeError. +def psiTCtteDT(double tau, + double dtn, + double dtnm1, + np.ndarray yn, + np.ndarray ypn, + np.ndarray ypnm1): + cdef double dtnp1 = 0.0 + cpsiTCtteDT(yn.size, + tau, + dtn, + dtnm1, + yn.data, + ypn.data, + ypnm1.data, + &dtnp1) + return dtnp1 -def applyDGlimitingP1Lagrange1d(int limiterFlag, - int nElements_global, - int nNodes_element, - int nElementBoundaries_element, - int nDOF_element, - np.ndarray elementNodesArray, - np.ndarray elementNeighborsArray, - np.ndarray nodeArray, - np.ndarray elementBarycentersArray, - np.ndarray l2g, - np.ndarray tag, - np.ndarray Uin, - np.ndarray Uout): +def applyDGlimitingP1Lagrange1d(np.ndarray elementNodesArray, + np.ndarray elementNeighborsArray, + np.ndarray nodeArray, + np.ndarray elementBarycentersArray, + np.ndarray l2g, + np.ndarray tag, + np.ndarray Uin, + np.ndarray Uout, + int limiterFlag=1): capplyDGlimitingP1Lagrange1d(limiterFlag, - nElements_global, - nNodes_element, - nElementBoundaries_element, - nDOF_element, - elementNodesArray.data, - elementNeighborsArray.data, - nodeArray.data, - elementBarycentersArray.data, - l2g.data, - tag.data, - Uin.data, - Uout.data) + elementNodesArray.shape[0], + elementNodesArray.shape[1], + elementNeighborsArray.shape[1], + l2g.shape[1], + elementNodesArray.data, + elementNeighborsArray.data, + nodeArray.data, + elementBarycentersArray.data, + l2g.data, + tag.data, + Uin.data, + Uout.data) def applyDGlimitingP1Lagrange1d_withVacuumTol(int enforcePositivity, - double vacuumTol, - int nElements_global, - int nNodes_element, - int nElementBoundaries_element, - int nDOF_element, - np.ndarray elementNodesArray, - np.ndarray elementNeighborsArray, - np.ndarray nodeArray, - np.ndarray elementBarycentersArray, - np.ndarray l2g, - np.ndarray tag, - np.ndarray Uin, - np.ndarray Uout): + double vacuumTol, + np.ndarray elementNodesArray, + np.ndarray elementNeighborsArray, + np.ndarray nodeArray, + np.ndarray elementBarycentersArray, + np.ndarray l2g, + np.ndarray tag, + np.ndarray Uin, + np.ndarray Uout): capplyDGlimitingP1Lagrange1d_withVacuumTol(enforcePositivity, - vacuumTol, - nElements_global, - nNodes_element, - nElementBoundaries_element, - nDOF_element, - elementNodesArray.data, - elementNeighborsArray.data, - nodeArray.data, - elementBarycentersArray.data, - l2g.data, - tag.data, - Uin.data, - Uout.data) - + vacuumTol, + elementNodesArray.shape[0], + elementNodesArray.shape[1], + elementNeighborsArray.shape[1], + l2g.shape[1], + elementNodesArray.data, + elementNeighborsArray.data, + nodeArray.data, + elementBarycentersArray.data, + l2g.data, + tag.data, + Uin.data, + Uout.data) def computeElementNeighborShapeGradients(np.ndarray elementBoundariesArray, np.ndarray elementNeighborsArray, @@ -271,128 +267,113 @@ def applyCockburnDGlimiterP1Lagrange2d(double nu, Uout.data) def applyDurlofskyDGlimiterP1Lagrange2d(int killExtrema, - int allowMinWithUndershoot, - int nElements_global, - int nElementBoundaries_element, - int nNodes_element, - int nSpace, - int nDOF_element, - np.ndarray elementNeighborsArray, - np.ndarray elementBoundariesArray, - np.ndarray elementNodesArray, - np.ndarray nodeArray, - np.ndarray elementBarycentersArray, - np.ndarray elementBoundaryBarycentersArray, - np.ndarray elementNeighborShapeGradients, - np.ndarray l2g, - np.ndarray grad_v0, - np.ndarray elementAverages, - np.ndarray tag, - np.ndarray Uin, - np.ndarray Uout): + int allowMinWithUndershoot, + np.ndarray elementNeighborsArray, + np.ndarray elementBoundariesArray, + np.ndarray elementNodesArray, + np.ndarray nodeArray, + np.ndarray elementBarycentersArray, + np.ndarray elementBoundaryBarycentersArray, + np.ndarray elementNeighborShapeGradients, + np.ndarray l2g, + np.ndarray grad_v0, + np.ndarray elementAverages, + np.ndarray tag, + np.ndarray Uin, + np.ndarray Uout): capplyDurlofskyDGlimiterP1Lagrange2d(killExtrema, - allowMinWithUndershoot, - nElements_global, - nElementBoundaries_element, - nNodes_element, - nSpace, - nDOF_element, - elementNeighborsArray.data, - elementBoundariesArray.data, - elementNodesArray.data, - nodeArray.data, - elementBarycentersArray.data, - elementBoundaryBarycentersArray.data, - elementNeighborShapeGradients.data, - l2g.data, - grad_v0.data, - elementAverages.data, - tag.data, - Uin.data, - Uout.data) + allowMinWithUndershoot, + elementNeighborsArray.shape[0], + elementNeighborsArray.shape[1], + elementNodesArray.shape[1], + grad_v0.shape[3], + l2g.shape[1], + elementNeighborsArray.data, + elementBoundariesArray.data, + elementNodesArray.data, + nodeArray.data, + elementBarycentersArray.data, + elementBoundaryBarycentersArray.data, + elementNeighborShapeGradients.data, + l2g.data, + grad_v0.data, + elementAverages.data, + tag.data, + Uin.data, + Uout.data) def applyDurlofskyDGlimiterP1Lagrange3d(int killExtrema, - int allowMinWithUndershoot, - int nElements_global, - int nElementBoundaries_element, - int nNodes_element, - int nSpace, - int nDOF_element, - np.ndarray elementNeighborsArray, - np.ndarray elementBoundariesArray, - np.ndarray elementNodesArray, - np.ndarray nodeArray, - np.ndarray elementBarycentersArray, - np.ndarray elementBoundaryBarycentersArray, - np.ndarray elementNeighborShapeGradients, - np.ndarray l2g, - np.ndarray grad_v0, - np.ndarray elementAverages, - np.ndarray tag, - np.ndarray Uin, - np.ndarray Uout): - capplyDurlofskyDGlimiterP1Lagrange3d(killExtrema, - allowMinWithUndershoot, - nElements_global, - nElementBoundaries_element, - nNodes_element, - nSpace, - nDOF_element, - elementNeighborsArray.data, - elementBoundariesArray.data, - elementNodesArray.data, - nodeArray.data, - elementBarycentersArray.data, - elementBoundaryBarycentersArray.data, - elementNeighborShapeGradients.data, - l2g.data, - grad_v0.data, - elementAverages.data, - tag.data, - Uin.data, - Uout.data) + int allowMinWithUndershoot, + np.ndarray elementNeighborsArray, + np.ndarray elementBoundariesArray, + np.ndarray elementNodesArray, + np.ndarray nodeArray, + np.ndarray elementBarycentersArray, + np.ndarray elementBoundaryBarycentersArray, + np.ndarray elementNeighborShapeGradients, + np.ndarray l2g, + np.ndarray grad_v0, + np.ndarray elementAverages, + np.ndarray tag, + np.ndarray Uin, + np.ndarray Uout): + capplyDurlofskyDGlimiterP1Lagrange3d(killExtrema, + allowMinWithUndershoot, + elementNeighborsArray.shape[0], + elementNeighborsArray.shape[1], + elementNodesArray.shape[1], + grad_v0.shape[3], + l2g.shape[1], + elementNeighborsArray.data, + elementBoundariesArray.data, + elementNodesArray.data, + nodeArray.data, + elementBarycentersArray.data, + elementBoundaryBarycentersArray.data, + elementNeighborShapeGradients.data, + l2g.data, + grad_v0.data, + elementAverages.data, + tag.data, + Uin.data, + Uout.data) def applyDurlofskyDGlimiterP1Lagrange2d_withVacuumTol(int killExtrema, - int allowMinWithUndershoot, - int enforcePositivity, - double vacuumTol, - int nElements_global, - int nElementBoundaries_element, - int nNodes_element, - int nSpace, - int nDOF_element, - np.ndarray elementNeighborsArray, - np.ndarray elementBoundariesArray, - np.ndarray elementNodesArray, - np.ndarray nodeArray, - np.ndarray elementBarycentersArray, - np.ndarray elementBoundaryBarycentersArray, - np.ndarray elementNeighborShapeGradients, - np.ndarray l2g, - np.ndarray grad_v0, - np.ndarray elementAverages, - np.ndarray tag, - np.ndarray Uin, - np.ndarray Uout): + int allowMinWithUndershoot, + int enforcePositivity, + double vacuumTol, + np.ndarray elementNeighborsArray, + np.ndarray elementBoundariesArray, + np.ndarray elementNodesArray, + np.ndarray nodeArray, + np.ndarray elementBarycentersArray, + np.ndarray elementBoundaryBarycentersArray, + np.ndarray elementNeighborShapeGradients, + np.ndarray l2g, + np.ndarray grad_v0, + np.ndarray elementAverages, + np.ndarray tag, + np.ndarray Uin, + np.ndarray Uout): capplyDurlofskyDGlimiterP1Lagrange2d_withVacuumTol(killExtrema, - allowMinWithUndershoot, - enforcePositivity, - vacuumTol, - nElements_global, - nElementBoundaries_element, - nNodes_element, - nSpace, - nDOF_element, - elementNeighborsArray.data, - elementBoundariesArray.data, - elementNodesArray.data, - nodeArray.data, - elementBarycentersArray.data, - elementBoundaryBarycentersArray.data, - elementNeighborShapeGradients.data, - l2g.data, - grad_v0.data, - elementAverages.data, - tag.data, - Uin.data, - Uout.data) + allowMinWithUndershoot, + enforcePositivity, + vacuumTol, + elementNeighborsArray.shape[0], + elementNeighborsArray.shape[1], + elementNodesArray.shape[1], + grad_v0.shape[3], + l2g.shape[1], + elementNeighborsArray.data, + elementBoundariesArray.data, + elementNodesArray.data, + nodeArray.data, + elementBarycentersArray.data, + elementBoundaryBarycentersArray.data, + elementNeighborShapeGradients.data, + l2g.data, + grad_v0.data, + elementAverages.data, + tag.data, + Uin.data, + Uout.data) diff --git a/proteus/mprans/RANS2P_IB.h b/proteus/mprans/RANS2P_IB.h index aa184af76c..e67d03cc58 100644 --- a/proteus/mprans/RANS2P_IB.h +++ b/proteus/mprans/RANS2P_IB.h @@ -5102,6 +5102,9 @@ namespace proteus proteus::pyarray& nodeDiametersArray = args.array("nodeDiametersArray"); double hFactor = args.scalar("hFactor"); int nElements_global = args.scalar("nElements_global"); + // the beam loads (q1-q3), vel_avg and netBeamDrag are summed over ranks by the caller, so only owned + // elements add to them; every local element still gets q_dragBeam*, the fluid's source + int nElements_owned = args.scalar("nElements_owned"); int nElementBoundaries_owned = args.scalar("nElementBoundaries_owned"); double useRBLES = args.scalar("useRBLES"); double useMetrics = args.scalar("useMetrics"); @@ -5257,6 +5260,8 @@ namespace proteus double globalConservationError=0.0; for(int eN=0;eN 128 else None + try: + name = signal.Signals(sig).name if sig else None + except ValueError: + name = None + if returncode < 0: + return "killed by %s" % (name or "signal %d" % -returncode) + return "exit status %d%s" % (returncode, " (%s?)" % name if name else "") + + +def report(cwd, what, argv, stdout, stderr): + """The assertion message for a failed launch; writes the full output under ``cwd``.""" + stdout, stderr = _text(stdout), _text(stderr) + paths = [] + for name, text in (("launch.stdout.txt", stdout), ("launch.stderr.txt", stderr)): + path = os.path.join(str(cwd), name) + try: + with open(path, "w") as f: + f.write(text) + paths.append(path) + except OSError as e: + paths.append("%s (not written: %s)" % (path, e)) + return ("%s: %s\nfull output: %s\n--- stderr, last %d of %d characters ---\n%s\n" + "--- stdout, last %d of %d characters ---\n%s" + % (" ".join(str(a) for a in argv), what, ", ".join(paths), STDERR_TAIL, len(stderr), + stderr[-STDERR_TAIL:], STDOUT_TAIL, len(stdout), stdout[-STDOUT_TAIL:])) + + +def run(argv, cwd, env, timeout): + """subprocess.run(argv, capture_output=True, text=True, ...) that raises AssertionError with ``report``'s message + on a nonzero exit status or a timeout, and returns the CompletedProcess otherwise.""" + try: + out = subprocess.run(argv, capture_output=True, text=True, timeout=timeout, env=env, cwd=str(cwd)) + except subprocess.TimeoutExpired as e: + raise AssertionError(report(cwd, "timed out after %s s" % timeout, argv, e.stdout, e.stderr)) from None + if out.returncode != 0: + raise AssertionError(report(cwd, _status(out.returncode), argv, out.stdout, out.stderr)) + return out diff --git a/test/mcorr_partition_mpi.py b/test/mcorr_partition_mpi.py index 9683354d6f..f57a420d19 100644 --- a/test/mcorr_partition_mpi.py +++ b/test/mcorr_partition_mpi.py @@ -10,6 +10,7 @@ import os, shutil, sys import numpy as np from proteus import Comm, MeshTools as MT, default_s +import rank_staging # every rank stages the same files: atomically SRC = os.path.join(os.path.dirname(os.path.abspath(__file__)), "LS_with_edgeBased_EV", "MCorr") os.makedirs("mcpkg", exist_ok=True) @@ -18,7 +19,7 @@ text = open(os.path.join(SRC, name)).read() if name == "thelper_cons_ls.py": text = text.replace("parallel = False", "parallel = " + os.environ.get("PAR", "True")) - open(os.path.join("mcpkg", name), "w").write(text) + rank_staging.write(os.path.join("mcpkg", name), text) sys.path.insert(0, os.getcwd()) import proteus diff --git a/test/ncls_edge_partition_mpi.py b/test/ncls_edge_partition_mpi.py index 6dde9da60e..655edd8e8a 100644 --- a/test/ncls_edge_partition_mpi.py +++ b/test/ncls_edge_partition_mpi.py @@ -13,6 +13,7 @@ import numpy as np from proteus import Comm, MeshTools as MT from mpi4py import MPI +import rank_staging # every rank stages the same files: atomically SRC = os.path.join(os.path.dirname(os.path.abspath(__file__)), "LS_with_edgeBased_EV", "NCLS") sys.path.insert(0, os.getcwd()) @@ -23,7 +24,7 @@ if name == "ncls_n.py": # converge the mass matrix solves far enough for a comparison to 1e-8 text += "\nnl_atol_res = 1.0e-13; l_atol_res = 1.0e-15; linTolFac = 0.0; tolFac = 0.0\n" - open(name, "w").write(text) + rank_staging.write(name, text) _partition = MT.Mesh.partitionMesh def partitionMesh(self, nLayersOfOverlap=1, parallelPartitioningType=MT.MeshParallelPartitioningTypes.node, diff --git a/test/rank_staging.py b/test/rank_staging.py new file mode 100644 index 0000000000..6ff5897ac0 --- /dev/null +++ b/test/rank_staging.py @@ -0,0 +1,32 @@ +"""Stage a deck's files into the working directory of a test run under mpiexec, from every rank. + +The ``*_mpi.py`` helpers run on every rank and each rank copies (and edits) the deck's files into the shared working +directory before importing them. Written in place (``open(path, "w").write(text)``, ``shutil.copy``, an append), one +rank truncates a file another rank is importing: that rank star-imports a truncated module and stops with NameError, +and the other ranks hang in their first collective (test_mcorr_partition_independence on the GitHub aarch64 runner, +2026-11-06: ``name 'RD_model' is not defined``, then a 2809 s timeout). Here every rank writes the same bytes to a name +of its own and renames it into place (``os.replace`` is atomic on POSIX), so a reader sees a whole file, the previous +one or this one, and the files are the same on every rank. No collective, so it can run before ``Comm.init``. +""" +import os +import time + + +def write(path, data, _pause=0.0): + """Put ``data`` (str or bytes) at ``path`` without a moment at which another process can read part of it. + ``_pause`` (seconds, tests of this module only) holds the written temporary file before the rename.""" + tmp = "%s.%d.staging" % (path, os.getpid()) + with open(tmp, "wb") as f: + f.write(data.encode() if isinstance(data, str) else data) + if _pause: + time.sleep(_pause) + os.replace(tmp, path) + + +def copy(src, dst): + """``shutil.copy(src, dst)`` (dst a file or a directory), staged as ``write`` does.""" + if os.path.isdir(dst): + dst = os.path.join(dst, os.path.basename(src)) + with open(src, "rb") as f: + write(dst, f.read()) + return dst diff --git a/test/rans2p_ib/beam_channel_mpi.py b/test/rans2p_ib/beam_channel_mpi.py new file mode 100644 index 0000000000..870416ca1d --- /dev/null +++ b/test/rans2p_ib/beam_channel_mpi.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python +"""Run the RANS2P_IB beam channel and print what a serial run must match. + +Run by test_rans2p_ib.py, under mpiexec for several ranks. Prints the sums over owned nodes of p, u, v, w (and of +their squares) at T, the beams' drag and the beam loads q1, q2, q3, summed. Direct solves (SuperLU / SuperLU_DIST) and +a 1e-10 Newton tolerance, so the comparison is not limited by a Krylov tolerance. Environment: PARTITION (node or +element, default node), NLAYERS (default 0), OPTS (the deck's context options, e.g. beamRigid=False). +""" +import os +import sys + +from proteus import Comm + +comm = Comm.init() + +import numpy as np +from mpi4py import MPI +from petsc4py import PETSc + +from proteus import Context, MeshTools, NumericalSolution, default_s, default_so, defaults, runtime + +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +Context.contextOptionsString = os.environ.get("OPTS", "") +p = defaults.load_physics("beam_channel_p", HERE) +n = defaults.load_numerics("beam_channel_n", HERE) +n.parallelPartitioningType = getattr(MeshTools.MeshParallelPartitioningTypes, os.environ.get("PARTITION", "node")) +n.nLayersOfOverlapForParallel = int(os.environ.get("NLAYERS", "0")) +options = PETSc.Options() +options.setValue("rans2p_ksp_type", "preonly") +options.setValue("rans2p_pc_type", "lu") +options.setValue("rans2p_pc_factor_mat_solver_type", "superlu_dist" if comm.size() > 1 else "superlu") +so = default_so +so.name = p.name +so.pnList = [(p, n)] +so.sList = [default_s] +so.tnList = [0.0] + [n.DT * (i + 1) for i in range(n.nDTout)] +ns = NumericalSolution.NS_base(so, [p], [n], so.sList, runtime.RunOptions()) +ns.calculateSolution(p.name) +model = ns.modelList[0].levelModelList[-1] +c = comm.comm.tompi4py() +owned = model.mesh.nNodes_owned +result = {} +for i, name in enumerate("puvw"): + dof = model.u[i].dof[:owned] + result[name + "2"] = c.allreduce(float((dof ** 2).sum()), op=MPI.SUM) + result[name + "s"] = c.allreduce(float(dof.sum()), op=MPI.SUM) +coefficients = model.coefficients +for i in range(3): + result["drag%d" % i] = float(coefficients.beamDrag[i]) +for name in ("q1", "q2", "q3", "xv", "yv", "zv"): + result[name] = float(getattr(coefficients, name).sum()) +if comm.rank() == 0: + print("RESULT ranks=%d " % comm.size() + " ".join("%s=%.17e" % kv for kv in result.items()), flush=True) diff --git a/test/rans2p_ib/beam_channel_n.py b/test/rans2p_ib/beam_channel_n.py new file mode 100644 index 0000000000..9667d5e540 --- /dev/null +++ b/test/rans2p_ib/beam_channel_n.py @@ -0,0 +1,42 @@ +from proteus import * +from proteus.default_n import * +from beam_channel_p import * + +timeIntegration = BackwardEuler_cfl +stepController = FixedStep +runCFL = 0.9 +DT = 0.01 +nDTout = 4 + +femSpaces = {0: C0_AffineLinearOnSimplexWithNodalBasis, + 1: C0_AffineLinearOnSimplexWithNodalBasis, + 2: C0_AffineLinearOnSimplexWithNodalBasis, + 3: C0_AffineLinearOnSimplexWithNodalBasis} +elementQuadrature = SimplexGaussQuadrature(nd, 3) +elementBoundaryQuadrature = SimplexGaussQuadrature(nd - 1, 3) + +nnx, nny, nnz = 9, 14, 5 +nLevels = 1 + +numericalFluxType = RANS2P_IB.RANS2P.NumericalFlux +subgridError = RANS2P_IB.RANS2P.SubgridError(coefficients, nd, lag=False, hFactor=1.0) +shockCapturing = RANS2P_IB.RANS2P.ShockCapturing(coefficients, nd, 0.0, lag=False) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton +nonlinearSmoother = None +linearSmoother = None +matrix = SparseMatrix +multilevelLinearSolver = KSP_petsc4py +levelLinearSolver = KSP_petsc4py +linear_solver_options_prefix = "rans2p_" +levelNonlinearSolverConvergenceTest = "r" +linearSolverConvergenceTest = "r-true" +tolFac = 0.0 +nl_atol_res = 1.0e-10 +l_atol_res = 1.0e-12 +maxNonlinearIts = 20 +maxLineSearches = 0 +parallelPartitioningType = MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 diff --git a/test/rans2p_ib/beam_channel_p.py b/test/rans2p_ib/beam_channel_p.py new file mode 100644 index 0000000000..afcedec991 --- /dev/null +++ b/test/rans2p_ib/beam_channel_p.py @@ -0,0 +1,93 @@ +"""One rigid beam (a vegetation stem) in a 3D water channel: the RANS2P_IB deck. + +RANS2P_IB is RANS2P plus the drag of immersed Euler-Bernoulli beams (beamFEM), spread onto the mesh with a smoothed +delta function. Single phase, gravity along -z. Starts from rest: inflow U at x = 0, p = 0 at x = Lx, no slip on the +bottom and the side walls, free slip on top. The beam stands at (0.2, 0.345) from z = 0, inside the drag-averaging box +y in [0.1725, 0.5175] that RANS2P_IB hard-codes. Context option beamRigid=False makes it a soft flexible beam (it +bends about 1 cm in 0.04 s). +""" +from proteus import * +from proteus.default_p import * +from proteus import Context +from proteus.mprans import RANS2P_IB + +opts = Context.Options([("beamRigid", True, "a rigid beam, or a flexible one (EI = GJ = 1e-4)")]) + +nd = 3 +L = (0.4, 0.69, 0.2) +domain = None +name = "beam_channel" +T = 0.04 +U = 0.1 +rho_0, nu_0 = 998.2, 1.0e-3 +g = [0.0, 0.0, -9.8] +eps = 1.0e-8 + +LevelModelType = RANS2P_IB.LevelModel +coefficients = RANS2P_IB.Coefficients(epsFact=1.5, rho_0=rho_0, nu_0=nu_0, rho_1=rho_0, nu_1=nu_0, g=g, nd=nd, + LS_model=None, forceStrongDirichlet=False, useMetrics=1.0, + beamLocation=[(0.2, 0.345)], beamLength=[0.1], beamRadius=[0.005], + EI=[1.0e3 if opts.beamRigid else 1.0e-4], + GJ=[1.0e3 if opts.beamRigid else 1.0e-4], nBeamElements=4, beam_quadOrder=3, + beamRigid=opts.beamRigid, beamFilename="Beams") + + +def _inflow(x): + return x[0] < eps + + +def _outflow(x): + return x[0] > L[0] - eps + + +def _wall(x): + return x[1] < eps or x[1] > L[1] - eps or x[2] < eps + + +def getDBC_p(x, flag): + if _outflow(x): + return lambda x, t: 0.0 + + +def getDBC_u(x, flag): + if _inflow(x): + return lambda x, t: U + if _wall(x) and not _outflow(x): + return lambda x, t: 0.0 + + +def getDBC_vw(x, flag): + if _inflow(x) or (_wall(x) and not _outflow(x)): + return lambda x, t: 0.0 + + +dirichletConditions = {0: getDBC_p, 1: getDBC_u, 2: getDBC_vw, 3: getDBC_vw} + + +def getAFBC_p(x, flag): + if _inflow(x): + return lambda x, t: -U + if not _outflow(x): + return lambda x, t: 0.0 + + +def getAFBC_vel(x, flag): + if not (_inflow(x) or _outflow(x)): + return lambda x, t: 0.0 + + +def getDFBC_vel(x, flag): + if not (_inflow(x) or _wall(x)): + return lambda x, t: 0.0 + + +advectiveFluxBoundaryConditions = {0: getAFBC_p, 1: getAFBC_vel, 2: getAFBC_vel, 3: getAFBC_vel} +diffusiveFluxBoundaryConditions = {0: {}, 1: {1: getDFBC_vel}, 2: {2: getDFBC_vel}, 3: {3: getDFBC_vel}} + + +class Zero(object): + def uOfXT(self, x, t): + return 0.0 + + +initialConditions = {0: Zero(), 1: Zero(), 2: Zero(), 3: Zero()} diff --git a/test/rans2p_ib/beam_channel_so.py b/test/rans2p_ib/beam_channel_so.py new file mode 100644 index 0000000000..03def98a83 --- /dev/null +++ b/test/rans2p_ib/beam_channel_so.py @@ -0,0 +1,12 @@ +"""The beam channel with a second model after RANS2P_IB: a dye (VOF, carried by the flow), so a failure of the dye +model comes after RANS2P_IB's postStep has moved a flexible beam (test_step_retry's multi-model RANS2P_IB cases).""" +from proteus.default_so import * +from proteus import Context +import beam_channel_p + +pnList = [("beam_channel_p", "beam_channel_n"), ("dye_p", "dye_n")] +name = "beam_channel_dye" +systemStepControllerType = Sequential_MinModelStep +needEBQ_GLOBAL = False +needEBQ = False +tnList = [0.0, 0.01, 0.02, 0.03, 0.04] diff --git a/test/rans2p_ib/dye_n.py b/test/rans2p_ib/dye_n.py new file mode 100644 index 0000000000..5bbbf499ac --- /dev/null +++ b/test/rans2p_ib/dye_n.py @@ -0,0 +1,37 @@ +from proteus.default_n import * +from proteus.mprans import VOF +import dye_p as physics +from beam_channel_n import nnx, nny, nnz, nLevels, DT, nDTout, femSpaces as _flow_spaces + +timeIntegration = BackwardEuler_cfl +stepController = FixedStep +runCFL = 0.9 + +femSpaces = {0: C0_AffineLinearOnSimplexWithNodalBasis} +elementQuadrature = SimplexGaussQuadrature(physics.nd, 3) +elementBoundaryQuadrature = SimplexGaussQuadrature(physics.nd - 1, 3) + +numericalFluxType = VOF.NumericalFlux +subgridError = VOF.SubgridError(coefficients=physics.coefficients, nd=physics.nd) +shockCapturing = VOF.ShockCapturing(coefficients=physics.coefficients, nd=physics.nd, shockCapturingFactor=0.2, + lag=True) + +fullNewtonFlag = True +multilevelNonlinearSolver = Newton +levelNonlinearSolver = Newton +nonlinearSmoother = None +linearSmoother = None +matrix = SparseMatrix +multilevelLinearSolver = KSP_petsc4py +levelLinearSolver = KSP_petsc4py +linear_solver_options_prefix = "dye_" +levelNonlinearSolverConvergenceTest = "rits" +linearSolverConvergenceTest = "r-true" +tolFac = 0.0 +nl_atol_res = 1.0e-10 +linTolFac = 0.0 +l_atol_res = 1.0e-12 +maxNonlinearIts = 20 +maxLineSearches = 0 +parallelPartitioningType = MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 diff --git a/test/rans2p_ib/dye_p.py b/test/rans2p_ib/dye_p.py new file mode 100644 index 0000000000..3f24f902c5 --- /dev/null +++ b/test/rans2p_ib/dye_p.py @@ -0,0 +1,35 @@ +"""A dye in the beam channel (VOF with no level set, carried by RANS2P_IB's velocity): dye 1 enters at x = 0 over +the lower half of the inflow, the channel starts clean.""" +from proteus.default_p import * +from proteus.mprans import VOF +from beam_channel_p import nd, L, T, eps, _inflow, _outflow + +name = "dye" +LevelModelType = VOF.LevelModel +coefficients = VOF.Coefficients(LS_model=None, V_model=0, RD_model=None, ME_model=1, checkMass=False, useMetrics=1.0, + epsFact=0.0, sc_uref=1.0, sc_beta=1.0) + + +def getDBC(x, flag): + if _inflow(x): + return lambda x, t: 1.0 if x[2] < 0.5 * L[2] else 0.0 + + +dirichletConditions = {0: getDBC} + + +def getAFBC(x, flag): + if not (_inflow(x) or _outflow(x)): + return lambda x, t: 0.0 + + +advectiveFluxBoundaryConditions = {0: getAFBC} +diffusiveFluxBoundaryConditions = {0: {}} + + +class Clean(object): + def uOfXT(self, x, t): + return 0.0 + + +initialConditions = {0: Clean()} diff --git a/test/rans2p_ib/test_rans2p_ib.py b/test/rans2p_ib/test_rans2p_ib.py new file mode 100644 index 0000000000..53f566b448 --- /dev/null +++ b/test/rans2p_ib/test_rans2p_ib.py @@ -0,0 +1,104 @@ +"""RANS2P_IB (RANS2P with immersed Euler-Bernoulli beams): the beam channel deck, serial against several ranks. + +Until 2026-11-05 no deck ran RANS2P_IB, and it could not have: its Jacobian raised NameError (``elf.coefficients``), +its Coefficients passed RANS2P the default arguments instead of its own (any rho, nu, g, nd or LS_model was dropped: +every RANS2P_IB run was 3D single-phase water), and the beam loads were integrated over ghost elements too before the +sum over ranks, so on 2 ranks the drag on a rigid beam was 40% too large and a flexible beam bent 70% further. +""" +import os +import shutil +import subprocess +import sys + +import pytest + +from proteus.mprans import RANS2P_IB + +HERE = os.path.dirname(os.path.abspath(__file__)) +HELPER = os.path.join(HERE, "beam_channel_mpi.py") +MPIEXEC = shutil.which("mpiexec") + + +def test_coefficients_pass_their_arguments_to_rans2p(): + c = RANS2P_IB.Coefficients(nu_0=1.0e-3, rho_0=1000.0, nd=2, g=[0.0, -9.81, 0.0], LS_model=2, useMetrics=1.0, + beamLocation=[], beamLength=[], beamRadius=[], EI=[], GJ=[]) + assert (c.nu_0, c.rho_0, c.nd, list(c.g), c.LS_model, c.useMetrics) == (1.0e-3, 1000.0, 2, [0.0, -9.81, 0.0], 2, + 1.0) + + +def _run(launch_policy, tmp_path, ranks, opts): + env = dict(os.environ, OPTS=opts, OMP_NUM_THREADS="1") + env.setdefault("PETSC_OPTIONS", "-mat_partitioning_type current") + out = launch_policy.retrying_run(subprocess.run)([MPIEXEC, "-n", str(ranks), sys.executable, HELPER], + cwd=str(tmp_path), env=env, capture_output=True, text=True, + timeout=600) + assert out.returncode == 0, (out.stdout + out.stderr)[-4000:] + line = [l for l in out.stdout.splitlines() if l.startswith("RESULT")][-1] + return {k: float(v) for k, v in (item.split("=") for item in line.split()[2:])} + + +@pytest.mark.skipif(MPIEXEC is None, reason="mpiexec not available") +@pytest.mark.parametrize("opts", ["beamRigid=True", "beamRigid=False"]) +def test_beam_channel_on_two_ranks_is_the_serial_answer(tmp_path, launch_policy, opts): + """Sums of p, u, v, w (and squares) over owned nodes, the beam drag, the beam loads q1-q3 and the beam's vertex + coordinates. Measured: equal to serial within 4e-16 relative on 2 and 3 ranks (drag1, itself 4e-6, within 4e-14); + before the owned-element fix the drag was 1.12e-2 against 7.99e-3 serial (rigid), and the flexible beam's summed x + coordinates 1.0331 against 1.0195.""" + (tmp_path / "1").mkdir() + (tmp_path / "2").mkdir() + serial = _run(launch_policy, tmp_path / "1", 1, opts) + parallel = _run(launch_policy, tmp_path / "2", 2, opts) + assert serial["drag0"] < -1e-3 and serial["q1"] > 0.5 # the beam is loaded + if opts == "beamRigid=False": + assert serial["xv"] > 1.01 # and bends (1.0 is upright) + worst = {k: abs(serial[k] - parallel[k]) / abs(serial[k]) for k in serial} + assert max(worst.values()) < 1e-11, worst + + +STEP_RETRY = os.path.join(os.path.dirname(HERE), "step_retry_mpi.py") +_LU = " ".join("-%s_ksp_type preonly -%s_pc_type lu -%s_pc_factor_mat_solver_type superlu" % (p, p, p) + for p in ("rans2p", "dye")) + + +def _retry_run(cwd, **overrides): + import numpy + cwd.mkdir() + env = dict(os.environ, OMP_NUM_THREADS="1", DIR="rans2p_ib", SO="beam_channel_so.py", PETSC_OPTIONS=_LU) + env.setdefault("FI_PROVIDER", "tcp") + env["PYTHONPATH"] = os.pathsep.join([p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) + env.update(overrides) + out = subprocess.run([sys.executable, STEP_RETRY], capture_output=True, text=True, timeout=600, env=env, + cwd=str(cwd)) + assert out.returncode == 0, (out.stdout + out.stderr)[-4000:] + return numpy.loadtxt(str(cwd / "u.csv"), delimiter=",", skiprows=1), str(cwd / "accepted.txt") + + +@pytest.mark.parametrize("opts", ["beamRigid=True", "beamRigid=False"]) +def test_a_later_models_failure_is_retried_from_the_step_start(tmp_path, opts): + """The beam channel with a dye after RANS2P_IB (beam_channel_so.py, Sequential_MinModelStep): the dye's third + solve is reported failed, after RANS2P_IB has taken its step, so the system step is cut and retried. The retry + must give the run that takes the same steps without failing (step_retry_mpi.py's REPLAY), in both models. + + Measured, largest relative difference over each model's dofs (rigid and flexible alike): + main: 9.2e-2. The retry skipped RANS2P_IB (Sequential_MinModelStep had already moved its last step's end to + the failed attempt's end), which then ran one step behind the dye and never reached T. + RANS2P_IB stepped again, without its state(): 6.4e-6, the lagged beam drag of the failed attempt. + With RANS2P_IB.Coefficients.state(): 0.0. (Restoring the drag but not the beam: 6.6e-7, and the flexible + beam's deflection at the retried step 1.8e-3 off.) + A flexible beam's deflection history has no row from the failed attempt and is the replay's.""" + import numpy + forced, accepted = _retry_run(tmp_path / "forced", OPTS=opts, FAILMODEL="dye", FAILAT="3") + replay, _ = _retry_run(tmp_path / "replay", OPTS=opts, REPLAY=accepted) + assert int(open(str(tmp_path / "forced" / "retries.txt")).read()) == 1 + times = numpy.loadtxt(accepted) + assert 0.025 in numpy.round(times, 12) and abs(times[-1] - 0.04) < 1e-12, times + assert forced.shape == replay.shape and numpy.abs(forced[:, :4] - replay[:, :4]).max() == 0.0 + for key in numpy.unique(forced[:, 0]): + k = forced[:, 0] == key + worst = numpy.abs(forced[k, 4] - replay[k, 4]).max() / numpy.abs(replay[k, 4]).max() + assert worst < 1e-12, (key, worst) + deflection = tmp_path / "forced" / "deflectionHistory.txt" + if opts == "beamRigid=False": + numpy.testing.assert_array_equal(numpy.loadtxt(str(deflection)), + numpy.loadtxt(str(tmp_path / "replay" / "deflectionHistory.txt"))) + numpy.testing.assert_array_equal(numpy.loadtxt(str(deflection))[:, 0], times) diff --git a/test/shallow_water/dambreak_1d/sw_dam_break_1d_p.py b/test/shallow_water/dambreak_1d/sw_dam_break_1d_p.py index e3c96d0dd7..466e4fba14 100644 --- a/test/shallow_water/dambreak_1d/sw_dam_break_1d_p.py +++ b/test/shallow_water/dambreak_1d/sw_dam_break_1d_p.py @@ -1,7 +1,7 @@ from proteus import * from proteus.default_p import * -""" +r""" This is a 1D dam break (Riemann) problem on the flat domain [0,L] and time interval [0,T]. .. math :: diff --git a/test/shallow_water/dambreak_1d/test_dambreak_1d.py b/test/shallow_water/dambreak_1d/test_dambreak_1d.py new file mode 100644 index 0000000000..8645fe3da1 --- /dev/null +++ b/test/shallow_water/dambreak_1d/test_dambreak_1d.py @@ -0,0 +1,99 @@ +"""test/shallow_water/dambreak_1d runs: RKDG P1 shallow water, Rusanov flux, P1 limiter, a dry-bed dam break between +two walls. + +Until 2026-11-05 it stopped in RusanovNumericalFlux's constructor (a missing periodic argument, since 2014) and then +in the 1D limiter (ctimeIntegration wrappers that did not take TimeIntegration.py's arguments, since 2019). +""" +import inspect +import os +import shutil +import subprocess +import sys + +import h5py +import numpy +import pytest + +from proteus import TimeIntegration + +# Cython 3 binds self to a compiled function stored as a class attribute; TimeIntegration's limiter classes need +# staticmethod for it (cekees/proteus#204, C47). Until that is on the branch the limiter call fails before this deck +# can run, so the test waits for it. +pytestmark = pytest.mark.skipif( + not isinstance(inspect.getattr_static(TimeIntegration.DGlimiterP1Lagrange1d, "applyDGlimitingP1Lagrange1d"), + staticmethod), + reason="needs the staticmethod limiter attributes of cekees/proteus#204") + +_here = os.path.dirname(os.path.abspath(__file__)) +_parun = os.path.join(_here, os.pardir, os.pardir, os.pardir, "scripts", "parun") +_P, _N = "sw_dam_break_1d_p.py", "sw_dam_break_1d_dgpk_n.py" +MPIEXEC = shutil.which("mpiexec") + +# The DG flux needs the face neighbour of every owned element (Transport refuses a node partition), and LU is serial. +_PARALLEL_N = """ +parallelPartitioningType = MeshParallelPartitioningTypes.element +nLayersOfOverlapForParallel = 1 +multilevelLinearSolver = KSP_petsc4py +levelLinearSolver = KSP_petsc4py +linear_solver_options_prefix = "sw_" +""" +_PARALLEL_PETSC = "-mat_partitioning_type current -sw_ksp_type preonly -sw_pc_type bjacobi -sw_sub_pc_type lu" + + +def _run(launch_policy, workdir, ranks=1, parallel=False): + """Run the deck in workdir; parallel: with the element partition and the PETSc solve (any number of ranks).""" + workdir.mkdir() + for name in (_P, _N): + shutil.copy(os.path.join(_here, name), workdir) + env = dict(os.environ) + cmd = [sys.executable, _parun if os.path.isfile(_parun) else shutil.which("parun"), "-l", "1", "-D", "data", _P, _N] + if parallel: + with open(workdir / _N, "a") as f: + f.write(_PARALLEL_N) + env["PETSC_OPTIONS"] = _PARALLEL_PETSC + cmd = [MPIEXEC, "-n", str(ranks)] + cmd + out = launch_policy.retrying_run(subprocess.run)(cmd, cwd=workdir, env=env, capture_output=True, text=True, + timeout=300) + assert out.returncode == 0, (out.stdout + out.stderr)[-4000:] + return _read(workdir / "data" / "sw_dam_break_1d.h5") + + +def _read(path): + """[(element midpoint, x, h, hu) per DG dof, sorted] for each of the 101 outputs.""" + steps = [] + with h5py.File(path, "r") as f: + for i in range(101): + x = f["nodes_dgp1_Lagrange%d" % i][:, 0] + elements = f["elements_dgp1_Lagrange%d" % i][:] + h, hu = f["h_t%d" % i][:], f["hu_t%d" % i][:] + rows = {(0.5 * (x[a] + x[b]), x[d], h[d], hu[d]) for a, b in elements for d in (a, b)} + steps.append(numpy.array(sorted(rows))) + return steps + + +def _mass(rows): + left, right = rows[0::2], rows[1::2] + return numpy.sum(0.5 * (left[:, 2] + right[:, 2]) * (right[:, 1] - left[:, 1])) + + +def test_dambreak_1d_conserves_mass_and_stays_nonnegative(tmp_path, launch_policy): + """Measured: mass 29/6 (the P1 interpolation of the step) at every output to 1.2e-14; h below zero at 3 of 101 + outputs, by at most 5.0e-8 (t = 0.6, the wall cell x = 0 as the front arrives: the P1 limiter does not preserve + positivity); the water settling towards the rest depth 29/60 (h in [0.3645, 0.4985] at T = 20).""" + steps = _run(launch_policy, tmp_path / "serial") + assert all(len(rows) == 60 for rows in steps) + masses = numpy.array([_mass(rows) for rows in steps]) + assert abs(masses[0] - 29.0 / 6.0) < 1e-12 + assert numpy.abs(masses - masses[0]).max() < 1e-12 + assert min(rows[:, 2].min() for rows in steps) > -1e-7 + assert numpy.ptp(steps[-1][:, 2]) < 0.15 + + +@pytest.mark.skipif(MPIEXEC is None, reason="mpiexec not available") +def test_dambreak_1d_on_two_ranks_is_the_serial_answer(tmp_path, launch_policy): + """Measured: identical to the last bit on 2, 3 and 4 ranks (explicit SSPRK, element-local work, a block-diagonal + mass matrix solved exactly per rank).""" + serial = _run(launch_policy, tmp_path / "1", 1, parallel=True) + parallel = _run(launch_policy, tmp_path / "2", 2, parallel=True) + for s, p in zip(serial, parallel): + assert numpy.array_equal(s, p) diff --git a/test/subdomain_mesh_size_mpi.py b/test/subdomain_mesh_size_mpi.py index b804ea70e4..349b41427d 100644 --- a/test/subdomain_mesh_size_mpi.py +++ b/test/subdomain_mesh_size_mpi.py @@ -10,6 +10,7 @@ import os, shutil, sys, runpy import numpy as np from proteus import Comm, MeshTools as MT +import rank_staging # every rank stages the same files: atomically SRC = os.path.join(os.path.dirname(os.path.abspath(__file__)), "LS_with_edgeBased_EV", "NCLS") sys.path.insert(0, os.getcwd()) @@ -19,7 +20,7 @@ text = text.replace("parallel = False", "parallel = True") if name == "ncls_n.py": text += "\nnl_atol_res = 1.0e-13; l_atol_res = 1.0e-15; linTolFac = 0.0; tolFac = 0.0\n" - open(name, "w").write(text) + rank_staging.write(name, text) _partition = MT.Mesh.partitionMesh def partitionMesh(self, nLayersOfOverlap=1, parallelPartitioningType=MT.MeshParallelPartitioningTypes.node, diff --git a/test/sw2dcv_partition_mpi.py b/test/sw2dcv_partition_mpi.py index ddc5c21dc7..a6a248d8c6 100644 --- a/test/sw2dcv_partition_mpi.py +++ b/test/sw2dcv_partition_mpi.py @@ -11,11 +11,13 @@ import os, shutil, sys, runpy import numpy as np from proteus import Comm +import rank_staging -shutil.copy(os.path.join(os.path.dirname(os.path.abspath(__file__)), "SWFlow", "circular_damBreak.py"), - "circular_damBreak.py") -open("__init__.py", "w").close() -open("circular_damBreak.py", "a").write("\nmySWFlowProblem.swe_parameters['LUMPED_MASS_MATRIX'] = %d\n" % (0 if os.environ.get("CONSISTENT") else 1)) +# every rank stages the same files: atomically (rank_staging) +rank_staging.write("circular_damBreak.py", open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "SWFlow", + "circular_damBreak.py")).read() + + "\nmySWFlowProblem.swe_parameters['LUMPED_MASS_MATRIX'] = %d\n" % (0 if os.environ.get("CONSISTENT") else 1)) +rank_staging.write("__init__.py", "") from proteus.mprans import SW2DCV if os.environ.get("NOSCAT"): diff --git a/test/tadr_partition_mpi.py b/test/tadr_partition_mpi.py index a064143c07..89133e4c68 100644 --- a/test/tadr_partition_mpi.py +++ b/test/tadr_partition_mpi.py @@ -13,6 +13,7 @@ import os, shutil, sys, runpy import numpy as np from proteus import Comm, MeshTools as MT +import rank_staging # every rank stages the same files: atomically SRC = os.path.join(os.path.dirname(os.path.abspath(__file__)), "TADR") sys.path.insert(0, os.getcwd()) @@ -27,7 +28,7 @@ text += ("\nif ct.STABILIZATION_TYPE == 5:\n levelNonlinearSolver = Newton\n fullNewtonFlag = True\n" " updateJacobian = True\n timeIntegration = BackwardEuler_cfl\n nl_atol_res = 1.0e-10\n" " l_atol_res = 1.0e-12\n maxNonlinearIts = 100\n") - open(name, "w").write(text) + rank_staging.write(name, text) _partition = MT.Mesh.partitionMesh def partitionMesh(self, nLayersOfOverlap=1, parallelPartitioningType=MT.MeshParallelPartitioningTypes.node, diff --git a/test/test_av_partition_independence.py b/test/test_av_partition_independence.py index cd564221cc..34197b5b51 100644 --- a/test/test_av_partition_independence.py +++ b/test/test_av_partition_independence.py @@ -6,15 +6,19 @@ t=0.1 with the cut at x=1.6. The indicator is now computed in Python on the owned dofs and scattered to the ghosts; the control runs without that scatter and differs. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "av_partition_mpi.py") @@ -39,8 +43,7 @@ def _pressure(ranks, cwd, **overrides): [p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "pressure.csv"), delimiter=",", skiprows=1) diff --git a/test/test_gn_sw2dcv_partition_independence.py b/test/test_gn_sw2dcv_partition_independence.py index a0bec2f408..507706fa66 100644 --- a/test/test_gn_sw2dcv_partition_independence.py +++ b/test/test_gn_sw2dcv_partition_independence.py @@ -9,15 +9,19 @@ The time step is chosen from the edge based cfl of the owned dofs only, since the cfl of a ghost dof comes from a truncated row and can exceed the global maximum; the ghost-cfl control puts garbage there, which the step must ignore. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "gn_sw2dcv_partition_mpi.py") @@ -42,8 +46,7 @@ def _field(ranks, cwd, **overrides): [p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) diff --git a/test/test_launch_output.py b/test/test_launch_output.py new file mode 100644 index 0000000000..d556002ef7 --- /dev/null +++ b/test/test_launch_output.py @@ -0,0 +1,46 @@ +"""launch_output: a failed helper launch keeps its whole output (the partition-independence tests kept the last 2000 +characters of each stream, which lost a 2-rank segfault's backtrace, and nothing at all of a hang).""" +import importlib.util +import os +import sys + +import pytest + +HERE = os.path.dirname(os.path.abspath(__file__)) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) + +LONG = "frame %06d\n" # 13 characters a line: 20000 lines is 260000 characters, far past any tail + + +def test_a_crash_keeps_its_whole_output_and_names_the_signal(tmp_path): + code = ("import os, signal, sys\n" + "sys.stderr.write(''.join(%r %% i for i in range(20000)))\nsys.stderr.flush()\n" + "print('stdout line', flush=True)\nos.kill(os.getpid(), signal.SIGSEGV)\n" % LONG) + with pytest.raises(AssertionError) as e: + launch_output.run([sys.executable, "-c", code], cwd=tmp_path, env=dict(os.environ), timeout=60) + message = str(e.value) + assert "SIGSEGV" in message and "frame 019999" in message and "stdout line" in message + assert (tmp_path / "launch.stderr.txt").read_text().startswith("frame 000000\n") + assert "frame 019999\n" in (tmp_path / "launch.stderr.txt").read_text() # then faulthandler's traceback, if on + assert str(tmp_path / "launch.stderr.txt") in message + + +def test_a_hang_keeps_what_it_printed(tmp_path): + code = "import time\nprint('started the solve', flush=True)\ntime.sleep(60)\n" + with pytest.raises(AssertionError) as e: + launch_output.run([sys.executable, "-c", code], cwd=tmp_path, env=dict(os.environ), timeout=3) + assert "timed out after 3 s" in str(e.value) and "started the solve" in str(e.value) + assert "started the solve" in (tmp_path / "launch.stdout.txt").read_text() + + +def test_a_success_returns_the_completed_process(tmp_path): + out = launch_output.run([sys.executable, "-c", "print('ok')"], cwd=tmp_path, env=dict(os.environ), timeout=60) + assert out.stdout == "ok\n" and not (tmp_path / "launch.stdout.txt").exists() + + +@pytest.mark.parametrize("code, status", [(1, "exit status 1"), (139, "exit status 139 (SIGSEGV?)"), + (-11, "killed by SIGSEGV")]) +def test_status(code, status): + assert launch_output._status(code) == status diff --git a/test/test_mcorr_partition_independence.py b/test/test_mcorr_partition_independence.py index d03f255230..dfd2328d5d 100644 --- a/test/test_mcorr_partition_independence.py +++ b/test/test_mcorr_partition_independence.py @@ -13,15 +13,19 @@ The run couples NCLS, VOF and MCorr. 2 ranks match serial to 1e-8 in all three; the control without the MCorr scatters differs; the KSP path agrees with the LU path in serial. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "mcorr_partition_mpi.py") @@ -47,8 +51,7 @@ def _field(ranks, cwd, **overrides): env["PETSC_OPTIONS"] = "-ncls_ksp_rtol 1e-14 -ncls_ksp_atol 1e-40 -vof_ksp_rtol 1e-14 -vof_ksp_atol 1e-40 -mcorr_ksp_rtol 1e-14 -mcorr_ksp_atol 1e-40 -mcorr_ksp_max_it 5000" env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) result = numpy.load(os.path.join(str(cwd), "final.npz")) return {key: result[key][:, 2] for key in result.files} diff --git a/test/test_ncls_edge_partition_independence.py b/test/test_ncls_edge_partition_independence.py index 4bc1ad60ff..5817393a31 100644 --- a/test/test_ncls_edge_partition_independence.py +++ b/test/test_ncls_edge_partition_independence.py @@ -14,15 +14,19 @@ 9e-16 on macOS arm64, the same on every run of a host, so it is round-off through that branch, not a partition effect. The control runs without the scatter of psi and of the entropy residual and differs. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "ncls_edge_partition_mpi.py") @@ -48,8 +52,7 @@ def _field(ranks, cwd, **overrides): env["PETSC_OPTIONS"] = "-ncls_ksp_rtol 1e-14 -ncls_ksp_atol 1e-40 -ncls_ksp_max_it 5000" env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) diff --git a/test/test_numerical_flux_constructors.py b/test/test_numerical_flux_constructors.py new file mode 100644 index 0000000000..62c17f8e39 --- /dev/null +++ b/test/test_numerical_flux_constructors.py @@ -0,0 +1,114 @@ +"""Every numerical flux in proteus/NumericalFlux.py can be built the two ways OneLevelTransport builds one. + +Transport.py (and the mprans models) call ``numericalFluxType(vt, getPointwiseBoundaryConditions, +getAdvectiveFluxBoundaryConditions, getDiffusiveFluxBoundaryConditions)``, and append +``options.periodicDirichletConditions`` when the p module sets periodic conditions. Until 2026-11-05 +RusanovNumericalFlux_Diagonal and ConvexOneSonicPointNumericalFlux required the periodic argument, so the ordinary +call raised TypeError, and RusanovNumericalFlux called RusanovNumericalFlux_Diagonal.__init__ without it, so every +deck using RusanovNumericalFlux (test/shallow_water/dambreak_1d) failed in its constructor (since b8912a02, 2014). +""" +import ast +import inspect + +from proteus import NumericalFlux + +_FLUXES = {name: cls for name, cls in vars(NumericalFlux).items() + if inspect.isclass(cls) and cls.__module__ == NumericalFlux.__name__ and issubclass(cls, NumericalFlux.NF_base)} + +# These take no periodic argument, so a periodic deck using them stops with TypeError in the constructor. That is the +# loud refusal, kept: most of them have no periodic face pairing in their exterior fluxes, so accepting the argument +# would only hide it. A class leaving this set must also pair periodic faces (NF_base.periodicBoundaryConditionsDictList). +NO_PERIODIC_ARGUMENT = { + "Advection_DiagonalUpwind_Diffusion_IIPG", "Advection_DiagonalUpwind_Diffusion_LDG", + "Advection_DiagonalUpwind_Diffusion_NIPG", "Advection_DiagonalUpwind_Diffusion_NIPG_exterior", + "Advection_DiagonalUpwind_Diffusion_SIPG", "Advection_DiagonalUpwind_Diffusion_SIPG_exterior", + "Advection_Diagonal_average", "ConstantAdvection_Diffusion_IIPG_exterior", + "ConstantAdvection_Diffusion_NIPG_exterior", "ConstantAdvection_Diffusion_SIPG_exterior", + "ConstantAdvection_exterior", "Curvature_exterior", "DarcyFCFF_IIPG_exterior", "DarcyFCPP_IIPG_exterior", + "DarcyFC_IIPG_exterior", "DarcySplitPressure_IIPG_exterior", "Diffusion_LDG", "MixedDarcy_exterior", + "Richards_IIPG_exterior", "Richards_SIPG_exterior", "RusanovLDG", "RusanovNumericalFlux_Diagonal_Diffusion_IIPG", + "StokesP_Advection_DiagonalUpwind_Diffusion_IIPG_exterior", +} + + +def _binds(cls, nargs): + try: + inspect.signature(cls.__init__).bind(*range(1 + nargs)) # self + arguments + return True + except TypeError: + return False + + +def test_every_flux_takes_the_ordinary_call(): + refused = sorted(name for name, cls in _FLUXES.items() if not _binds(cls, 4)) + assert refused == [], refused + + +def test_the_fluxes_that_take_no_periodic_argument_are_the_known_set(): + refused = {name for name, cls in _FLUXES.items() if not _binds(cls, 5)} + assert refused == NO_PERIODIC_ARGUMENT, (sorted(refused - NO_PERIODIC_ARGUMENT), + sorted(NO_PERIODIC_ARGUMENT - refused)) + + +def _base_init_calls(): + """(class, line, base, parameters of the calling __init__, call) for every Base.__init__(self, ...) in a flux's + __init__.""" + tree = ast.parse(inspect.getsource(NumericalFlux)) + for cdef in (node for node in tree.body if isinstance(node, ast.ClassDef) and node.name in _FLUXES): + for fdef in cdef.body: + if not (isinstance(fdef, ast.FunctionDef) and fdef.name == "__init__"): + continue + params = [a.arg for a in fdef.args.args] + for call in ast.walk(fdef): + if (isinstance(call, ast.Call) and isinstance(call.func, ast.Attribute) + and call.func.attr == "__init__" and isinstance(call.func.value, ast.Name) + and call.func.value.id in _FLUXES): + yield cdef.name, call.lineno, _FLUXES[call.func.value.id], params, call + + +def test_every_base_constructor_call_binds_and_forwards_the_periodic_conditions(): + wrong = [] + for name, line, base, params, call in _base_init_calls(): + names = [arg.id if isinstance(arg, ast.Name) else ast.unparse(arg) for arg in call.args] + try: + inspect.signature(base.__init__).bind(*names, **{kw.arg: kw.arg for kw in call.keywords}) + except TypeError as e: + wrong.append("%s:%d %s.__init__: %s" % (name, line, base.__name__, e)) + continue + periodic = [p for p in params if p.startswith("getPeriodic")] + if periodic and periodic[0] not in names + [kw.arg for kw in call.keywords]: + wrong.append("%s:%d drops %s" % (name, line, periodic[0])) + assert wrong == [], "\n".join(wrong) + + +# Seven fluxes accept the periodic argument and pair no periodic faces. Measured on a C0 P1 periodic Poisson problem +# (NumericalFlux.PERIODIC_SEAM_REFUSED says how): the four diffusion ones are wrong with periodic conditions and refuse +# them; the other three are not changed by them and keep accepting them: +# StrongDirichlet: 4.6e-2 without the periodic pair, 4.7e-2 with it, 1 Newton iteration (ADRProblem's periodic +# problems use it through StrongDirichletFactory: test_adr_problem's periodic Poisson, channel, several ranks); +# DoNothing: no boundary term at all, the right one for a model without boundary conditions on a periodic C0 space +# (the MCorr and RDLS models of test/periodic and test/hotstart_measure use it); +# NoFlux: the same answer with and without the periodic pair (it imposes no Dirichlet condition either way). +PERIODIC_ACCEPTED_WITHOUT_PAIRING = {"StrongDirichlet", "DoNothing", "NoFlux"} + + +def _periodic(): + return {0: lambda x, flag: None} + + +def test_the_diffusion_exterior_fluxes_refuse_periodic_conditions(): + import pytest + for name in sorted(NumericalFlux.PERIODIC_SEAM_REFUSED): + with pytest.raises(ValueError, match="periodic"): + _FLUXES[name](None, {}, {}, {}, _periodic()) + + +def test_the_other_unpaired_fluxes_still_take_periodic_conditions(): + # they get past the refusal and fail later on the missing transport (None), or not at all (DoNothing) + for name in sorted(PERIODIC_ACCEPTED_WITHOUT_PAIRING): + try: + _FLUXES[name](None, {}, {}, {}, _periodic()) + except ValueError as e: + raise AssertionError("%s refused: %s" % (name, e)) + except (AttributeError, TypeError): + pass diff --git a/test/test_rank_staging.py b/test/test_rank_staging.py new file mode 100644 index 0000000000..b67bc18f76 --- /dev/null +++ b/test/test_rank_staging.py @@ -0,0 +1,88 @@ +"""rank_staging: the ``*_mpi.py`` helpers stage a deck's files from every rank without a moment at which another rank +can import part of one (test_mcorr_partition_independence hung the GitHub aarch64 leg on 2026-11-06: one rank +star-imported a module another rank had just truncated, NameError, and the other rank waited 2809 s).""" +import ast +import glob +import importlib.util +import os +import threading +import time + +HERE = os.path.dirname(os.path.abspath(__file__)) +# by path: CI's --import-mode=importlib does not put test/ on sys.path (the helpers, run as scripts, have it there) +_spec = importlib.util.spec_from_file_location("rank_staging", os.path.join(HERE, "rank_staging.py")) +rank_staging = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(rank_staging) +A, B = "x = 1\n" * 20000, "y = 2\n" * 30000 + + +def _read_while_writing(tmp_path, writer): + path = str(tmp_path / "m.py") + rank_staging.write(path, A) + seen, stop = set(), threading.Event() + + def reader(): + while not stop.is_set(): + with open(path) as f: + seen.add(len(f.read())) + t = threading.Thread(target=reader) + t.start() + for text in (B, A, B, A): + writer(path, text) + stop.set() + t.join() + return seen + + +def test_a_reader_sees_only_whole_files(tmp_path): + seen = _read_while_writing(tmp_path, lambda path, text: rank_staging.write(path, text, _pause=0.02)) + assert seen <= {len(A), len(B)}, sorted(seen) + + +def test_the_in_place_write_it_replaces_is_seen_part_written(tmp_path): + """The control: open(path, "w") then write, held 20 ms in between as a descheduled rank would be.""" + def in_place(path, text): + with open(path, "w") as f: + time.sleep(0.02) + f.write(text) + seen = _read_while_writing(tmp_path, in_place) + assert seen - {len(A), len(B)}, sorted(seen) + + +def test_copy_into_a_directory(tmp_path): + src = tmp_path / "src.py" + src.write_text(A) + (tmp_path / "d").mkdir() + assert rank_staging.copy(str(src), str(tmp_path / "d")) == str(tmp_path / "d" / "src.py") + assert (tmp_path / "d" / "src.py").read_text() == A + assert os.listdir(str(tmp_path / "d")) == ["src.py"] + + +# serial only (test_step_retry runs it without mpiexec) +SERIAL_HELPERS = {"step_retry_mpi.py"} +#: what a helper writes as its result (per rank or from rank 0), not a file another rank imports +OUTPUTS = (".json", ".csv", ".txt", ".npz", ".npy", ".h5") + + +def test_no_mpi_helper_stages_files_in_place(): + """Every rank runs a helper's module level: none may copy, write or append a deck file there except through + rank_staging (results, by extension, writes under ``if rank ...`` and inside functions are not looked at).""" + wrong = [] + for path in sorted(glob.glob(os.path.join(HERE, "**", "*_mpi.py"), recursive=True)): + if os.path.basename(path) in SERIAL_HELPERS: + continue + tree = ast.parse(open(path).read()) + stack = [n for n in tree.body if not isinstance(n, (ast.FunctionDef, ast.ClassDef, ast.If))] + stack += [n for n in tree.body if isinstance(n, ast.If) and "rank" not in ast.unparse(n.test)] + for top in stack: + for node in ast.walk(top): + if not isinstance(node, ast.Call): + continue + name = ast.unparse(node.func) + mode = [a for a in node.args[1:2] if isinstance(a, ast.Constant)] + target = ast.unparse(node.args[0]) if node.args else "" + if name in ("shutil.copy", "shutil.copyfile", "shutil.copy2") or ( + name == "open" and mode and str(mode[0].value)[:1] in "wa" + and not any(ext in target for ext in OUTPUTS)): + wrong.append("%s:%d %s" % (os.path.relpath(path, HERE), node.lineno, ast.unparse(node)[:80])) + assert wrong == [], "\n".join(wrong) diff --git a/test/test_step_retry.py b/test/test_step_retry.py index 15a4c5f1c6..b9e7366230 100644 --- a/test/test_step_retry.py +++ b/test/test_step_retry.py @@ -56,6 +56,8 @@ def _worst(a, b): return worst +_LU = "-rans2p_ksp_type preonly -rans2p_pc_type lu -rans2p_pc_factor_mat_solver_type superlu" + CASES = { "fillingTank-clsvof": dict(DECK="fillingTank.py", ATOL=1e-9, FAILMODEL="clsvof", FAILAT=10, BIGAUGES="1,2,3,4"), "fillingTank-pressureIncrement": dict(DECK="fillingTank.py", ATOL=1e-9, FAILMODEL="pressureIncrement", FAILAT=10), @@ -85,6 +87,13 @@ def _worst(a, b): ATOL=1e-10, MAXITS=1000, FAILMODEL="kappa", FAILAT=5), "damBreakLevelSet-rans2p-clsvof": dict(DECK="damBreakLevelSet.py", OPTS="he=0.1 interface=1", FT=0.1, ATOL=1e-10, MAXITS=1000, FAILMODEL="clsvof", FAILAT=6), + # RANS2P_IB with a rigid and a flexible beam (the beam bends in postStep, after an accepted step only): measured + # 0.0, the failed attempt leaving no row in the deflection history + "beam_channel-rans2p_ib": dict(DIR="rans2p_ib", SO="beam_channel_p.py beam_channel_n.py", OPTS="beamRigid=True", + PETSC_OPTIONS=_LU, FAILMODEL="beam_channel", FAILAT=3), + "beam_channel-rans2p_ib-flexible": dict(DIR="rans2p_ib", SO="beam_channel_p.py beam_channel_n.py", + OPTS="beamRigid=False", PETSC_OPTIONS=_LU, FAILMODEL="beam_channel", + FAILAT=3), "circular_damBreak-sw2dcv": dict(DIR="SWFlow", SO="--SWEs circular_damBreak.py", OPTS="sw_model=0 final_time=0.1 dt_output=0.1", FAILMODEL="sw", FAILAT=6), "circular_damBreak-gn_sw2dcv": dict(DIR="SWFlow", SO="--SWEs circular_damBreak.py", @@ -104,6 +113,12 @@ def test_a_retried_step_is_the_step_taken_without_failing(tmp_path, case): replay, _ = _run(tmp_path / "replay", REPLAY=accepted, **_replay_options(options)) assert int(open(os.path.join(str(tmp_path / "forced"), "retries.txt")).read()) >= 1, "the failure was not retried" assert _worst(forced, replay) < 1e-10 + # the failed attempt leaves no row in the beam deflection history (RANS2P_IB, a flexible beam) + deflection = os.path.join(str(tmp_path / "forced"), "deflectionHistory.txt") + if os.path.exists(deflection): + numpy.testing.assert_array_equal(numpy.loadtxt(deflection), + numpy.loadtxt(os.path.join(str(tmp_path / "replay"), "deflectionHistory.txt"))) + numpy.testing.assert_array_equal(numpy.loadtxt(deflection)[:, 0], numpy.loadtxt(accepted)) # the failed attempt leaves no row in the histories RANS2P/RANS3PF write history = os.path.join(str(tmp_path / "forced"), "timeHistory.txt") if os.path.exists(history): diff --git a/test/test_subdomain_mesh_size.py b/test/test_subdomain_mesh_size.py index 7864ed63d6..869debdc1f 100644 --- a/test/test_subdomain_mesh_size.py +++ b/test/test_subdomain_mesh_size.py @@ -5,15 +5,19 @@ rank 1 had h = 0.035294 against 0.035538 serially, and 2 ranks differed from serial by 5.1e-7 at t=0.08 (3 ranks 5.1e-7). With the sizes reduced over the ranks: 4.2e-11 and 1.6e-11. The control keeps each rank's own sizes. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "subdomain_mesh_size_mpi.py") @@ -39,8 +43,7 @@ def _run(ranks, cwd, **overrides): env["PETSC_OPTIONS"] = "-ncls_ksp_rtol 1e-14 -ncls_ksp_atol 1e-40 -ncls_ksp_max_it 5000" env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) u = numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) sizes = numpy.atleast_2d(numpy.loadtxt(os.path.join(str(cwd), "h.csv"), delimiter=",", skiprows=1)) return u, sizes diff --git a/test/test_sw2dcv_partition_independence.py b/test/test_sw2dcv_partition_independence.py index be7098ad53..f1283a6d05 100644 --- a/test/test_sw2dcv_partition_independence.py +++ b/test/test_sw2dcv_partition_independence.py @@ -10,15 +10,19 @@ The time step is chosen from the edge based cfl of the owned dofs only, since the cfl of a ghost dof comes from a truncated row and can exceed the global maximum; the ghost-cfl control puts garbage there, which the step must ignore. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "sw2dcv_partition_mpi.py") @@ -43,8 +47,7 @@ def _field(ranks, cwd, **overrides): [p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) diff --git a/test/test_tadr_partition_independence.py b/test/test_tadr_partition_independence.py index 6e9291bd62..18ecb88ea6 100644 --- a/test/test_tadr_partition_independence.py +++ b/test/test_tadr_partition_independence.py @@ -17,15 +17,19 @@ and the limited solution, 0.49 and 0.8% of the mass (the defect reported on a branch without them: 0.077 and 3.8% on 4 ranks). """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "tadr_partition_mpi.py") @@ -51,8 +55,7 @@ def _field(ranks, cwd, **overrides): env["PETSC_OPTIONS"] = "-tadr_ksp_rtol 1e-14 -tadr_ksp_atol 1e-40 -tadr_ksp_max_it 5000" env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) diff --git a/test/test_time_integration_wrappers.py b/test/test_time_integration_wrappers.py new file mode 100644 index 0000000000..73d4ad6f11 --- /dev/null +++ b/test/test_time_integration_wrappers.py @@ -0,0 +1,79 @@ +"""The ctimeIntegration wrappers take the arguments TimeIntegration.py passes. + +The 2019 Cython rewrite of the C module (#879) gave seven wrappers new signatures (the sizes as leading arguments, +the limiter flag first, psiTCtteDT writing dt into an argument instead of returning it) while TimeIntegration.py kept +calling them the old way, so every 1D DG limiter (DGlimiterP1Lagrange1d, its _Sw and P2/P3 variants, the vacuum +tolerance ones), the Durlofsky limiters and pseudo-transient continuation (PsiTCtte) raised TypeError on their first +call. test/shallow_water/dambreak_1d stopped there. +""" +import ast +import os +import re + +import numpy + +from proteus import ctimeIntegration + +ROOT = os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir, "proteus") + + +def _signatures(): + """name -> (number of parameters, number with defaults) for every def in ctimeIntegration.pyx.""" + pyx = open(os.path.join(ROOT, "ctimeIntegration.pyx")).read() + sigs = {} + for m in re.finditer(r"^def (\w+)\((.*?)\):", pyx, re.S | re.M): + params = [p for p in m.group(2).replace("\n", " ").split(",") if p.strip()] + sigs[m.group(1)] = (len(params), sum("=" in p for p in params)) + return sigs + + +def test_every_call_from_time_integration_binds(): + sigs = _signatures() + tree = ast.parse(open(os.path.join(ROOT, "TimeIntegration.py")).read()) + calls, wrong = 0, [] + for call in ast.walk(tree): + if isinstance(call, ast.Call) and isinstance(call.func, ast.Attribute) and call.func.attr in sigs: + calls += 1 + nparams, ndefaults = sigs[call.func.attr] + if not nparams - ndefaults <= len(call.args) <= nparams or call.keywords: + wrong.append("TimeIntegration.py:%d %s: %d arguments, the wrapper takes %d to %d" + % (call.lineno, call.func.attr, len(call.args), nparams - ndefaults, nparams)) + assert calls >= 15 + assert wrong == [], "\n".join(wrong) + + +def _line(u): + """Four P1 DG elements on [0, 4] with dof values u (two per element).""" + nodes = numpy.zeros((5, 3)) + nodes[:, 0] = numpy.arange(5.0) + elementNodes = numpy.array([[0, 1], [1, 2], [2, 3], [3, 4]], "i") + neighbours = numpy.array([[1, -1], [2, 0], [3, 1], [-1, 2]], "i") # across the face opposite local node 0, 1 + barycentres = numpy.zeros((4, 3)) + barycentres[:, 0] = numpy.arange(4) + 0.5 + l2g = numpy.arange(8, dtype="i").reshape(4, 2) + return elementNodes, neighbours, nodes, barycentres, l2g, numpy.zeros(4, "i"), numpy.array(u, dtype=float) + + +def test_p1_1d_limiter_keeps_linear_data_and_flattens_extrema(): + for flag in (0, 1, 2): + *mesh, tag, u = _line([0, 1, 1, 2, 2, 3, 3, 4]) + out = numpy.zeros(8) + ctimeIntegration.applyDGlimitingP1Lagrange1d(*mesh, tag, u, out, flag) + assert numpy.array_equal(out, u) and tag.tolist() == [1, 1, 1, 1] + # element 1 (mean 3) is above both neighbours' means (0.5, 2.5): flattened; element 2 (mean 2.5, slope 1) + # sees a falling left difference: flattened + *mesh, tag, u = _line([0, 1, 1, 5, 2, 3, 3, 4]) + ctimeIntegration.applyDGlimitingP1Lagrange1d(*mesh, tag, u, out, flag) + assert out.tolist() == [0, 1, 3, 3, 2.5, 2.5, 3, 4] and tag.tolist() == [1, 0, 0, 1] + *mesh, tag, u = _line([0, 1, 1, 5, 2, 3, 3, 4]) + out = numpy.zeros(8) + ctimeIntegration.applyDGlimitingP1Lagrange1d_withVacuumTol(0, 0.0, *mesh, tag, u, out) + assert out.tolist() == [0, 1, 3, 3, 2.5, 2.5, 3, 4] + + +def test_psi_tc_returns_the_step_from_every_point(): + yn, ypn = numpy.zeros((2, 3)), numpy.ones((2, 3)) + ypnm1 = ypn.copy() + ypnm1[1, 2] -= 0.4 # the largest change at the last of the 6 points + dt = ctimeIntegration.psiTCtteDT(1.0e-2, 0.1, 0.2, yn, ypn, ypnm1) + assert dt == numpy.sqrt(1.0e-2 / (0.2 + 1.0e-8)) * 0.9 diff --git a/test/test_twophaseflow_partition_independence.py b/test/test_twophaseflow_partition_independence.py index 79e8d57924..a191341342 100644 --- a/test/test_twophaseflow_partition_independence.py +++ b/test/test_twophaseflow_partition_independence.py @@ -42,15 +42,19 @@ so recovery is partition-correct; but the retried serial run differs from the unforced one by 16% in velocity at t=0.05, which is a time-stepping question, not a partition one. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "twophaseflow_partition_mpi.py") @@ -78,8 +82,7 @@ def _field(ranks, cwd, **overrides): [p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) names = dict(line.strip().split(",") for line in open(os.path.join(str(cwd), "models.csv"))) u = numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) return u, {int(k): v for k, v in names.items()} diff --git a/test/test_vof_partition_independence.py b/test/test_vof_partition_independence.py index e9f4174cf2..d346c8a776 100644 --- a/test/test_vof_partition_independence.py +++ b/test/test_vof_partition_independence.py @@ -10,15 +10,19 @@ the cfl of a ghost is built from a truncated row and can exceed the global maximum; the ghost-cfl control puts garbage there, which the step must ignore. """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "vof_fct_partition_mpi.py") @@ -43,8 +47,7 @@ def _field(ranks, cwd, **overrides): [p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) diff --git a/test/test_vos_partition_independence.py b/test/test_vos_partition_independence.py index c0622b2fc0..2dda713f8d 100644 --- a/test/test_vos_partition_independence.py +++ b/test/test_vos_partition_independence.py @@ -7,15 +7,19 @@ that made the unit-square run differ by 1e-4): the ghost-cfl control puts garbage there, which the step must ignore. (Scattering the limiter coefficients and the limited solution as well, as VOF needs, changed nothing here.) """ +import importlib.util import os import shutil -import subprocess import sys import numpy import pytest HERE = os.path.dirname(os.path.abspath(__file__)) +# a failed launch keeps its whole output (loaded by path: --import-mode=importlib) +_spec = importlib.util.spec_from_file_location("launch_output", os.path.join(HERE, "launch_output.py")) +launch_output = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(launch_output) HELPER = os.path.join(HERE, "vos_partition_mpi.py") @@ -40,8 +44,7 @@ def _field(ranks, cwd, **overrides): [p for p in sys.path if p] + [env.get("PYTHONPATH", "")]).strip(os.pathsep) env.update({k: str(v) for k, v in overrides.items()}) argv = [sys.executable, HELPER] if ranks == 1 else [MPIEXEC, "-n", str(ranks), sys.executable, HELPER] - out = subprocess.run(argv, capture_output=True, text=True, timeout=1800, env=env, cwd=str(cwd)) - assert out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:] + out = launch_output.run(argv, cwd=cwd, env=env, timeout=1800) return numpy.loadtxt(os.path.join(str(cwd), "u.csv"), delimiter=",", skiprows=1) diff --git a/test/twophaseflow_partition_mpi.py b/test/twophaseflow_partition_mpi.py index cd35db44a8..fed43721e8 100644 --- a/test/twophaseflow_partition_mpi.py +++ b/test/twophaseflow_partition_mpi.py @@ -12,13 +12,14 @@ """ import glob, os, runpy, shutil, sys import numpy as np +import rank_staging # every rank stages the same files: atomically SO = os.environ.get("SO") SRC = os.path.join(os.path.dirname(os.path.abspath(__file__)), os.environ.get("DIR", "TwoPhaseFlow")) DECK = os.environ.get("DECK", "damBreak.py") for f in (glob.glob(os.path.join(SRC, "*")) if SO else [os.path.join(SRC, DECK)] + glob.glob(os.path.join(SRC, "mesh*"))): if os.path.isfile(f): - shutil.copy(f, ".") + rank_staging.copy(f, ".") from proteus import Comm from proteus import NumericalSolution as NS diff --git a/test/vos_partition_mpi.py b/test/vos_partition_mpi.py index 2315846a08..cc2c8c9ba8 100644 --- a/test/vos_partition_mpi.py +++ b/test/vos_partition_mpi.py @@ -10,6 +10,7 @@ import os, shutil, sys, runpy import numpy as np from proteus import Comm, MeshTools as MT +import rank_staging # every rank stages the same files: atomically SRC = os.path.join(os.path.dirname(os.path.abspath(__file__)), "VOS") for name in ("vos.py", "vos_p.py", "vos_n.py"): @@ -18,7 +19,7 @@ text = text.replace("ct.T=1.0\n ct.nDTout=20", "ct.T=0.4\n ct.nDTout=8") text = text.replace("unstructured=False", "unstructured=True").replace("nny=int((nnx-1)/10+1)", "nny=nnx") text = text.replace("0.1 if ct.problem==0 else 1.0", "1.0") - open(name, "w").write(text) + rank_staging.write(name, text) _partition = MT.Mesh.partitionMesh def partitionMesh(self, nLayersOfOverlap=1, parallelPartitioningType=MT.MeshParallelPartitioningTypes.node,