Merge branch 'IEPC2025' into 'development'

Merge IEPC2025

See merge request JorgeGonz/fpakc!54

data/collisions/IO_e-Kr.dat

@@ -0,0 +1,52 @@
+# D108525 "refs": {"B56": {"note": "CLM-R294 (1989)"}}
+# Relative energy (eV) cross section (m^2)
+1.40E+01 0
+1.62E+01 7.249E-21
+1.88E+01 1.199E-20
+2.18E+01 1.644E-20
+2.53E+01 2.1E-20
+2.94E+01 2.542E-20
+3.41E+01 2.937E-20
+3.95E+01 3.26E-20
+4.58E+01 3.499E-20
+5.32E+01 3.653E-20
+6.17E+01 3.726E-20
+7.15E+01 3.728E-20
+8.29E+01 3.671E-20
+9.62E+01 3.566E-20
+1.12E+02 3.426E-20
+1.29E+02 3.259E-20
+1.50E+02 3.075E-20
+1.74E+02 2.881E-20
+2.02E+02 2.682E-20
+2.34E+02 2.484E-20
+2.72E+02 2.289E-20
+3.15E+02 2.101E-20
+3.65E+02 1.922E-20
+4.24E+02 1.751E-20
+4.91E+02 1.592E-20
+5.70E+02 1.443E-20
+6.61E+02 1.305E-20
+7.67E+02 1.177E-20
+8.89E+02 1.06E-20
+1.03E+03 9.526E-21
+1.20E+03 8.547E-21
+1.39E+03 7.658E-21
+1.61E+03 6.851E-21
+1.87E+03 6.121E-21
+2.16E+03 5.462E-21
+2.51E+03 4.868E-21
+2.91E+03 4.334E-21
+3.38E+03 3.855E-21
+3.92E+03 3.426E-21
+4.54E+03 3.041E-21
+5.27E+03 2.698E-21
+6.11E+03 2.391E-21
+7.09E+03 2.118E-21
+8.22E+03 1.875E-21
+9.53E+03 1.658E-21
+1.11E+04 1.466E-21
+1.28E+04 1.295E-21
+1.49E+04 1.143E-21
+1.72E+04 1.009E-21
+2.00E+04 8.898E-22

data/collisions/IO_e-Xe.dat

@@ -0,0 +1,52 @@
+# EL cross sections extracted from PROGRAM MAGBOLTZ, VERSION 7.1 JUNE 2004 www.lxcat.net/Biagi-v7.1
+# Relative energy (eV) cross section (m^2)
+1.21E+01 0
+1.41E+01 3.923E-21
+1.64E+01 1.194E-20
+1.91E+01 2.1E-20
+2.22E+01 2.946E-20
+2.58E+01 3.65E-20
+3.00E+01 4.185E-20
+3.49E+01 4.552E-20
+4.06E+01 4.766E-20
+4.72E+01 4.85E-20
+5.49E+01 4.828E-20
+6.39E+01 5.031E-20
+7.43E+01 5.1E-20
+8.64E+01 5.1E-20
+1.01E+02 5.032E-20
+1.17E+02 4.906E-20
+1.36E+02 4.732E-20
+1.58E+02 4.521E-20
+1.84E+02 4.283E-20
+2.14E+02 4.029E-20
+2.49E+02 3.764E-20
+2.90E+02 3.497E-20
+3.37E+02 3.233E-20
+3.92E+02 2.975E-20
+4.56E+02 2.726E-20
+5.31E+02 2.489E-20
+6.17E+02 2.266E-20
+7.18E+02 2.056E-20
+8.35E+02 1.861E-20
+9.72E+02 1.68E-20
+1.13E+03 1.514E-20
+1.32E+03 1.361E-20
+1.53E+03 1.221E-20
+1.78E+03 1.094E-20
+2.07E+03 9.781E-21
+2.41E+03 8.735E-21
+2.80E+03 7.789E-21
+3.26E+03 6.938E-21
+3.79E+03 6.171E-21
+4.41E+03 5.484E-21
+5.13E+03 4.868E-21
+5.97E+03 4.316E-21
+6.94E+03 3.824E-21
+8.07E+03 3.385E-21
+9.39E+03 2.994E-21
+1.09E+04 2.646E-21
+1.27E+04 2.336E-21
+1.48E+04 2.062E-21
+1.72E+04 1.818E-21
+2.00E+04 1.602E-21

