Merge branch 'issue/collisionWeight' into 'development'

Adjusting weights for collisions

See merge request JorgeGonz/fpakc!32

src/modules/init/moduleInput.f90

@@ -612,7 +612,7 @@ MODULE moduleInput
       nPartOld = 0
 
       !Initialize the lock for the non-analogue (NA) list of particles
-      CALL OMP_INIT_LOCK(lockWScheme)
+      CALL OMP_INIT_LOCK(partWScheme%lock)
 
     END SUBROUTINE readSpecies
 
@@ -676,7 +676,7 @@ MODULE moduleInput
           CALL config%get('interactions.folderCollisions', pathCollisions, found)
 
           !Inits lock for list of particles
-          CALL OMP_INIT_LOCK(lockCollisions)
+          CALL OMP_INIT_LOCK(partCollisions%lock)
 
           CALL config%info('interactions.collisions', found, n_children = nPairs)
           DO i = 1, nPairs
@@ -770,6 +770,7 @@ MODULE moduleInput
       USE moduleErrors
       USE moduleSpecies
       USE moduleRefParam
+      USE moduleList, ONLY: partSurfaces
       USE json_module
       IMPLICIT NONE
 
@@ -862,6 +863,9 @@ MODULE moduleInput
 
       END DO
 
+      !Init the list of particles from surfaces
+      CALL OMP_INIT_LOCK(partSurfaces%lock)
+
     END SUBROUTINE readBoundary
 
     !Read the geometry (mesh) for the case

src/modules/mesh/moduleMeshBoundary.f90

@@ -169,10 +169,10 @@ MODULE moduleMeshBoundary
             newIon%n_in      = .TRUE.
 
             !Add particles to list
-            CALL OMP_SET_LOCK(lockSurfaces)
+            CALL partSurfaces%setLock()
             CALL partSurfaces%add(newElectron)
             CALL partSurfaces%add(newIon)
-            CALL OMP_UNSET_LOCK(lockSurfaces)
+            CALL partSurfaces%unsetLock()
 
             !Electron loses energy due to ionization
             eRel = eRel - bound%eThreshold

src/modules/moduleCollisions.f90

@@ -305,7 +305,7 @@ MODULE moduleCollisions
       REAL(8):: rMass, eRel
       TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
       REAL(8), DIMENSION(1:3):: vChange
-      TYPE(particle), POINTER:: newElectron
+      TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL()
 
       rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
       eRel = rMass*vRel**2
@@ -324,6 +324,38 @@ MODULE moduleCollisions
 
         END IF
 
+        !Exchange of velocity between particles
+        vChange = self%deltaV*randomDirectionVHS()
+        
+        !Energy is loss by the primary electron
+        electron%v = electron%v - vChange
+
+        !Creates a new electron from ionization
+        ALLOCATE(newElectron)
+
+        !Copy basic information from primary electron
+        newElectron = electron
+
+        !Secondary electorn gains energy from ionization
+        newElectron%v       =  vChange
+
+        !Correct the weight of the particles
+        IF (electron%weight >= neutral%weight) THEN
+          !If primary electron is hevier than neutral, reduce weight to secondary electron
+          newElectron%weight = neutral%weight
+
+        ELSEIF (electron%weight < neutral%weight) THEN
+          !If primary electron is ligther than neutral, change weight of neutral and create new neutral
+          ALLOCATE(remainingNeutral)
+
+          remainingNeutral = neutral
+
+          remainingNeutral%weight = neutral%weight - electron%weight
+
+          neutral%weight = electron%weight
+
+        END IF
+
         !Ionize neutral particle
         SELECT TYPE(sp => neutral%species)
         TYPE IS(speciesNeutral)
@@ -335,28 +367,14 @@ MODULE moduleCollisions
 
         END SELECT
 
-        !Exchange of velocity between particles
+        !Adds new particles to the list
-        vChange = self%deltaV*randomDirectionVHS()
+        CALL partCollisions%setLock()
-        
-        !Energy is loss by the primary electron
-        electron%v = electron%v - vChange
-
-        !Creates a new electron from ionization
-        ALLOCATE(newElectron)
-        newElectron%species => electron%species
-        !Secondary electorn gains energy from ionization
-        newElectron%v       =  vChange
-        newElectron%r       =  neutral%r
-        newElectron%xi      =  neutral%xi
-        newElectron%n_in    =  .TRUE.
-        newElectron%vol     =  neutral%vol
-        newElectron%volColl =  neutral%volColl
-        newElectron%weight  =  neutral%weight
-
-        !Adds new electron to list of new particles from collisions
-        CALL OMP_SET_LOCK(lockCollisions)
         CALL partCollisions%add(newElectron)
-        CALL OMP_UNSET_LOCK(lockCollisions)
+        IF (ASSOCIATED(remainingNeutral)) THEN
+          CALL partCollisions%add(remainingNeutral)
+
+        END IF
+        CALL partCollisions%unsetLock()
 
       END IF
 
