diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 2523b3b..ac40c3b 100644 --- a/src/modules/init/moduleInput.f90 +++ b/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 diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index c5f3fc3..2e15c5c 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/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 diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 3c241c1..ba1edfd 100644 --- a/src/modules/moduleCollisions.f90 +++ b/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 - vChange = self%deltaV*randomDirectionVHS() - - !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) + !Adds new particles to the list + CALL partCollisions%setLock() 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,8 +460,27 @@ MODULE moduleCollisions !TODO: This energy should be transformed into photons vp_i = ion%v* (1.D0 - (vRel + self%deltaV)/NORM2(ion%v)) - !Remove electron from simulation - electron%n_in = .FALSE. + 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) @@ -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 diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index cf27388..b040c80 100644 --- a/src/modules/moduleList.f90 +++ b/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) - CLASS(listNode):: self + PURE SUBROUTINE eraseList(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 diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 88f7bf8..135b08f 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/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