First implementation of multiple pushers for different species

runs/cylFlow/input.json

@@ -3,12 +3,12 @@
     "path": "./runs/cylFlow/",
 	"trigger": 10,
 	"cpuTime": true,
-	"numColl": true
+	"numColl": false
   },
   "geometry": {
     "type":     "2DCyl",
     "meshType": "gmsh",
-    "meshFile": "meshTria.msh"
+    "meshFile": "mesh.msh"
   },
   "species": [
     {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 5.0e7}
@@ -30,9 +30,9 @@
     "radius": 1.88e-10
   },
   "case": {
-    "tau":  6.0e-7,
+    "tau":  [6.0e-7],
     "time":  3.0e-3,
-    "solver": "2DCylNeutral"
+    "pusher": ["2DCylNeutral"]
   },
   "interactions": {
     "folderCollisions": "./data/collisions/",

src/modules/moduleEM.f95

@@ -58,7 +58,7 @@ MODULE moduleEM
 
       xi  = part%xi
 
-      elField = mesh%vols(part%e_p)%obj%gatherEF(xi)
+      elField = mesh%vols(part%vol)%obj%gatherEF(xi)
 
     END FUNCTION gatherElecField
 

src/modules/moduleInject.f95

@@ -10,7 +10,7 @@ MODULE moduleInject
     REAL(8):: vTh(1:3) !Thermal Velocity
     REAL(8):: n(1:3) !Direction of injection
     INTEGER:: nParticles !Number of particles to introduce each time step
-    INTEGER:: pt !Species of injection
+    INTEGER:: sp !Species of injection
     INTEGER:: nEdges
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
     REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection
@@ -26,7 +26,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, T, flow, units, pt, physicalSurface)
+    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
       USE moduleMesh
       USE moduleMeshCyl
       USE moduleRefParam
@@ -38,7 +38,7 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
       REAL(8), INTENT(in):: v, n(1:3), T(1:3)
-      INTEGER, INTENT(in):: pt, physicalSurface
+      INTEGER, INTENT(in):: sp, physicalSurface
       REAL(8), INTENT(in):: flow
       CHARACTER(:), ALLOCATABLE, INTENT(in):: units
       INTEGER:: e, et
@@ -51,21 +51,21 @@ MODULE moduleInject
       SELECT CASE(units)
       CASE ("sccm")
         !Standard cubic centimeter per minute
-        self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
+        self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(sp)%obj%weight)
 
       CASE ("A")
         !Input current in Ampers
-        self%nParticles = INT(flow*tau*ti_ref/(qe*species(pt)%obj%weight))
+        self%nParticles = INT(flow*tau*ti_ref/(qe*species(sp)%obj%weight))
 
       CASE DEFAULT
         CALL criticalError("No support for units: " // units, 'initInject')
         
       END SELECT
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
-      self%pt  = pt
+      self%sp  = sp
 
       self%v   = self%v * self%n
-      self%vTh = DSQRT(self%T/species(self%pt)%obj%m)
+      self%vTh = DSQRT(self%T/species(self%sp)%obj%m)
 
       !Gets the edge elements from which particles are injected
       !TODO: Improve this A LOT
@@ -145,13 +145,13 @@ MODULE moduleInject
       nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
       nMax = nMin + self%nParticles - 1
       !Assign particle type
-      partInj(nMin:nMax)%pt = self%pt
+      partInj(nMin:nMax)%sp = self%sp
       !Assign weight to particle.
-      partInj(nMin:nMax)%weight = species(self%pt)%obj%weight
+      partInj(nMin:nMax)%weight = species(self%sp)%obj%weight
       !Particle is considered to be outside the domain
       partInj(nMin:nMax)%n_in = .FALSE.
       !Assign charge/mass to charged particle.
-      SELECT TYPE(sp => species(self%pt)%obj)
+      SELECT TYPE(sp => species(self%sp)%obj)
       TYPE IS(speciesCharged)
         partInj(nMin:nMax)%qm = sp%qm
 
@@ -173,10 +173,10 @@ MODULE moduleInject
         partInj(n)%r = randomEdge%randPos()
         !Volume associated to the edge:
         IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN
-          partInj(n)%e_p = randomEdge%e1%n
+          partInj(n)%vol = randomEdge%e1%n
         
         ELSE
-          partInj(n)%e_p = randomEdge%e2%n
+          partInj(n)%vol = randomEdge%e2%n
 
         END IF
 
@@ -185,9 +185,9 @@ MODULE moduleInject
                           vBC(self%v(3), self%vTh(3)) /)
 
         !Push new particle
