First version with possibility for charged particles to be included.

Now, the solver needs to be an input parameter of the case, to select if
it is for charged or neutral particles.

Resolution of Poisson equation with Dirichlet boundary conditions is
possible. The source vector is the charge density. This resolution is
done in two steps to save computational time:
  1. When reading the mesh, the PLU factorization of the K matrix is
  computed.
  2. In each iteration, the system K*u = f is solved, in which f is the
  source vector (charge density) and u is the solution (potential) in
  each node.

No case has been added to the repository. This will be done in next
commit.

The 'non-analog' scheme has been commented. It still needs to split
the particle to avoid 'overweight' particles.

.gitignore

@@ -1,4 +1,4 @@
-DSMC_Neutrals
+PPartiC
 mod/
 obj/
 doc/user_manual/

makefile

@@ -12,20 +12,24 @@ JSONINC := $(JSONDIR)/include
 INCDIR := 
 INCDIR += $(JSONINC)
 
+VAR := ""
+
 # compile flags
-FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -fno-stack-arrays -g
+FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -g
 
 #Output file
-OUTPUT = DSMC_Neutrals
+OUTPUT = PPartiC
 
 export
 
-$(OUTPUT): src.o 
+all: $(OUTPUT)
+
+$(OUTPUT):  src.o 
 
 src.o:
 	@mkdir -p $(MODDIR)
 	@mkdir -p $(OBJDIR)
-	$(MAKE) -C src all
+	$(MAKE) -C src $(OUTPUT)
 
 clean:
 	rm -f $(OUTPUT)

runs/cylFlow/input.json

