From c0c447faf944fcb5b1d616bad2d0cdff8e2eaad4 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Mon, 8 Sep 2025 10:38:01 +0200 Subject: [PATCH 1/3] add checking for zero mass of particle and remove properly --- avaframe/com1DFA/DFAfunctionsCython.pyx | 22 ++++++++++++++++++---- avaframe/com1DFA/com1DFACfg.ini | 2 +- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 4afcad126..d48e6b190 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -761,6 +761,7 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): cdef int Lx0, Ly0, LxNew0, LyNew0, iCell, iCellNew cdef double w[4] cdef double wNew[4] + # loop on particles for k in range(nPart): m = mass[k] @@ -784,6 +785,19 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): uMag = DFAtlsC.norm(ux, uy, uz) uMagt0 = DFAtlsC.norm(ux, uy, uz) + # check if particle's mass is zero then remove particle + if m == 0: + # if the particle's mass is zero remove particle but keep info in separate stoppedParticles dict + xStoppedArray = np.append(xStoppedArray, xArray[k]) + yStoppedArray = np.append(yStoppedArray, yArray[k]) + mStoppedArray = np.append(mStoppedArray, mass[k]) + idStoppedArray = np.append(idStoppedArray, ID[k]) + uMagStoppedArray = np.append(uMagStoppedArray, uMagNew) + massStopped = massStopped + mass[k] + notStopParticle[k] = 0 # particle is deleted + nStop = nStop + 1 + continue # this particle will be removed, skip what is below and directly go to the next particle + # procede to time integration # operator splitting # estimate new velocity due to driving force @@ -804,9 +818,8 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): massEntrained = massEntrained + dM[k] massDetrained = massDetrained + dMDet[k] - # Stopped particles with velocity = 0 or mass = 0 - if delStoppedParticles == 1: - if uMagNew == 0 or mNew == 0: + # Stopped particles with velocity = 0 + if delStoppedParticles == 1 and uMagNew == 0: xStoppedArray = np.append(xStoppedArray, xArray[k]) yStoppedArray = np.append(yStoppedArray, yArray[k]) mStoppedArray = np.append(mStoppedArray, mass[k]) @@ -1052,7 +1065,7 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): particles = particleTools.removePart(particles, mask, nRemove, 'because they exited the domain', snowSlide=snowSlide) # remove particles that have mass = 0 or velocity = 0 - if nStop > 0 and delStoppedParticles == 1: + if nStop > 0: indRemoveParticle = np.array([], dtype=np.int64) for k in range(len(keepParticle)): # consider particles that are removed because they exit the domain @@ -1063,6 +1076,7 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): notStopParticle = np.delete(notStopParticle, indRemoveParticle) mask = np.array(np.asarray(notStopParticle), dtype=bool) particles = particleTools.removePart(particles, mask, nStop, 'because their mass or velocity is zero', snowSlide=snowSlide) + return particles diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index d1dceab93..b581311c0 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -386,7 +386,7 @@ entShearResistance = 0 entDefResistance = 0 #++++++++++++ Deposition, Erosion and adaptive surface -# delete stopped particles (velocity = 0, mass = 0 or stopping criterion) but save in a additional dictionary +# delete stopped particles (velocity = 0 or stopping criterion) but save in an additional dictionary # activate with 1 delStoppedParticles = 0 # adapt the topography every time step From 06606ec0189eedb9a78ff3890b2c7e8a0cc08a6c Mon Sep 17 00:00:00 2001 From: Anna Wirbel <68687423+awirb@users.noreply.github.com> Date: Mon, 8 Sep 2025 13:23:07 +0200 Subject: [PATCH 2/3] Update avaframe/com1DFA/DFAfunctionsCython.pyx Co-authored-by: Paula Spannring <95042192+PaulaSp3@users.noreply.github.com> --- avaframe/com1DFA/DFAfunctionsCython.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index d48e6b190..15b116e23 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -792,7 +792,7 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): yStoppedArray = np.append(yStoppedArray, yArray[k]) mStoppedArray = np.append(mStoppedArray, mass[k]) idStoppedArray = np.append(idStoppedArray, ID[k]) - uMagStoppedArray = np.append(uMagStoppedArray, uMagNew) + uMagStoppedArray = np.append(uMagStoppedArray, 0.0) massStopped = massStopped + mass[k] notStopParticle[k] = 0 # particle is deleted nStop = nStop + 1 From f46156313065c1c6c1593ff62b02767a20e69ca8 Mon Sep 17 00:00:00 2001 From: Anna Wirbel Date: Tue, 9 Sep 2025 14:16:55 +0200 Subject: [PATCH 3/3] remove from stopped part dict the detrained ones stopped particles are always and only deleted and stored in stoppedParticles when surface is adapted due to stopping little comment Update avaframe/com1DFA/DFAfunctionsCython.pyx Co-authored-by: Anna Wirbel <68687423+awirb@users.noreply.github.com> --- avaframe/com1DFA/DFAfunctionsCython.pyx | 22 +++++++--------------- avaframe/com1DFA/com1DFACfg.ini | 5 ++--- avaframe/tests/test_DFAfunctionsCython.py | 3 +-- docs/com1DFAAlgorithm.rst | 2 +- docs/theoryCom1DFA.rst | 3 +++ 5 files changed, 14 insertions(+), 21 deletions(-) diff --git a/avaframe/com1DFA/DFAfunctionsCython.pyx b/avaframe/com1DFA/DFAfunctionsCython.pyx index 15b116e23..0877190d7 100644 --- a/avaframe/com1DFA/DFAfunctionsCython.pyx +++ b/avaframe/com1DFA/DFAfunctionsCython.pyx @@ -660,7 +660,6 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): cdef double thresholdProjection = cfg.getfloat('thresholdProjection') cdef double centeredPosition = cfg.getfloat('centeredPosition') cdef int snowSlide = cfg.getint('snowSlide') - cdef int delStoppedParticles = cfg.getint('delStoppedParticles') cdef int adaptSfcStopped = cfg.getint('adaptSfcStopped') cdef int adaptSfcDetrainment = cfg.getint('adaptSfcDetrainment') cdef int adaptSfcEntrainment = cfg.getint('adaptSfcEntrainment') @@ -787,15 +786,8 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): # check if particle's mass is zero then remove particle if m == 0: - # if the particle's mass is zero remove particle but keep info in separate stoppedParticles dict - xStoppedArray = np.append(xStoppedArray, xArray[k]) - yStoppedArray = np.append(yStoppedArray, yArray[k]) - mStoppedArray = np.append(mStoppedArray, mass[k]) - idStoppedArray = np.append(idStoppedArray, ID[k]) - uMagStoppedArray = np.append(uMagStoppedArray, 0.0) - massStopped = massStopped + mass[k] - notStopParticle[k] = 0 # particle is deleted - nStop = nStop + 1 + keepParticle[k] = 0 + nRemove = nRemove + 1 continue # this particle will be removed, skip what is below and directly go to the next particle # procede to time integration @@ -818,8 +810,8 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): massEntrained = massEntrained + dM[k] massDetrained = massDetrained + dMDet[k] - # Stopped particles with velocity = 0 - if delStoppedParticles == 1 and uMagNew == 0: + # Stop particles with velocity = 0 + if adaptSfcStopped == 1 and uMagNew == 0: xStoppedArray = np.append(xStoppedArray, xArray[k]) yStoppedArray = np.append(yStoppedArray, yArray[k]) mStoppedArray = np.append(mStoppedArray, mass[k]) @@ -1059,13 +1051,13 @@ def updatePositionC(cfg, particles, dem, force, fields, int typeStop=0): particles['massStopped'] = - massStopped - # remove particles that are not located on the mesh any more + # remove particles that are not located on the mesh any more OR have zero mass if nRemove > 0: mask = np.array(np.asarray(keepParticle), dtype=bool) particles = particleTools.removePart(particles, mask, nRemove, 'because they exited the domain', snowSlide=snowSlide) - # remove particles that have mass = 0 or velocity = 0 - if nStop > 0: + # remove particles that have velocity = 0 (only when surface is adapted due to stopping) + if nStop > 0 and adaptSfcStopped == 1: indRemoveParticle = np.array([], dtype=np.int64) for k in range(len(keepParticle)): # consider particles that are removed because they exit the domain diff --git a/avaframe/com1DFA/com1DFACfg.ini b/avaframe/com1DFA/com1DFACfg.ini index b581311c0..58391766b 100644 --- a/avaframe/com1DFA/com1DFACfg.ini +++ b/avaframe/com1DFA/com1DFACfg.ini @@ -386,11 +386,10 @@ entShearResistance = 0 entDefResistance = 0 #++++++++++++ Deposition, Erosion and adaptive surface -# delete stopped particles (velocity = 0 or stopping criterion) but save in an additional dictionary -# activate with 1 -delStoppedParticles = 0 # adapt the topography every time step # activate with 1 +# if adaptSfcStopped = 1: stopped particles (velocity = 0 or stopping criterion is reached) +# are deleted and saved in an additional dictionary (particles["stoppedParticles"]) adaptSfcStopped = 0 adaptSfcDetrainment = 0 adaptSfcEntrainment = 0 diff --git a/avaframe/tests/test_DFAfunctionsCython.py b/avaframe/tests/test_DFAfunctionsCython.py index 5567d80dc..ede4dc84e 100644 --- a/avaframe/tests/test_DFAfunctionsCython.py +++ b/avaframe/tests/test_DFAfunctionsCython.py @@ -211,7 +211,6 @@ def test_updatePositionC(): "dissDam": "1", "snowSlide": "1", "wetSnow": "1", - "delStoppedParticles": "0", "adaptSfcStopped": "0", "adaptSfcDetrainment": "0", "adaptSfcEntrainment": "0", @@ -537,7 +536,7 @@ def test_updatePositionC(): "forceSPHZ": np.asarray([0.0, 0.0, 0.0]), } - cfg["GENERAL"]["delStoppedParticles"] = "1" + cfg["GENERAL"]["adaptSfcStopped"] = "1" cfg["GENERAL"]["explicitFriction"] = "1" cfg["GENERAL"]["snowSlide"] = "0" typeStop = 1 diff --git a/docs/com1DFAAlgorithm.rst b/docs/com1DFAAlgorithm.rst index cfe1f20d4..4c2321094 100644 --- a/docs/com1DFAAlgorithm.rst +++ b/docs/com1DFAAlgorithm.rst @@ -168,7 +168,7 @@ Other particles properties are also initialized here: - ``dmEnt`` - entrained mass of particle [kg] - - ``stoppedParticles`` - dictionary with particles (containing x-, y-coordinates, mass and ID) that are stopped (mass or velocity is zero) and deleted from the particles in each time step + - ``stoppedParticles`` - dictionary with particles (containing x-, y-coordinates, mass and ID) that are stopped (velocity is zero) and deleted from the particles in each time step (only if ``adaptSfcStopped`` is set to ``1`` in the configuration file) For more details, see :py:func:`com1DFA.com1DFA.initializeParticles`. diff --git a/docs/theoryCom1DFA.rst b/docs/theoryCom1DFA.rst index 7406cae8b..c25129259 100644 --- a/docs/theoryCom1DFA.rst +++ b/docs/theoryCom1DFA.rst @@ -254,6 +254,9 @@ The adapted surface at a specific time step :math:`z(t)` is computed as (with :m In every time step the surface is adapted, the :ref:`DFAnumerics:Cell normals` and :ref:`DFAnumerics:Cell area` are also adapted. +If the surface is adapted due to stopped particles (velocity = 0 and ``adaptSfcStopped = 1``), the stopped particles +are deleted and stored in a separate dictionary ``stoppedParticles`` (see :ref:`com1DFAAlgorithm:Particle properties`). + .. Note:: This feature requires more detailed testing and may cause numerical problems if the change in snow