-        CALL solver%push(partInj(n))
+        CALL solver%pusher(self%sp)%pushParticle(partInj(n))
         !Assign cell to new particle
-        CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n))
+        CALL solver%updateParticleCell(partInj(n))
 
       END DO
       !$OMP END DO

src/modules/moduleInput.f95

@@ -25,15 +25,15 @@ MODULE moduleInput
       !Reads reference parameters
       CALL readReference(config)
 
-      !Reads case parameters
-      CALL readCase(config)
-
       !Reads output parameters
       CALL readOutput(config)
 
       !Read species
       CALL readSpecies(config)
 
+      !Reads case parameters
+      CALL readCase(config)
+
       !Read interactions between species
       CALL readInteractions(config)
 
@@ -108,44 +108,69 @@ MODULE moduleInput
       USE moduleErrors
       USE moduleCaseParam
       USE moduleSolver
+      USE moduleSpecies
       USE json_module
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
       LOGICAL:: found
+      REAL(8), ALLOCATABLE:: tauSpecies(:)
       CHARACTER(:), ALLOCATABLE:: object
       REAL(8):: time !simulation time in [t]
-      CHARACTER(:), ALLOCATABLE:: solverType
+      CHARACTER(:), ALLOCATABLE:: pusherType, EMType, NAType
+      INTEGER:: nTau, nSolver
+      INTEGER:: i
+      CHARACTER(2):: iString
 
       object = 'case'