src/modules/common/moduleRandom.f90

@@ -47,24 +47,42 @@ MODULE moduleRandom
 
     END FUNCTION randomIntAB
 
-    !Returns a random number in a Maxwellian distribution of mean 0 and width 1
+    !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method
     FUNCTION randomMaxwellian() RESULT(rnd)
       USE moduleConstParam, ONLY: PI
       IMPLICIT NONE
 
       REAL(8):: rnd
-      REAL(8):: x, y
+      REAL(8):: v1, v2, Rsquare
+
+      Rsquare = 1.D0
+      do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0)
+        v1 = 2.D0 * random() - 1.D0
+        v2 = 2.D0 * random() - 1.D0
+        Rsquare = v1**2 + v2**2
+
+      end do
+
+      rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare)
+
+    END FUNCTION randomMaxwellian
+
+    !Returns a random number in a Maxwellian distribution of mean 0 and width 1
+    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+      REAL(8):: x
 
       rnd = 0.D0
       x = 0.D0
       DO WHILE (x == 0.D0)
         CALL RANDOM_NUMBER(x)
       END DO
-      CALL RANDOM_NUMBER(y)
 
-      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
+      rnd = DSQRT(-DLOG(x))
 
-    END FUNCTION randomMaxwellian
+    END FUNCTION randomHalfMaxwellian
 
     !Returns a random number weighted with the cumWeight array
     FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd)

src/modules/init/moduleInput.f90

@@ -259,6 +259,10 @@ MODULE moduleInput
           !Read BC
           CALL readEMBoundary(config)
 
+        CASE("ElectrostaticBoltzmann")
+          !Read BC
+          CALL readEMBoundary(config)
+
         CASE("ConstantB")
           !Read BC
           CALL readEMBoundary(config)
@@ -1327,6 +1331,7 @@ MODULE moduleInput
       USE moduleCaseParam, ONLY: tauMin
       USE moduleMesh, ONLY: mesh
       USE moduleSpecies, ONLY: nSpecies
+      USE moduleRefParam, ONLY: ti_ref
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
@@ -1340,8 +1345,11 @@ MODULE moduleInput
         CALL config%get('average.startTime', tStart, found)
 
         IF (found) THEN
-          tAverageStart = INT(tStart / tauMin)
+          tAverageStart = INT(tStart / ti_ref / tauMin)
 
+        ELSE
+          tAverageStart = 0
+          
         END IF
 
         ALLOCATE(averageScheme(1:mesh%numNodes))

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -495,34 +495,36 @@ MODULE moduleMesh2DCart
 
     END FUNCTION insideQuad
 
-    !Transform physical coordinates to element coordinates
+    !Transform physical coordinates to element coordinates with a Taylor series
     PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
       IMPLICIT NONE
 
       CLASS(meshCell2DCartQuad), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
+      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
       REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
       REAL(8):: pDer(1:3, 1:3)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
-      conv = 1.D0
+      conv  = 1.D0
-      XiO  = 0.D0
+      Xi0   = 0.D0
+      Xi(3) = 0.D0
 
-      f(3) = 0.D0
       DO WHILE(conv > 1.D-4)
-        dPsi = self%dPsi(XiO, 4)
+        fPsi           = self%fPsi(Xi0, 4)
-        pDer = self%partialDer(4, dPsi)
+        x0(1)          = dot_product(fPsi, self%x)
-        detJ = self%detJac(pDer)
+        x0(2)          = dot_product(fPsi, self%y)
-        invJ = self%invJac(pDer)
+        deltaR         = r(1:2) - x0
-        fPsi = self%fPsi(XiO, 4)
+        dPsi           = self%dPsi(Xi0, 4)
-        f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), &
+        pDer           = self%partialDer(4, dPsi)
-                    DOT_PRODUCT(fPsi,self%y) /) - r(1:2)
+        detJ           = self%detJac(pDer)
-        Xi = XiO - MATMUL(invJ, f)/detJ
+        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-        conv = MAXVAL(DABS(Xi-XiO),1)
+        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
-        XiO = Xi
+        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
+        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
+        Xi0(1:2)       = Xi(1:2)
 
       END DO
 
