Implementation of Boltzmann electrons

Still not working, just saving code.

src/modules/init/moduleInput.f90

@@ -273,6 +273,10 @@ MODULE moduleInput
           !Read BC
           CALL readEMBoundary(config)
 
+        CASE("ElectrostaticBoltzmann")
+          !Read BC
+          CALL readEMBoundary(config)
+
         CASE("ConstantB")
           !Read BC
           CALL readEMBoundary(config)

src/modules/solver/electromagnetic/moduleEM.f90

@@ -49,7 +49,7 @@ MODULE moduleEM
     END SUBROUTINE
 
     !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
@@ -58,15 +58,16 @@ 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
       CLASS(meshNode), POINTER:: node
 
-      !$OMP SINGLE
+      ! !$OMP SINGLE
       vectorF = 0.D0
-      !$OMP END SINGLE
+      ! !$OMP END SINGLE
 
-      !$OMP DO REDUCTION(+:vectorF)
+      ! !$OMP DO REDUCTION(+:vectorF)
       DO e = 1, mesh%numCells
         nNodes = mesh%cells(e)%obj%nNodes
         nodes  = mesh%cells(e)%obj%getNodes(nNodes)
@@ -77,6 +78,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
 
@@ -94,10 +99,10 @@ MODULE moduleEM
         DEALLOCATE(nodes, rho)
 
       END DO
-      !$OMP END DO
+      ! !$OMP END DO
 
       !Apply boundary conditions
-      !$OMP DO
+      ! !$OMP DO
       DO i = 1, mesh%numNodes
         node => mesh%nodes(i)%obj
 
@@ -108,7 +113,7 @@ MODULE moduleEM
         END SELECT
 
       END DO
-      !$OMP END DO
+      ! !$OMP END DO
 
     END SUBROUTINE assembleSourceVector
 
@@ -156,4 +161,80 @@ 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 = -520.0D0, T_e = 11604.0
+
+      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:: n_e(:), phi_old(:)
+      INTEGER:: k
+      EXTERNAL:: dgetrs
+
+      !$OMP SINGLE
+      ALLOCATE(tempF(1:mesh%numNodes))
+      ALLOCATE(n_e(1:mesh%numNodes))
+      ALLOCATE(phi_old(1:mesh%numNodes))
+      DO n = 1, mesh%numNodes
+        tempF(n) = mesh%nodes(n)%obj%emData%phi
+
+      END DO
+      n_e = BoltzmannElectron(tempF, mesh%numNodes)
+      !$OMP END SINGLE
+
+      CALL assembleSourceVector(tempF, n_e)
+
+      !$OMP SINGLE
+      DO k = 1, 5
+        phi_old = tempF
+        CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
+                                mesh%IPIV, tempF,  mesh%numNodes, info)
+
+        PRINT*, k, "diff = ", MAXVAL(ABS(tempF - phi_old))
+        n_e = BoltzmannElectron(tempF, mesh%numNodes)
+        CALL assembleSourceVector(tempF, n_e)
+
+      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 = tempF(n)
+
+        END DO
+        !$OMP END DO
+
+      ELSE
+        !$OMP SINGLE
+        CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann')
+        !$OMP END SINGLE
+
+      END IF
+
+      !$OMP SINGLE
+      DEALLOCATE(tempF)
+      !$OMP END SINGLE
+
+    END SUBROUTINE solveElecFieldBoltzmann
+
 END MODULE moduleEM

src/modules/solver/moduleSolver.f90

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