Avoid very low macroparticle weight

doc/user-manual/fpakc_UserManual.tex

@@ -401,6 +401,9 @@ make
           \item \textbf{bounraiesParticle}: Logical.
                                             Determines if the properties of the boundaries are written.
                                             \textit{Warning}: Not all particle models produce this output, only those whose properties change during the simulation.
+          \item \textbf{injects}: Logical.
+                                  Determines if the properties of the injection of particles are written.
+                                  \textit{Warning}: Not all injection types produce this output, only those whose properties change during the simulation.
           \item \textbf{probes}: Array of objects.
                                  Defines the probes employed for obtaining the distribution function at specific positions.
                                  See Sec.~\ref{sec:probing} for more information.
@@ -688,6 +691,14 @@ make
                             Units of $\unit{K}$
                             Array dimension $3$.
                             Temperature in each direction.
+          \item \textbf{type}: Character.
+                               Type of injection.
+                               Available values are:
+                               \begin{itemize}
+                               \item \textbf{constant}: The flow does not change.
+                                 This is the default value if none is provided.
+                               \item \textbf{quasiNeutral}: The weight of the injected species' particles is changed over time to keep a zero electric field at the boundary.
+                               \end{itemize}
           \item \textbf{physicalSurface}: Integer.
                                           Identification of the edge in the mesh file.
           \item \textbf{particlesPerEdge}: Integer.

src/fpakc.f90

@@ -1,15 +1,39 @@
 !  FPAKC main program
 PROGRAM fpakc
-  USE moduleCompTime
-  USE moduleCaseParam
-  USE moduleInput
-  USE moduleInject
-  USE moduleSolver
-  USE moduleMesh
-  USE moduleProbe
-  USE moduleErrors
-  USE OMP_LIB
-  IMPLICIT NONE
+  use moduleCompTime,  only: tStep,    &
+                             tEMField, &
+                             tCoul,    &
+                             tColl,    &
+                             tPush,    &
+                             tReset,   &
+                             tWeight
+  use omp_lib,         only: omp_get_wtime
+  use moduleErrors,    only: criticalError, &
+                             verboseError
+  use moduleInput,     only: readCOnfig, &
+                             initOutput
+  use moduleCaseParam, only: timeStep, &
+                             tInitial, &
+                             tFinal
+  use moduleProbe,     only: resetProbes
+  use moduleSolver,    only: doScatter, &
+                             doEMField, &
+                             doOutput,  &
+                             solver,    &
+                             doPushes,  &
+                             doReset,   &
+                             doAverage, &
+                             doInjects
+  use moduleMesh,      only: boundariesEM_update,       &
+                             boundariesEM_update,       &
+                             boundariesParticle_update, &
+                             mesh,                      &
+                             meshForMCC,                &
+                             doMCCollisions,            &
+                             doCollisions,              &
+                             doCoulombScattering
+  use moduleInject,    only: updateInjects
+  implicit none
 
   ! arg1 = Input argument 1 (input file)
   CHARACTER(200):: arg1
@@ -71,6 +95,9 @@ PROGRAM fpakc
     ! Update EM boundary models
     call boundariesEM_update()
 
+    ! Update injects
+    call updateInjects()
+
     !Checks if a species needs to me moved in this iteration
     CALL solver%updatePushSpecies()
 

src/modules/common/moduleRandom.f90

@@ -53,7 +53,7 @@ MODULE moduleRandom
       implicit none
 
       real(8):: rnd
-      real(8):: v1, v2, Rsquare
+      real(8):: v1, v2
 
       v1 = 0.d0
       do while (v1 <= 0.d0)
@@ -67,19 +67,20 @@ MODULE moduleRandom
     end function randomMaxwellian
 
     !Returns a random number in a Maxwellian distribution of mean 0 and width 1
-    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+    function randomHalfMaxwellian() result(rnd)
       IMPLICIT NONE
 
-      REAL(8):: rnd
-      REAL(8):: x
+      real(8):: rnd
+      real(8):: v1
 
       rnd = 0.D0
-      x = 0.D0
-      DO WHILE (x == 0.D0)
-        CALL RANDOM_NUMBER(x)
-      END DO
+      v1 = 0.D0
+      do while (v1 == 0.D0)
+        v1 = random()
+      
+      end do
 
-      rnd = DSQRT(-DLOG(x))
+      rnd = sqrt(-2.d0*log(v1))
 
     END FUNCTION randomHalfMaxwellian
 

src/modules/common/velocityDistribution.f90

@@ -53,12 +53,11 @@ module velocityDistribution
     function randomVelMaxwellian(self) result (v)
       use moduleRandom
       implicit none
-
       class(velDistMaxwellian), intent(in):: self
       real(8):: v
       v = 0.D0
 
-      v = self%vTh*randomMaxwellian()/sqrt(2.d0)
+      v = self%vTh*randomMaxwellian()
 
     end function randomVelMaxwellian
 

src/modules/init/moduleInput.f90

