From 7d95b91377cf344dfc4acfbe7970d2fd181eca6e Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 21:58:05 +0100 Subject: [PATCH 01/10] NumericalFlux: Rusanov and convex-sonic fluxes take the ordinary 4-argument call OneLevelTransport builds a flux as numericalFluxType(vt, dirichlet, advective, diffusive) and appends the periodic conditions only when the p module sets them. RusanovNumericalFlux_Diagonal and ConvexOneSonicPointNumericalFlux required the periodic argument, so the ordinary call raised TypeError, and RusanovNumericalFlux called RusanovNumericalFlux_Diagonal.__init__ with four arguments, so it failed even when given five. Every deck using RusanovNumericalFlux stopped in its constructor (test/shallow_water/dambreak_1d at NumericalFlux.py:4726), since b8912a02 (2014) added the periodic argument. The periodic argument now defaults to None in all three and RusanovNumericalFlux forwards it; the Rusanov exterior flux already pairs periodic faces. Sweep of the 50 NF_base subclasses (test/test_numerical_flux_constructors.py, which imports the module and binds every constructor and every Base.__init__(self, ...) call): ordinary call refused: main 2 (ConvexOneSonicPoint, Rusanov_Diagonal) -> 0 base __init__ call unbound: main 1 (RusanovNumericalFlux:4726) -> 0 periodic argument dropped by a subclass that takes one: 0 on main 23 classes take no periodic argument at all; a periodic deck with one of them stops with TypeError. The test pins that set: most of them have no periodic face pairing, so accepting the argument would hide the gap, not close it. Seven more accept it and pair no faces (Diffusion_IIPG/SIPG_exterior, Stokes_*_IIPG/SIPG_exterior, NoFlux, StrongDirichlet, DoNothing): noted in the round handoff, not changed. dambreak_1d after this commit gets past the flux and stops in the DG limiter (ctimeIntegration, next commit). Tests: test_numerical_flux_constructors (3, all 3 fail on main), test_numerical_flux_arguments. Co-Authored-By: Claude Opus 5.5 --- proteus/NumericalFlux.py | 10 +-- test/test_numerical_flux_constructors.py | 81 ++++++++++++++++++++++++ 2 files changed, 87 insertions(+), 4 deletions(-) create mode 100644 test/test_numerical_flux_constructors.py diff --git a/proteus/NumericalFlux.py b/proteus/NumericalFlux.py index 2ed93b556d..1a942cc989 100644 --- a/proteus/NumericalFlux.py +++ b/proteus/NumericalFlux.py @@ -3026,7 +3026,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 +3270,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 +4722,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/test/test_numerical_flux_constructors.py b/test/test_numerical_flux_constructors.py new file mode 100644 index 0000000000..1fddda17ab --- /dev/null +++ b/test/test_numerical_flux_constructors.py @@ -0,0 +1,81 @@ +"""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) From bb126792840944970c0d359e1bf772fad83610b6 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 21:59:17 +0100 Subject: [PATCH 02/10] ctimeIntegration: wrappers take the arguments TimeIntegration.py passes The 2019 Cython rewrite of the ctimeIntegration C module (#879, f2771244) gave seven wrappers new signatures: the sizes as leading arguments (nElements_global, nNodes_element, ...), the 1D limiter's flag first instead of last, and psiTCtteDT writing dt into a seventh argument instead of returning it. TimeIntegration.py kept the old calls, so 11 of its 16 calls into the module raised TypeError on first use: every 1D DG limiter (DGlimiterP1Lagrange1d and the classes built on it, the _withVacuumTol ones), the Durlofsky 2D/3D limiters and PsiTCtte's dt choice. The three Cockburn/shape-gradient wrappers already took the old arguments and are unchanged. The wrappers again infer the sizes from the array shapes, exactly as the C module did (nSpace from grad_v's last axis, nDOF from l2g), take the limiter flag last with its old default 1 (MUSCL), and psiTCtteDT returns dt. Nothing else calls these wrappers. TimeIntegration.py calls that do not bind: main 11 of 16 -> 0 P1 1D limiter, 4 elements: linear data kept, an extremum cell and its falling neighbour flattened to their means (flags 0/1/2; vacuum-tolerance variant the same); psiTCtteDT reads all 6 points of a (2,3) array and returns sqrt(tau/(max dy + 1e-8))*0.9 exactly. test/test_time_integration_wrappers.py: 3 tests, all 3 fail on main. test/shallow_water/dambreak_1d (RKDG P1, SSPRK2, Rusanov, P1 limiter) then runs to T=20 in 2.9 s, given the Cython 3 staticmethod wrapping of the class-stored limiter functions in TimeIntegration.py (cekees/proteus#204, C47; without it Cython 3 passes self as an extra argument). Measured with #204's TimeIntegration.py hunks applied locally: mass 29/6 at all 101 outputs to 1.2e-14 (the P1 interpolation of the step), min h -5.0e-8 (3 outputs below 0, at the wall cell as the front arrives), and with an element partition (1 layer) serial = 2, 3 and 4 ranks to the last bit in h and hu at all 101 outputs. The deck test lands after #204. Co-Authored-By: Claude Opus 5.5 --- proteus/ctimeIntegration.pyx | 355 ++++++++++++------------- test/test_time_integration_wrappers.py | 79 ++++++ 2 files changed, 247 insertions(+), 187 deletions(-) create mode 100644 test/test_time_integration_wrappers.py 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/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 From f2938dc9ed3d713eb80a9c81baa5403ad18ebda7 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 22:19:54 +0100 Subject: [PATCH 03/10] RANS2P_IB: a deck (one beam in a 3D channel), and the four faults it found RANS2P_IB (RANS2P plus immersed Euler-Bernoulli beams, beamFEM) had no deck and could not have run one: 1. getJacobian read `elf.coefficients` (NameError on the first Newton iteration). 2. Coefficients passed RANS2P its default arguments instead of its own, so rho, nu, g, nd, LS_model, useMetrics, ... were dropped: every RANS2P_IB run was 3D single-phase water (nu_0=1e-3, nd=2, LS_model=2 came out as 1.004e-6, 3, None). It now forwards them; the defaults it had (dragAlpha 0.01, eb_penalty_constant 10) are unchanged. 3. calculateBeams integrated the beam loads q1-q3, vel_avg and netBeamDrag over ghost elements too, then summed over ranks, so elements shared by two ranks counted twice. They are now integrated over owned elements; every local element still gets q_dragBeam*, the fluid's source. (C0's kernel territory, mprans/*.h.) 4. postStep printed the forces to stdout, so the history files RANS2P opens stayed empty, and wrote the deflection row directly. It now queues RANS2P's rows (time, wetted area, forces, moments), the beam drag and the deflection with deferWrite, so a retried step leaves no row. (RANS2P's postStep itself does not fit: this LevelModel is an older copy without the lagged dV/active-element arrays, which its kernel does not read.) Deck test/rans2p_ib/beam_channel: 0.4 x 0.69 x 0.2 channel from rest, inflow 0.1 m/s, no slip on the bottom and side walls, one beam (radius 5 mm, height 0.1) at (0.2, 0.345); rigid (default) or flexible (beamRigid=False, EI=GJ=1e-4, bends ~1 cm by t=0.04). P1, Backward Euler, dt 0.01 to 0.04, direct solves, Newton atol 1e-10. Serial vs N ranks (node partition, no overlap, sums over owned nodes of p, u, v, w and squares, beam drag, q1-q3, beam vertex coordinates), worst relative difference: before 3. (1, 2 fixed) after rigid, 2 ranks drag 1.12e-2 vs 7.99e-3 (+40%), q1 +41% 4e-16 (drag1, 4e-6 itself, 4e-14) flexible, 2 ranks drag +42%, beam x-sum 1.0331 vs 1.0195, 4e-16 fluid u-sum 3.7e-8 rigid/flexible, 3 ranks 4e-16 (the fluid matched to 1e-15 before 3. for a rigid beam: the loads feed back only through a flexible beam's shape) Retry (test_step_retry, two new cases): RANS2P_IB failing its 3rd solve, retried at half the step, against a failure-free run on the same times: 0.0 in p, u, v, w, rigid and flexible; the deflection and time histories hold one row per accepted step (before 4.: timeHistory.txt empty). NORESTORE gives 0.0 too on this single-model deck: postStep runs only after an accepted step, so nothing of the beam's is ahead of the retry. In a multi-model deck a later model's failure would leave a flexible beam moved (Beam_Solver Phi, xv/yv/zv, xq/yq/zq): RANS2P_IB needs a state() for that, noted in the round handoff. Tests: test/rans2p_ib/test_rans2p_ib.py (3; on main the deck stops at 1., and the coefficient test fails), test_step_retry.py beam_channel cases (2). Co-Authored-By: Claude Opus 5.5 --- proteus/mprans/RANS2P_IB.h | 8 ++- proteus/mprans/RANS2P_IB.py | 108 ++++++++++++++--------------- test/rans2p_ib/beam_channel_mpi.py | 54 +++++++++++++++ test/rans2p_ib/beam_channel_n.py | 42 +++++++++++ test/rans2p_ib/beam_channel_p.py | 93 +++++++++++++++++++++++++ test/rans2p_ib/test_rans2p_ib.py | 55 +++++++++++++++ test/test_step_retry.py | 15 ++++ 7 files changed, 319 insertions(+), 56 deletions(-) create mode 100644 test/rans2p_ib/beam_channel_mpi.py create mode 100644 test/rans2p_ib/beam_channel_n.py create mode 100644 test/rans2p_ib/beam_channel_p.py create mode 100644 test/rans2p_ib/test_rans2p_ib.py 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 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/test_rans2p_ib.py b/test/rans2p_ib/test_rans2p_ib.py new file mode 100644 index 0000000000..43e0abcf75 --- /dev/null +++ b/test/rans2p_ib/test_rans2p_ib.py @@ -0,0 +1,55 @@ +"""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 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): From 64124044699d2b06f0d1259cab63b5302e6b4fce Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 22:23:52 +0100 Subject: [PATCH 04/10] test/shallow_water/dambreak_1d: run the deck, serial and on 2 ranks The deck (RKDG P1 shallow water, Rusanov flux, P1 limiter, SSPRK2, a dry-bed dam break between two walls, 30 elements to T=20) stopped in RusanovNumericalFlux's constructor (7d95b913) and then in the 1D limiter wrappers (bb126792). Measured (macOS arm64, with cekees/proteus#204's TimeIntegration.py, which Cython 3 needs for the limiter class attributes): mass 29/6 (the P1 interpolation of the step) at all 101 outputs to 1.2e-14; h below zero at 3 of 101 outputs, by at most 5.0e-8 (t=0.6, the wall cell as the front arrives: the P1 limiter does not preserve positivity); h in [0.3645, 0.4985] at T=20, settling to the rest depth 29/60; element partition, 1 layer, bjacobi/LU mass solve: serial = 2, 3 and 4 ranks to the last bit in h and hu at every output (2 ranks in the test). The test skips, saying why, until #204's staticmethod attributes are on the branch, and runs from then on. The p module's docstring is now raw (it held "\in", a SyntaxWarning). Co-Authored-By: Claude Opus 5.5 --- .../dambreak_1d/sw_dam_break_1d_p.py | 2 +- .../dambreak_1d/test_dambreak_1d.py | 99 +++++++++++++++++++ 2 files changed, 100 insertions(+), 1 deletion(-) create mode 100644 test/shallow_water/dambreak_1d/test_dambreak_1d.py 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) From 7c2a5d52fd60970fd5a5e55e0122c45125912274 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 23:23:36 +0100 Subject: [PATCH 05/10] Coverage floors: NCLS3P.py 76.2 -> 75.6 (round 5 moved 15 covered statements to Transport) The floor was recorded on main ce43b736 (33f30b1a): 522 of 685 statements, 76.2%. cekees/proteus#199 (C0 round 5, 76f3e125) replaced NCLS3P's own scatterDOFArray (15 statements, all covered) by the one on OneLevelTransport, so main c09280cb has 507 of 670, 75.67%: the same lines untested, a smaller file. Every pull request on c09280cb without #204 (whose staticmethod lines add covered statements to the file) is below the floor: this one's coverage job and #212's both report NCLS3P 75.7%. Line by line, the NCLS3P lines hit in this branch's run and in #204's run (mapped through #204's edits) are the same set. The new floor is the measured value rounded down to a tenth. Co-Authored-By: Claude Opus 5.5 --- .github/coverage/floors.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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, From 83ac62c9150328b02acb9c31853d1f135a345993 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 6 Nov 2026 08:51:06 +0100 Subject: [PATCH 06/10] Step retry: a model that finished in the failed attempt steps again SO_base.modelStepTaken (Sequential_MinModelStep, default_so's default, 11 so files under test/; and Sequential_tnList) moves a model's t_model_last to the end of the system step as soon as that model's step is taken. When a later model of the same system step failed, restoreSystemStepState put the dofs and the lag state back but not t_model_last, so the retry, stepping to an earlier time, found the model already past it and skipped it; from then on it ran one step behind the others and never reached T. The controllers that override modelStepTaken (Sequential_MinAdaptiveModelStep, the TwoPhaseFlow decks, and the FLCBDF ones) do not move it there, which is why test_step_retry's cases never saw it. saveSystemStepState now keeps each model's t_model_last and restoreSystemStepState puts it back. The case: test/rans2p_ib/beam_channel_so.py, the beam channel with a dye (VOF, carried by RANS2P_IB's velocity) after RANS2P_IB; the dye's third solve is reported failed (step_retry_mpi.py FAILMODEL=dye FAILAT=3), the step 0.02 -> 0.03 is cut to 0.025, against a failure-free run on the same steps (REPLAY). Largest relative difference over each model's dofs, macOS arm64, serial, SuperLU: | | RANS2P_IB p, u, v, w | dye | |---|---|---| | main c09280cb (rigid and flexible beam) | 4.8e-2, 8.0e-2, 5.2e-2, 9.2e-2 | 3.3e-4 | | this commit | 1.2e-6, 2.8e-6, 6.4e-6, 3.3e-6 | 2.7e-8 | Model solves of the forced run on main: beam_channel at 0.01, 0.02, 0.03 | retry: dye 0.025, beam_channel 0.02 (dt 0.005), dye 0.03, beam_channel 0.025, ... beam_channel last at 0.035. With the fix both models solve at 0.025, 0.03, 0.035, 0.04. What remains is RANS2P_IB's lagged beam drag, which the retry does not restore yet (next commit); the new test bounds it at 1e-5 until then. test_step_retry.py (all cases) passes unchanged. Co-Authored-By: Claude Opus 5.5 --- proteus/NumericalSolution.py | 8 ++++++ test/rans2p_ib/beam_channel_so.py | 12 +++++++++ test/rans2p_ib/dye_n.py | 37 ++++++++++++++++++++++++++++ test/rans2p_ib/dye_p.py | 35 ++++++++++++++++++++++++++ test/rans2p_ib/test_rans2p_ib.py | 41 +++++++++++++++++++++++++++++++ 5 files changed, 133 insertions(+) create mode 100644 test/rans2p_ib/beam_channel_so.py create mode 100644 test/rans2p_ib/dye_n.py create mode 100644 test/rans2p_ib/dye_p.py 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/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 index 43e0abcf75..e33187f159 100644 --- a/test/rans2p_ib/test_rans2p_ib.py +++ b/test/rans2p_ib/test_rans2p_ib.py @@ -53,3 +53,44 @@ def test_beam_channel_on_two_ranks_is_the_serial_answer(tmp_path, launch_policy, 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.""" + 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-5, (key, worst) # the lagged drag: next commit From 07f4d6a0583534ce1371d9db9d2736d05b63c3c7 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 6 Nov 2026 08:52:59 +0100 Subject: [PATCH 07/10] RANS2P_IB: a state() for the lagged beam drag and the beams RANS2P_IB's beam drag is lagged: after each model step its LevelModel's beamStep spreads the drag from the beams onto the mesh (Coefficients.q_dragBeam1-3, ebqe_dragBeam1-3, which the next step's residual reads) and integrates the beam loads (q1-q3); postStep then moves a flexible beam (xv..zq, each beamFEM solver's Phi and x, y, z) and numbers the beam file (tStep). Both run before the later models of the system step. Its Coefficients had no state() (RANS2P.Coefficients.state returns None for this older LevelModel copy), so when a later model failed, the retry restored the dofs and kept the failed attempt's drag and beam. Coefficients.state()/set_state() now carry them (vel_avg is set in attachModels so state() has it before the first beamStep). The same state is archived for a hot start (solver_state//coefficients). test/rans2p_ib: the dye after RANS2P_IB fails its third solve (beam_channel_so.py, previous commit), the retry against a failure-free run on the same steps. Largest relative difference over each model's dofs, macOS arm64, serial, SuperLU: | | RANS2P_IB p, u, v, w | dye | flexible beam deflection at the retried step | |---|---|---|---| | previous commit, rigid | 1.2e-6, 2.8e-6, 6.4e-6, 3.3e-6 | 2.7e-8 | | | previous commit, flexible | 1.3e-6, 3.2e-6, 6.4e-6, 3.8e-6 | 2.7e-8 | | | drag restored, beam not (flexible) | 8.8e-8, 1.3e-7, 6.6e-7, 1.8e-7 | 7.7e-10 | 8.2941e-3 vs 8.2956e-3 | | this commit, rigid and flexible | 0.0 | 0.0 | equal, no row from the failed attempt | The test's bound goes from 1e-5 to 1e-12, and the flexible case checks the deflection history. Ran: test/rans2p_ib (5, incl. serial = 2 ranks), test_step_retry beam_channel cases (2). Not checked: a hot start of the deck, which stops before reading any state (the restart reads regular_9x14x5_grid.node/.ele, which an in-memory box mesh never writes). Co-Authored-By: Claude Opus 5.5 --- proteus/mprans/RANS2P_IB.py | 42 ++++++++++++++++++++++++++++++++ test/rans2p_ib/test_rans2p_ib.py | 12 +++++++-- 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/proteus/mprans/RANS2P_IB.py b/proteus/mprans/RANS2P_IB.py index c9a88b320e..5f386b7846 100644 --- a/proteus/mprans/RANS2P_IB.py +++ b/proteus/mprans/RANS2P_IB.py @@ -118,6 +118,7 @@ def attachModels(self, modelList): RANS2P.Coefficients.attachModels(self, modelList) self.initializeBeams() self.netBeamDrag = np.array([0.0]) + self.vel_avg = np.array([0.0, 0.0, 0.0]) # set by beamStep after each step; here so state() has it from the start def initializeElementQuadrature(self, t, cq): RANS2P.Coefficients.initializeElementQuadrature(self, t, cq) @@ -141,6 +142,47 @@ def initializeMesh(self, mesh): if self.beamRigid == False: self.deflectionHistory = open("deflectionHistory.txt", "w") + #: what a step leaves for the next one in the coefficients. The drag the residual reads is lagged: the + #: LevelModel's beamStep, after each model step, spreads it from the beams' quadrature points (q_dragBeam*, + #: ebqe_dragBeam*) and integrates the beam loads (q1-q3, with vel_avg and netBeamDrag); postStep then moves a + #: flexible beam (xv..zq; updateBeams falls back to the last step's coordinates when a deflected beam would cross + #: z = 0) + BEAM_STATE = ("q_dragBeam1", "q_dragBeam2", "q_dragBeam3", "ebqe_dragBeam1", "ebqe_dragBeam2", "ebqe_dragBeam3", + "q1", "q2", "q3", "vel_avg", "netBeamDrag", "xv", "yv", "zv", "xq", "yq", "zq") + + def state(self): + """What a step leaves for the next that the solution does not determine: the lagged beam drag and loads, the + beams (their coordinates, each beam solver's rotations Phi and vertex coordinates) and the beam file counter, + plus what RANS2P.Coefficients keeps, if anything (its LevelModel arrays are not in this older copy, so nothing + today). + + beamStep and postStep run after this model's step, before the later models of the system step: without this + a later model's failure left the drag of the failed attempt (and a flexible beam moved), and the retry started + from them.""" + if not hasattr(self, "Beam_Solver"): + return None + st = dict(super().state() or {}) + for name in self.BEAM_STATE: + st["beam_" + name] = getattr(self, name).copy() + for I, b in enumerate(self.Beam_Solver): + for name in ("Phi", "x", "y", "z"): + st["beam%d_%s" % (I, name)] = np.array(getattr(b, name), dtype=float) + st["beam_tStep"] = int(self.tStep) + return st + + def set_state(self, state): + if any(not k.startswith("beam") for k in state): + super().set_state({k: v for k, v in state.items() if not k.startswith("beam")}) + for name in self.BEAM_STATE: + # vel_avg is rebound by beamStep (a new array): restore by value into whatever is bound now + getattr(self, name)[...] = numpy.reshape(state["beam_" + name], getattr(self, name).shape) + for I, b in enumerate(self.Beam_Solver): + # updateSolution changes Phi in place and the fallback replaces it with a copy: replace it here too + b.Phi = np.array(state["beam%d_Phi" % I], dtype=float) + for name in ("x", "y", "z"): + getattr(b, name)[:] = state["beam%d_%s" % (I, name)] + self.tStep = int(state["beam_tStep"]) + def postStep(self, t, firstStep=False): # Not RANS2P.Coefficients.postStep: this LevelModel is an older copy of RANS2P's, without the lagged dV and # active-element arrays that one updates (and the RANS2P_IB kernel does not read them). diff --git a/test/rans2p_ib/test_rans2p_ib.py b/test/rans2p_ib/test_rans2p_ib.py index e33187f159..53f566b448 100644 --- a/test/rans2p_ib/test_rans2p_ib.py +++ b/test/rans2p_ib/test_rans2p_ib.py @@ -82,7 +82,10 @@ def test_a_later_models_failure_is_retried_from_the_step_start(tmp_path, opts): 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.""" + 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) @@ -93,4 +96,9 @@ def test_a_later_models_failure_is_retried_from_the_step_start(tmp_path, opts): 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-5, (key, worst) # the lagged drag: next commit + 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) From 8f919e2825dcf831a110a891bcb73afc0b22cc9a Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 6 Nov 2026 09:01:28 +0100 Subject: [PATCH 08/10] NumericalFlux: the diffusion exterior fluxes refuse periodic conditions Round 6's sweep found seven fluxes that accept the periodic argument and pair no periodic faces (Diffusion_IIPG/SIPG_exterior, Stokes_Advection_DiagonalUpwind_Diffusion_IIPG/SIPG_exterior, NoFlux, StrongDirichlet, DoNothing). All seven are boundary-only, so they serve C0 spaces, where the periodic DOF map already joins the seam; whether skipping the pairing matters was measured, not assumed. NF_base marks the face points of a periodic pair as DOF boundary with no value. The diffusion fluxes then apply their weak-Dirichlet flux there with the trace equal to the interior value: the residual keeps a one-sided -a grad(u).n on each side of the seam (it does not cancel, grad(u) jumps across) and the Jacobian a penalty the residual lacks. The other three put nothing at the seam that depends on it. C0 P1, -Laplace u = 5 pi^2 sin(2 pi x) sin(pi y) on the unit square, u = 0 at y = 0, 1, 8x8 cells, L2 error and Newton iterations; "without" makes x = 0, 1 Dirichlet (the exact solution vanishes there too): | flux | without the periodic pair | periodic in x | |---|---|---| | Diffusion_IIPG_exterior | 3.2e-2, 1 | 3.5e-1, not converged in 30 (10 retries) | | Diffusion_SIPG_exterior | 4.3e-2, 1 | 4.4e-1, not converged | | StrongDirichlet | 4.6e-2, 1 | 4.7e-2, 1 | | NoFlux | 4.1 (imposes no Dirichlet) | 4.1, 1 | | DoNothing | 1.5 (imposes no Dirichlet) | 1.5, 1 | So the four diffusion fluxes (Stokes: the same diffusive flux at the same points, not run, no periodic Stokes deck exists) now raise ValueError when built with periodic conditions, naming StrongDirichlet. NF_base checks the exact class name (PERIODIC_SEAM_REFUSED); subclasses with their own exterior flux (the HamiltonJacobi Lesaint-Raviart SIPG pair) are not covered. StrongDirichlet, DoNothing and NoFlux keep accepting periodic conditions. Who builds them with periodic conditions (Transport hook over test/periodic, TADR, LS_with_edgeBased_EV, POD, test_adr_problem, test_pro_cli; plus a grep of every p/n module that sets periodicDirichletConditions): StrongDirichlet: ADRProblem's StrongDirichletFactory (test_adr_problem's periodic Poisson, periodic channel, periodic on 2-4 ranks); DoNothing: test/periodic ls_consrv_n, redist_n; la_periodicGauss_n, hotstart_measure/rkev_gauss_n (MCorr, RDLS); the four refused classes: nobody. The Stokes ones are used without periodic conditions only (solver_tests stokesDrivenCavity_2d_n, mesh_tests stokes_2d_n), unaffected. Not changed, for the next round: Advection_DiagonalUpwind_Diffusion_IIPG_exterior pairs the periodic faces and has the same seam defect on the same problem (1.06e-1 vs 3.2e-2, Newton not converged). The mprans fluxes built on it (VOF, NCLS, CLSVOF, Kappa, Dissipation) run periodic decks in test/periodic; VOF and NCLS have no diffusion there. Refusing it would break them; the fix is to drop the exterior diffusive flux at C0 periodic face points. Ran: test_numerical_flux_constructors (5), test_numerical_flux_arguments, test/periodic, test/TADR, test/LS_with_edgeBased_EV, test/POD, test_adr_problem, solver_tests/test_stokesDC_TH_2D, mesh_tests/test_quads: 127 passed, 6 skipped, 1 failed (test_ymfrun_knows_which_domains_fix_the_mesh: the env's ymf has no discretizations.mesh_from_domain; fails the same without this commit). Co-Authored-By: Claude Opus 5.5 --- proteus/NumericalFlux.py | 24 +++++++++++++++++ test/test_numerical_flux_constructors.py | 33 ++++++++++++++++++++++++ 2 files changed, 57 insertions(+) diff --git a/proteus/NumericalFlux.py b/proteus/NumericalFlux.py index 1a942cc989..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() diff --git a/test/test_numerical_flux_constructors.py b/test/test_numerical_flux_constructors.py index 1fddda17ab..62c17f8e39 100644 --- a/test/test_numerical_flux_constructors.py +++ b/test/test_numerical_flux_constructors.py @@ -79,3 +79,36 @@ def test_every_base_constructor_call_binds_and_forwards_the_periodic_conditions( 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 From 61f08b510937bf10d4d3d4fa527d09d7f4fe2121 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 6 Nov 2026 14:57:00 +0100 Subject: [PATCH 09/10] Partition-independence tests: a failed helper launch keeps its whole output The ten tests that launch a *_mpi.py helper (mcorr, av, ncls_edge, gn_sw2dcv, twophaseflow, sw2dcv, vof, tadr, vos partition independence; subdomain mesh size) asserted `out.returncode == 0, out.stdout[-2000:] + out.stderr[-2000:]`. #208's CI had a 2-rank risingBubble segfault (signal 11, both ranks) on ceds-model-4 in the twophaseflow fixture, and the 2000-character tail had cut the backtrace frames away; C17/C14 could not reproduce it in 300 runs. A timeout lost everything: subprocess.run raises TimeoutExpired, and the mcorr staging race (next commit) ended that way, 2809 s. test/launch_output.py runs the launch; on a nonzero exit status or a timeout it writes the whole stdout and stderr to launch.stdout.txt / launch.stderr.txt in the launch's directory (under the test's tmp_path) and raises AssertionError naming them, with the exit status and signal ("killed by SIGSEGV", "exit status 139 (SIGSEGV?)"), the command, and the last 8000 characters of stderr and 3000 of stdout. The tests load it by path: CI's --import-mode=importlib does not put test/ on sys.path. test_launch_output.py: a process writing 260000 characters to stderr then SIGSEGV -> the message has the signal name and the last frame, the file has the first (before: the last 2000 characters only); a hang (3 s timeout) -> the message and file have what it printed (before: TimeoutExpired, no output). Ran (pytest from outside test/, CI's import mode): the ten modules collect (43 tests); test_launch_output (6), test_subdomain_mesh_size, test_vos_partition_independence pass. Co-Authored-By: Claude Opus 5.5 --- test/launch_output.py | 63 +++++++++++++++++++ test/test_av_partition_independence.py | 9 ++- test/test_gn_sw2dcv_partition_independence.py | 9 ++- test/test_launch_output.py | 46 ++++++++++++++ test/test_mcorr_partition_independence.py | 9 ++- test/test_ncls_edge_partition_independence.py | 9 ++- test/test_subdomain_mesh_size.py | 9 ++- test/test_sw2dcv_partition_independence.py | 9 ++- test/test_tadr_partition_independence.py | 9 ++- ...est_twophaseflow_partition_independence.py | 9 ++- test/test_vof_partition_independence.py | 9 ++- test/test_vos_partition_independence.py | 9 ++- 12 files changed, 169 insertions(+), 30 deletions(-) create mode 100644 test/launch_output.py create mode 100644 test/test_launch_output.py diff --git a/test/launch_output.py b/test/launch_output.py new file mode 100644 index 0000000000..5cb34e701a --- /dev/null +++ b/test/launch_output.py @@ -0,0 +1,63 @@ +"""Run a helper (under mpiexec or not) for a test and, when it fails, keep everything it printed. + +The partition-independence tests used to assert ``out.stdout[-2000:] + out.stderr[-2000:]``: a crash's backtrace (a +2-rank segfault in the twophaseflow fixture on ceds-model-4, 2026-11-05) was cut away, and a hang lost all of its +output (subprocess.TimeoutExpired). Here a failed launch writes the whole stdout and stderr to ``launch.stdout.txt`` and +``launch.stderr.txt`` in its working directory (under the test's tmp_path) and the assertion message names them, with +the exit status (and signal), the command and the last 8000 characters of stderr and 3000 of stdout. + +Load it by path (CI runs pytest with --import-mode=importlib, which does not put test/ on sys.path). +""" +import os +import signal +import subprocess + +STDERR_TAIL, STDOUT_TAIL = 8000, 3000 + + +def _text(data): + if data is None: + return "" + return data.decode(errors="replace") if isinstance(data, bytes) else data + + +def _status(returncode): + """'exit status 139 (SIGSEGV)', 'killed by SIGSEGV', 'exit status 1'.""" + sig = -returncode if returncode < 0 else returncode - 128 if returncode > 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/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_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_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) From e2aa8b11979126fe32d7c5da937640883c086126 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 6 Nov 2026 15:01:32 +0100 Subject: [PATCH 10/10] MPI test helpers: stage the deck's files atomically from every rank proteustoolkit/proteus#1343's linux-aarch64 leg (run 42184914225, GitHub-hosted) hung for 2809 s in test_mcorr_partition_independence's 2-rank launch: one rank stopped with `mcpkg/thelper_ncls_p.py, line 32: RD_model=RD_model ... NameError: name 'RD_model' is not defined`, the other waited in its first collective. test/mcorr_partition_mpi.py (from f2b0ab95, already on main) runs its staging on every rank: each rank rewrites mcpkg/*.py in the shared working directory with open(path, "w").write(text) and then imports them, so a rank can star-import a module another rank has just truncated. Nine helpers stage that way at module level (every rank): mcorr, ncls_edge, tadr, vos, subdomain_mesh_size (open "w"), av, twophaseflow (shutil.copy, also truncate-then-write), sw2dcv, gn_sw2dcv (copy, then every rank appends a line to circular_damBreak.py and truncates __init__.py). test/rank_staging.py writes the bytes to a per-process name and os.replace()s it into place (atomic on POSIX): a reader sees a whole file, and every rank writes the same bytes. No collective, so it runs before Comm.init. The two shallow-water helpers now write the file with its appended line in one go. The race, measured (serial-free 2-rank launches of mcorr_partition_mpi.py's staging and its first imports, macOS arm64, MPICH): it needs one rank to start late (proteus imports on a loaded runner) and a writer to be descheduled between truncate and write. With neither, 0 failures in 480 rank-imports (4 and 8 ranks). With rank 1 starting up to 1.5 s late and 50 ms between truncate and write (after: between writing the temporary file and the rename), Mac at load 36-52: | | launches | rank-imports failed | hung launches | |---|---|---|---| | main's staging | 30 + 30 | 4/60 + 2/58 (NameError: RD_model 4, NCLS 2, VOF 1 incl. an aborted batch) | 0 + 1 | | rank_staging | 30 + 30 | 0/60 + 0/60 | 0 | The same errors as on CI. A third, aborted batch of main's staging also hung one launch (3 h 50 min, mpiexec waiting with both ranks gone) after 5 NameErrors in its first ~8 launches. test/test_rank_staging.py: a reader thread never sees a partial file through rank_staging, and does through the in-place write it replaces (the control); copy() into a directory; and no *_mpi.py helper copies, writes or appends a deck file at module level outside rank_staging (results by extension, writes under `if rank` and inside functions excepted; step_retry_mpi.py is serial). On main the last test lists all 14 sites in the nine helpers. It loads rank_staging by path (CI's --import-mode=importlib); the helpers import it by name, as scripts, with test/ on sys.path. Ran: the nine helpers' tests (mcorr, ncls_edge, tadr, vos, av, sw2dcv, gn_sw2dcv, twophaseflow partition independence, subdomain mesh size) and test_rank_staging: 39 passed, -n 2; test_rank_staging (4) from outside test/ with CI's import mode. Co-Authored-By: Claude Opus 5.5 --- test/av_partition_mpi.py | 5 +- test/gn_sw2dcv_partition_mpi.py | 10 ++-- test/mcorr_partition_mpi.py | 3 +- test/ncls_edge_partition_mpi.py | 3 +- test/rank_staging.py | 32 +++++++++++ test/subdomain_mesh_size_mpi.py | 3 +- test/sw2dcv_partition_mpi.py | 10 ++-- test/tadr_partition_mpi.py | 3 +- test/test_rank_staging.py | 88 ++++++++++++++++++++++++++++++ test/twophaseflow_partition_mpi.py | 3 +- test/vos_partition_mpi.py | 3 +- 11 files changed, 147 insertions(+), 16 deletions(-) create mode 100644 test/rank_staging.py create mode 100644 test/test_rank_staging.py diff --git a/test/av_partition_mpi.py b/test/av_partition_mpi.py index 4955831b1b..e526a324fa 100644 --- a/test/av_partition_mpi.py +++ b/test/av_partition_mpi.py @@ -9,17 +9,18 @@ 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 HERE = os.path.dirname(os.path.abspath(__file__)) SRC = os.path.join(HERE, "TwoPhaseFlow") for name in os.listdir(SRC): if name.startswith("meshDambreak") or name == "__init__.py": - shutil.copy(os.path.join(SRC, name), name) + rank_staging.copy(os.path.join(SRC, name), name) case = open(os.path.join(SRC, "damBreak.py")).read() case += ("\n[(setattr(mm.n, 'nl_atol_res', 1.0e-10), setattr(mm.n, 'l_atol_res', 1.0e-14), " "(setattr(mm.n, 'linTolFac', 1.0e-14) if hasattr(mm.n, 'linTolFac') else None)) " "for mm in myTpFlowProblem.SystemPhysics.modelDict.values()]\n") -open("damBreak.py", "w").write(case) +rank_staging.write("damBreak.py", case) if os.environ.get("NOFIX"): from proteus.mprans import RANS3PF diff --git a/test/gn_sw2dcv_partition_mpi.py b/test/gn_sw2dcv_partition_mpi.py index 121ccb91da..80b3a7d2b6 100644 --- a/test/gn_sw2dcv_partition_mpi.py +++ b/test/gn_sw2dcv_partition_mpi.py @@ -8,11 +8,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.GN_swe_parameters['LUMPED_MASS_MATRIX'] = 1\n") +# 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.GN_swe_parameters['LUMPED_MASS_MATRIX'] = 1\n") +rank_staging.write("__init__.py", "") from proteus.mprans import GN_SW2DCV as SW2DCV if os.environ.get("NOSCAT"): 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/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_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/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,