From 98ee3e9c9c38ec6d461c50268a44a99912c8eb70 Mon Sep 17 00:00:00 2001
From: Jorge Gonzalez <jorge.gonzalez@upm.es>
Date: Mon, 30 Sep 2024 17:06:25 +0200
Subject: [PATCH] Still not working, but will be fixed

I have the solution in the plasma expansion code, but I need to do other
stuff.
---
 src/modules/moduleInject.f90                  |  3 +-
 .../solver/electromagnetic/moduleEM.f90       | 50 +++++++++----------
 2 files changed, 26 insertions(+), 27 deletions(-)

diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 4e57083..496ea6a 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -108,7 +108,8 @@ MODULE moduleInject
 
       CASE ("A")
         !Input current in Ampers
-        self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
+        ! TODO: Make this only available for charge species
+        self%nParticles = INT(flow*tauInject*ti_ref/(qe*abs(species(sp)%obj%qm*species(sp)%obj%m)*species(sp)%obj%weight))
 
       CASE ("part/s")
         !Input current in Ampers
diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90
index 0f02bed..a8a5323 100644
--- a/src/modules/solver/electromagnetic/moduleEM.f90
+++ b/src/modules/solver/electromagnetic/moduleEM.f90
@@ -46,7 +46,7 @@ MODULE moduleEM
 
       END DO
 
-    END SUBROUTINE
+    END SUBROUTINE apply
 
     !Assemble the source vector based on the charge density to solve Poisson's equation
     SUBROUTINE assembleSourceVector(vectorF, n_e)
@@ -169,25 +169,16 @@ MODULE moduleEM
       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
+      REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0
       INTEGER:: i
 
-      DO i =1, n 
-        IF (phi(i)*Volt_ref >= phi_0) THEN
-          n_e(i) = n_e0 / n_ref * EXP(-qe * (phi(i)*Volt_ref - phi_0) / (kb * T_e))
-
-        ELSE
-          n_e(i) = 0.D0
-
-        END IF
-
-      END DO
+      n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e))
 
       RETURN
 
     END FUNCTION BoltzmannElectron
 
-    SUBROUTINE solveElecFieldBoltzmann
+    SUBROUTINE solveElecFieldBoltzmann()
       USE moduleMesh
       USE moduleErrors
       IMPLICIT NONE
@@ -195,7 +186,7 @@ MODULE moduleEM
       INTEGER, SAVE:: INFO
       INTEGER:: n
       REAL(8), ALLOCATABLE, SAVE:: tempF(:)
-      REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:)
+      REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:)
       INTEGER:: k
       EXTERNAL:: dgetrs
 
@@ -203,22 +194,29 @@ MODULE moduleEM
       ALLOCATE(tempF(1:mesh%numNodes))
       ALLOCATE(n_e(1:mesh%numNodes))
       ALLOCATE(phi_old(1:mesh%numNodes))
-      n_e = 0.D0
-      CALL assembleSourceVector(tempF, n_e)
-
+      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, 10
-        phi_old = tempF
+      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(tempF), MINVAL(tempF)
-        PRINT*, k, "diff = ", MAXVAL(ABS(tempF - phi_old))
-        n_e = BoltzmannElectron(tempF, mesh%numNodes)
-        CALL assembleSourceVector(tempF, n_e)
+        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
@@ -227,7 +225,7 @@ MODULE moduleEM
         !Suscessful resolution of Poission equation
         !$OMP DO
         DO n = 1, mesh%numNodes
-          mesh%nodes(n)%obj%emData%phi = tempF(n)
+          mesh%nodes(n)%obj%emData%phi = phi_old(n)
 
         END DO
         !$OMP END DO
@@ -240,7 +238,7 @@ MODULE moduleEM
       END IF
 
       !$OMP SINGLE
-      DEALLOCATE(tempF, n_e, phi_old)
+      DEALLOCATE(tempF, n_e, phi_old, phi)
       !$OMP END SINGLE
 
     END SUBROUTINE solveElecFieldBoltzmann