@@ -680,8 +682,8 @@ MODULE moduleMesh2DCart
 
       dPsi = 0.D0
 
-      dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
+      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
 
     END FUNCTION dPsiTria
 
@@ -826,19 +828,19 @@ MODULE moduleMesh2DCart
       CLASS(meshCell2DCartTria), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: deltaR(1:3)
+      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
-      REAL(8):: dPsi(1:3, 1:3)
+      REAL(8):: dPsi(1:3,1:4)
       REAL(8):: pDer(1:3, 1:3)
-      REAL(8):: invJ(1:3,  1:3), detJ
 
       !Direct method to convert coordinates
-      Xi     = 0.D0
+      Xi(3)          = 0.D0
-      deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /)
+      deltaR         = (/ r(1) - self%x(1), r(2) - self%y(1) /)
-      dPsi   = self%dPsi(Xi, 3)
+      dPsi           = self%dPsi(Xi, 3)
-      pDer   = self%partialDer(3, dPsi)
+      pDer           = self%partialDer(3, dPsi)
-      invJ   = self%invJac(pDer)
+      detJ           = self%detJac(pDer)
-      detJ   = self%detJac(pDer)
+      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-      Xi     = MATMUL(invJ,deltaR)/detJ
+      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ
 
     END FUNCTION phy2logTria
 
@@ -913,8 +915,8 @@ MODULE moduleMesh2DCart
 
       invJ = 0.D0
 
-      invJ(1, 1:2) = (/  pDer(2,2), -pDer(1,2) /)
+      invJ(1, 1:2) = (/  pDer(2,2), -pDer(2,1) /)
-      invJ(2, 1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
       invJ(3, 3)   = 1.D0
 
     END FUNCTION invJ2DCart

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -510,34 +510,36 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION insideQuad
 
-    !Transform physical coordinates to element coordinates
+    !Transform physical coordinates to element coordinates with a Taylor series
     PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
       IMPLICIT NONE
 
       CLASS(meshCell2DCylQuad), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
+      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
       REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
       REAL(8):: pDer(1:3, 1:3)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
-      conv = 1.D0
+      conv  = 1.D0
-      XiO  = 0.D0
+      Xi0   = 0.D0
+      Xi(3) = 0.D0
 
-      f(3) = 0.D0
       DO WHILE(conv > 1.D-4)
-        dPsi = self%dPsi(XiO, 4)
+        fPsi           = self%fPsi(Xi0, 4)
-        pDer = self%partialDer(4, dPsi)
+        x0(1)          = dot_product(fPsi, self%z)
-        detJ = self%detJac(pDer)
+        x0(2)          = dot_product(fPsi, self%r)
-        invJ = self%invJac(pDer)
+        deltaR         = r(1:2) - x0
-        fPsi = self%fPsi(XiO, 4)
+        dPsi           = self%dPsi(Xi0, 4)
-        f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), &
+        pDer           = self%partialDer(4, dPsi)
-                    DOT_PRODUCT(fPsi,self%r) /) - r(1:2)
+        detJ           = self%detJac(pDer)
-        Xi = XiO - MATMUL(invJ, f)/detJ
+        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-        conv = MAXVAL(DABS(Xi-XiO),1)
+        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
-        XiO = Xi
+        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
+        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
+        Xi0(1:2)       = Xi(1:2)
 
       END DO
 
@@ -705,8 +707,8 @@ MODULE moduleMesh2DCyl
 
       dPsi = 0.D0
 
-      dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
+      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
 
     END FUNCTION dPsiTria
 
@@ -858,19 +860,19 @@ MODULE moduleMesh2DCyl
       CLASS(meshCell2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: deltaR(1:3)
+      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
-      REAL(8):: dPsi(1:3, 1:3)
+      REAL(8):: dPsi(1:3,1:4)
       REAL(8):: pDer(1:3, 1:3)
-      REAL(8):: invJ(1:3,  1:3), detJ
 
       !Direct method to convert coordinates
-      Xi     = 0.D0
+      Xi(3)          = 0.D0
-      deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /)
+      deltaR         = (/ r(1) - self%z(1), r(2) - self%r(1) /)
-      dPsi   = self%dPsi(Xi, 3)
+      dPsi           = self%dPsi(Xi, 3)
-      pDer   = self%partialDer(3, dPsi)
+      pDer           = self%partialDer(3, dPsi)
-      invJ   = self%invJac(pDer)
+      detJ           = self%detJac(pDer)
-      detJ   = self%detJac(pDer)
+      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-      Xi     = MATMUL(invJ,deltaR)/detJ
+      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ
 
     END FUNCTION phy2logTria
 
