[com1] add checking for zero mass of particle and remove properly

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')
@@ -761,6 +760,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 +784,12 @@ 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:
+      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
     # operator splitting
     # estimate new velocity due to driving force
@@ -804,9 +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 or mass = 0
-    if delStoppedParticles == 1:
-      if uMagNew == 0 or mNew == 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])
@@ -1046,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 and delStoppedParticles == 1:
+  # 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
@@ -1063,6 +1068,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
 
 

avaframe/com1DFA/com1DFACfg.ini

@@ -386,11 +386,10 @@ 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
-# 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

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

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`.
 

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