Implementation of different time steps per species.

runs/ALPHIE_Grid/input.json

@@ -39,9 +39,10 @@
     "temperature": 2500.0
   },
   "case": {
-    "tau":  1.0e-11,
+    "tau":  [1.0e-11, 1.0e-11],
     "time": 2.0e-6,
-    "solver": "2DCylCharged"
+    "pusher": ["2DCylCharged", "2DCylCharged"],
+    "EMSolver": "Electrostatic"
   },
   "parallel": {
     "OpenMP":{

src/fpakc.f95

@@ -40,9 +40,13 @@ PROGRAM fpakc
     !Insert new particles and push them
     !$OMP SINGLE
     tStep = omp_get_wtime()
+
+    !Checks if a species needs to me moved in this iteration
+    CALL solver%updatePushSpecies(t)
     tPush = omp_get_wtime()
     !$OMP END SINGLE
 
+    !Injects new particles
     CALL doInjects()
 
     !Push old particles
@@ -50,6 +54,7 @@ PROGRAM fpakc
 
     !$OMP SINGLE
     tPush = omp_get_wtime() - tPush
+
     !Collisions
     tColl = omp_get_wtime()
     !$OMP END SINGLE
@@ -58,6 +63,7 @@ PROGRAM fpakc
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl
+
     !Reset particles
     tReset = omp_get_wtime()
     !$OMP END SINGLE
@@ -75,6 +81,7 @@ PROGRAM fpakc
 
     !$OMP SINGLE
     tWeight = omp_get_wtime() - tWeight
+
     tEMField = omp_get_wtime()
     !$OMP END SINGLE
 
@@ -82,7 +89,9 @@ PROGRAM fpakc
 
     !$OMP SINGLE
     tEMField = omp_get_wtime() - tEMField
+
     tStep = omp_get_wtime() - tStep
+
     !Output data
     CALL doOutput(t)
     !$OMP END SINGLE

src/modules/moduleCaseParam.f95

@@ -2,7 +2,8 @@
 MODULE moduleCaseParam
   !Maximum number of iterations and number of species
   INTEGER:: tmax
-  REAL(8):: tau
+  REAL(8), ALLOCATABLE:: tau(:)
+  REAL(8):: tauMin
 
 END MODULE moduleCaseParam
 

src/modules/moduleInject.f95

@@ -32,6 +32,7 @@ MODULE moduleInject
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
+      USE moduleSolver
       USE moduleErrors
       IMPLICIT NONE
 
@@ -51,17 +52,19 @@ MODULE moduleInject
       SELECT CASE(units)
       CASE ("sccm")
         !Standard cubic centimeter per minute
-        self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(sp)%obj%weight)
+        self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight)
 
       CASE ("A")
         !Input current in Ampers
-        self%nParticles = INT(flow*tau*ti_ref/(qe*species(sp)%obj%weight))
+        self%nParticles = INT(flow*tauMin*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%nParticles = self%nParticles * solver%pusher(sp)%every
       self%sp  = sp
 
       self%v   = self%v * self%n
@@ -92,19 +95,30 @@ MODULE moduleInject
 
       END DO
       self%sumWeight = SUM(self%weight)
-      nPartInj = nPartInj + self%nParticles
 
     END SUBROUTINE
 
     !Injection of particles
     SUBROUTINE doInjects()
-    
+      USE moduleSpecies
+      USE moduleSolver
       IMPLICIT NONE
 
       INTEGER:: i
+      LOGICAL:: mask(1:nInject)
+      
+      !$OMP SINGLE
+      nPartInj = 0
+      DO i = 1, nInject
+        IF (solver%pusher(inject(i)%sp)%pushSpecies) nPartInj = nPartInj + inject(i)%nParticles
+
+      END DO
+      IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
+      ALLOCATE(partInj(1:nPartInj))
+      !$OMP END SINGLE
 
       DO i=1, nInject
-        CALL inject(i)%addParticles()
+        IF (solver%pusher(inject(i)%sp)%pushSpecies) CALL inject(i)%addParticles()
       END DO
 
     END SUBROUTINE doInjects
@@ -135,14 +149,19 @@ MODULE moduleInject
 
       CLASS(injectGeneric), INTENT(in):: self
       REAL(8):: randomX
-      INTEGER:: j
+      INTEGER:: i, j
       INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
       CLASS(meshEdge), POINTER:: randomEdge
 
       !Insert particles
       !$OMP SINGLE
-      nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
+      nMin = 0
+      DO i = 1, self%id - 1
+        IF (solver%pusher(inject(i)%sp)%pushSpecies) nMin = nMin + inject(i)%nParticles
+
+      END DO
+      nMin = nMin + 1
       nMax = nMin + self%nParticles - 1
       !Assign particle type
       partInj(nMin:nMax)%sp = self%sp
@@ -184,8 +203,11 @@ MODULE moduleInject
                           vBC(self%v(2), self%vTh(2)), &
                           vBC(self%v(3), self%vTh(3)) /)
 