@@ -24,14 +24,15 @@
     {"name": "Nozzle", "species": "Argon", "flow": 10.0, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], "n": [1, 0, 0], "physicalSurface": 1}
   ],
   "reference": {
-	"density": 1.0e20,
+    "density": 1.0e20,
-	"mass": 6.633e-26,
+    "mass": 6.633e-26,
-	"temperature": 300.0,
+    "temperature": 300.0,
-	"radius": 1.88e-10
+    "radius": 1.88e-10
   },
   "case": {
-    "tau":  6.e-7,
+    "tau":  6.0e-7,
-    "time": 3.e-3
+    "time": 3.0e-3,
+    "solver": "2DCylNeutral"
   },
   "interactions": {
     "folderCollisions": "./data/collisions/",

src/PPartiC.f95

@@ -1,5 +1,5 @@
-!  PPartiC neutral DSMC main program.
+!  PPartiC neutral PIC main program.
-PROGRAM DSMC_Neutrals
+PROGRAM PPartiC
   USE moduleInput
   USE moduleErrors
   USE OMP_LIB
@@ -10,28 +10,32 @@ PROGRAM DSMC_Neutrals
   USE moduleMesh
   IMPLICIT NONE
 
-  ! t = time step; n = number of particle; e = number of element; i = number of inject
+  ! t = time step
   INTEGER:: t = 0
+  ! arg1 = Input argument 1 (input file)
   CHARACTER(200):: arg1
+  ! inputFile = path+name of input file
   CHARACTER(:), ALLOCATABLE:: inputFile
 
-
   tStep = omp_get_wtime()
   !Gets the input file
   CALL getarg(1, arg1)
   inputFile = TRIM(arg1)
   !If no input file is declared, a critical error is called
-  IF (inputFile == "") CALL criticalError("No input file provided", "DSMC_Neutrals")
+  IF (inputFile == "") CALL criticalError("No input file provided", "PPartiC")
 
   !Reads the json configuration file
   CALL readConfig(inputFile)
 
+  CALL verboseError("Calculating initial EM field...")
+  CALL doEMField()
+
   tStep = omp_get_wtime() - tStep
 
   !Output initial state
   CALL doOutput(t)
   CALL verboseError('Starting main loop...')
-  !$OMP PARALLEL PRIVATE(t) DEFAULT(SHARED) 
+  !$OMP PARALLEL DEFAULT(SHARED) 
   DO t = 1, tmax
     !Insert new particles and push them
     !$OMP SINGLE
@@ -46,18 +50,18 @@ PROGRAM DSMC_Neutrals
 
     !$OMP SINGLE
     tPush = omp_get_wtime() - tPush
+    !Collisions
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
-    !Collisions
     CALL doCollisions()
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl
+    !Reset particles
     tReset = omp_get_wtime()
     !$OMP END SINGLE
 
-    !Reset particles
     CALL doReset()
 
     !$OMP SINGLE
@@ -71,6 +75,13 @@ PROGRAM DSMC_Neutrals
 
     !$OMP SINGLE
     tWeight = omp_get_wtime() - tWeight
+    tEMField = omp_get_wtime()
+    !$OMP END SINGLE
+
+    CALL doEMField()
+
+    !$OMP SINGLE
+    tEMField = omp_get_wtime() - tEMField
     tStep = omp_get_wtime() - tStep
     !Output data
     CALL doOutput(t)
@@ -79,5 +90,5 @@ PROGRAM DSMC_Neutrals
   END DO
   !$OMP END PARALLEL
 
-END PROGRAM DSMC_Neutrals
+END PROGRAM PPartiC
 

src/makefile

@@ -1,17 +1,16 @@
-OBJECTS = $(OBJDIR)/$(OUTPUT).o \
+OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
-	        $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
 			  	$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
-					$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
+			  	$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
-					$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
+			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
-					$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o
+			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
-
+			  	$(OBJDIR)/moduleEM.o
 
 all: $(OUTPUT)
-	$(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(JSONLIB)
 
 $(OUTPUT): modules.o $(OUTPUT).f95
 	$(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f95
+	$(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(OBJDIR)/$(OUTPUT).o $(JSONLIB) -L/usr/local/lib -llapack -lopenblas
 
 modules.o:
 	$(MAKE) -C modules all

src/modules/makefile

@@ -1,9 +1,9 @@
 
 OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o\
-			 moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\
+		   moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\
-			 moduleMeshCylRead.o moduleOutput.o moduleInput.o \
+		   moduleMeshCylRead.o moduleOutput.o moduleInput.o \
-			 moduleSolver.o moduleCollisions.o moduleTable.o \
+		   moduleSolver.o moduleCollisions.o moduleTable.o \
-			 moduleParallel.o
+		   moduleParallel.o moduleEM.o
 
 all: $(OBJS)
 
@@ -40,11 +40,14 @@ moduleRefParam.o: moduleConstParam.o moduleRefParam.f95
 moduleSpecies.o: moduleCaseParam.o moduleSpecies.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
-moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
+moduleSolver.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
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
+moduleEM.o: moduleSpecies.o moduleMesh.o moduleEM.f95
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
+
 %.o: %.f95
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

src/modules/moduleCompTime.f95

@@ -5,6 +5,7 @@ MODULE moduleCompTime
   REAL(8):: tReset=0.D0
   REAL(8):: tColl=0.D0
   REAL(8):: tWeight=0.D0
+  REAL(8):: tEMField=0.D0
 
 END MODULE moduleCompTime
 

src/modules/moduleConstParam.f95

@@ -6,9 +6,10 @@ MODULE moduleConstParam
 
   REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi
   REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1
-  REAL(8), PARAMETER:: eV2J = 1.60217662D-19 !Electron volt to Joule (elementary charge)
+  REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge
   REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI
-  ! REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0 (Not used in Neutrals)
+  REAL(8), PARAMETER:: eV2J = qe !Electron volt to Joule conversion
+  REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0
 
 END MODULE moduleConstParam
 

src/modules/moduleEM.f95

@@ -0,0 +1,133 @@
+!Module to solve the electromagnetic field
+MODULE moduleEM
+  IMPLICIT NONE
+
+  TYPE:: boundaryEM
+    CHARACTER(:), ALLOCATABLE:: typeEM
+    INTEGER:: physicalSurface
+    REAL(8):: potential
+
+    CONTAINS
+      PROCEDURE, PASS:: apply
+
+  END TYPE boundaryEM
+
+  INTEGER:: nBoundaryEM
+  TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
+
+  !Information of charge and reference parameters for rho vector
+  REAL(8), ALLOCATABLE:: qSpecies(:)
+
+  CONTAINS
+    SUBROUTINE apply(self, edge)
+      USE moduleMesh
+      IMPLICIT NONE
+
+      CLASS(boundaryEM), INTENT(in):: self
+      CLASS(meshEdge):: edge
+      INTEGER, ALLOCATABLE:: nodes(:)
+      INTEGER:: n
+
+      nodes = edge%getNodes()
+
+      DO n = 1, SIZE(nodes)
+        SELECT CASE(self%typeEM)
+        CASE ("dirichlet")
+          mesh%K(nodes(n), :)        = 0.D0
+          mesh%K(nodes(n), nodes(n)) = 1.D0
+        
+          mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
+          mesh%nodes(nodes(n))%obj%emData%phi  = self%potential
+
+        END SELECT
+
+      END DO
+
+    END SUBROUTINE
+
+    PURE FUNCTION gatherElecField(part) RESULT(elField)
+      USE moduleSpecies
+      USE moduleMesh
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(in):: part
+      REAl(8):: xi(1:3) !Logical coordinates of particle in element
+      REAL(8):: elField(1:3)
+
+      elField = 0.D0
+
+      xi  = part%xLog
+
+      elField = mesh%vols(part%e_p)%obj%gatherEF(xi)
+
+    END FUNCTION gatherElecField
+
+    SUBROUTINE assembleSourceVector(vectorF)
+      USE moduleMesh
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      REAL(8), INTENT(out):: vectorF(1:mesh%numNodes)
+      REAL(8), ALLOCATABLE:: localF(:)
+      INTEGER, ALLOCATABLE:: nodes(:)
+      REAL(8), ALLOCATABLE:: rho(:)
+      INTEGER:: nNodes
+      INTEGER:: e, i, ni
+      CLASS(meshNode), POINTER:: node
+
+      !$OMP SINGLE
+      vectorF = 0.D0
+      !$OMP END SINGLE
+
+      !$OMP SINGLE
+      DO e = 1, mesh%numVols
+        nodes  = mesh%vols(e)%obj%getNodes()
+        nNodes = SIZE(nodes)
+        !Calculates charge density (rho) in element nodes
+        ALLOCATE(rho(1:nNodes))
+        rho = 0.D0
+        DO i = 1, nNodes
+          ni = nodes(i)
+          node => mesh%nodes(ni)%obj
+          rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))
+
+        END DO
+
+        !If there is charge in the nodes, compute the local F vector
+        IF (ANY(rho > 0.D0)) THEN
+          !Calculates local F vector
+          localF = mesh%vols(e)%obj%elemF(rho)
+
+          !Assign local F to global F
+          DO i = 1, nNodes
+            ni = nodes(i)
+            vectorF(ni) = vectorF(ni) + localF(i)
+
+          END DO
+
+          DEALLOCATE(localF)
+
+        END IF
+
+        DEALLOCATE(nodes, rho)
+
+      END DO
+      !$OMP END SINGLE
+
+      !Apply boundary conditions
+      !$OMP DO
+      DO i = 1, mesh%numNodes
+        node => mesh%nodes(i)%obj
+
+        SELECT CASE(node%emData%type)
+        CASE ("dirichlet")
+          vectorF(i) = node%emData%phi
+
+        END SELECT
+
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE assembleSourceVector
+
+END MODULE moduleEM

src/modules/moduleErrors.f95

@@ -37,7 +37,6 @@ MODULE moduleErrors
       errorMsg  = msg
 
       WRITE (*, '(A)') errorMsg
-      WRITE (*, *)
 
     END SUBROUTINE verboseError
 

src/modules/moduleInject.f95

@@ -25,24 +25,14 @@ MODULE moduleInject
   TYPE(injectGeneric), ALLOCATABLE:: inject(:)
 
   CONTAINS
-    !Injection of particles
+    !Initialize an injection of particles
-    SUBROUTINE doInjects()
+    SUBROUTINE initInject(self, i, v, n, T, flow, units, pt, physicalSurface)
-      IMPLICIT NONE
-
-      INTEGER:: i
-
-      DO i=1, nInject
-        CALL inject(i)%addParticles()
-      END DO
-
-    END SUBROUTINE doInjects
-
-    SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface)
       USE moduleMesh
       USE moduleMeshCyl
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
+      USE moduleErrors
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(inout):: self
@@ -50,6 +40,7 @@ MODULE moduleInject
       REAL(8), INTENT(in):: v, n(1:3), T(1:3)
       INTEGER, INTENT(in):: pt, physicalSurface
       REAL(8), INTENT(in):: flow
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: units
       INTEGER:: e, et
       INTEGER:: phSurface(1:mesh%numEdges)
 
@@ -57,8 +48,21 @@ MODULE moduleInject
       self%vMod       = v/v_ref
       self%n          = n
       self%T          = T/T_ref
-      self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
+      SELECT CASE(units)
-      self%pt         = pt
+      CASE ("sccm")
+        !Standard cubic centimeter per minute
+        self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
+
+      CASE ("A")
+        !Input current in Ampers
+        self%nParticles = INT(flow*tau*ti_ref/(qe*species(pt)%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%v   = self%v * self%n
       self%vTh = DSQRT(self%T/species(self%pt)%obj%m)
@@ -92,6 +96,19 @@ MODULE moduleInject
 
     END SUBROUTINE
 
+    !Injection of particles
+    SUBROUTINE doInjects()
+    
+      IMPLICIT NONE
+
+      INTEGER:: i
+
+      DO i=1, nInject
+        CALL inject(i)%addParticles()
+      END DO
+
+    END SUBROUTINE doInjects
+
     !Random velocity from maxwellian distribution
     REAL(8) FUNCTION vBC(u, vth)
       USE moduleConstParam, ONLY: PI
@@ -119,18 +136,28 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(in):: self
       REAL(8):: randomX
       INTEGER:: j
-      REAL(8):: randomPos
+      INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
-      INTEGER:: 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
-      nMax = nMin + self%nParticles
+      nMax = nMin + self%nParticles - 1
       !Assign particle type
       partInj(nMin:nMax)%pt = self%pt
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = species(self%pt)%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)
+      TYPE IS(speciesCharged)
+        partInj(nMin:nMax)%qm = sp%qm
+
+      END SELECT
+      !$OMP END SINGLE
+
       !$OMP DO
       DO n = nMin, nMax
         !Select edge randomly from which inject particle
@@ -143,8 +170,7 @@ MODULE moduleInject
 
         randomEdge => mesh%edges(self%edges(j))%obj
         !Random position in edge
-        randomPos = RAND()
+        partInj(n)%r = randomEdge%randPos()
-        partInj(n)%r = randomEdge%randPos(randomPos)
         !Volume associated to the edge:
         IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN
           partInj(n)%e_p = randomEdge%e1%n
@@ -159,7 +185,7 @@ MODULE moduleInject
                           vBC(self%v(3), self%vTh(3)) /)
 
         !Push new particle
-        CALL push(partInj(n))
+        CALL solver%push(partInj(n))
         !Assign cell to new particle
         CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n))
 

src/modules/moduleInput.f95

@@ -44,6 +44,9 @@ MODULE moduleInput
       CALL verboseError('Reading Geometry...')
       CALL readGeometry(config)
 