@@ -493,6 +493,7 @@ MODULE moduleInput
       CALL config%get(object // '.numColl', collOutput, found)
       CALL config%get(object // '.EMField', emOutput, found)
       call config%get(object // '.boundariesParticle', boundaryParticleOutput, found)
+      call config%get(object // '.injects', injectOutput, found)
 
       CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found)
       IF (.NOT. found) THEN
@@ -628,6 +629,7 @@ MODULE moduleInput
       USE moduleMesh
       USE moduleCaseParam
       USE moduleRefParam
+      use moduleOutput, only: collOutput
       USE OMP_LIB
       USE json_module
       IMPLICIT NONE
@@ -1426,53 +1428,19 @@ MODULE moduleInput
       CHARACTER(2):: iString
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
-      CHARACTER(:), ALLOCATABLE:: speciesName
-      CHARACTER(:), ALLOCATABLE:: name
-      REAL(8):: v
-      REAL(8), ALLOCATABLE:: temperature(:), normal(:)
-      REAL(8):: flow
-      CHARACTER(:), ALLOCATABLE:: units
-      INTEGER:: physicalSurface
-      INTEGER:: particlesPerEdge
-      INTEGER:: sp
 
       CALL config%info('inject', found, n_children = nInject)
-      ALLOCATE(inject(1:nInject))
-      nPartInj = 0
+      ALLOCATE(injects(1:nInject))
       DO i = 1, nInject
         WRITE(iString, '(i2)') i
         object = 'inject(' // trim(iString) // ')'
 
-        !Find species
-        CALL config%get(object // '.species', speciesName, found)
-        sp = speciesName2Index(speciesName)
+        CALL initInject(injects(i)%obj, i, object, config)
 
-        CALL config%get(object // '.name',            name,            found)
-        CALL config%get(object // '.v',               v,               found)
-        CALL config%get(object // '.T',               temperature,     found)
-        CALL config%get(object // '.n',               normal,          found)
-        IF (.NOT. found) THEN
-          ALLOCATE(normal(1:3))
-          normal = 0.D0
-        END IF
-        CALL config%get(object // '.flow',            flow,            found)
-        CALL config%get(object // '.units',           units,           found)
-        CALL config%get(object // '.physicalSurface', physicalSurface, found)
-        particlesPerEdge = 0
-        CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
-
-        CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
-
-        CALL readVelDistr(config, inject(i), object)
+        CALL readVelDistr(config, injects(i)%obj, object)
         
       END DO
 
-      !Allocate array for injected particles
-      IF (nPartInj > 0) THEN
-        ALLOCATE(partInj(1:nPartInj))
-
-      END IF
-      
     END SUBROUTINE readInject
 
     SUBROUTINE readAverage(config)
@@ -1522,7 +1490,7 @@ MODULE moduleInput
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
-      TYPE(injectGeneric), INTENT(inout):: inj
+      class(injectGeneric), INTENT(inout):: inj
       CHARACTER(:), ALLOCATABLE, INTENT(in):: object
       INTEGER:: i
       CHARACTER(2):: iString

src/modules/makefile

@@ -14,7 +14,7 @@ output.o: moduleSpecies.o common.o
 mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o 
 	$(MAKE) -C mesh all
 
-solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
+solver.o: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o
 	$(MAKE) -C solver all
 
 init.o: common.o solver.o moduleInject.o

src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90

@@ -1,4 +1,5 @@
 MODULE moduleMeshOutputGmsh2
+  use moduleOutput, only: prefix, informFileCreation, formatFileName, generateFilePath
 
   CONTAINS
     !Header for mesh format
@@ -16,7 +17,7 @@ MODULE moduleMeshOutputGmsh2
     !Node data subroutines
     !Header
     SUBROUTINE writeGmsh2HeaderNodeData(fileID, title, iteration, time, dimensions, nNodes)
-      use moduleOutput
+      use moduleOutput, only: fmtInt, fmtReal
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: fileID
@@ -50,7 +51,7 @@ MODULE moduleMeshOutputGmsh2
     !Element data subroutines
     !Header
     SUBROUTINE writeGmsh2HeaderElementData(fileID, title, iteration, time, dimensions, nVols)
-      use moduleOutput
+      use moduleOutput, only: fmtInt, fmtReal
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: fileID
@@ -86,7 +87,7 @@ MODULE moduleMeshOutputGmsh2
       USE moduleMesh
       USE moduleRefParam
       USE moduleSpecies
-      USE moduleOutput
+      use moduleOutput, only: outputFormat, calculateOutput, fmtReal
       USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
@@ -141,7 +142,7 @@ MODULE moduleMeshOutputGmsh2
       USE moduleRefParam
       USE moduleCaseParam
       USE moduleCollisions
-      USE moduleOutput
+      use moduleOutput, only: collOutput, fmtInt
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
@@ -199,7 +200,7 @@ MODULE moduleMeshOutputGmsh2
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
-      USE moduleOutput
+      use moduleOutput, only: emOutput
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
@@ -248,7 +249,6 @@ MODULE moduleMeshOutputGmsh2
       USE moduleMesh
       USE moduleRefParam
       USE moduleSpecies
-      USE moduleOutput
       USE moduleAverage
       IMPLICIT NONE
 

src/modules/mesh/inout/text/moduleMeshOutputText.f90

@@ -1,9 +1,14 @@
 module moduleMeshOutputText
+  use moduleOutput, only: prefix,             &
+                          informFileCreation, &
+                          formatFileName,     &
+                          generateFilePath
+
   contains
 
     subroutine writeSpeciesOutput(self, fileID, speciesIndex)
       use moduleMesh
-      use moduleOutput
+      use moduleOutput, only: calculateOutput, outputFormat
       use moduleRefParam, only: L_ref
       implicit none
 
@@ -88,7 +93,6 @@ module moduleMeshOutputText
                                   speciesIndex)
 
       use moduleMesh
-      use moduleOutput
       use moduleAverage
       use moduleRefParam, only: L_ref
       implicit none
@@ -151,8 +155,8 @@ module moduleMeshOutputText
 
     subroutine printCollText(self)
       use moduleMesh
-      use moduleOutput
       use moduleCaseParam, only: timeStep
+      use moduleOutput, only: collOutput
       implicit none
       
       class(meshGeneric), intent(in):: self
@@ -190,6 +194,7 @@ module moduleMeshOutputText
     subroutine printEMText(self)
       use moduleMesh
       use moduleCaseParam, only: timeStep
+      use moduleOutput, only: emOutput
       implicit none
       
       class(meshParticles), intent(in):: self
@@ -236,7 +241,7 @@ module moduleMeshOutputText
 
         write(fileIDMean, '(5(A,","),A)') '"Position (m)"', &
                                           '"Density, mean (m^-3)"', &
-                                          '"Velocity, mean (m s^-1):0"', '"Velocity (m s^-1):1"', '"Velocity (m s^-1):2"', &
+                                          '"Velocity, mean (m s^-1):0"', '"Velocity, mean (m s^-1):1"', '"Velocity, mean (m s^-1):2"', &
                                           '"Temperature, mean (K)"'
 
         fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv')
@@ -245,7 +250,7 @@ module moduleMeshOutputText
 
         write(fileIDDeviation, '(5(A,","),A)') '"Position (m)"', &
                                                '"Density, deviation (m^-3)"', &
-                                               '"Velocity, deviation (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (m s^-1):2"', &
+                                               '"Velocity, deviation (m s^-1):0"', 'Velocity, deviation (m s^-1):1"', 'Velocity, deviation (m s^-1):2"', &
                                                '"Temperature, deviation (K)"'
 
         call writeAverage(self, fileIDMean, fileIDDeviation, s)

src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90

@@ -1,9 +1,11 @@
 MODULE moduleMeshOutputVTU
+  use moduleOutput, only: prefix, informFileCreation, formatFileName, generateFilePath
 
   CONTAINS
 
     SUBROUTINE writeHeader(nNodes, nCells, fileID)
       USE moduleMesh
+      use moduleOutput, only: fmtInt
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: nNodes, nCells
@@ -77,6 +79,7 @@ MODULE moduleMeshOutputVTU
     SUBROUTINE writeMesh(self, fileID)
       USE moduleMesh
       USE moduleRefParam
+      use moduleOutput, only: fmtReal, fmtInt
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
@@ -159,6 +162,7 @@ MODULE moduleMeshOutputVTU
     SUBROUTINE writeCollOutput(self,fileID)
       USE moduleMesh
       USE moduleCollisions
+      use moduleOutput, only: fmtInt
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
@@ -185,6 +189,7 @@ MODULE moduleMeshOutputVTU
     SUBROUTINE writeEM(self, fileID)
       USE moduleMesh
       USE moduleRefParam
+      use moduleOutput, only: fmtReal
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
@@ -381,6 +386,7 @@ MODULE moduleMeshOutputVTU
     SUBROUTINE printEMVTU(self)
       USE moduleMesh
       USE moduleCaseParam, ONLY: timeStep
+      use moduleOutput, only: emOutput
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
@@ -414,7 +420,6 @@ MODULE moduleMeshOutputVTU
 
     SUBROUTINE printAverageVTU(self)
       USE moduleMesh
-      use moduleOutput
       USE moduleSpecies
       IMPLICIT NONE
 

src/modules/mesh/moduleMesh.f90

@@ -1,7 +1,7 @@
 !moduleMesh: General module for Finite Element mesh
 MODULE moduleMesh
   USE moduleList
-  USE moduleOutput
+  USE moduleOutput, only: outputNode, emNode
   USE moduleCollisions
   use moduleSpecies, only: nSpecies
 
@@ -663,50 +663,6 @@ MODULE moduleMesh
 
   end type boundaryParticleGeneric
 
-  interface
-    ! Init interfaces
-    ! wallTemeprature
-    module subroutine wallTemperature_init(boundary, T, c)
-      CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary
-      REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
-
-    end subroutine wallTemperature_init
-
-    module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
-      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
-      real(8), intent(in):: mImpact
-      real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties
-      integer, intent(in):: ion
-      real(8), intent(in):: effTime
-      character(:), allocatable, intent(in):: crossSection
-      real(8), intent(in):: eThreshold
-      integer, optional, intent(in):: electronSecondary
-
-    end subroutine ionization_init
-
-    ! quasiNeutrality
-    module subroutine quasiNeutrality_init(boundary, s_incident)
-      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
-      integer, intent(in):: s_incident
-
-    end subroutine quasiNeutrality_init
-
-    ! outflowAdaptive
-    module subroutine outflowAdaptive_init(boundary, s_incident)
-      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
-      integer, intent(in):: s_incident
-
-    end subroutine outflowAdaptive_init
-
-    ! Get the index of the species from its name
-    module function boundaryParticleName_to_Index(boundaryName) result(bp)
-      character(:), allocatable:: boundaryName
-      integer:: bp
-
-    end function boundaryParticleName_to_Index
-
-  end interface
-
   abstract interface
     ! Apply the effects of the boundary to the particle
     subroutine applyParticle_interface(self, edge, part)
@@ -745,6 +701,18 @@ MODULE moduleMesh
 
   END TYPE boundaryReflection
 
+  interface
+    module subroutine reflection_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryReflection), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine reflection_apply
+
+  end interface
+
   !Absorption boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAbsorption
     CONTAINS
@@ -752,6 +720,18 @@ MODULE moduleMesh
 
   END TYPE boundaryAbsorption
 
+  interface
+    module subroutine absorption_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryAbsorption), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine absorption_apply
+
+  end interface
+
   !Transparent boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryTransparent
     CONTAINS
@@ -759,6 +739,18 @@ MODULE moduleMesh
 
   END TYPE boundaryTransparent
 
+  interface
+    module subroutine transparent_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryTransparent), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine transparent_apply
+
+  end interface
+
   !Symmetry axis
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAxis
     CONTAINS
@@ -766,6 +758,18 @@ MODULE moduleMesh
 
   END TYPE boundaryAxis
 
+  interface
+    module subroutine symmetryAxis_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryAxis), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine symmetryAxis_apply
+
+  end interface
+
   !Wall Temperature boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryWallTemperature
     !Thermal velocity of the wall: square root(Wall temperature X specific heat)
@@ -775,6 +779,24 @@ MODULE moduleMesh
 
   END TYPE boundaryWallTemperature
 
+  interface
+    module subroutine wallTemperature_init(boundary, T, c)
+      CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary
+      REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
+
+    end subroutine wallTemperature_init
+
+    module subroutine wallTemperature_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryWallTemperature), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine wallTemperature_apply
+
+  end interface
+
   !Ionization boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryIonization
     REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
@@ -791,6 +813,30 @@ MODULE moduleMesh
 
   END TYPE boundaryIonization
 
+  interface
+    module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      real(8), intent(in):: mImpact
+      real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties
+      integer, intent(in):: ion
+      real(8), intent(in):: effTime
+      character(:), allocatable, intent(in):: crossSection
+      real(8), intent(in):: eThreshold
+      integer, optional, intent(in):: electronSecondary
+
+    end subroutine ionization_init
+
+    module subroutine ionization_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryIonization), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine ionization_apply
+
+  end interface
+
   ! Ensures quasi-neutrality by changing the reflection coefficient
   type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality
     real(8), allocatable:: alpha(:) ! Reflection parameter
@@ -801,6 +847,24 @@ MODULE moduleMesh
 
   end type boundaryQuasiNeutrality
 
+  interface
+    module subroutine quasiNeutrality_init(boundary, s_incident)
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
+
+    end subroutine quasiNeutrality_init
+
+    module subroutine quasiNeutrality_apply(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryQuasiNeutrality), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine quasiNeutrality_apply
+
+  end interface
+
   !Boundary for quasi-neutral outflow adjusting reflection coefficient
   type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive
     real(8), allocatable:: velocity_shift(:)
@@ -812,68 +876,11 @@ MODULE moduleMesh
   end type boundaryOutflowAdaptive
 
   interface
-    module subroutine reflection_apply(self, edge, part)
-      use moduleSpecies
+    module subroutine outflowAdaptive_init(boundary, s_incident)
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
 
-      class(boundaryReflection), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine reflection_apply
-
-    module subroutine absorption_apply(self, edge, part)
-      use moduleSpecies
-
-      class(boundaryAbsorption), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine absorption_apply
-
-    module subroutine transparent_apply(self, edge, part)
-      use moduleSpecies
-
-      class(boundaryTransparent), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine transparent_apply
-
-    module subroutine symmetryAxis_apply(self, edge, part)
-      use moduleSpecies
-
-      class(boundaryAxis), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine symmetryAxis_apply
-
-    module subroutine wallTemperature_apply(self, edge, part)
-      use moduleSpecies
-
-      class(boundaryWallTemperature), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine wallTemperature_apply
-
-    module subroutine ionization_apply(self, edge, part)
-      use moduleSpecies
-
-      class(boundaryIonization), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine ionization_apply
-
-    module subroutine quasiNeutrality_apply(self, edge, part)
-      use moduleSpecies
-
-      class(boundaryQuasiNeutrality), intent(inout):: self
-      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout):: part
-
-    end subroutine quasiNeutrality_apply
+    end subroutine outflowAdaptive_init
 
     module subroutine outflowAdaptive_apply(self, edge, part)
       use moduleSpecies
@@ -884,7 +891,10 @@ MODULE moduleMesh
 
     end subroutine outflowAdaptive_apply
 
-    ! Generic basic boundary conditions to use internally in the code
+  end interface
+
+  ! Generic basic boundary conditions to use internally in the code
+  interface
     module subroutine genericReflection(edge, part) 
       use moduleSpecies
 
@@ -927,6 +937,13 @@ MODULE moduleMesh
   type(boundaryParticleCont), allocatable, target:: boundariesParticle(:)
 
   interface 
+    ! Get the index of the species from its name
+    module function boundaryParticleName_to_Index(boundaryName) result(bp)
+      character(:), allocatable:: boundaryName
+      integer:: bp
+
+    end function boundaryParticleName_to_Index
+
     ! Update the particle boundary models
     module subroutine boundariesParticle_update()
 
@@ -952,33 +969,6 @@ MODULE moduleMesh
 
   end type boundaryEMGeneric
 
-  interface
-    module subroutine initDirichlet(self, config, object)
-      use json_module
-
-      class(boundaryEMGeneric), allocatable, intent(inout):: self
-      type(json_file), intent(inout):: config
-      character(:), allocatable, intent(in):: object
-
-    end subroutine initDirichlet
-
-    module subroutine initDirichletTime(self, config, object)
-      use json_module
-
-      class(boundaryEMGeneric), allocatable, intent(inout):: self
-      type(json_file), intent(inout):: config
-      character(:), allocatable, intent(in):: object
-
-    end subroutine initDirichletTime
-
-    module function boundaryEMName_to_Index(boundaryName) result(bp)
-      character(:), allocatable:: boundaryName
-      integer:: bp
-
-    end function boundaryEMName_to_Index
-
-  end interface
-
   abstract interface
     ! Apply boundary condition to the load vector for the Poission equation
     subroutine applyEM_interface(self, vectorF)
@@ -1008,6 +998,24 @@ MODULE moduleMesh
 
   END TYPE boundaryEMDirichlet
 
+  interface
+    module subroutine initDirichlet(self, config, object)
+      use json_module
+
+      class(boundaryEMGeneric), allocatable, intent(inout):: self
+      type(json_file), intent(inout):: config
+      character(:), allocatable, intent(in):: object
+
+    end subroutine initDirichlet
+
+    module subroutine applyDirichlet(self, vectorF)
+      class(boundaryEMDirichlet), intent(in):: self
+      real(8), intent(inout):: vectorF(:)
+
+    end subroutine applyDirichlet
+
+  end interface
+
   TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
     real(8):: potential
     real(8):: timeFactor
@@ -1019,11 +1027,14 @@ MODULE moduleMesh
   END TYPE boundaryEMDirichletTime
 
   interface
-    module subroutine applyDirichlet(self, vectorF)
-      class(boundaryEMDirichlet), intent(in):: self
-      real(8), intent(inout):: vectorF(:)
+    module subroutine initDirichletTime(self, config, object)
+      use json_module
 
-    end subroutine applyDirichlet
+      class(boundaryEMGeneric), allocatable, intent(inout):: self
+      type(json_file), intent(inout):: config
+      character(:), allocatable, intent(in):: object
+
+    end subroutine initDirichletTime
 
     module subroutine applyDirichletTime(self, vectorF)
       class(boundaryEMDirichletTime), intent(in):: self
@@ -1047,6 +1058,12 @@ MODULE moduleMesh
 
   ! Update the EM boundary models
   interface 
+    module function boundaryEMName_to_Index(boundaryName) result(bp)
+      character(:), allocatable:: boundaryName
+      integer:: bp
+
+    end function boundaryEMName_to_Index
+
     module subroutine boundariesEM_update()
 
     end subroutine boundariesEM_update

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -56,6 +56,7 @@ submodule(moduleMesh) boundaryEM
       use json_module
       use moduleRefParam, ONLY: Volt_ref, ti_ref
       use moduleErrors
+      use moduleOutput, only: path
       implicit none
 
       class(boundaryEMGeneric), allocatable, intent(inout):: self

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -342,6 +342,7 @@ submodule(moduleMesh) boundaryParticle
 
     ! Print
     subroutine ionization_print(self, fileID)
+      use moduleOutput, only: fmtInt
       implicit none
 
       class(boundaryParticleGeneric), intent(inout):: self
@@ -473,6 +474,7 @@ submodule(moduleMesh) boundaryParticle
 
     ! Print output
     subroutine quasiNeutrality_print(self, fileID)
+      use moduleOutput, only: fmtColInt, fmtColReal
       implicit none
 
       class(boundaryParticleGeneric), intent(inout):: self
@@ -613,6 +615,7 @@ submodule(moduleMesh) boundaryParticle
 
     ! Print
     subroutine outflowAdaptive_print(self, fileID)
+      use moduleOutput, only: fmtColInt, fmtColReal
       implicit none
 
       class(boundaryParticleGeneric), intent(inout):: self
@@ -691,7 +694,12 @@ submodule(moduleMesh) boundaryParticle
     ! Writes the output into the Step_XXXXX_boundaryParticles file when the %print procedure is associated
     module subroutine boundariesParticle_write()
       use moduleCaseparam, only: timeStep
-      use moduleOutput, only:fileID_boundaryParticle, formatFileName, informFileCreation
+      use moduleOutput, only: fileID_boundaryParticle, &
+                              formatFileName, &
+                              informFileCreation,&
+                              boundaryParticleOutput,&
+                              generateFilePath, &
+                              prefix
       implicit none
 
       integer:: b
@@ -700,7 +708,7 @@ submodule(moduleMesh) boundaryParticle
       if (boundaryParticleOutput) then
         fileName = formatFileName(prefix, 'boundariesParticle', 'csv', timeStep)
         call informFileCreation(fileName)
-        open(fileID_boundaryParticle, file = path // folder // '/' //  fileName)
+        open(fileID_boundaryParticle, file = generateFilePath(fileName))
 
         do b = 1, nBoundariesParticle
           if (associated(boundariesParticle(b)%obj%print)) then

src/modules/moduleInject.f90

@@ -5,14 +5,14 @@ MODULE moduleInject
   use moduleMesh, only: meshEdgePointer
 
   !Generic injection of particles
-  TYPE:: injectGeneric
+  TYPE, abstract:: injectGeneric
     INTEGER:: id
     CHARACTER(:), ALLOCATABLE:: name
     REAL(8):: vMod !Velocity (module)
     REAL(8):: temperature(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
     LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
-    INTEGER:: nParticles !Number of particles to introduce each time step
+    INTEGER:: nParticles !Number of particles to inject each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     type(meshEdgePointer), allocatable:: edges(:)
@@ -20,38 +20,125 @@ MODULE moduleInject
     REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
     REAL(8):: surface ! Total surface of injection
     TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
+    procedure(updateInject_interface), pointer, pass:: update => null()
+    procedure(printInject_interface),  pointer, pass:: print  => null()
     CONTAINS
-      PROCEDURE, PASS:: init         => initInject
       PROCEDURE, PASS:: addParticles
 
   END TYPE injectGeneric
 
+  abstract interface
+    ! Update the values of the particle boundary model
+    subroutine updateInject_interface(self)
+      import injectGeneric
+
+      class(injectGeneric), intent(inout):: self
+
+    end subroutine updateInject_interface
+
+    ! Write the values of the particle boundary model
+    subroutine printInject_interface(self, fileID)
+      import injectGeneric
+
+      class(injectGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+
+    end subroutine printInject_interface
+
+  end interface
+
+  ! Default type, constant inject of particles
+  type, extends(injectGeneric):: injectConstant
+
+  end type injectconstant
+
+  type, extends(injectGeneric):: injectQuasiNeutral
+    real(8):: qSign_incident
+
+  end type injectQuasiNeutral
+
+  ! Wrapper for injects
+  type:: injectWrapper
+    class(injectGeneric), allocatable:: obj
+
+  end type injectWrapper
+
   INTEGER:: nInject
-  TYPE(injectGeneric), ALLOCATABLE:: inject(:)
+  TYPE(injectWrapper), ALLOCATABLE, target:: injects(:)
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, ps, particlesPerEdge)
+    SUBROUTINE initInject(self, i, object, config)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
-      USE moduleSolver
       USE moduleErrors
+      use json_module
       IMPLICIT NONE
 
-      CLASS(injectGeneric), INTENT(inout):: self
+      CLASS(injectGeneric), allocatable, INTENT(inout):: self
       INTEGER, INTENT(in):: i
-      REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
-      INTEGER, INTENT(in):: sp, ps, particlesPerEdge
+      character(:), allocatable, intent(in):: object
+      type(json_file), intent(inout):: config
+      character(:), allocatable:: type
+      logical:: found
+      REAL(8):: v
+      real(8), allocatable:: n(:), temperature(:)
+      INTEGER:: sp, ps, particlesPerEdge
+      character(:), allocatable:: speciesName
       integer:: ps_index
       REAL(8):: tauInject
-      REAL(8), INTENT(in):: flow
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: units
+      REAL(8):: flow
+      CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: e
       REAL(8):: fluxPerStep = 0.D0
 
-      self%id          =  i
+      !Type of injection
+      call config%get(object // '.type',            type,        found)
+      if (.not. found) then
+        ! If no type is found, assume constant injection of particles
+        type = 'constant'
+      end if
+
+      ! Assign type to inject and allocates it
+      select case(type)
+      case ('constant')
+        allocate(injectConstant:: self)
+
+      case ('quasiNeutral')
+        allocate(injectQuasiNeutral:: self)
+
+        self%update => updateQuasiNeutral
+        self%print  => printQuasiNeutral
+
+      case default
+        call criticalError('No injection type ' // type // ' defined', 'initInject')
+      end select
+
+      self%id =  i
+      call config%get(object // '.name',             self%name,        found)
+
+      !Find species
+      CALL config%get(object // '.species', speciesName, found)
+      sp = speciesName2Index(speciesName)
+
+      ! Velocity and direction of injection
+      call config%get(object // '.v',                v,                found)
+      call config%get(object // '.n',                n,                found)
+      if (.not. found) then
+        allocate(n(1:3))
+        n = 0.D0
+      end if
+      ! Temperature of injection
+      call config%get(object // '.T',                temperature,      found)
+      ! Flow
+      call config%get(object // '.flow',             flow,             found)
+      call config%get(object // '.units',            units,            found)
+      CALL config%get(object // '.physicalSurface',  ps,               found)
+      particlesPerEdge = 0
+      CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
+
       self%vMod        =  v / v_ref
       self%n           =  n / NORM2(n)
       self%temperature =  temperature / T_ref
@@ -124,8 +211,6 @@ MODULE moduleInject
 
         END DO
 
-        self%nParticles = SUM(self%particlesPerEdge)
-
       ELSE
         ! No particles assigned per edge, use the species weight
         self%weightPerEdge = self%species%weight
@@ -133,56 +218,235 @@ MODULE moduleInject
           self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight))
         END DO
 
-        self%nParticles = SUM(self%particlesPerEdge)
-
         !Rescale weight to match flux
         self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles))
 
       END IF
 
+      ! Total number of particles to inject
+      self%nParticles = SUM(self%particlesPerEdge)
+
       !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
+      ! Specific parameters per type
+      select type(self)
+      type is(injectQuasiNeutral)
+        ! Get the sign of the charged species
+        self%qSign_incident = sign(1.d0, self%species%qm)
+
+      end select
+
     END SUBROUTINE initInject
 
-    !Injection of particles
-    SUBROUTINE doInjects()
+    ! Update the value of injects
+    subroutine updateInjects()
+      implicit none
+      
+      integer:: i
+
+      do i = 1, nInject
+        if (associated(injects(i)%obj%update)) then
+          call injects(i)%obj%update()
+
+        end if
+
+      end do
+
+    end subroutine updateInjects
+
+    ! Updates the flow in the injection to maintain quasineutrality
+    subroutine updateQuasiNeutral(self)
+      use moduleMesh, only: meshEdge, meshCell
+      implicit none
+
+      class(injectGeneric), intent(inout):: self
+      integer:: e
+      class(meshEdge), pointer:: edge
+      class(meshCell), pointer:: cell
+      real(8):: Xi(1:3)
+      real(8):: EF_normal
+      real(8):: alpha
+
+      select type(self)
+      type is (injectQuasiNeutral)
+        do e = 1, self%nEdges
+          edge => self%edges(e)%obj
+
+          Xi = edge%centerXi()
+
+          if (associated(edge%e1)) then
+            cell => edge%e1
+
+          else
+            cell => edge%e2
+
+          end if
+
+          ! Projection of EF on the edge normal vector
+          EF_normal = dot_product(cell%gatherElectricField(Xi), edge%normal)
+
+          ! Correction for this time step
+          alpha = self%qSign_incident * EF_normal
+
+          ! Adjust the weight of particles to match the new current
+          self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha
+
+          ! Limit the weight of a macroparticle to 1.0
+          self%weightPerEdge(e) = max(self%weightPerEdge(e), 1.d0)
+
+        end do
+
+      end select
+
+    end subroutine updateQuasiNeutral
+
+    !Add particles for the injection
+    SUBROUTINE addParticles(self)
       USE moduleSpecies
-      USE moduleSolver
+      USE moduleMesh
+      USE moduleRandom
+      USE moduleErrors
+      use moduleList, only: partInj
       IMPLICIT NONE
 
-      INTEGER:: i
+      CLASS(injectGeneric), INTENT(in):: self
+      type(particle), pointer:: part
+      INTEGER:: i, e
+      INTEGER:: sp
+      CLASS(meshEdge), POINTER:: edge
+
+      REAL(8):: direction(1:3)
+
+      !Insert particles
+      !$OMP DO
+      DO e = 1, self%nEdges
+        ! Select edge for injection
+        edge => self%edges(e)%obj
+        ! Inject particles in edge
+        DO i = 1, self%particlesPerEdge(e)
+          allocate(part)
+          ! Index in the partInj array
+          !Particle is considered to be inside the domain
+          part%n_in = .true.
+          !Random position in edge
+          part%r = edge%randPos()
+          !Assign weight to particle.
+          part%weight = self%weightPerEdge(e)
+          !Volume associated to the edge:
+          IF (ASSOCIATED(edge%e1)) THEN
+            part%cell = edge%e1%n
+
+          ELSEIF (ASSOCIATED(edge%e2)) THEN
+            part%cell = edge%e2%n
+
+          ELSE
+            CALL criticalError("No Volume associated to edge", 'addParticles')
+
+          END IF
+          part%cellColl = edge%eColl%n
+          sp = self%species%n
+
+          !Assign particle type
+          part%species => self%species
+
+          if (all(self%n == 0.D0)) then
+            direction =  edge%normal
+
+          else
+            direction = self%n
+
+          end if
+
+          part%v = 0.D0
+
+          ! do while(dot_product(partInj(n)%v, direction) <= 0.d0)
+            part%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+                                                     self%v(2)%obj%randomVel(), &
+                                                     self%v(3)%obj%randomVel() /)
+            !If injecting a no-drift distribution and velocity is negative, reflect
+            if ((self%vMod == 0.D0) .and. &
+                (dot_product(part%v, direction) <= 0.D0)) then
+              part%v = - part%v
+
+            end if
+
+          ! end do
+
+          !Obtain natural coordinates of particle in cell
+          part%Xi = mesh%cells(part%cell)%obj%phy2log(part%r)
+
+          ! Add particle to global list
+          call partInj%setLock()
+          call partInj%add(part)
+          call partInj%unsetLock()
+
+        END DO
+
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE addParticles
+
+    ! Writes the value of injects
+    subroutine writeInjects()
+      use moduleCaseparam, only: timeStep
+      use moduleOutput, only: fileID_inject, &
+                              formatFileName, &
+                              informFileCreation,&
+                              injectOutput,&
+                              generateFilePath, &
+                              prefix
+      implicit none
       
-      !$OMP SINGLE
-      nPartInj = 0
-      DO i = 1, nInject
-        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          nPartInj = nPartInj + inject(i)%nParticles
+      integer:: i
+      character(:), allocatable:: fileName
 
-        END IF
+      if (injectOutput) then
+        fileName = formatFileName(prefix, 'injects', 'csv', timeStep)
+        call informFileCreation(fileName)
+        open(fileID_inject, file = generateFilePath(fileName))
 
-      END DO
+        do i = 1, nInject
+          if (associated(injects(i)%obj%print)) then
+            call injects(i)%obj%print(fileID_inject)
 
-      IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
-      ALLOCATE(partInj(1:nPartInj))
-      !$OMP END SINGLE
+          end if
 
-      DO i=1, nInject
-        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          CALL inject(i)%addParticles()
+        end do
 
-        END IF
-      END DO
+        close(fileID_inject)
 
-    END SUBROUTINE doInjects
+      end if
 
+    end subroutine writeInjects
+
+    ! Write information about quasiNeutral inject
+    subroutine printQuasiNeutral(self, fileID)
+      use moduleOutput, only: fmtColInt, fmtColReal
+      implicit none
+
+      class(injectGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+      integer:: e
+
+      write(fileID, '(A)') self%name
+      write(fileID, '(A,",",A)') '"Edge id"', '"weight of particles"'
+      do e = 1, self%nEdges
+        write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%weightPerEdge(e)
+
+      end do
+
+    end subroutine printQuasiNeutral
+
+    ! Inits for different velocity distribution functions
     SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
@@ -192,7 +456,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
 
     END SUBROUTINE initVelDistHalfMaxwellian
 
@@ -205,101 +469,4 @@ MODULE moduleInject
 
     END SUBROUTINE initVelDistDelta
 
-    !Add particles for the injection
-    SUBROUTINE addParticles(self)
-      USE moduleSpecies
-      USE moduleSolver
-      USE moduleMesh
-      USE moduleRandom
-      USE moduleErrors
-      IMPLICIT NONE
-
-      CLASS(injectGeneric), INTENT(in):: self
-      INTEGER, SAVE:: nMin
-      INTEGER:: i, e
-      INTEGER:: n, sp
-      CLASS(meshEdge), POINTER:: randomEdge
-      REAL(8):: direction(1:3)
-
-      !Insert particles
-      !$OMP SINGLE
-      nMin = 0
-      DO i = 1, self%id -1
-        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          nMin = nMin + inject(i)%nParticles
-
-        END IF
-
-      END DO
-      nMin = nMin + 1
-      !$OMP END SINGLE
-
-      !$OMP DO
-      DO e = 1, self%nEdges
-        ! Select edge for injection
-        randomEdge => self%edges(e)%obj
-        ! Inject particles in edge
-        DO i = 1, self%particlesPerEdge(e)
-          ! Index in the global partInj array
-          n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i
-          !Particle is considered to be outside the domain
-          partInj(n)%n_in = .FALSE.
-          !Random position in edge
-          partInj(n)%r = randomEdge%randPos()
-          !Assign weight to particle.
-          partInj(n)%weight = self%weightPerEdge(e)
-          !Volume associated to the edge:
-          IF (ASSOCIATED(randomEdge%e1)) THEN
-            partInj(n)%cell = randomEdge%e1%n
-
-          ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
-            partInj(n)%cell = randomEdge%e2%n
-
-          ELSE
-            CALL criticalError("No Volume associated to edge", 'addParticles')
-
-          END IF
-          partInj(n)%cellColl = randomEdge%eColl%n
-          sp = self%species%n
-
-          !Assign particle type
-          partInj(n)%species => self%species
-
-          if (all(self%n == 0.D0)) then
-            direction =  randomEdge%normal
-
-          else
-            direction = self%n
-
-          end if
-
-          partInj(n)%v = 0.D0
-
-          do while(dot_product(partInj(n)%v, direction) <= 0.d0)
-            partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
-                                                     self%v(2)%obj%randomVel(), &
-                                                     self%v(3)%obj%randomVel() /)
-            !If injecting a no-drift distribution and velocity is negative, reflect
-            if ((self%vMod == 0.D0) .and. &
-                (dot_product(partInj(n)%v, direction) <= 0.D0)) then
-              partInj(n)%v = - partInj(n)%v
-
-            end if
-
-          end do
-
-          !Obtain natural coordinates of particle in cell
-          partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
-          !Push new particle with the minimum time step
-          CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
-          !Assign cell to new particle
-          CALL solver%updateParticleCell(partInj(n))
-
-        END DO
-
-      END DO
-      !$OMP END DO
-
-    END SUBROUTINE addParticles
-
 END MODULE moduleInject

src/modules/moduleList.f90

@@ -23,6 +23,7 @@ MODULE moduleList
 
   END TYPE listNode
 
+  TYPE(listNode):: partInj !Particles comming from injections
   TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme
   TYPE(listNode):: partCollisions !Particles created in collisional process
   TYPE(listNode):: partSurfaces !Particles created in surface interactions

src/modules/moduleSpecies.f90

@@ -48,11 +48,8 @@ MODULE moduleSpecies
 
   !Number of old particles
   INTEGER:: nPartOld
-  !Number of injected particles
-  INTEGER:: nPartInj
   !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
 
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(sp)

src/modules/output/moduleOutput.f90

@@ -14,6 +14,7 @@ MODULE moduleOutput
   LOGICAL:: collOutput = .FALSE.
   LOGICAL:: emOutput   = .FALSE.
   logical:: boundaryParticleOutput = .false.
+  logical:: injectOutput = .false.
 
   ! Prefix for iteration files
   character(len=*), parameter:: prefix = 'Step'
@@ -33,6 +34,7 @@ MODULE moduleOutput
   integer, parameter:: fileID_output = 20           ! Base id for species/collisions/EM output
   integer, parameter:: fileID_boundaryParticle = 30 ! Particle boundaries
   integer, parameter:: fileID_boundaryEM = 31       ! EM boundaries
+  integer, parameter:: fileID_inject = 32          ! Injects
   integer, parameter:: fileID_reference = 40        ! Reference values
   integer, parameter:: fileID_time =50              ! Computation time
 

src/modules/solver/moduleSolver.f90

@@ -235,23 +235,21 @@ MODULE moduleSolver
 
       INTEGER:: nn, n, e
       INTEGER, SAVE:: nPartNew
-      INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
+      INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
       INTEGER:: s
 
       !$OMP SECTIONS
       !$OMP SECTION
-      nInjIn = 0
-      IF (ALLOCATED(partInj)) THEN
-        nInjIn = COUNT(partInj%n_in)
-
-      END IF
-      !$OMP SECTION
       nOldIn = 0
       IF (ALLOCATED(partOld)) THEN
         nOldIn = COUNT(partOld%n_in)
 
       END IF
+
+      !$OMP SECTION
+      nInj = partInj%amount
+
       !$OMP SECTION
       nWScheme = partWScheme%amount
 
@@ -267,30 +265,14 @@ MODULE moduleSolver
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces
+      nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces
       ALLOCATE(partOld(1:nPartNew))
       !$OMP END SINGLE
 
       !$OMP SECTIONS
-      !$OMP SECTION
-      !Reset particles from injection
-      nn = 0
-      DO n = 1, nPartInj
-        IF (partInj(n)%n_in) THEN
-          nn = nn + 1
-          partOld(nn) = partInj(n)
-          IF (nProbes > 0) THEN
-            CALL doProbes(partOld(nn))
-
-          END IF
-
-        END IF
-
-      END DO
-
       !$OMP SECTION
       !Reset particles from previous iteration
-      nn = nInjIn
+      nn = 0
       DO n = 1, nPartOld
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
@@ -304,19 +286,24 @@ MODULE moduleSolver
 
       END DO
 
+      !$OMP SECTION
+      !Reset particles from injection
+      nn = nOldIn
+      call resetList(partInj, partOld, nn)
+
       !$OMP SECTION
       !Reset particles from weighting scheme
-      nn = nInjIn + nOldIn
+      nn = nOldIn + nInj
       CALL resetList(partWScheme, partOld, nn)
 
       !$OMP SECTION
       !Reset particles from collisional process
-      nn = nInjIn + nOldIn + nWScheme
+      nn = nOldIn + nInj + nWScheme
       CALL resetList(partCollisions, partOld, nn)
 
       !$OMP SECTION
       !Reset particles from surface process
-      nn = nInjIn + nOldIn + nWScheme + nCollisions
+      nn = nOldIn + nInj + nWScheme + nCollisions
       CALL resetList(partSurfaces, partOld, nn)
 
       !$OMP SECTION
@@ -473,6 +460,28 @@ MODULE moduleSolver
 
     END SUBROUTINE splitParticle
 
+    !Injection of particles
+    SUBROUTINE doInjects()
+      USE moduleSpecies
+      use moduleInject
+      IMPLICIT NONE
+
+      INTEGER:: i, sp
+      
+      DO i=1, nInject
+        associate(inject => injects(i)%obj)
+        sp = inject%species%n
+        IF (solver%pusher(sp)%pushSpecies) THEN 
+
+          CALL inject%addParticles()
+
+        END IF
+        end associate
+
+      END DO
+
+    END SUBROUTINE doInjects
+
     SUBROUTINE updateParticleCell(self, part)
       USE moduleSpecies
       USE moduleMesh
@@ -520,6 +529,7 @@ MODULE moduleSolver
       USE moduleSpecies
       USE moduleCompTime
       USE moduleProbe
+      use moduleInject, only: writeInjects
       USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
@@ -539,6 +549,9 @@ MODULE moduleSolver
         ! Output of boundaries
         call boundariesParticle_write()
 
+        ! Output of inejcts
+        call writeInjects()
+
         WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
         IF (timeStep == 0) THEN