+        IF (solver%pusher(self%sp)%pushSpecies) THEN
           !Push new particle
           CALL solver%pusher(self%sp)%pushParticle(partInj(n))
+
+        END IF
         !Assign cell to new particle
         CALL solver%updateParticleCell(partInj(n))
 

src/modules/moduleInput.f95

@@ -114,7 +114,6 @@ MODULE moduleInput
 
       TYPE(json_file), INTENT(inout):: config
       LOGICAL:: found
-      REAL(8), ALLOCATABLE:: tauSpecies(:)
       CHARACTER(:), ALLOCATABLE:: object
       REAL(8):: time !simulation time in [t]
       CHARACTER(:), ALLOCATABLE:: pusherType, EMType, NAType
@@ -127,25 +126,24 @@ MODULE moduleInput
       !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))
+      ALLOCATE(tau(1:nSpecies))
       DO i = 1, nTau
         WRITE(iString, '(I2)') i
-        CALL config%get(object // '.tau(' // TRIM(iString) // ')', tauSpecies(i), found)
+        CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(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))
+        tau(nTau+1:nSpecies) = MINVAL(tau(1:nTau))
 
       END IF
-      !Selects the minimum tau as reference value
+      tauMin = MINVAL(tau)
-      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')
       !Convert simulation time to number of iterations
-      tmax = INT(time/tau)
+      tmax = INT(time/tauMin)
 
       !Gest the pusher for each species
       CALL config%info(object // '.pusher', found, n_children = nSolver)
@@ -158,7 +156,7 @@ MODULE moduleInput
         CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)
 
         !Associate the type of solver
-        CALL solver%pusher(i)%init(pusherType, tau, tauSpecies(i))
+        CALL solver%pusher(i)%init(pusherType, tauMin, tau(i))
 
       END DO
       
@@ -173,6 +171,7 @@ MODULE moduleInput
 
       !Makes tau non-dimensional
       tau = tau / ti_ref
+      tauMin = tauMin / ti_ref
 
     END SUBROUTINE readCase
 

src/modules/moduleMeshCyl.f95

@@ -1102,7 +1102,7 @@ MODULE moduleMeshCyl
       self%nColl = 0
       nPart = self%listPart_in%amount
       IF (nPart > 1) THEN
-        pMax = self%totalWeight*self%sigmaVrelMax*tau/self%volume
+        pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
         self%nColl = INT(REAL(nPart)*pMax*0.5D0)
 
         !Converts the list of particles to an array for easy access

src/modules/moduleMeshCylRead.f95

@@ -602,7 +602,7 @@ MODULE moduleMeshCylRead
       CHARACTER(:), ALLOCATABLE:: fileName
       CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
 
-      time = DBLE(t)*tau*ti_ref
+      time = DBLE(t)*tauMin*ti_ref
 
       DO i = 1, nSpecies
         WRITE(tstring, '(I6.6)') t
