From ab231116fdad24ec76f40344e380a81393283d98 Mon Sep 17 00:00:00 2001 From: Felix Oesterle Date: Tue, 12 May 2026 14:45:27 +0200 Subject: [PATCH 1/5] fix(com1DFA): ensure arrays returned by `updateFieldsC` are owned copies (BUG COM5) - Updated multiple `fields` and `particles` assignments to ensure `.copy()` is used when converting to `np.asarray`, preventing potential issues with memory views. - Adjusted handling of `bondDist` to avoid creating problematic views during inplace modifications. - Added tests to verify that arrays returned by `updateFieldsC` have `owndata=True` and no base reference. --- avaframe/com1DFA/DFAfunctionsCython.pyx | 52 ++++++++++++----------- avaframe/tests/test_DFAfunctionsCython.py | 51 ++++++++++++++++++++++ 2 files changed, 78 insertions(+), 25 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 919707a3c..b0b8c5d56 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -481,7 +481,7 @@ def computeForceC(cfg, particles, fields, dem, int frictType, int resistanceType if entrMassCell < 0: entrMassCell = 0 entrMassRaster[indCellY, indCellX] = entrMassCell - fields['entrMassRaster'] = np.asarray(entrMassRaster) + fields['entrMassRaster'] = np.asarray(entrMassRaster).copy() return particles, force, fields @@ -1398,27 +1398,27 @@ def updateFieldsC(cfg, particles, dem, fields): FTStopBilinear[j, i] = - MassStopBilinear[j, i] / (areaRaster[j, i] * rho) # / m * FTBilinear[j, i] FTEntBilinear[j, i] = - MassEntBilinear[j, i] / (areaRaster[j, i] * rhoEnt) - fields['FM'] = np.asarray(MassBilinear) - fields['FV'] = np.asarray(VBilinear) - fields['Vx'] = np.asarray(VXBilinear) - fields['Vy'] = np.asarray(VYBilinear) - fields['Vz'] = np.asarray(VZBilinear) - fields['FT'] = np.asarray(FTBilinear) - fields['pfv'] = np.asarray(PFV) - fields['pft'] = np.asarray(PFT) - fields['dmDet'] = np.asarray(DMDet) - fields['FTStop'] = np.asarray(FTStopBilinear) - fields['FTDet'] = np.asarray(FTDetBilinear) - fields['FTEnt'] = np.asarray(FTEntBilinear) + fields['FM'] = np.asarray(MassBilinear).copy() + fields['FV'] = np.asarray(VBilinear).copy() + fields['Vx'] = np.asarray(VXBilinear).copy() + fields['Vy'] = np.asarray(VYBilinear).copy() + fields['Vz'] = np.asarray(VZBilinear).copy() + fields['FT'] = np.asarray(FTBilinear).copy() + fields['pfv'] = np.asarray(PFV).copy() + fields['pft'] = np.asarray(PFT).copy() + fields['dmDet'] = np.asarray(DMDet).copy() + fields['FTStop'] = np.asarray(FTStopBilinear).copy() + fields['FTDet'] = np.asarray(FTDetBilinear).copy() + fields['FTEnt'] = np.asarray(FTEntBilinear).copy() if computeP: - fields['ppr'] = np.asarray(PP) - fields['P'] = np.asarray(PBilinear) + fields['ppr'] = np.asarray(PP).copy() + fields['P'] = np.asarray(PBilinear).copy() if computeTA: - fields['TA'] = np.asarray(travelAngleField) - fields['pta'] = np.asarray(PTA) + fields['TA'] = np.asarray(travelAngleField).copy() + fields['pta'] = np.asarray(PTA).copy() if computeKE: - fields['pke'] = np.asarray(PKE) + fields['pke'] = np.asarray(PKE).copy() for k in range(nPart): @@ -1674,7 +1674,9 @@ def computeCohesionForceC(cfg, particles, force): force['forceSPHX'] = np.asarray(forceSPHX) force['forceSPHY'] = np.asarray(forceSPHY) force['forceSPHZ'] = np.asarray(forceSPHZ) - particles['bondDist'] = np.asarray(bondDist) + # bondDist is modified in-place via the memoryview; do not reassign + # particles['bondDist'] here, otherwise numpy will create a VIEW with + # base=memoryview which can hang on deallocation on Windows. return force, particles @@ -1775,9 +1777,9 @@ def initializeBondsC(particles, triangles): bondDist[bondStart2[p2+1]-1] = distanceIni bondStart2[p2+1] = bondStart2[p2+1] - 1 - particles['bondStart'] = np.asarray(bondStart) - particles['bondPart'] = np.asarray(bondPart) - particles['bondDist'] = np.asarray(bondDist) + particles['bondStart'] = np.asarray(bondStart).copy() + particles['bondPart'] = np.asarray(bondPart).copy() + particles['bondDist'] = np.asarray(bondDist).copy() return particles def countRemovedBonds(particles, mask, nRemove): @@ -1870,9 +1872,9 @@ def removedBonds(particles, mask, nRemove, nBondRemove): countBond = countBond + 1 - particles['bondStart'] = np.asarray(bondStartNew) - particles['bondPart'] = np.asarray(bondPartNew) - particles['bondDist'] = np.asarray(bondDistNew) + particles['bondStart'] = np.asarray(bondStartNew).copy() + particles['bondPart'] = np.asarray(bondPartNew).copy() + particles['bondDist'] = np.asarray(bondDistNew).copy() return particles diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index b7375191d..97e9c3b97 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -619,6 +619,12 @@ def test_initializeBondsC(): bondDist = particles["bondDist"] bondPart = particles["bondPart"] assert np.array_equal(bondStart, np.asarray([0, 2, 4, 6])) + assert particles["bondStart"].flags.owndata + assert particles["bondStart"].base is None + assert particles["bondPart"].flags.owndata + assert particles["bondPart"].base is None + assert particles["bondDist"].flags.owndata + assert particles["bondDist"].base is None for k in range(nPart): # loop on all bonded particles neighbors = list() @@ -664,6 +670,12 @@ def test_removeBondsC(): bondDist = particles["bondDist"] bondPart = particles["bondPart"] assert np.array_equal(bondStart, np.asarray([0, 1, 1, 2])) + assert particles["bondStart"].flags.owndata + assert particles["bondStart"].base is None + assert particles["bondPart"].flags.owndata + assert particles["bondPart"].base is None + assert particles["bondDist"].flags.owndata + assert particles["bondDist"].base is None for k in range(nPart): # loop on all bonded particles neighbors = list() @@ -682,6 +694,45 @@ def test_removeBondsC(): assert np.array_equal(bondDist, np.asarray([1, 1])) +def test_updateFieldsC_returns_owned_arrays(): + cfg = configparser.ConfigParser() + cfg["GENERAL"] = {"rho": "200.", "rhoEnt": "100.", "interpOption": "2"} + header = {"nrows": 5, "ncols": 5, "cellsize": 1} + dem = {"header": header, "areaRaster": 25 * np.ones((header["nrows"], header["ncols"]))} + + particles = { + "m": np.array([100.0, 100.0, 100.0, 100.0]), + "dmDet": np.array([1.0, 1.0, 1.0, 1.0]), + "dmEnt": np.array([0.0, 0.0, 0.0, 0.0]), + "x": np.array([2.0, 1.0, 2.0, 1.0]), + "y": np.array([2.0, 2.0, 1.0, 1.0]), + "ux": np.array([10.0, 10.0, 10.0, 10.0]), + "uy": np.array([10.0, 10.0, 10.0, 10.0]), + "uz": np.array([10.0, 10.0, 10.0, 10.0]), + "trajectoryAngle": np.array([10.0, 10.0, 10.0, 10.0]), + "stoppedParticles": {"m": np.empty(0), "x": np.empty(0), "y": np.empty(0)}, + "nPart": 4, + } + + fields = { + "computeTA": False, + "computeKE": False, + "computeP": False, + "pfv": np.zeros((1, 1)), + "ppr": np.zeros((1, 1)), + "pft": np.zeros((1, 1)), + "pta": np.zeros((1, 1)), + "pke": np.zeros((1, 1)), + "dmDet": np.zeros((header["nrows"], header["ncols"])), + } + + particles, fields = DFAfunC.updateFieldsC(cfg["GENERAL"], particles, dem, fields) + + for key in ["FM", "FV", "Vx", "Vy", "Vz", "FT", "pfv", "pft", "dmDet", "FTStop", "FTDet", "FTEnt"]: + assert fields[key].flags.owndata + assert fields[key].base is None + + def test_computeCohesionForceC(): cfg = configparser.ConfigParser() cfg["GENERAL"] = { From 2903c04858f8a75d5e7b05f7bdbdb0a0146e8fee Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Tue, 12 May 2026 15:20:22 +0200 Subject: [PATCH 2/5] fix(com1DFA): convergence check and division-by-zero in reprojection samosProjectionIteratrive: fix Lxi/Lyi never updated inside loop, causing convergence check to always compare against initial cell (perf-only, no correctness change). distConservProjectionIteratrive: add distn==0 guard before dist/distn to prevent NaN for stopped particles. --- avaframe/com1DFA/DFAToolsCython.pyx | 7 +- avaframe/tests/test_DFAToolsCython.py | 99 +++++++++++++++++++++++++++ 2 files changed, 104 insertions(+), 2 deletions(-) diff --git a/avaframe/com1DFA/DFAToolsCython.pyx b/avaframe/com1DFA/DFAToolsCython.pyx index 63708e754..809e12d1c 100644 --- a/avaframe/com1DFA/DFAToolsCython.pyx +++ b/avaframe/com1DFA/DFAToolsCython.pyx @@ -573,8 +573,8 @@ cpdef (double, double, int, int, int, double, double, double, double) samosProje # are we in the same cell? if Lxi==Lxj and Lyi==Lyj: return xNew, yNew, iCell, Lx0, Ly0, w[0], w[1], w[2], w[3] - Lxj = Lxi - Lyj = Lyi + Lxi = Lxj + Lyi = Lyj return xNew, yNew, iCell, Lx0, Ly0, w[0], w[1], w[2], w[3] @@ -675,6 +675,9 @@ cpdef (double, double, double, int, int, int, double, double, double, double) di # second step: conserve the distance dist (correction ov the position) # measure distance between this new point and the previous time step position distn = norm(xNew-xPrev, yNew-yPrev, zNew-zPrev) + # guard against stopped particle (distn == 0) + if distn == 0: + break # adjust position on the Xprev, Xnew line xNew = xPrev + (xNew-xPrev) * dist / distn yNew = yPrev + (yNew-yPrev) * dist / distn diff --git a/avaframe/tests/test_DFAToolsCython.py b/avaframe/tests/test_DFAToolsCython.py index cbcdf0c34..54c4cd1d4 100644 --- a/avaframe/tests/test_DFAToolsCython.py +++ b/avaframe/tests/test_DFAToolsCython.py @@ -310,6 +310,105 @@ def test_reprojectionC(capfd): assert abs(dist-dist0) <= threshold*(dist0 + csz) +def test_samosProjectionConvergence(capfd): + """Test samosProjectionIteratrive convergence: result is stable across iteration counts.""" + ncols = 100 + nrows = 50 + csz = 5 + header = {} + header['ncols'] = ncols + header['nrows'] = nrows + header['cellsize'] = csz + dem = {} + dem['header'] = header + x0 = 300 + z0 = 500 + X = np.linspace(0, csz * (ncols - 1), ncols) + Y = np.linspace(0, csz * (nrows - 1), nrows) + XX, YY = np.meshgrid(X, Y) + ZZ = z0 / (x0 * x0) * (XX - x0) * (XX - x0) + dem['rasterData'] = ZZ + num = 1 + interpOption = 2 + dem = gT.getNormalMesh(dem, num=num) + Nx = dem['Nx'] + Ny = dem['Ny'] + Nz = dem['Nz'] + + # Point above the parabola surface, offset to trigger cell change during reprojection + xIn = 350 + yIn = 125 + Lx0, Ly0, iCell, w0, w1, w2, w3 = DFAfunC.getCellAndWeights( + xIn, yIn, ncols, nrows, csz, interpOption) + zSurf = DFAfunC.getScalar(Lx0, Ly0, w0, w1, w2, w3, ZZ) + zIn = zSurf + 20 + + # Reference: many iterations + xRef, yRef, iCellRef, LxR0, LyR0, wR0, wR1, wR2, wR3 = DFAfunC.samosProjectionIteratrive( + xIn, yIn, zIn, ZZ, Nx, Ny, Nz, csz, ncols, nrows, interpOption, 50) + zRef = DFAfunC.getScalar(LxR0, LyR0, wR0, wR1, wR2, wR3, ZZ) + + # With enough iterations, result matches the reference (non-regression) + for nIter in [29, 30, 50]: + xN, yN, iCellN, LxN0, LyN0, wN0, wN1, wN2, wN3 = DFAfunC.samosProjectionIteratrive( + xIn, yIn, zIn, ZZ, Nx, Ny, Nz, csz, ncols, nrows, interpOption, nIter) + zN = DFAfunC.getScalar(LxN0, LyN0, wN0, wN1, wN2, wN3, ZZ) + assert xN == pytest.approx(xRef, rel=1e-4) + assert yN == pytest.approx(yRef, rel=1e-4) + assert zN == pytest.approx(zRef, rel=1e-4) + + +def test_distConservStoppedParticle(capfd): + """distConservProjectionIteratrive handles distn=0 inside loop without NaN. + + prev is on a flat surface, current is above the surface at the same (x,y). + Initial distn is non-zero (current is above surface), so the loop enters. + After orthogonal reprojection, the particle lands on the surface and distn + becomes zero, triggering the division by zero in the distance conservation step.""" + ncols = 20 + nrows = 20 + csz = 5 + header = {} + header['ncols'] = ncols + header['nrows'] = nrows + header['cellsize'] = csz + dem = {} + dem['header'] = header + X = np.linspace(0, csz * (ncols - 1), ncols) + Y = np.linspace(0, csz * (nrows - 1), nrows) + XX, YY = np.meshgrid(X, Y) + # Flat surface so normal is straight up, no lateral shift + ZZ = np.zeros_like(XX) + dem['rasterData'] = ZZ + num = 1 + interpOption = 2 + dem = gT.getNormalMesh(dem, num=num) + Nx = dem['Nx'] + Ny = dem['Ny'] + Nz = dem['Nz'] + + xPrev = 50.0 + yPrev = 50.0 + zPrev = 0.0 # on flat surface + xCur = xPrev # same x, y — no lateral movement + yCur = yPrev + zCur = 20.0 # above surface + threshold = 0.001 + + xNew, yNew, zNew, iCell, Lx0, Ly0, w0, w1, w2, w3 = DFAfunC.distConservProjectionIteratrive( + xPrev, yPrev, zPrev, ZZ, Nx, Ny, Nz, + xCur, yCur, zCur, + csz, ncols, nrows, interpOption, 10, threshold) + + # Must not be NaN or inf (would be if division by zero is unguarded) + assert not np.isnan(xNew) + assert not np.isnan(yNew) + assert not np.isnan(zNew) + assert np.isfinite(xNew) + assert np.isfinite(yNew) + assert np.isfinite(zNew) + + def test_SamosATfric(capfd): """ Test the account4FrictionForce function""" tau0 = 0 From ff7594d8de791c600bcec40e1b4d498d183eb075 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Wed, 13 May 2026 08:12:13 +0200 Subject: [PATCH 3/5] fix(com1DFA): remap bond indices after particle removal removedBonds now builds old-to-new particle index mapping and remaps both bondStart (reindexed by new particle indices) and bondPart (target indices remapped). removePart skips bond keys since they are already correctly rebuilt by removedBonds. Previously bondPart stored original indices that became stale after np.extract compacted particle arrays, causing silent data corruption or out-of-bounds access when particles exited the domain. --- avaframe/com1DFA/DFAfunctionsCython.pyx | 28 ++++++++++++++++++++++- avaframe/com1DFA/particleTools.py | 3 ++- avaframe/tests/test_DFAfunctionsCython.py | 13 ++++++----- 3 files changed, 36 insertions(+), 8 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index b0b8c5d56..a7d5ccdfa 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -1872,7 +1872,33 @@ def removedBonds(particles, mask, nRemove, nBondRemove): countBond = countBond + 1 - particles['bondStart'] = np.asarray(bondStartNew).copy() + # Remap bondStart and bondPart to use compacted (new) particle indices + # After removePart compacts particle arrays via np.extract, old indices become stale + cdef int[:] oldToNew = np.full(nPart, -1, dtype=np.int32) + cdef int newIdx = 0 + cdef int kk + for kk in range(nPart): + if keepParticle[kk] == 1: + oldToNew[kk] = newIdx + newIdx = newIdx + 1 + cdef int nPartNew = newIdx + + # Rebuild bondStart indexed by new particle indices + cdef int[:] bondStartFinal = np.zeros(nPartNew + 1, dtype=np.int32) + cdef int bondCumul = 0 + newIdx = 0 + for kk in range(nPart): + if keepParticle[kk] == 1: + bondCumul = bondCumul + (bondStartNew[kk + 1] - bondStartNew[kk]) + bondStartFinal[newIdx + 1] = bondCumul + newIdx = newIdx + 1 + + # Remap bondPartNew target indices from old to new + cdef int ibb + for ibb in range(countBondNew): + bondPartNew[ibb] = oldToNew[bondPartNew[ibb]] + + particles['bondStart'] = np.asarray(bondStartFinal).copy() particles['bondPart'] = np.asarray(bondPartNew).copy() particles['bondDist'] = np.asarray(bondDistNew).copy() return particles diff --git a/avaframe/com1DFA/particleTools.py b/avaframe/com1DFA/particleTools.py index a8b13fbed..6676407f0 100644 --- a/avaframe/com1DFA/particleTools.py +++ b/avaframe/com1DFA/particleTools.py @@ -229,7 +229,8 @@ def removePart(particles, mask, nRemove, reasonString="", snowSlide=0): for key in particles: if key == "nPart": particles["nPart"] = particles["nPart"] - nRemove - elif key == "stoppedParticles": + elif key == "stoppedParticles" or key == "bondStart" or key == "bondPart" or key == "bondDist": + # bond arrays are already rebuilt by removedBonds with new-indexed data continue # for all keys in particles that are arrays of size nPart do: elif type(particles[key]).__module__ == np.__name__: diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 97e9c3b97..b17678115 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -669,15 +669,18 @@ def test_removeBondsC(): bondStart = particles["bondStart"] bondDist = particles["bondDist"] bondPart = particles["bondPart"] - assert np.array_equal(bondStart, np.asarray([0, 1, 1, 2])) + # After removing particle 1 (old index), indices are remapped: + # old 0 -> new 0, removed, old 2 -> new 1 + assert np.array_equal(bondStart, np.asarray([0, 1, 2])) assert particles["bondStart"].flags.owndata assert particles["bondStart"].base is None assert particles["bondPart"].flags.owndata assert particles["bondPart"].base is None assert particles["bondDist"].flags.owndata assert particles["bondDist"].base is None - for k in range(nPart): - # loop on all bonded particles + # Only 2 particles remain; bondStart is new-indexed + nPartNew = nPart - nRemove + for k in range(nPartNew): neighbors = list() for ib in range(bondStart[k], bondStart[k + 1]): l = bondPart[ib] @@ -685,10 +688,8 @@ def test_removeBondsC(): neighbors.sort() if k == 0: - assert neighbors == [2] + assert neighbors == [1] if k == 1: - assert neighbors == [] - if k == 2: assert neighbors == [0] bondDist.sort() assert np.array_equal(bondDist, np.asarray([1, 1])) From 1793088be5029117e0884167edc4a564ac205a01 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Wed, 13 May 2026 08:23:42 +0200 Subject: [PATCH 4/5] fix(com1DFA): remove circular reference in stoppedParticles particles["stoppedParticles"] = particles created a self-referencing dict that overwrote the correct stopped-particle arrays (x, y, m, ID) with the full flowing-particle data. Downstream in updateFieldsC this caused stopped mass rasters to use flowing instead of stopped masses. Also prevented garbage collection and broke serialization. Removed the offending line. stoppedParticles is already correctly populated earlier in the function, and massStopped captures the accumulated total. --- avaframe/com1DFA/DFAfunctionsCython.pyx | 1 - avaframe/tests/test_DFAfunctionsCython.py | 106 ++++++++++++++++++++++ 2 files changed, 106 insertions(+), 1 deletion(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index a7d5ccdfa..568eb1268 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -1115,7 +1115,6 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): if stop: particles['iterate'] = False - particles['stoppedParticles'] = particles massStopped = massStopped + massFlowing if typeStop == 1: diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index b17678115..7ff816dfb 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -582,6 +582,112 @@ def test_updatePositionC(): assert particles["massStopped"] == -np.sum(particles["stoppedParticles"]["m"]) +def test_updatePositionC_no_circular_reference(): + """updatePositionC must not create circular reference in stoppedParticles. + + Uses stopCrit=1.0 with typeStop=0 and stopCritType=kinEnergy so that the stop + criterion is met on the first call. This exercises the `if stop:` code path + where the circular reference was previously created.""" + cfg = configparser.ConfigParser() + cfg["GENERAL"] = { + "stopCrit": "1.01", + "stopCritIni": "0.1", + "stopCritIniSmall": "1.001", + "stopCritType": "kinEnergy", + "uFlowingThreshold": "0.1", + "gravAcc": "9.81", + "velMagMin": "1.e-6", + "rho": "100.", + "interpOption": "2", + "explicitFriction": "0", + "centeredPosition": "1", + "reprojMethodPosition": "2", + "reprojectionIterations": "5", + "thresholdProjection": "0.001", + "dissDam": "1", + "snowSlide": "1", + "wetSnow": "1", + "adaptSfcStopped": "0", + "adaptSfcDetrainment": "0", + "adaptSfcEntrainment": "0", + } + + particles = { + "dt": 1.0, + "m": np.asarray([10.0, 10.0, 10.0], dtype=np.float64), + "idFixed": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "trajectoryLengthXY": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "trajectoryLengthXYCor": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "trajectoryLengthXYZ": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "x": np.asarray([0.0, 1.0, 2.0], dtype=np.float64), + "y": np.asarray([2.0, 3.0, 4.0], dtype=np.float64), + "z": np.asarray([1.0, 1.0, 1.0], dtype=np.float64), + "ux": np.asarray([1.0, 1.0, 1.0], dtype=np.float64), + "uy": np.asarray([1.0, 1.0, 1.0], dtype=np.float64), + "uz": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "uAcc": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "kineticEne": 0.0, + "peakKinEne": 0.0, + "peakForceSPH": 0.0, + "forceSPHIni": 0.0, + "nPart": 3, + "peakMassFlowing": 0.0, + "iterate": True, + "totalEnthalpy": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "velocityMag": np.asarray([1.0, 1.0, 1.0], dtype=np.float64), + "h": np.asarray([1.0, 1.0, 1.0], dtype=np.float64), + "ID": np.array([0, 1, 2], dtype=np.int64), + "stoppedParticles": {}, + "indXDEM": np.array([0], dtype=np.int32), + "indYDEM": np.array([0], dtype=np.int32), + } + particles["potentialEne"] = np.sum(9.81 * particles["z"] * particles["m"]) + + demHeader = {} + demHeader["xllcenter"] = 0.0 + demHeader["yllcenter"] = 0.0 + demHeader["cellsize"] = 5.0 + demHeader["nodata_value"] = -9999 + demHeader["nrows"] = 10 + demHeader["ncols"] = 10 + dem = {"header": demHeader} + dem["rasterData"] = np.ones((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64) + dem["outOfDEM"] = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64).flatten() + dem["Nx"] = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64) + dem["Ny"] = np.zeros((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64) + dem["Nz"] = np.ones((demHeader["nrows"], demHeader["ncols"]), dtype=np.float64) + + force = { + "forceZ": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "forceFrict": np.asarray([10.0, 10.0, 10.0], dtype=np.float64), + "dM": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "dMDet": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + "forceX": np.asarray([50.0, 50.0, 50.0], dtype=np.float64), + "forceY": np.asarray([50.0, 50.0, 50.0], dtype=np.float64), + "forceSPHX": np.asarray([50.0, 50.0, 50.0], dtype=np.float64), + "forceSPHY": np.asarray([50.0, 50.0, 50.0], dtype=np.float64), + "forceSPHZ": np.asarray([0.0, 0.0, 0.0], dtype=np.float64), + } + fields = {"FT": np.zeros((2, 2))} + wallLineDict = { + "dam": 0, + "nIterDam": 1, + "cellsCrossed": np.zeros((dem["header"]["ncols"] * dem["header"]["nrows"]), dtype=np.int64), + } + for key in ["x", "y", "z", "xCrown", "yCrown", "zCrown", + "xTangent", "yTangent", "zTangent"]: + wallLineDict[key] = np.ones((1)) * 1.0 + for key in ["nPoints", "height", "slope", "restitutionCoefficient"]: + wallLineDict[key] = 0 + dem["damLine"] = wallLineDict + + particles = DFAfunC.updatePositionC(cfg["GENERAL"], particles, dem, force, fields, typeStop=0) + + # Verify no circular reference was created + assert particles["stoppedParticles"] is not particles + assert "stoppedParticles" not in particles["stoppedParticles"] + + def test_computeTrajectoryAngle(): # first compute travel angle for each particle # get parent Id in order to get the first z position From 3933b0458a434645a36de19c2ce29a2aafbbf140 Mon Sep 17 00:00:00 2001 From: Felix Oesterle <6945681+fso42@users.noreply.github.com> Date: Wed, 13 May 2026 09:05:54 +0200 Subject: [PATCH 5/5] fix(tests): correct array initialization dimensions in `test_DFAfunctionsCython` Updated array initializations (`pfv`, `ppr`, `pft`, `pta`, `pke`) to use dynamic dimensions based on `header["nrows"]` and `header["ncols"]` instead of static `(1, 1)` arrays. This ensures compatibility with test inputs of varying sizes. --- avaframe/tests/test_DFAfunctionsCython.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 7ff816dfb..123ede57c 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -825,11 +825,11 @@ def test_updateFieldsC_returns_owned_arrays(): "computeTA": False, "computeKE": False, "computeP": False, - "pfv": np.zeros((1, 1)), - "ppr": np.zeros((1, 1)), - "pft": np.zeros((1, 1)), - "pta": np.zeros((1, 1)), - "pke": np.zeros((1, 1)), + "pfv": np.zeros((header["nrows"], header["ncols"])), + "ppr": np.zeros((header["nrows"], header["ncols"])), + "pft": np.zeros((header["nrows"], header["ncols"])), + "pta": np.zeros((header["nrows"], header["ncols"])), + "pke": np.zeros((header["nrows"], header["ncols"])), "dmDet": np.zeros((header["nrows"], header["ncols"])), }