@@ -948,8 +950,8 @@ MODULE moduleMesh2DCyl
 
       invJ = 0.D0
 
-      invJ(1, 1:2) = (/  pDer(2,2), -pDer(1,2) /)
+      invJ(1, 1:2) = (/  pDer(2,2), -pDer(2,1) /)
-      invJ(2, 1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
       invJ(3, 3)   = 1.D0
 
     END FUNCTION invJ2DCyl

src/modules/moduleInject.f90

@@ -199,17 +199,21 @@ MODULE moduleInject
 
         END DO
 
+        self%nParticles = SUM(self%particlesPerEdge)
+
       ELSE
         ! No particles assigned per edge, use the species weight
         self%weightPerEdge = self%species%weight
         DO et = 1, self%nEdges
-          self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight)
+          self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight))
-
         END DO
 
-      END IF
+        self%nParticles = SUM(self%particlesPerEdge)
 
-      self%nParticles = SUM(self%particlesPerEdge)
+        !Rescale weight to match flux
+        self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles))
+
+      END IF
 
       !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
@@ -285,7 +289,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = self%vTh*randomMaxwellian()
+      v = sqrt(2.0)*self%vTh*randomMaxwellian()
 
     END FUNCTION randomVelMaxwellian
 
@@ -298,10 +302,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      DO WHILE (v <= 0.D0)
+      v = sqrt(2.0)*self%vTh*randomHalfMaxwellian()
-        v = self%vTh*randomMaxwellian()
-
-      END DO
 
     END FUNCTION randomVelHalfMaxwellian
 
@@ -376,7 +377,13 @@ MODULE moduleInject
           !Assign particle type
           partInj(n)%species => self%species
 
-          direction = self%n
+          if (all(self%n == 0.D0)) then
+            direction =  randomEdge%normal
+
+          else
+            direction = self%n
+
+          end if
 
           partInj(n)%v = 0.D0
 
@@ -384,11 +391,12 @@ MODULE moduleInject
                                                    self%v(2)%obj%randomVel(), &
                                                    self%v(3)%obj%randomVel() /)
 
-          !If velocity is not in the right direction, invert it
+          !If injecting a no-drift distribution and velocity is negative, reflect
-          IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
+          if ((self%vMod == 0.D0) .and. &
+              (dot_product(direction, partInj(n)%v) < 0.D0)) then
             partInj(n)%v = - partInj(n)%v
 
-          END IF
+          end if
 
           !Obtain natural coordinates of particle in cell
           partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)

src/modules/output/moduleOutput.f90

@@ -22,7 +22,6 @@ MODULE moduleOutput
 
   !Type for EM data in node
   TYPE emNode
-    CHARACTER(:), ALLOCATABLE:: type
     REAL(8):: phi
     REAL(8):: B(1:3)
 

src/modules/solver/electromagnetic/moduleEM.f90

@@ -11,7 +11,6 @@ MODULE moduleEM
 
     CONTAINS
       PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
-      !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better.
 
   END TYPE boundaryEMGeneric
 
@@ -168,8 +167,8 @@ MODULE moduleEM
       INTEGER:: n, ni
 
       DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi = self%potential
+        self%nodes(n)%obj%emData%phi  = self%potential
-        vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
+        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
 
       END DO
 
@@ -189,15 +188,15 @@ MODULE moduleEM
       timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
 
       DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi = self%potential * timeFactor
+        self%nodes(n)%obj%emData%phi  = self%potential * timeFactor
-        vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
+        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
 
       END DO
 
     END SUBROUTINE applyDirichletTime
 
     !Assemble the source vector based on the charge density to solve Poisson's equation
-    SUBROUTINE assembleSourceVector(vectorF)
+    SUBROUTINE assembleSourceVector(vectorF, n_e)
       USE moduleMesh
       USE moduleRefParam
       IMPLICIT NONE
