From 0677684f8527bb1bec24cfb9a369039dec88ccd9 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 31 Dec 2022 10:46:25 +0100 Subject: [PATCH 1/2] Adjusting weights for collisions Ionization and recombination collisions have been modified to have the right products accounting for the possibility that primary electron and target particle have different weight. --- src/modules/moduleCollisions.f90 | 93 ++++++++++++++++++++++++-------- 1 file changed, 70 insertions(+), 23 deletions(-) diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 3c241c1..702b54a 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,39 @@ 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) + PRINT *, "ionize" + + 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,27 +368,13 @@ 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 + !Adds new particles to the list CALL OMP_SET_LOCK(lockCollisions) CALL partCollisions%add(newElectron) + IF (ASSOCIATED(remainingNeutral)) THEN + CALL partCollisions%add(remainingNeutral) + + END IF CALL OMP_UNSET_LOCK(lockCollisions) END IF @@ -399,7 +418,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 +442,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 +461,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 +493,13 @@ MODULE moduleCollisions END SELECT + !Adds new particles to the list + IF (ASSOCIATED(remainingIon)) THEN + CALL OMP_SET_LOCK(lockCollisions) + CALL partCollisions%add(remainingIon) + CALL OMP_UNSET_LOCK(lockCollisions) + END IF + END SUBROUTINE collideBinaryRecombination !RESONANT CHARGE EXCHANGE From 8199a228c8380c30db73c303fd183f32e4bc87a1 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 31 Dec 2022 11:22:02 +0100 Subject: [PATCH 2/2] Locks for particle lists are now inside the type. The lock of a particle list is no longer an external variable, it is now part of the type. New procedures have been added to set and unset the lock. --- src/modules/init/moduleInput.f90 | 8 +++++-- src/modules/mesh/moduleMeshBoundary.f90 | 4 ++-- src/modules/moduleCollisions.f90 | 9 ++++---- src/modules/moduleList.f90 | 30 ++++++++++++++++++++----- src/modules/solver/moduleSolver.f90 | 4 ++-- 5 files changed, 39 insertions(+), 16 deletions(-) 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 702b54a..ba1edfd 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -347,7 +347,6 @@ MODULE moduleCollisions ELSEIF (electron%weight < neutral%weight) THEN !If primary electron is ligther than neutral, change weight of neutral and create new neutral ALLOCATE(remainingNeutral) - PRINT *, "ionize" remainingNeutral = neutral @@ -369,13 +368,13 @@ MODULE moduleCollisions END SELECT !Adds new particles to the list - CALL OMP_SET_LOCK(lockCollisions) + CALL partCollisions%setLock() CALL partCollisions%add(newElectron) IF (ASSOCIATED(remainingNeutral)) THEN CALL partCollisions%add(remainingNeutral) END IF - CALL OMP_UNSET_LOCK(lockCollisions) + CALL partCollisions%unsetLock() END IF @@ -495,9 +494,9 @@ MODULE moduleCollisions !Adds new particles to the list IF (ASSOCIATED(remainingIon)) THEN - CALL OMP_SET_LOCK(lockCollisions) + CALL partCollisions%setLock() CALL partCollisions%add(remainingIon) - CALL OMP_UNSET_LOCK(lockCollisions) + CALL partCollisions%unsetLock() END IF END SUBROUTINE collideBinaryRecombination 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