-      !Time parameters
-      CALL config%get(object // '.tau', tau, found)
-      IF (.NOT. found) CALL criticalError('Required parameter tau not found','readCase')
 
+      !Time parameters
+      CALL config%info(object // '.tau', found, n_children = nTau)
+      IF (.NOT. found .OR. nTau == 0) CALL criticalError('Required parameter tau not found','readCase')
+      ALLOCATE(tauSpecies(1:nSpecies))
+      DO i = 1, nTau
+        WRITE(iString, '(I2)') i
+        CALL config%get(object // '.tau(' // TRIM(iString) // ')', tauSpecies(i), found)
+
+      END DO
+      IF (nTau < nSpecies) THEN
+        CALL warningError('Using minimum time step for some species')
+        tauSpecies(nTau+1:nSpecies) = MINVAL(tauSpecies(1:nTau))
+
+      END IF
+      !Selects the minimum tau as reference value
+      tau = MINVAL(tauSpecies)
+
+      !Gets the simulation time
       CALL config%get(object // '.time', time, found)
       IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
-
-      CALL config%get(object // '.solver', solverType, found)
-      IF (.NOT. found) CALL criticalError('Required parameter solver not found','readCase')
-
-      !Allocate the type of solver
-      SELECT CASE(solverType)
-      CASE('2DCylNeutral')
-        !Allocate the solver for neutral pusher 2D Cylindrical
-        ALLOCATE(solver, source=solverCylNeutral())
-
-      CASE('2DCylCharged')
-        !Allocate the solver for charged pusher 2D Cylindrical
-        ALLOCATE(solver, source=solverCylCharged())
-
-      CASE DEFAULT
-        CALL criticalError('Solver ' // solverType // ' not found','readCase')
-
-      END SELECT
-
       !Convert simulation time to number of iterations
       tmax = INT(time/tau)
 
+      !Gest the pusher for each species
+      CALL config%info(object // '.pusher', found, n_children = nSolver)
+      IF (.NOT. found .OR. nSolver /= nSpecies) CALL criticalError('Required parameter pusher not found','readCase')
+      !Allocates all the pushers for particles
+      ALLOCATE(solver%pusher(1:nSpecies))
+      !Initialize pushers
+      DO i = 1, nSolver
+        WRITE(iString, '(I2)') i
+        CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)
+
+        !Associate the type of solver
+        CALL solver%pusher(i)%init(pusherType, tau, tauSpecies(i))
+
+      END DO
+      
+      !Gets the solver for the electromagnetic field
+      CALL config%get(object // '.EMSolver', EMType, found)
+      CALL solver%initEM(EMType)
+
+      !Gest the non-analogue scheme
+      CALL config%get(object // '.NAScheme', NAType, found)
+      CALL solver%initNA(NAType)
+
+
       !Makes tau non-dimensional
       tau = tau / ti_ref
 
@@ -234,7 +259,7 @@ MODULE moduleInput
         !Assign shared parameters for all species
         CALL config%get(object // '.name',   species(i)%obj%name,   found)
         CALL config%get(object // '.weight', species(i)%obj%weight, found)
-        species(i)%obj%pt = i
+        species(i)%obj%sp = i
         species(i)%obj%m  = mass
 
       END DO
@@ -462,7 +487,7 @@ MODULE moduleInput
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
-      INTEGER:: pt
+      INTEGER:: sp
 
       CALL config%info('inject', found, n_children = nInject)
       ALLOCATE(inject(1:nInject))
@@ -473,7 +498,7 @@ MODULE moduleInput
 
         !Find species
         CALL config%get(object // '.species', speciesName, found)
-        pt = speciesName2Index(speciesName)
+        sp = speciesName2Index(speciesName)
 
         CALL config%get(object // '.name',            name,            found)
         CALL config%get(object // '.v',               v,               found)
@@ -483,7 +508,7 @@ MODULE moduleInput
         CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, units, pt, physicalSurface)
+        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
 
       END DO
 

src/modules/moduleMeshCyl.f95

@@ -463,22 +463,22 @@ MODULE moduleMeshCyl
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%pt)
+      vertex => self%n1%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%pt)
+      vertex => self%n2%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%pt)
+      vertex => self%n3%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
 
-      vertex => self%n4%output(part%pt)
+      vertex => self%n4%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(4)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
@@ -852,17 +852,17 @@ MODULE moduleMeshCyl
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%pt)
+      vertex => self%n1%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%pt)
+      vertex => self%n2%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%pt)
+      vertex => self%n3%output(part%sp)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
@@ -1041,7 +1041,7 @@ MODULE moduleMeshCyl
         !   part%weight = part%weight*oldCell%volume/self%volume
         !
         ! END IF
-        part%e_p  = self%n
+        part%vol  = self%n
         part%xi   = xi
         part%n_in = .TRUE.
         !Assign particle to listPart_in
@@ -1117,7 +1117,7 @@ MODULE moduleMeshCyl
           part_i => partTemp(rnd)%part
           rnd = 1 + FLOOR(nPart*RAND())
           part_j => partTemp(rnd)%part
-          ij = interactionIndex(part_i%pt, part_j%pt)
+          ij = interactionIndex(part_i%sp, part_j%sp)
           sigmaVrelMaxNew = 0.D0
           DO k = 1, interactionMatrix(ij)%amount
             CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)

src/modules/moduleSolver.f95

@@ -1,13 +1,31 @@
 MODULE moduleSolver
 
-  TYPE, PUBLIC, ABSTRACT:: solverGeneric
+  !Generic type for pusher of particles
+  TYPE, PUBLIC:: pusherGeneric
+    PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL()
+    !Boolean to indicate if the species is moved in the iteration
+    LOGICAL:: pushSpecies
+    !How many interations between advancing the species
+    INTEGER:: every
     CONTAINS
-      PROCEDURE(push_interafece), DEFERRED, NOPASS:: push
+      PROCEDURE, PASS:: init => initPusher
-      PROCEDURE(solveEM_interface), DEFERRED, NOPASS:: solveEMField
+
+  END TYPE pusherGeneric
+
+  !Generic type for solver
+  TYPE, PUBLIC:: solverGeneric
+    TYPE(pusherGeneric), ALLOCATABLE:: pusher(:)
+    PROCEDURE(solveEM_interface), POINTER, NOPASS:: solveEM => NULL()
+    PROCEDURE(nonAnalogue_interface), POINTER, NOPASS:: nonAnalogue => NULL()
+    CONTAINS
+      PROCEDURE, PASS:: initEM
+      PROCEDURE, PASS:: initNA
+      PROCEDURE, PASS:: updateParticleCell
+      PROCEDURE, PASS:: updatePushSpecies
 
   END TYPE solverGeneric
 
-  ABSTRACT INTERFACE
+  INTERFACE
     !Push a particle
     PURE SUBROUTINE push_interafece(part)
       USE moduleSpecies
@@ -21,38 +39,94 @@ MODULE moduleSolver
     
     END SUBROUTINE solveEM_interface
 
+    !Apply nonAnalogue scheme to a particle
+    SUBROUTINE nonAnalogue_interface(part)
+      USE moduleSpecies
+
+      TYPE(particle), INTENT(inout):: part
+
+    END SUBROUTINE nonAnalogue_interface
+
   END INTERFACE
 
-  TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylNeutral
+  TYPE(solverGeneric):: solver
-    CONTAINS
-      PROCEDURE, NOPASS:: push => pushCylNeutral
-      PROCEDURE, NOPASS:: solveEMField => noEMField
-
-  END TYPE solverCylNeutral
-
-  TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylCharged
-    CONTAINS
-      PROCEDURE, NOPASS:: push => pushCylCharged
-      PROCEDURE, NOPASS:: solveEMField => solveElecField
-
-  END TYPE solverCylCharged
-
-  CLASS(solverGeneric), ALLOCATABLE:: solver
 
   CONTAINS
+    !Init Pusher
+    SUBROUTINE initPusher(self, pusherType, tau, tauSp)
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(pusherGeneric), INTENT(out):: self
+      CHARACTER(:), ALLOCATABLE:: pusherType
+      REAL(8):: tau, tauSp
+
+      SELECT CASE(pusherType)
+      CASE('2DCylNeutral')
+        self%pushParticle => pushCylNeutral
+
+      CASE('2DCylCharged')
+        self%pushParticle => pushCylCharged
+
+      CASE DEFAULT
+        CALL criticalError('Solver ' // pusherType // ' not found','readCase')
+
+      END SELECT
+
+      self%pushSpecies = .FALSE.
+      self%every = INT(tau/tauSp)
+
+    END SUBROUTINE initPusher
+
+    SUBROUTINE initEM(self, EMType)
+      IMPLICIT NONE
+
+      CLASS(solverGeneric), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE:: EMType
+
+      SELECT CASE(EMType)
+      CASE('Electrostatic')
+        self%solveEM => solveElecField
+
+      CASE DEFAULT
+        self%solveEM => noEMField
+
+      END SELECT
+
+    END SUBROUTINE initEM
+
+    !Initialize the non-analogue scheme
+    SUBROUTINE initNA(self, NAType)
+      IMPLICIT NONE
+
+      CLASS(solverGeneric), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE:: NAType
+
+      SELECT CASE(NAType)
+      CASE DEFAULT
+        self%nonAnalogue => noNAScheme
+
+      END SELECT
+
+    END SUBROUTINE initNA
+    
+    !Do all pushes for particles
     SUBROUTINE doPushes()
       USE moduleSpecies
       USE moduleMesh
       IMPLICIT NONE
 
       INTEGER:: n
+      INTEGER:: sp
 
       !$OMP DO
       DO n=1, nPartOld
+        !Select species type
+        sp = partOld(n)%sp
         !Push particle
-        CALL solver%push(partOld(n))
+        CALL solver%pusher(sp)%pushParticle(partOld(n))
         !Find cell in wich particle reside
-        CALL mesh%vols(partOld(n)%e_p)%obj%findCell(partOld(n))
+        CALL solver%updateParticleCell(partOld(n))
 
       END DO
       !$OMP END DO
@@ -89,8 +163,8 @@ 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%pt  = part%pt
+      part_temp%sp  = part%sp
-      part_temp%e_p = part%e_p
+      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
@@ -224,7 +298,7 @@ MODULE moduleSolver
       !Loops over the particles to scatter them
       !$OMP DO
       DO n=1, nPartOld
-        CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n))
+        CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n))
 
       END DO
       !$OMP END DO
@@ -234,7 +308,7 @@ MODULE moduleSolver
     SUBROUTINE doEMField()
       IMPLICIT NONE
 
-      CALL solver%solveEMField()
+      CALL solver%solveEM()
 
     END SUBROUTINE doEMField
 
@@ -283,12 +357,59 @@ MODULE moduleSolver
 
     END SUBROUTINE solveElecField
 
-    !Empty module that does no computation of EM field for neutral cases
+    !Empty procedure that does no computation of EM field for neutral cases
     SUBROUTINE noEMField()
       IMPLICIT NONE
 
     END SUBROUTINE noEMField
 
+    !Empty procedure that does no computation of EM field for neutral cases
+    SUBROUTINE noNAScheme(part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+
+    END SUBROUTINE noNAScheme
+
+    SUBROUTINE updateParticleCell(self, part)
+      USE moduleSpecies
+      USE moduleMesh
+      IMPLICIT NONE
+
+      CLASS(solverGeneric), INTENT(in):: self
+      TYPE(particle), INTENT(inout):: part
+      CLASS(meshVol), POINTER:: vol
+      INTEGER:: volOld !Old cell for particle
+
+      volOld = part%vol
+      vol => mesh%vols(volOld)%obj
+      CALL vol%findCell(part)
+      !If particle has change cell, call nonAnalogue scheme
+      IF (volOld /= part%vol) THEN
+        CALL self%nonAnalogue(part)
+
+      END IF
+
+    END SUBROUTINE updateParticleCell
+
+    !Update the information about if a species needs to be moved this iteration
+    SUBROUTINE updatePushSpecies(self, t)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(solverGeneric), INTENT(inout):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: s
+
+      DO s=1, nSpecies
+        self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0
+
+      END DO
+
+    END SUBROUTINE updatePushSpecies
+
+    !Output the different data and information
     SUBROUTINE doOutput(t)
       USE moduleMesh
       USE moduleOutput

src/modules/moduleSpecies.f95

@@ -6,7 +6,7 @@ MODULE moduleSpecies
   TYPE, ABSTRACT:: speciesGeneric
     CHARACTER(:), ALLOCATABLE:: name
     REAL(8):: m=0.D0, weight=0.D0
-    INTEGER:: pt=0
+    INTEGER:: sp=0
   END TYPE speciesGeneric
 
   TYPE, EXTENDS(speciesGeneric):: speciesNeutral
@@ -29,8 +29,8 @@ MODULE moduleSpecies
   TYPE particle
     REAL(8):: r(1:3) !Position
     REAL(8):: v(1:3) !Velocity
-    INTEGER:: pt !Particle species id
+    INTEGER:: sp !Particle species id
-    INTEGER:: e_p !Index of element in which the particle is located
+    INTEGER:: vol !Index of element in which the particle is located
     REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p.
     LOGICAL:: n_in !Flag that indicates if a particle is in the domain
     REAL(8):: weight=0.D0 !weight of particle
@@ -46,17 +46,17 @@ MODULE moduleSpecies
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
 
   CONTAINS
-    FUNCTION speciesName2Index(speciesName) RESULT(pt)
+    FUNCTION speciesName2Index(speciesName) RESULT(sp)
       IMPLICIT NONE
 
       CHARACTER(:), ALLOCATABLE:: speciesName
-      INTEGER:: pt 
+      INTEGER:: sp 
       INTEGER:: n
 
-      pt = 0
+      sp = 0
       DO n = 1, nSpecies
         IF (speciesName == species(n)%obj%name) THEN
-          pt = species(n)%obj%pt
+          sp = species(n)%obj%sp
           EXIT
         END IF
       END DO