@@ -399,7 +417,8 @@ MODULE moduleCollisions
           collision%electron => sp
 
         CLASS DEFAULT
-          CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization')
+          CALL criticalError("Species " // sp%name // " chosen for recombination is not a charged species", &
+                             'initBinaryRecombination')
 
         END SELECT
 
@@ -422,6 +441,7 @@ MODULE moduleCollisions
       TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
       REAL(8):: sigmaVrel
       REAL(8), DIMENSION(1:3):: vp_i
+      TYPE(particle), POINTER:: remainingIon => NULL()
 
       IF    (ASSOCIATED(part_i%species, self%electron)) THEN
         electron => part_i
@@ -440,9 +460,28 @@ MODULE moduleCollisions
       !TODO: This energy should be transformed into photons
       vp_i = ion%v*     (1.D0 -  (vRel + self%deltaV)/NORM2(ion%v))
 
+      IF (electron%weight > ion%weight) THEN
+        !Reduce weight of primary electron but particle continues
+        electron%weight = electron%weight - ion%weight
+
+      ELSE
         !Remove electron from simulation
         electron%n_in = .FALSE.
 
+        !There is some ion remaining
+        IF (electron%weight < ion%weight) THEN
+          ALLOCATE(remainingIon)
+
+          remainingIon = ion
+
+          remainingIon%weight = ion%weight - electron%weight
+
+          ion%weight = electron%weight
+
+        END IF
+
+      END IF
+
       !Neutralize ion particle
       SELECT TYPE(sp => ion%species)
       TYPE IS(speciesCharged)
@@ -453,6 +492,13 @@ MODULE moduleCollisions
 
       END SELECT
 
+      !Adds new particles to the list
+      IF (ASSOCIATED(remainingIon)) THEN
+        CALL partCollisions%setLock()
+        CALL partCollisions%add(remainingIon)
+        CALL partCollisions%unsetLock()
+      END IF
+
     END SUBROUTINE collideBinaryRecombination
 
     !RESONANT CHARGE EXCHANGE

src/modules/moduleList.f90

@@ -13,19 +13,19 @@ MODULE moduleList
     INTEGER:: amount = 0
     TYPE(lNode),POINTER:: head => NULL()
     TYPE(lNode),POINTER:: tail => NULL()
+    INTEGER(KIND=OMP_LOCK_KIND):: lock
     CONTAINS
       PROCEDURE,PASS:: add => addToList
       PROCEDURE,PASS:: convert2Array
       PROCEDURE,PASS:: erase => eraseList
+      PROCEDURE,PASS:: setLock
+      PROCEDURE,PASS:: unsetLock
 
   END TYPE listNode
 
   TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme
-  INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles
   TYPE(listNode):: partCollisions !Particles created in collisional process
-  INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles
   TYPE(listNode):: partSurfaces !Particles created in surface interactions
-  INTEGER(KIND=OMP_LOCK_KIND):: lockSurfaces !Lock for the NA list of particles
   TYPE(listNode):: partInitial !Initial distribution of particles
 
   TYPE pointerArray
@@ -80,8 +80,10 @@ MODULE moduleList
     END FUNCTION convert2Array
 
     !Erase list
-    SUBROUTINE eraseList(self)
+    PURE SUBROUTINE eraseList(self)
-      CLASS(listNode):: self
+      IMPLICIT NONE
+
+      CLASS(listNode), INTENT(inout):: self
       TYPE(lNode),POINTER:: current, next
       
       current => self%head
@@ -95,4 +97,22 @@ MODULE moduleList
       self%amount = 0
     END SUBROUTINE eraseList
 
+    SUBROUTINE setLock(self)
+      USE OMP_LIB
+      IMPLICIT NONE
+
+      CLASS(listNode):: self
+      CALL OMP_SET_LOCK(self%lock)
+
+    END SUBROUTINE setLock
+
+    SUBROUTINE unsetLock(self)
+      USE OMP_LIB
+      IMPLICIT NONE
+
+      CLASS(listNode):: self
+      CALL OMP_UNSET_LOCK(self%lock)
+
+    END SUBROUTINE unsetLock
+
 END MODULE moduleList

src/modules/solver/moduleSolver.f90

@@ -434,9 +434,9 @@ MODULE moduleSolver
         !Copy data from original particle
         newPart = part
         !Add particle to list of new particles from weighting scheme
-        CALL OMP_SET_LOCK(lockWScheme)
+        CALL partWScheme%setLock()
         CALL partWScheme%add(newPart)
-        CALL OMP_UNSET_LOCK(lockWScheme)
+        CALL partWScheme%unsetLock()
         !Add particle to cell list
         CALL OMP_SET_lock(vol%lock)
         sp = part%species%n