+      !Read boundary for EM field
+      CALL readEMBoundary(config)
+
       !Read injection of particles
       CALL verboseError('Reading Interactions between species...')
       CALL readInject(config)
@@ -63,7 +66,7 @@ MODULE moduleInput
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
-      LOGICAL:: found
+      LOGICAL:: found, found_r
       CHARACTER(:), ALLOCATABLE:: object
 
       object = 'reference'
@@ -77,15 +80,25 @@ MODULE moduleInput
       CALL config%get(object // '.temperature', T_ref, found)
       IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
 
-      CALL config%get(object // '.radius', r_ref, found)
+      CALL config%get(object // '.radius', r_ref, found_r)
-      IF (.NOT. found) CALL criticalError('Reference radius not found','readReference')
+      
-
       !Derived parameters
       v_ref     = DSQRT(kb*T_ref/m_ref)  !reference velocity
-      sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
+      !TODO: Make this solver dependent
-      L_ref     = 1.D0/(sigma_ref*n_ref) !mean free path
+      IF (found_r) THEN
+        L_ref = 1.D0/(sigma_ref*n_ref)     !mean free path
+        sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
+
+      ELSE
+        L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length
+        !TODO: Obtain right sigma_ref for PIC case
+        sigma_ref = PI*(4.D-10)**2        !reference cross section
+
+      END IF
       ti_ref    = L_ref/v_ref            !reference time
       Vol_ref   = L_ref**3               !reference volume
+      EF_ref    = qe*n_ref*L_ref/eps_0   !reference electric field
+      Volt_ref  = EF_ref*L_ref           !reference voltage
 
     END SUBROUTINE readReference
 
@@ -94,6 +107,7 @@ MODULE moduleInput
       USE moduleRefParam
       USE moduleErrors
       USE moduleCaseParam
+      USE moduleSolver
       USE json_module
       IMPLICIT NONE
 
@@ -101,6 +115,7 @@ MODULE moduleInput
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: object
       REAL(8):: time !simulation time in [t]
+      CHARACTER(:), ALLOCATABLE:: solverType
 
       object = 'case'
       !Time parameters
@@ -110,6 +125,21 @@ MODULE moduleInput
       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())
+
+      END SELECT
+
       !Convert simulation time to number of iterations
       tmax = INT(time/tau)
 
@@ -149,6 +179,7 @@ MODULE moduleInput
 
       CALL config%get(object // '.cpuTime', timeOutput, found)
       CALL config%get(object // '.numColl', collOutput, found)
+      CALL config%get(object // '.EMField', emOutput,   found)
 
     END SUBROUTINE readOutput
 
@@ -164,12 +195,12 @@ MODULE moduleInput
       CHARACTER(2):: iString
       CHARACTER(:), ALLOCATABLE:: object
       CHARACTER(:), ALLOCATABLE:: speciesType
-      REAL(8):: mass
+      REAL(8):: mass, charge
       LOGICAL:: found
       INTEGER:: i
 
       !Gets the number of species
-      CALL config%info('species', n_children = nSpecies)
+      CALL config%info('species', found, n_children = nSpecies)
       !Zero species means critical error
       IF (nSpecies == 0) CALL criticalError("No species found", "configRead")
 
@@ -180,21 +211,28 @@ MODULE moduleInput
         WRITE(iString, '(I2)') i
         object = 'species(' // TRIM(iString) // ')'
         CALL config%get(object // '.type', speciesType, found)
+        CALL config%get(object // '.mass', mass,        found)
+        mass = mass/m_ref
 
+        !Allocate species depending on type and assign specific parameters
         SELECT CASE(speciesType)
         CASE ("neutral")
           ALLOCATE(species(i)%obj, source=speciesNeutral())
-          !TODO: move to subroutine
+
-          CALL config%get(object // '.name',   species(i)%obj%name,        found)
+        CASE ("charged")
-          CALL config%get(object // '.mass',   mass,                       found)
+          CALL config%get(object // '.charge', charge,      found)
-          CALL config%get(object // '.weight', species(i)%obj%weight,      found)
+          ALLOCATE(species(i)%obj, source=speciesCharged(q  = charge, &
-          species(i)%obj%pt = i
+                                                         qm = charge/mass))
-          species(i)%obj%m  = mass/m_ref
 
         CASE DEFAULT
           CALL warningError("Species " // speciesType // " not supported yet")
 
         END SELECT
+        !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%m  = mass
 
       END DO
 
@@ -217,29 +255,30 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: species_i, species_j
       CHARACTER(:), ALLOCATABLE:: crossSecFile
       CHARACTER(:), ALLOCATABLE:: crossSecFilePath
+      LOGICAL:: found
       INTEGER:: nInteractions, nCollisions
       INTEGER:: i, k, ij
       INTEGER:: pt_i, pt_j
 
       CALL initInteractionMatrix(interactionMatrix)
 
-      CALL config%get('interactions.folderCollisions', pathCollisions)
+      CALL config%get('interactions.folderCollisions', pathCollisions, found)
 
-      CALL config%info('interactions.collisions', n_children = nInteractions)
+      CALL config%info('interactions.collisions', found, n_children = nInteractions)
       DO i = 1, nInteractions
         WRITE(iString, '(I2)') i
         object = 'interactions.collisions(' // TRIM(iString) // ')'
-        CALL config%get(object // '.species_i', species_i)
+        CALL config%get(object // '.species_i', species_i, found)
         pt_i = speciesName2Index(species_i)
-        CALL config%get(object // '.species_j', species_j)
+        CALL config%get(object // '.species_j', species_j, found)
         pt_j = speciesName2Index(species_j)
-        CALL config%info(object // '.crossSections', n_children = nCollisions)
+        CALL config%info(object // '.crossSections', found, n_children = nCollisions)
         ij = interactionIndex(pt_i,pt_j)
         CALL interactionMatrix(ij)%init(nCollisions)
 
         DO k = 1, nCollisions
           WRITE (kString, '(I2)') k
-          CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile)
+          CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile, found)
           crossSecFilePath = pathCollisions // crossSecFile
           CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
@@ -247,7 +286,6 @@ MODULE moduleInput
 
       END DO
 
-
     END SUBROUTINE readInteractions
 
     !Reads boundary conditions for the mesh
@@ -263,7 +301,7 @@ MODULE moduleInput
       LOGICAL:: found
       INTEGER:: i
 
-      CALL config%info('boundary', n_children = nBoundary)
+      CALL config%info('boundary', found, n_children = nBoundary)
       ALLOCATE(boundary(1:nBoundary))
       DO i = 1, nBoundary
         WRITE(istring, '(i2)') i
@@ -322,6 +360,85 @@ MODULE moduleInput
 
     END SUBROUTINE readGeometry
 
+    SUBROUTINE readEMBoundary(config)
+      USE moduleMesh
+      USE moduleOutput
+      USE moduleErrors
+      USE moduleEM
+      USE moduleRefParam
+      USE moduleSpecies
+      USE json_module
+      IMPLICIT NONE
+
+      TYPE(json_file), INTENT(inout):: config
+      CHARACTER(:), ALLOCATABLE:: object
+      LOGICAL:: found
+      CHARACTER(2):: istring
+      INTEGER:: i, e, s
+      INTEGER:: info
+      EXTERNAL:: dgetrf
+
+      CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
+      
+      IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
+      
+      DO i = 1, nBoundaryEM
+        WRITE(istring, '(i2)') i
+        object = 'boundaryEM(' // trim(istring) // ')'
+
+        CALL config%get(object // '.type', boundEM(i)%typeEM, found)
+
+        SELECT CASE(boundEM(i)%typeEM)
+        CASE ("dirichlet")
+          CALL config%get(object // '.potential', boundEM(i)%potential, found)
+          IF (.NOT. found) &
+            CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
+          boundEM(i)%potential = boundEM(i)%potential/Volt_ref
+
+          CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
+          IF (.NOT. found) &
+            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
+
+        CASE DEFAULT
+          CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
+
+        END SELECT
+
+      END DO
+
+      ALLOCATE(qSpecies(1:nSpecies))
+      DO s = 1, nSpecies
+        SELECT TYPE(sp => species(s)%obj)
+        TYPE IS (speciesCharged)
+          qSpecies(s) = sp%q
+
+        CLASS DEFAULT
+          qSpecies(s) = 0.D0
+
+        END SELECT
+
+      END DO
+
+      !TODO: Improve this
+      IF (ALLOCATED(boundEM)) THEN
+        DO e = 1, mesh%numEdges
+          IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
+            DO i = 1, nBoundaryEM
+              IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
+                CALL boundEM(i)%apply(mesh%edges(e)%obj)
+
+              END IF
+            END DO
+          END IF
+        END DO
+      END IF
+
+      !Compute the PLU factorization of K once boundary conditions have been read
+      CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
+      IF (info /= 0) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
+
+    END SUBROUTINE readEMBoundary
+
     !Reads the injection of particles from the boundaries
     SUBROUTINE readInject(config)
       USE moduleSpecies
@@ -340,10 +457,11 @@ MODULE moduleInput
       REAL(8):: v
       REAL(8), ALLOCATABLE:: T(:), normal(:)
       REAL(8):: flow
+      CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
       INTEGER:: pt
 
-      CALL config%info('inject', n_children = nInject)
+      CALL config%info('inject', found, n_children = nInject)
       ALLOCATE(inject(1:nInject))
       nPartInj = 0
       DO i = 1, nInject
@@ -359,11 +477,13 @@ MODULE moduleInput
         CALL config%get(object // '.T',               T,               found)
         CALL config%get(object // '.n',               normal,          found)
         CALL config%get(object // '.flow',            flow,            found)
+        CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, pt, physicalSurface)
+        CALL inject(i)%init(i, v, normal, T, flow, units, pt, physicalSurface)
 
       END DO
+      PRINT *, inject(:)%nParticles
 
       !Allocate array for injected particles
       IF (nPartInj > 0) THEN

src/modules/moduleMesh.f95

@@ -2,7 +2,6 @@
 MODULE moduleMesh
   USE moduleList
   USE moduleOutput
-  USE OMP_LIB
   IMPLICIT NONE
 
   !Parent of Node element
@@ -13,6 +12,7 @@ MODULE moduleMesh
     REAL(8):: v = 0.D0
     !Output values
     TYPE(outputNode), ALLOCATABLE:: output(:)
+    TYPE(emNode):: emData
     CONTAINS
       PROCEDURE(initNode_interface), DEFERRED, PASS:: init
       PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
@@ -58,15 +58,17 @@ MODULE moduleMesh
     !id for boundary condition
     INTEGER:: bt = 0
     CONTAINS
-      PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
+      PROCEDURE(initEdge_interface),     DEFERRED, PASS:: init
-      PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary
+      PROCEDURE(boundary_interface),     DEFERRED, PASS:: fBoundary
-      PROCEDURE(randPos_interface), DEFERRED, PASS::  randPos
+      PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
+      PROCEDURE(randPos_interface),      DEFERRED, PASS:: randPos
 
   END TYPE meshEdge
 
   ABSTRACT INTERFACE
     SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
       IMPORT:: meshEdge
+
       CLASS(meshEdge), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
@@ -84,10 +86,16 @@ MODULE moduleMesh
 
     END SUBROUTINE
 
-    PURE FUNCTION randPos_interface(self, rnd) RESULT(r)
+    PURE FUNCTION getNodesEdge_interface(self) RESULT(n)
+      IMPORT:: meshEdge
+      CLASS(meshEdge), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+    END FUNCTION
+
+    FUNCTION randPos_interface(self) RESULT(r)
       IMPORT:: meshEdge
       CLASS(meshEdge), INTENT(in):: self
-      REAL(8), INTENT(in):: rnd
       REAL(8):: r(1:3)
 
     END FUNCTION randPos_interface
@@ -119,6 +127,9 @@ MODULE moduleMesh
     CONTAINS
       PROCEDURE(initVol_interface),     DEFERRED, PASS:: init
       PROCEDURE(scatter_interface),     DEFERRED, PASS:: scatter
+      PROCEDURE(gatherEF_interface),    DEFERRED, PASS:: gatherEF
+      PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
+      PROCEDURE(elemF_interface),       DEFERRED, PASS:: elemF
       PROCEDURE(collision_interface),   DEFERRED, PASS:: collision
       PROCEDURE(findCell_interface),    DEFERRED, PASS:: findCell
       PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput
@@ -143,6 +154,32 @@ MODULE moduleMesh
 
     END SUBROUTINE scatter_interface
 
+    PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF)
+      
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: EF(1:3)
+
+    END FUNCTION gatherEF_interface
+
+    PURE FUNCTION getNodesVol_interface(self) RESULT(n)
+      
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+    END FUNCTION getNodesVol_interface
+
+    PURE FUNCTION elemF_interface(self, source) RESULT(localF)
+      
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8), INTENT(in):: source(1:)
+      REAL(8), ALLOCATABLE:: localF(:)
+
+    END FUNCTION elemF_interface
+
     SUBROUTINE collision_interface(self)
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(inout):: self
@@ -184,13 +221,14 @@ MODULE moduleMesh
     TYPE(meshVolCont),  ALLOCATABLE:: vols(:)
     !Global stiffness matrix
     REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
-    !Global load vector
+    !Permutation matrix for P L U factorization
-    REAL(8), ALLOCATABLE, DIMENSION(:):: F 
+    INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
 
     CONTAINS
       PROCEDURE(readMesh_interface),    PASS, DEFERRED:: readMesh
       PROCEDURE(printOutput_interface), PASS, DEFERRED:: printOutput
       PROCEDURE(printColl_interface),   PASS, DEFERRED:: printColl
+      PROCEDURE(printEM_interface),     PASS, DEFERRED:: printEM
 
   END TYPE meshGeneric
 
@@ -213,7 +251,7 @@ MODULE moduleMesh
 
     END SUBROUTINE printOutput_interface
 
-    !Prints unmber of collisions
+    !Prints number of collisions
     SUBROUTINE printColl_interface(self, t)
       IMPORT meshGeneric
 
@@ -222,6 +260,14 @@ MODULE moduleMesh
 
     END SUBROUTINE printColl_interface
 
+    !Prints EM info
+    SUBROUTINE printEM_interface(self, t)
+      IMPORT meshGeneric
+
+      CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+
+    END SUBROUTINE printEM_interface
   END INTERFACE
 
   !Generic mesh

src/modules/moduleMeshCyl.f95

@@ -25,8 +25,9 @@ MODULE moduleMeshCyl
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
     CONTAINS
-      PROCEDURE, PASS:: init    => initEdgeCyl
+      PROCEDURE, PASS:: init     => initEdgeCyl
-      PROCEDURE, PASS:: randPos => randPosCyl
+      PROCEDURE, PASS:: getNodes => getNodesCyl
+      PROCEDURE, PASS:: randPos  => randPosCyl
 
   END TYPE meshEdgeCyl
 
@@ -43,22 +44,20 @@ MODULE moduleMeshCyl
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
     !Connectivity to adjacent elements
     CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
-    !Maximum collision rate
-    !TODO: Move to other more appropiate section
-    REAL(8):: phi(1:4) = 0.D0
-    REAL(8):: Ke(1:4,1:4) = 0.D0
     REAL(8):: arNodes(1:4) = 0.D0
 
     CONTAINS
       PROCEDURE, PASS::   init        => initVolQuadCyl
-      PROCEDURE, PASS::   locKe       => localKeQuad
+      PROCEDURE, PASS::   elemK       => elemKQuad
+      PROCEDURE, PASS::   elemF       => elemFQuad
       PROCEDURE, PASS::   detJac      => detJQuad
       PROCEDURE, PASS::   invJ        => invJQuad
       PROCEDURE, PASS::   area        => areaQuad
       PROCEDURE, NOPASS:: weight      => weightQuad
       PROCEDURE, NOPASS:: inside      => insideQuad
-      PROCEDURE, PASS::   dVal        => dValQuad
       PROCEDURE, PASS::   scatter     => scatterQuad
+      PROCEDURE, PASS::   gatherEF    => gatherEFQuad
+      PROCEDURE, PASS::   getNodes    => getNodesQuad
       PROCEDURE, PASS::   phy2log     => phy2logQuad
       PROCEDURE, PASS::   findCell    => findCellCylQuad
       PROCEDURE, PASS::   resetOutput => resetOutputQuad
@@ -78,8 +77,8 @@ MODULE moduleMeshCyl
       REAL(8), INTENT(in):: r(1:3)
 
       self%n = n
-      self%r = r(1)/L_ref
       self%z = r(2)/L_ref
+      self%r = r(1)/L_ref
       !Node volume, to be determined in mesh
       self%v = 0.D0
 
@@ -128,13 +127,25 @@ MODULE moduleMeshCyl
 
     END SUBROUTINE initEdgeCyl
 
+    PURE FUNCTION getNodesCyl(self) RESULT(n)
+      IMPLICIT NONE
 
-    PURE FUNCTION randPosCyl(self, rnd) RESULT(r)
       CLASS(meshEdgeCyl), INTENT(in):: self
-      REAL(8), INTENT(in):: rnd
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1:2))
+      n = (/self%n1%n, self%n2%n /)
+
+    END FUNCTION getNodesCyl
+
+    FUNCTION randPosCyl(self) RESULT(r)
+      CLASS(meshEdgeCyl), INTENT(in):: self
+      REAL(8):: rnd
       REAL(8):: r(1:3)
       REAL(8):: p1(1:2), p2(1:2)
 
+      rnd = RAND()
+
       p1 = (/self%z(1), self%r(1) /)
       p2 = (/self%z(2), self%r(2) /)
       r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
@@ -173,15 +184,14 @@ MODULE moduleMeshCyl
       self%n3%v = self%n3%v + self%arNodes(3)
       self%n4%v = self%n4%v + self%arNodes(4)
 
-      CALL self%locKe()
-
       self%sigmaVrelMax = sigma_ref/L_ref**2
 
       CALL OMP_INIT_LOCK(self%lock)
 
     END SUBROUTINE initVolQuadCyl
 
-    FUNCTION fpsi(xi1,xi2) RESULT(psi) 
+    !TODO: MAKE ELEMENT DEPENDENT
+    PURE FUNCTION fpsi(xi1,xi2) RESULT(psi) 
       IMPLICIT NONE
 
       REAL(8),INTENT(in):: xi1,xi2
@@ -195,7 +205,8 @@ MODULE moduleMeshCyl
     END FUNCTION fpsi
 
     !Derivative element function (xi1)
-    FUNCTION dpsiXi1(xi2) RESULT(psiXi1) 
+    !TODO: MAKE ELEMENT DEPENDENT
+    PURE FUNCTION dpsiXi1(xi2) RESULT(psiXi1) 
       IMPLICIT NONE
 
       REAL(8),INTENT(in):: xi2
@@ -209,7 +220,8 @@ MODULE moduleMeshCyl
     END FUNCTION dpsiXi1
 
     !Derivative element function (xi2)
-    FUNCTION dpsiXi2(xi1) RESULT(psiXi2) 
+    !TODO: MAKE ELEMENT DEPENDENT
+    PURE FUNCTION dpsiXi2(xi1) RESULT(psiXi2) 
       IMPLICIT NONE
 
       REAL(8),INTENT(in):: xi1
@@ -223,12 +235,13 @@ MODULE moduleMeshCyl
     END FUNCTION dpsiXi2
 
     !Transforms physical coordinates to element coordinates
-    FUNCTION phy2logQuad(this,r) RESULT(xN) 
+    PURE FUNCTION phy2logQuad(self,r) RESULT(xN) 
       IMPLICIT NONE
 
-      CLASS(meshVolCylQuad):: this
+      CLASS(meshVolCylQuad), INTENT(in):: self
-      REAL(8):: r(1:2)
+      REAL(8), INTENT(in):: r(1:3)
-      REAL(8):: xN(1:2), xO(1:2), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2)
+      REAL(8):: xN(1:3)
+      REAL(8):: xO(1:3), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
@@ -238,12 +251,12 @@ MODULE moduleMeshCyl
       DO WHILE(conv>1.D-4)
         dPsi(1,:) = dpsiXi1(xO(2))
         dPsi(2,:) = dpsiXi2(xO(1))
-        j    = this%invJ(xO(1),xO(2),dPsi)
+        j    = self%invJ(xO, dPsi)
-        psi = fpsi(xO(1), xO(2))
+        psi  = fpsi(xO(1), xO(2))
-        f(1) = DOT_PRODUCT(psi,this%z)-r(1)
+        f(1) = DOT_PRODUCT(psi,self%z)-r(1)
-        f(2) = DOT_PRODUCT(psi,this%r)-r(2)
+        f(2) = DOT_PRODUCT(psi,self%r)-r(2)
-        detJ = this%detJac(xO(1),xO(2),dPsi)
+        detJ = self%detJac(xO,dPsi)
-        xN=xO - MATMUL(j, f)/detJ
+        xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ
         conv=MAXVAL(DABS(xN-xO),1)
         xO=xN
 
@@ -252,106 +265,139 @@ MODULE moduleMeshCyl
     END FUNCTION phy2logQuad
 
     !Computes element local stiffness matrix
-    SUBROUTINE localKeQuad(this) 
+    PURE FUNCTION elemKQuad(self) RESULT(ke)
       USE moduleConstParam
       IMPLICIT NONE
 
-      CLASS(meshVolCylQuad):: this
+      CLASS(meshVolCylQuad), INTENT(in):: self
-      REAL(8):: r, xi1, xi2, dPsi(1:2,1:4)
+      REAL(8):: r, xi(1:3), dPsi(1:2,1:4)
+      REAL(8):: ke(1:4,1:4)
       REAL(8):: j(1:2,1:2), detJ
       INTEGER:: l, m
 
-      this%Ke=0.D0
+      ke=0.D0
       !Start 2D Gauss Quad Integral
       DO l=1, 3
+        xi(1)     = cor(l)
+        dPsi(2,:) = dpsiXi2(xi(1))
         DO m = 1, 3
-          xi1 = cor(l)
+          xi(2)     = cor(m)
-          xi2 = cor(m)
+          dPsi(1,:) = dpsiXi1(xi(2))
-          dPsi(1,:) = dpsiXi1(xi2)
+          detJ      = self%detJac(xi,dPsi)
-          dPsi(2,:) = dpsiXi2(xi1)
+          j         = self%invJ(xi,dPsi)
-          detJ = this%detJac(xi1,xi2,dPsi)
+          r         = DOT_PRODUCT(fpsi(xi(1),xi(2)),self%r)
-          j    = this%invJ(xi1,xi2,dPsi)
+          ke        = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ
-          r    = DOT_PRODUCT(fpsi(xi1,xi2),this%r)
-          this%Ke = this%Ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ
 
         END DO
       END DO
-      this%Ke = this%Ke*2.D0*PI
+      ke = ke*2.D0*PI
 
-    END SUBROUTINE localKeQuad
+    END FUNCTION elemKQuad
 
-    !Computes element Jacobian determinant
+    !Computes the local source vector for a force f
-    FUNCTION detJQuad(this,xi1,xi2,dPsi_in) RESULT(dJ) 
+    PURE FUNCTION elemFQuad(self, source) RESULT(localF)
+      USE moduleConstParam
       IMPLICIT NONE
 
-      REAL(8),OPTIONAL:: dPsi_in(1:2,1:4)
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: source(1:)
+      REAL(8), ALLOCATABLE:: localF(:)
+      REAL(8):: r, xi(1:3), psi(1:4)
+      REAL(8):: detJ, f
+      INTEGER:: l, m
+
+      ALLOCATE(localF(1:4))
+      localF = 0.D0
+      xi = 0.D0
+      DO l=1, 3
+        DO m = 1, 3
+          xi(1) = cor(l)
+          xi(2) = cor(m)
+          detJ = self%detJac(xi)
+          psi = fpsi(xi(1),xi(2))
+          r = DOT_PRODUCT(psi,self%r)
+          f = DOT_PRODUCT(psi,source)
+          localF = localF + r*f*psi*w(l)*w(m)*detJ
+
+         END DO
+      END DO
+      localF = localF*2.D0*PI
+
+    END FUNCTION elemFQuad
+
+    !Computes element Jacobian determinant
+    PURE FUNCTION detJQuad(self,xi,dPsi_in) RESULT(dJ) 
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:2,1:4)
+      REAL(8):: dJ
       REAL(8):: dPsi(1:2,1:4)
       REAL(8):: dz(1:2), dr(1:2)
-      REAL(8):: xi1,xi2, dJ
-      CLASS(meshVolCylQuad):: this
 
       IF(PRESENT(dPsi_in)) THEN
         dPsi=dPsi_in
 
       ELSE
-        dPsi(1,:) = dpsiXi1(xi2)
+        dPsi(1,:) = dpsiXi1(xi(2))
-        dPsi(2,:) = dpsiXi2(xi1)
+        dPsi(2,:) = dpsiXi2(xi(1))
 
       END IF
 
-      dz(1) = DOT_PRODUCT(dPsi(1,:),this%z)
+      dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
-      dz(2) = DOT_PRODUCT(dPsi(2,:),this%z)
+      dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
-      dr(1) = DOT_PRODUCT(dPsi(1,:),this%r)
+      dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
-      dr(2) = DOT_PRODUCT(dPsi(2,:),this%r)
+      dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
       dJ = dz(1)*dr(2)-dz(2)*dr(1)
 
     END FUNCTION detJQuad
 
-    FUNCTION invJQuad(this,xi1,xi2,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant)
+    PURE FUNCTION invJQuad(self,xi,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant)
       IMPLICIT NONE
 
-      REAL(8),OPTIONAL:: dpsi_in(1:2,1:4)
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), INTENT(in), OPTIONAL:: dpsi_in(1:2,1:4)
       REAL(8):: dPsi(1:2,1:4)
       REAL(8):: dz(1:2), dr(1:2)
-      REAL(8):: xi1,xi2, j(1:2,1:2)
+      REAL(8):: j(1:2,1:2)
-      CLASS(meshVolCylQuad):: this
 
       IF(PRESENT(dPsi_in)) THEN
         dPsi=dPsi_in
 
       ELSE
-        dPsi(1,:) = dpsiXi1(xi2)
+        dPsi(1,:) = dpsiXi1(xi(2))
-        dPsi(2,:) = dpsiXi2(xi1)
+        dPsi(2,:) = dpsiXi2(xi(1))
 
       END IF
 
-      dz(1) = DOT_PRODUCT(dPsi(1,:),this%z)
+      dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
-      dz(2) = DOT_PRODUCT(dPsi(2,:),this%z)
+      dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
-      dr(1) = DOT_PRODUCT(dPsi(1,:),this%r)
+      dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
-      dr(2) = DOT_PRODUCT(dPsi(2,:),this%r)
+      dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
       j(1,:) = (/  dr(2), -dz(2) /)
       j(2,:) = (/ -dr(1),  dz(1) /)
 
     END FUNCTION invJQuad
 
-    SUBROUTINE areaQuad(this)!Computes element area
+    SUBROUTINE areaQuad(self)!Computes element area
       USE moduleConstParam
       IMPLICIT NONE
 
-      CLASS(meshVolCylQuad):: this
+      CLASS(meshVolCylQuad):: self
-      REAL(8):: r, xi1, xi2
+      REAL(8):: r, xi(1:3)
       REAL(8):: detJ, fpsi0(1:4)
 
-      this%volume  = 0.D0
+      self%volume  = 0.D0
-      this%arNodes = 0.D0
+      self%arNodes = 0.D0
       !2D 1 point Gauss Quad Integral
-      xi1 = 0.D0
+      xi = 0.D0
-      xi2 = 0.D0
+      detJ = self%detJac(xi)*8.D0*PI !4*2*pi
-      detJ = this%detJac(xi1,xi2)*8.D0*PI !4*2*pi
+      fpsi0 = fpsi(xi(1),xi(2))
-      fpsi0 = fpsi(xi1,xi2)
+      r = DOT_PRODUCT(fpsi0,self%r)
-      r = DOT_PRODUCT(fpsi0,this%r)
+      self%volume  =       r*detJ
-      this%volume  =       r*detJ
+      self%arNodes = fpsi0*r*detJ
-      this%arNodes = fpsi0*r*detJ
 
     END SUBROUTINE areaQuad
 
@@ -376,24 +422,6 @@ MODULE moduleMeshCyl
 
     END FUNCTION insideQuad
 
-    FUNCTION dValQuad(this,xi1,xi2) RESULT(EF)!Computes electric field in xi1,xi2
-      IMPLICIT NONE
-
-      CLASS(meshVolCylQuad):: this
-      REAL(8):: xi1, xi2, dPsi(1:2,1:4)
-      REAL(8):: j(1:2,1:2), detJ
-      REAL(8):: EF(1:3)
-
-      dPsi(1,:) = dpsiXi1(xi2)
-      dPsi(2,:) = dpsiXi2(xi1)
-      detJ = this%detJac(xi1,xi2,dPsi)
-      j    = this%invJ(xi1,xi2,dPsi)
-      EF(1) = -DOT_PRODUCT(dPsi(1,:), this%phi)*j(2,2)/detJ
-      EF(2) = -DOT_PRODUCT(dPsi(2,:), this%phi)*j(1,1)/detJ
-      EF(3) =  0.D0
-
-    END FUNCTION dValQuad
-
     SUBROUTINE scatterQuad(self, part)
       USE moduleOutput
       USE moduleSpecies
@@ -428,29 +456,65 @@ MODULE moduleMeshCyl
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
 
-
     END SUBROUTINE scatterQuad
 
+    PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)!Computes electric field inside element
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: dPsi(1:2,1:4)
+      REAL(8):: j(1:2,1:2), detJ
+      REAL(8):: phi(1:4)
+      REAL(8):: EF(1:3)
+
+      phi = (/self%n1%emData%phi, &
+              self%n2%emData%phi, &
+              self%n3%emData%phi, &
+              self%n4%emData%phi /)
+
+      dPsi(1,:) = dpsiXi1(xi(2))
+      dPsi(2,:) = dpsiXi2(xi(1))
+      detJ = self%detJac(xi,dPsi)
+      j    = self%invJ(xi,dPsi)
+      EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ
+      EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ
+      EF(3) =  0.D0
+
+    END FUNCTION gatherEFQuad
+
+    PURE FUNCTION getNodesQuad(self) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1:4))
+      n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
+
+    END FUNCTION getNodesQuad
+
     RECURSIVE SUBROUTINE findCellCylQuad(self, part, oldCell)
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVolCylQuad), INTENT(inout):: self
       CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
       CLASS(particle), INTENT(inout), TARGET:: part
-      REAL(8):: xLog(1:2)
+      REAL(8):: xLog(1:3)
       REAL(8):: xLogArray(1:4)
       CLASS(*), POINTER:: nextElement
       INTEGER:: nextInt
 
-      xLog = self%phy2log(part%r(1:2))
+      xLog = self%phy2log(part%r)
+      !Checks if particle is inside of current cell
       IF (self%inside(xLog(1), xLog(2))) THEN
-        !Checks if particle is inside of current cell
+        ! IF (PRESENT(oldCell)) THEN
-        IF (PRESENT(oldCell)) THEN
+        !   !If oldCell, recalculate particle weight, as particle has entered a new cell.
-          !If oldCell, recalculate particle weight, as particle has entered a new cell.
+        !   part%weight = part%weight*oldCell%volume/self%volume
-          part%weight = part%weight*oldCell%volume/self%volume
+        !
-
+        ! END IF
-        END IF
         part%e_p  = self%n
         part%xLog = xLog
         part%n_in = .TRUE.

src/modules/moduleMeshCylBoundary.f95

@@ -74,8 +74,6 @@ MODULE moduleMeshCylBoundary
       CLASS(meshEdgeCylAxis), INTENT(inout):: self
       CLASS(particle), INTENT(inout):: part
 
-      !Do to the 2D Cyl pusher, this function should never be called.
-
     END SUBROUTINE symmetryAxis
 
 

src/modules/moduleMeshCylRead.f95

@@ -8,6 +8,7 @@ MODULE moduleMeshCylRead
       PROCEDURE, PASS:: readMesh => readMeshCyl
       PROCEDURE, PASS:: printOutput => printOutputCyl
       PROCEDURE, PASS:: printColl   => printCollisionsCyl
+      PROCEDURE, PASS:: printEM     => printEMCyl
 
   END TYPE
 
@@ -41,8 +42,10 @@ MODULE moduleMeshCylRead
       READ(10, *) self%numNodes
       !Allocate required matrices and vectors
       ALLOCATE(self%nodes(1:self%numNodes))
-      ! ALLOCATE(self%K(1:numNodes,1:numNodes))
+      ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
-      ! ALLOCATE(self%F(1:numNodes))
+      ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
+      self%K = 0.D0
+      self%IPIV = 0
       !Read nodes cartesian coordinates (x=z, y=r, z=null)
       DO e=1, self%numNodes
         READ(10, *) n, z, r
@@ -123,24 +126,10 @@ MODULE moduleMeshCylRead
 
         END DO
 
-      END DO
+        !Constructs the global K matrix
+        CALL constructGlobalK(self%K, self%vols(e)%obj)
 
-      ! !Compute global stiffness matrix
+      END DO
-      ! GlobalK=0.D0
-      ! DO e=1, numElem
-        ! DO i=1, 4
-          ! DO j=1, 4
-            ! GlobalK(elems(e)%p(i),elems(e)%p(j)) = GlobalK(elems(e)%p(i),elems(e)%p(j)) + elems(e)%Ke(i,j)
-          ! END DO
-        ! END DO
-      ! END DO
-      ! ! Apply Dirichlet boundary conditions to GlobalK
-      ! DO n=1, numNodes
-        ! IF (nodes(n)%bound == 1) THEN
-          ! GlobalK(n,:)=0.D0
-          ! GlobalK(n,n)=1.D0
-        ! END IF
-      ! END DO
 
     END SUBROUTINE
 
@@ -298,6 +287,41 @@ MODULE moduleMeshCylRead
 
     END SUBROUTINE connectedQuadEdge
 
+    SUBROUTINE constructGlobalK(K, elem)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(inout):: K(1:,1:)
+      CLASS(meshVol), INTENT(in):: elem
+      REAL(8), ALLOCATABLE:: localK(:,:)
+      INTEGER:: nNodes, i, j
+      INTEGER, ALLOCATABLE:: n(:)
+
+      SELECT TYPE(elem)
+      TYPE IS(meshVolCylQuad)
+        nNodes = 4
+        ALLOCATE(localK(1:nNodes,1:nNodes))
+        localK = elem%elemK()
+        ALLOCATE(n(1:nNodes))
+        n = (/ elem%n1%n, elem%n2%n, &
+               elem%n3%n, elem%n4%n /)
+
+      CLASS DEFAULT
+        nNodes = 0
+        ALLOCATE(localK(1:1, 1:1))
+        localK = 0.D0
+        ALLOCATE(n(1:1))
+        n = 0
+
+      END SELECT
+      
+      DO i = 1, nNodes
+        DO j = 1, nNodes
+          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
+        END DO
+      END DO
+
+    END SUBROUTINE constructGlobalK
+
     SUBROUTINE printOutputCyl(self, t)
       USE moduleRefParam
       USE moduleSpecies
@@ -425,4 +449,60 @@ MODULE moduleMeshCylRead
 
     END SUBROUTINE printCollisionsCyl
 
+    SUBROUTINE printEMCyl(self, t)
+      USE moduleRefParam
+      USE moduleCaseParam
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshCylGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, e
+      REAL(8):: time
+      CHARACTER(:), ALLOCATABLE:: fileName
+      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
+
+      IF (emOutput) THEN
+        time = DBLE(t)*tau*ti_ref
+        WRITE(tstring, '(I6.6)') t
+
+        fileName='OUTPUT_' // tstring// '_EMField.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (20, file = path // folder // '/' //  fileName)
+        WRITE(20, "(A)") '$MeshFormat'
+        WRITE(20, "(A)") '2.2 0 8'
+        WRITE(20, "(A)") '$EndMeshFormat'
+        WRITE(20, "(A)") '$NodeData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Potential (V)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 1
+        WRITE(20, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
+        END DO
+        WRITE(20, "(A)") '$EndNodeData'
+        
+        WRITE(20, "(A)") '$ElementData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Electric Field (V/m)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 3
+        WRITE(20, *) self%numVols
+        DO e=1, self%numVols
+          WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF((/0.D0, 0.D0, 0.D0/))*EF_ref
+        END DO
+        WRITE(20, "(A)") '$EndElementData'
+        CLOSE(20)
+
+      END IF
+
+    END SUBROUTINE printEMCyl
+  
 END MODULE moduleMeshCylRead

src/modules/moduleOutput.f95

@@ -7,6 +7,13 @@ MODULE moduleOutput
 
   END TYPE
 
+  !Type for EM data in node
+  TYPE emNode
+    CHARACTER(:), ALLOCATABLE:: type
+    REAL(8):: phi
+
+  END TYPE emNode
+
   !Output in dimensional units to print
   TYPE outputFormat
     REAL(8):: density, velocity(1:3), pressure, temperature
@@ -18,6 +25,7 @@ MODULE moduleOutput
   INTEGER:: triggerOutput, counterOutput
   LOGICAL:: timeOutput = .FALSE.
   LOGICAL:: collOutput = .FALSE.
+  LOGICAL:: emOutput   = .FALSE.
 
   CONTAINS
     FUNCTION outerProduct(a,b) RESULT(s)
@@ -42,6 +50,7 @@ MODULE moduleOutput
     END FUNCTION tensorTrace
 
     SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn)
+      USE moduleConstParam
       USE moduleRefParam
       USE moduleSpecies
       IMPLICIT NONE

src/modules/moduleRefParam.f95

@@ -1,10 +1,9 @@
 !Reference parameters
 MODULE moduleRefParam
-  USE moduleConstParam
   !Parameters that define the problem (inputs)
-  REAL(8):: n_ref, m_ref, T_ref, r_ref, sigma_ref!, q_ref=1.6022D-19
+  REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref
   !Reference parameters for non-dimensional problem
-  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref
+  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref
 
 END MODULE moduleRefParam
 

src/modules/moduleSolver.f95

@@ -1,5 +1,44 @@
 MODULE moduleSolver
 
+  TYPE, PUBLIC, ABSTRACT:: solverGeneric
+    CONTAINS
+      PROCEDURE(push_interafece), DEFERRED, NOPASS:: push
+      PROCEDURE(solveEM_interface), DEFERRED, NOPASS:: solveEMField
+
+  END TYPE solverGeneric
+
+  ABSTRACT INTERFACE
+    !Push a particle
+    PURE SUBROUTINE push_interafece(part)
+      USE moduleSpecies
+
+      TYPE(particle), INTENT(inout):: part
+
+    END SUBROUTINE push_interafece
+
+    !Solve the electromagnetic field
+    SUBROUTINE solveEM_interface()
+    
+    END SUBROUTINE solveEM_interface
+
+  END INTERFACE
+
+  TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylNeutral
+    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
     SUBROUTINE doPushes()
       USE moduleSpecies
@@ -10,7 +49,9 @@ MODULE moduleSolver
 
       !$OMP DO
       DO n=1, nPartOld
-        CALL push(partOld(n))
+        !Push particle
+        CALL solver%push(partOld(n))
+        !Find cell in wich particle reside
         CALL mesh%vols(partOld(n)%e_p)%obj%findCell(partOld(n))
 
       END DO
@@ -18,16 +59,15 @@ MODULE moduleSolver
 
     END SUBROUTINE doPushes
 
-    !Push one particle. Boris pusher for 2D Cyl
+    !Push one particle. Boris pusher for 2D Cyl Neutral particle
-    !TODO: make general for any pusher
+    PURE SUBROUTINE pushCylNeutral(part)
-    PURE SUBROUTINE push(part)
       USE moduleSpecies
-      USE moduleMesh
       IMPLICIT NONE
+
       TYPE(particle), INTENT(inout):: part
       REAL(8):: v_p_oh_star(2:3)
       TYPE(particle):: part_temp
-      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha, alpha
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
 
       part_temp = part
       !z
@@ -40,8 +80,6 @@ MODULE moduleSolver
       y_new =             v_p_oh_star(3)*tau
       r = DSQRT(x_new**2+y_new**2)
       part_temp%r(2) = r
-      alpha = 0.D0!ATAN2(y_new,x_new) !0 in 2D problem
-      part_temp%r(3) = part%r(3) + alpha
       IF (r > 0.D0) THEN
         sin_alpha = y_new/r
         cos_alpha = x_new/r
@@ -57,7 +95,47 @@ MODULE moduleSolver
       !Copy temporal particle to particle
       part=part_temp
 
-    END SUBROUTINE push
+    END SUBROUTINE pushCylNeutral
+
+    !Push one particle. Boris pusher for 2D Cyl Charged particle
+    PURE SUBROUTINE pushCylCharged(part)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      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
+
+      part_temp = part
+      !Get electric field at particle position
+      qmEFt = part_temp%qm*gatherElecField(part_temp)*tau
+      !z
+      part_temp%v(1) = part%v(1) + qmEFt(1)
+      part_temp%r(1) = part%r(1) + part_temp%v(1)*tau
+      !r,theta
+      v_p_oh_star(2) = part%v(2) + qmEFt(2)
+      x_new = part%r(2) + v_p_oh_star(2)*tau
+      v_p_oh_star(3) = part%v(3) + qmEFt(3)
+      y_new =             v_p_oh_star(3)*tau
+      r = DSQRT(x_new**2+y_new**2)
+      part_temp%r(2) = r
+      IF (r > 0.D0) THEN
+        sin_alpha = y_new/r
+        cos_alpha = x_new/r
+      ELSE
+        sin_alpha = 0.D0
+        cos_alpha = 1.D0
+      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%n_in = .FALSE. !Assume particle is outside until cell is found
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE pushCylCharged
 
     !Do the collisions in all the cells
     SUBROUTINE doCollisions()
@@ -153,6 +231,64 @@ MODULE moduleSolver
 
     END SUBROUTINE doScatter
 
+    SUBROUTINE doEMField()
+      IMPLICIT NONE
+
+      CALL solver%solveEMField()
+
+    END SUBROUTINE doEMField
+
+    !Solving the Poisson equation for electrostatic potential
+    SUBROUTINE solveElecField()
+      USE moduleEM
+      USE moduleMesh
+      USE moduleErrors
+      IMPLICIT NONE
+
+      INTEGER, SAVE:: INFO
+      INTEGER:: n
+      REAL(8), ALLOCATABLE, SAVE:: tempF(:)
+      EXTERNAL:: dgetrs
+
+      !$OMP SINGLE
+      ALLOCATE(tempF(1:mesh%numNodes))
+      !$OMP END SINGLE
+
+      CALL assembleSourceVector(tempF)
+
+      !$OMP SINGLE
+      CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
+                              mesh%IPIV, tempF,  mesh%numNodes, info)
+      !$OMP END SINGLE
+
+      IF (info == 0) THEN
+        !Suscessful resolution of Poission equation
+        !$OMP DO
+        DO n = 1, mesh%numNodes
+          mesh%nodes(n)%obj%emData%phi = tempF(n)
+
+        END DO
+        !$OMP END DO
+
+      ELSE
+        !$OMP SINGLE
+        CALL criticalError('Poisson equation failed', 'solveElecField')
+        !$OMP END SINGLE
+
+      END IF
+
+      !$OMP SINGLE
+      DEALLOCATE(tempF)
+      !$OMP END SINGLE
+
+    END SUBROUTINE solveElecField
+
+    !Empty module that does no computation of EM field for neutral cases
+    SUBROUTINE noEMField()
+      IMPLICIT NONE
+
+    END SUBROUTINE noEMField
+
     SUBROUTINE doOutput(t)
       USE moduleMesh
       USE moduleOutput
@@ -172,8 +308,9 @@ MODULE moduleSolver
         CALL mesh%printOutput(t)
         CALL printTime(t, t == 0)
         CALL mesh%printColl(t)
-        WRITE(*, "(5X,A21,I10,A1,I8)") "t/tmax: ",    t,       "/", tmax
+        CALL mesh%printEM(t)
-        WRITE(*, "(5X,A21,I10)")       "Particles: ", nPartOld
+        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tmax: ",    t,       "/", tmax
+        WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
         IF (t == 0) THEN
           WRITE(*, "(5X,A21,F8.1,A2)")   "     init time: ",    1.D3*tStep, "ms"
 
@@ -183,7 +320,7 @@ MODULE moduleSolver
         END IF
 
         IF (nPartOld > 0) THEN
-          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*(tPush+tReset+tColl+tWeight)/DBLE(nPartOld), "ns"
+          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns"
 
         END IF
         WRITE(*,*)

src/modules/moduleSpecies.f95

@@ -31,9 +31,10 @@ MODULE moduleSpecies
     REAL(8):: v(1:3) !Velocity
     INTEGER:: pt !Particle species id
     INTEGER:: e_p !Index of element in which the particle is located
-    REAL(8):: xLog(1:2) !Logical coordinates of particle in element e_p.
+    REAL(8):: xLog(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
+    REAL(8):: weight=0.D0 !weight of particle
+    REAL(8):: qm = 0.D0 !charge over mass
 
   END TYPE particle