@@ -206,6 +205,7 @@ MODULE moduleEM
       REAL(8), ALLOCATABLE:: localF(:)
       INTEGER, ALLOCATABLE:: nodes(:)
       REAL(8), ALLOCATABLE:: rho(:)
+      REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes)
       INTEGER:: nNodes
       INTEGER:: e, i, ni, b
       CLASS(meshNode), POINTER:: node
@@ -214,7 +214,6 @@ MODULE moduleEM
       vectorF = 0.D0
       !$OMP END SINGLE
 
-      ! Calculate charge density in each node
       !$OMP DO REDUCTION(+:vectorF)
       DO e = 1, mesh%numCells
         nNodes = mesh%cells(e)%obj%nNodes
@@ -226,6 +225,10 @@ MODULE moduleEM
           ni = nodes(i)
           node => mesh%nodes(ni)%obj
           rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))
+          IF (PRESENT(n_e)) THEN
+            rho(i) = rho(i) - n_e(i)
+
+          END IF
 
         END DO
 
@@ -247,10 +250,10 @@ MODULE moduleEM
 
       !Apply boundary conditions
       !$OMP SINGLE
-      DO b = 1, nBoundaryEM
+      do b = 1, nBoundaryEM
-        CALL boundaryEM(b)%obj%apply(vectorF)
+        call boundaryEM(b)%obj%apply(vectorF)
 
-      END DO
+      end do
       !$OMP END SINGLE
 
     END SUBROUTINE assembleSourceVector
@@ -299,4 +302,86 @@ MODULE moduleEM
 
     END SUBROUTINE solveElecField
 
+    FUNCTION BoltzmannElectron(phi, n) RESULT(n_e)
+      USE moduleRefParam
+      USE moduleConstParam
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: n
+      REAL(8), INTENT(in):: phi(1:n)
+      REAL(8):: n_e(1:n)
+      REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0
+      INTEGER:: i
+
+      n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e))
+
+      RETURN
+
+    END FUNCTION BoltzmannElectron
+
+    SUBROUTINE solveElecFieldBoltzmann()
+      USE moduleMesh
+      USE moduleErrors
+      IMPLICIT NONE
+
+      INTEGER, SAVE:: INFO
+      INTEGER:: n
+      REAL(8), ALLOCATABLE, SAVE:: tempF(:)
+      REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:)
+      INTEGER:: k
+      EXTERNAL:: dgetrs
+
+      !$OMP SINGLE
+      ALLOCATE(tempF(1:mesh%numNodes))
+      ALLOCATE(n_e(1:mesh%numNodes))
+      ALLOCATE(phi_old(1:mesh%numNodes))
+      ALLOCATE(phi(1:mesh%numNodes))
+      !$OMP END SINGLE
+
+      !$OMP DO
+      DO n = 1, mesh%numNodes
+        phi_old(n) = mesh%nodes(n)%obj%emData%phi
+
+      END DO
+      !$OMP END DO
+
+      !$OMP SINGLE
+      DO k = 1, 100
+        n_e = BoltzmannElectron(phi_old, mesh%numNodes)
+        CALL assembleSourceVector(tempF, n_e)
+
+        CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
+                                mesh%IPIV, tempF,  mesh%numNodes, info)
+        phi = tempF
+
+        PRINT *, MAXVAL(n_e), MINVAL(n_e)
+        PRINT *, MAXVAL(phi), MINVAL(phi)
+        PRINT*, k, "diff = ", MAXVAL(ABS(phi - phi_old))
+        phi_old = phi
+
+      END DO
+      !$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 = phi_old(n)
+
+        END DO
+        !$OMP END DO
+
+      ELSE
+        !$OMP SINGLE
+        CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann')
+        !$OMP END SINGLE
+
+      END IF
+
+      !$OMP SINGLE
+      DEALLOCATE(tempF, n_e, phi_old, phi)
+      !$OMP END SINGLE
+
+    END SUBROUTINE solveElecFieldBoltzmann
+
 END MODULE moduleEM

src/modules/solver/moduleSolver.f90

@@ -138,6 +138,9 @@ MODULE moduleSolver
       CASE('Electrostatic','ConstantB')
         self%solveEM => solveElecField
 
+      CASE('ElectrostaticBoltzmann')
+        self%solveEM => solveElecFieldBoltzmann
+
       END SELECT
 
     END SUBROUTINE initEM