First implementation of Non-Analogue Scheme using volume weighting. The

scheme to use is chosen in the input file. Additional schemes could be
added easily.

runs/cylFlow/input.json

@@ -11,7 +11,7 @@
     "meshFile": "mesh.msh"
   },
   "species": [
-    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 5.0e7}
+    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 5.0e8}
   ],
   "boundary": [ 
     {"name": "Injection",       "type": "absorption", "physicalSurface": 1},
@@ -32,7 +32,8 @@
   "case": {
     "tau":  [6.0e-7],
     "time":  3.0e-3,
-    "pusher": ["2DCylNeutral"]
+    "pusher": ["2DCylNeutral"],
+    "NAScheme": "Volume"
   },
   "interactions": {
     "folderCollisions": "./data/collisions/",

src/modules/makefile

@@ -37,10 +37,10 @@ moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95
 moduleRefParam.o: moduleConstParam.o moduleRefParam.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
-moduleSpecies.o: moduleCaseParam.o moduleSpecies.f95
+moduleSpecies.o: moduleErrors.o moduleCaseParam.o moduleSpecies.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
-moduleSolver.o: moduleEM.o moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
+moduleSolver.o: moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleTable.o: moduleErrors.o moduleTable.f95

src/modules/moduleInject.f95

@@ -13,8 +13,8 @@ MODULE moduleInject
     INTEGER:: sp !Species of injection
     INTEGER:: nEdges
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
-    REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection
+    ! REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection
-    REAL(8):: sumWeight
+    ! REAL(8):: sumWeight
     CONTAINS
       PROCEDURE, PASS:: init         => initInject
       PROCEDURE, PASS:: addParticles => addParticlesMaxwellian
@@ -79,22 +79,22 @@ MODULE moduleInject
 
       self%nEdges = COUNT(phSurface == physicalSurface)
       ALLOCATE(inject(i)%edges(1:self%nEdges))
-      ALLOCATE(inject(i)%weight(1:self%nEdges))
+      ! ALLOCATE(inject(i)%weight(1:self%nEdges))
       et = 0
       DO e=1, mesh%numEdges
         IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
           et = et + 1
           self%edges(et)  = mesh%edges(e)%obj%n
-          SELECT TYPE(edge => mesh%edges(e)%obj)
+          ! SELECT TYPE(edge => mesh%edges(e)%obj)
-          CLASS IS (meshEdgeCyl)
+          ! CLASS IS (meshEdgeCyl)
-            self%weight(et) = (edge%r(1)+edge%r(2))/2.D0
+          !   self%weight(et) = (edge%r(1)+edge%r(2))/2.D0
-
+          !
-          END SELECT
+          ! END SELECT
 
         END IF
 
       END DO
-      self%sumWeight = SUM(self%weight)
+      ! self%sumWeight = SUM(self%weight)
 
     END SUBROUTINE
 
@@ -105,7 +105,6 @@ MODULE moduleInject
       IMPLICIT NONE
 
       INTEGER:: i
-      LOGICAL:: mask(1:nInject)
       
       !$OMP SINGLE
       nPartInj = 0
@@ -148,8 +147,8 @@ MODULE moduleInject
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
-      REAL(8):: randomX
+      INTEGER:: randomX
-      INTEGER:: i, j
+      INTEGER:: i!, j
       INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
       CLASS(meshEdge), POINTER:: randomEdge
@@ -180,14 +179,15 @@ MODULE moduleInject
       !$OMP DO
       DO n = nMin, nMax
         !Select edge randomly from which inject particle
-        randomX = RAND()*self%sumWeight
+        ! randomX = RAND()*self%sumWeight
-        DO j = 1, self%nEdges
+        ! DO j = 1, self%nEdges
-          IF (randomX < self%weight(j)) EXIT
+        !   IF (randomX < self%weight(j)) EXIT
-          randomX = randomX - self%weight(j)
+        !   randomX = randomX - self%weight(j)
+        !
+        ! END DO
+        randomX = INT(DBLE(self%nEdges-1)*RAND()) + 1 
 
-        END DO
+        randomEdge => mesh%edges(self%edges(randomX))%obj
-
-        randomEdge => mesh%edges(self%edges(j))%obj
         !Random position in edge
         partInj(n)%r = randomEdge%randPos()
         !Volume associated to the edge:

src/modules/moduleInput.f95

@@ -31,9 +31,6 @@ MODULE moduleInput
       !Read species
       CALL readSpecies(config)
 
-      !Reads case parameters
-      CALL readCase(config)
-
       !Read interactions between species
       CALL readInteractions(config)
 
@@ -44,6 +41,10 @@ MODULE moduleInput
       CALL verboseError('Reading Geometry...')
       CALL readGeometry(config)
 
+      !Reads case parameters
+      CALL verboseError('Reading Case Parameters...')
+      CALL readCase(config)
+
       !Read boundary for EM field
       CALL readEMBoundary(config)
 
@@ -267,6 +268,9 @@ MODULE moduleInput
       !TODO: In a future, this should include the particles from init states
       nPartOld = 0
 
+      !Initialize the lock for the non-analogue (NA) list of particles
+      CALL OMP_INIT_LOCK(lockNAScheme)
+
     END SUBROUTINE readSpecies
 
     !Reads information about interactions between species

src/modules/moduleList.f95

@@ -1,7 +1,6 @@
 !Linked list of particles
 MODULE moduleList
   USE moduleSpecies
-  USE OMP_LIB
   IMPLICIT NONE
 
   TYPE lNode 
@@ -21,9 +20,7 @@ MODULE moduleList
 
   END TYPE listNode
 
-  INTEGER(KIND=OMP_LOCK_KIND):: resetLock
+  TYPE(listNode):: partNAScheme !Particles comming from the nonAnalogue scheme
-
-  TYPE(listNode):: partNew
 
   TYPE pointerArray
     TYPE(particle), POINTER:: part

src/modules/moduleMeshCyl.f95

@@ -1034,13 +1034,8 @@ MODULE moduleMeshCyl
       CLASS(*), POINTER:: nextElement
 
       xi = self%phy2log(part%r)
-      !Checks if particle is inside of current cell
+      !Checks if particle is inside 'self' cell
       IF (self%inside(xi)) THEN
-        ! IF (PRESENT(oldCell)) THEN
-        !   !If oldCell, recalculate particle weight, as particle has entered a new cell.
-        !   part%weight = part%weight*oldCell%volume/self%volume
-        !
-        ! END IF
         part%vol  = self%n
         part%xi   = xi
         part%n_in = .TRUE.

src/modules/moduleSolver.f95

@@ -40,10 +40,14 @@ MODULE moduleSolver
     END SUBROUTINE solveEM_interface
 
     !Apply nonAnalogue scheme to a particle
-    SUBROUTINE nonAnalogue_interface(part)
+    SUBROUTINE nonAnalogue_interface(part, volOld, volNew)
       USE moduleSpecies
+      USE moduleMesh
+      IMPLICIT NONE
 
       TYPE(particle), INTENT(inout):: part
+      CLASS(meshVol), POINTER, INTENT(in):: volOld
+      CLASS(meshVol), POINTER, INTENT(inout):: volNew
 
     END SUBROUTINE nonAnalogue_interface
 
@@ -97,12 +101,16 @@ MODULE moduleSolver
 
     !Initialize the non-analogue scheme
     SUBROUTINE initNA(self, NAType)
+      USE moduleMesh
       IMPLICIT NONE
 
       CLASS(solverGeneric), INTENT(inout):: self
       CHARACTER(:), ALLOCATABLE:: NAType
 
       SELECT CASE(NAType)
+      CASE ('Volume')
+        self%nonAnalogue => volumeNAScheme
+
       CASE DEFAULT
         self%nonAnalogue => noNAScheme
 
@@ -170,8 +178,6 @@ MODULE moduleSolver
       END IF
       part_temp%v(2) =  cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
       part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
-      part_temp%sp  = part%sp
-      part_temp%vol = part%vol
       part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
       !Copy temporal particle to particle
       part=part_temp
@@ -238,12 +244,14 @@ MODULE moduleSolver
     
     SUBROUTINE doReset()
       USE moduleSpecies
+      USE moduleList
       IMPLICIT NONE
 
       INTEGER:: nn, n
       INTEGER, SAVE:: nPartNew
-      INTEGER, SAVE:: nInjIn, nOldIn
+      INTEGER, SAVE:: nInjIn, nOldIn, nNAScheme
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
+      TYPE(lNode), POINTER:: partCurr, partNext
 
       !$OMP SECTIONS
       !$OMP SECTION
@@ -258,13 +266,15 @@ MODULE moduleSolver
         nOldIn = COUNT(partOld%n_in)
 
       END IF
+      !$OMP SECTION
+      nNAScheme = partNAScheme%amount
       !$OMP END SECTIONS
 
       !$OMP BARRIER
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      nPartNew = nInjIn + nOldIn
+      nPartNew = nInjIn + nOldIn + nNAScheme
       ALLOCATE(partOld(1:nPartNew))
       !$OMP END SINGLE
 
@@ -290,6 +300,20 @@ MODULE moduleSolver
         END IF
 
       END DO
+      !$OMP SECTION
+      nn = nInjIn + nOldIn
+      partCurr => partNAScheme%head
+      DO n = 1, nNAScheme
+        partNext => partCurr%next
+        partOld(nn+n) = partCurr%part
+        DEALLOCATE(partCurr)
+        partCurr => partNext
+
+      END DO
+      IF (ASSOCIATED(partNAScheme%head)) NULLIFY(partNAScheme%head)
+      IF (ASSOCIATED(partNAScheme%tail)) NULLIFY(partNAScheme%tail)
+      partNAScheme%amount = 0
+
       !$OMP END SECTIONS
 
       !$OMP SINGLE
@@ -373,12 +397,88 @@ MODULE moduleSolver
 
     END SUBROUTINE noEMField
 
-    !Empty procedure that does no computation of EM field for neutral cases
+    !Modify particle weight as a function of cell volume and splits particle
-    SUBROUTINE noNAScheme(part)
+    SUBROUTINE volumeNAScheme(part, volOld, volNew)
       USE moduleSpecies
+      USE moduleMesh
       IMPLICIT NONE
 
       TYPE(particle), INTENT(inout):: part
+      CLASS(meshVol), POINTER, INTENT(in):: volOld
+      CLASS(meshVol), POINTER, INTENT(inout):: volNew
+      REAL(8):: fractionVolume, fractionWeight
+      INTEGER:: nSplit
+
+      !If particle has change cell, call nonAnalogue scheme
+      IF (volOld%n /= volNew%n) THEN
+        fractionVolume = volOld%volume/volNew%volume
+
+        part%weight = part%weight * fractionVolume
+
+        fractionWeight = part%weight / species(part%sp)%obj%weight
+
+        IF (fractionWeight >= 2.D0) THEN
+          nSplit = FLOOR(fractionWeight)
+          CALL splitParticle(part, nSplit, volNew)
+
+        ELSEIF (part%weight < 1.D0) THEN
+          !Particle has lost statistical meaning and will be terminated
+          part%n_in = .FALSE.
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE volumeNAScheme
+
+    !Subroutine to split the particle 'part' into a number 'nSplit' of particles.
+    !'nSplit-1' particles are added to the partNAScheme list
+    SUBROUTINE splitParticle(part, nSplit, vol)
+      USE moduleSpecies
+      USE moduleList
+      USE moduleMesh
+      USE OMP_LIB
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      INTEGER, INTENT(in):: nSplit
+      CLASS(meshVol), INTENT(inout):: vol
+      REAL(8):: newWeight
+      TYPE(particle), POINTER:: newPart
+      INTEGER:: p
+
+      newWeight = part%weight / nSplit
+
+      !Assign new weight to original particle
+      part%weight = newWeight
+      
+      !Add new particles to list of NA particles
+      DO p = 2, nSplit
+        !Allocate the pointer for the new particles
+        ALLOCATE(newPart)
+        !Copy data from original particle
+        newPart = part
+        CALL OMP_SET_LOCK(lockNAScheme)
+        CALL partNAScheme%add(newPart)
+        CALL OMP_UNSET_LOCK(lockNASCheme)
+        !Add particle to cell list
+        CALL OMP_SET_lock(vol%lock)
+        CALL vol%listPart_in%add(newPart)
+        CALL OMP_UNSET_lock(vol%lock)
+
+      END DO
+
+    END SUBROUTINE splitParticle
+
+    !Empty procedure that does no computation of EM field for neutral cases
+    SUBROUTINE noNAScheme(part, volOld, volNew)
+      USE moduleSpecies
+      USE moduleMesh
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      CLASS(meshVol), POINTER, INTENT(in):: volOld
+      CLASS(meshVol), POINTER, INTENT(inout):: volNew
 
     END SUBROUTINE noNAScheme
 
@@ -389,17 +489,13 @@ MODULE moduleSolver
 
       CLASS(solverGeneric), INTENT(in):: self
       TYPE(particle), INTENT(inout):: part
-      CLASS(meshVol), POINTER:: vol
+      CLASS(meshVol), POINTER:: volOld, volNew
-      INTEGER:: volOld !Old cell for particle
 
-      volOld = part%vol
+      volOld => mesh%vols(part%vol)%obj
-      vol => mesh%vols(volOld)%obj
+      CALL volOld%findCell(part)
-      CALL vol%findCell(part)
+      volNew => mesh%vols(part%vol)%obj
-      !If particle has change cell, call nonAnalogue scheme
+      !Call the NA shcme
-      IF (volOld /= part%vol) THEN
+      CALL self%nonAnalogue(part, volOld, volNew)
-        CALL self%nonAnalogue(part)
-
-      END IF
 
     END SUBROUTINE updateParticleCell
 

src/modules/moduleSpecies.f95

@@ -1,6 +1,7 @@
 !Contains the information about species (particles)
 MODULE moduleSpecies
   USE moduleCaseParam
+  USE OMP_LIB
   IMPLICIT NONE
   
   TYPE, ABSTRACT:: speciesGeneric
@@ -44,9 +45,11 @@ MODULE moduleSpecies
   !Arrays that contain the particles
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
+  INTEGER(KIND=OMP_LOCK_KIND):: lockNAScheme !Lock for the NA list of particles
 
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(sp)
+      USE moduleErrors
       IMPLICIT NONE
 
       CHARACTER(:), ALLOCATABLE:: speciesName
@@ -60,7 +63,8 @@ MODULE moduleSpecies
           EXIT
         END IF
       END DO
-      !TODO: add error handling when species not found
+      !If no species is found, call a critical error
+      IF (sp == 0) CALL criticalError('Species ' // speciesName // ' not found.', 'speciesName2Index')
 
     END FUNCTION speciesName2Index