@@ -686,7 +686,7 @@ MODULE moduleMeshCylRead
 
 
       IF (collOutput) THEN
-        time = DBLE(t)*tau*ti_ref
+        time = DBLE(t)*tauMin*ti_ref
         WRITE(tstring, '(I6.6)') t
 
         fileName='OUTPUT_' // tstring// '_Collisions.msh'
@@ -730,7 +730,7 @@ MODULE moduleMeshCylRead
       REAL(8):: xi(1:3)
 
       IF (emOutput) THEN
-        time = DBLE(t)*tau*ti_ref
+        time = DBLE(t)*tauMin*ti_ref
         WRITE(tstring, '(I6.6)') t
 
         fileName='OUTPUT_' // tstring// '_EMField.msh'

src/modules/moduleSolver.f95

@@ -74,7 +74,7 @@ MODULE moduleSolver
       END SELECT
 
       self%pushSpecies = .FALSE.
-      self%every = INT(tau/tauSp)
+      self%every = INT(tauSp/tau)
 
     END SUBROUTINE initPusher
 
@@ -123,11 +123,15 @@ MODULE moduleSolver
       DO n=1, nPartOld
         !Select species type
         sp = partOld(n)%sp
+        !Checks if the species sp is update this iteration
+        IF (solver%pusher(sp)%pushSpecies) THEN
           !Push particle
           CALL solver%pusher(sp)%pushParticle(partOld(n))
           !Find cell in wich particle reside
           CALL solver%updateParticleCell(partOld(n))
 
+        END IF
+
       END DO
       !$OMP END DO
 
@@ -139,19 +143,22 @@ MODULE moduleSolver
       IMPLICIT NONE
 
       TYPE(particle), INTENT(inout):: part
-      REAL(8):: v_p_oh_star(2:3)
       TYPE(particle):: part_temp
+      REAL(8):: tauSp
       REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
+      REAL(8):: v_p_oh_star(2:3)
 
       part_temp = part
+      !Time step for the species
+      tauSp = tau(part_temp%sp)
       !z
       part_temp%v(1) = part%v(1)
-      part_temp%r(1) = part%r(1) + part_temp%v(1)*tau
+      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp
       !r,theta
       v_p_oh_star(2) = part%v(2)
-      x_new = part%r(2) + v_p_oh_star(2)*tau
+      x_new = part%r(2) + v_p_oh_star(2)*tauSp
       v_p_oh_star(3) = part%v(3)
-      y_new =             v_p_oh_star(3)*tau
+      y_new =             v_p_oh_star(3)*tauSp
       r = DSQRT(x_new**2+y_new**2)
       part_temp%r(2) = r
       IF (r > 0.D0) THEN
@@ -181,19 +188,22 @@ MODULE moduleSolver
       REAL(8):: v_p_oh_star(2:3)
       TYPE(particle):: part_temp
       REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
-      REAL(8):: qmEFt(1:3)!charge*tau*EF/mass
+      REAL(8):: tauSp
+      REAL(8):: qmEFt(1:3)!charge*tauSp*EF/mass
 
       part_temp = part
+      !Time step for the species
+      tauSp = tau(part_temp%sp)
       !Get electric field at particle position
-      qmEFt = part_temp%qm*gatherElecField(part_temp)*tau
+      qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp
       !z
       part_temp%v(1) = part%v(1) + qmEFt(1)
-      part_temp%r(1) = part%r(1) + part_temp%v(1)*tau
+      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp
       !r,theta
       v_p_oh_star(2) = part%v(2) + qmEFt(2)
-      x_new = part%r(2) + v_p_oh_star(2)*tau
+      x_new = part%r(2) + v_p_oh_star(2)*tauSp
       v_p_oh_star(3) = part%v(3) + qmEFt(3)
-      y_new =             v_p_oh_star(3)*tau
+      y_new =             v_p_oh_star(3)*tauSp
       r = DSQRT(x_new**2+y_new**2)
       part_temp%r(2) = r
       IF (r > 0.D0) THEN