From 9da03ee0ea8e44dca7929ed8a45bb41a42a2ce65 Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Tue, 1 Oct 2024 14:51:46 +0200
Subject: [PATCH] Is it working?

Okay, so after going to Cartesian coordinantes, going back to spherical
coordinates, reviewing Poisson and Vlasov equations in spherical
coordinates (be careful on how you write Vlasov grad.(v f)).
---
 plasmaExpansion.f90 | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index e15f053..b05dc7e 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -395,7 +395,7 @@ program plasmaExpansion
   t0 = 0.0_dp
   tf = 1.0e-6_dp / t_ref
   ! tf = 1.0e1_dp * (rf - r0) / c_s
-  dt = 1.0e-2_dp*dr/(10.0_dp*c_s)
+  dt = 1.0e-1_dp*dr/(10.0_dp*c_s)
   nt = nint((tf - t0) / dt) 
   dt =      (tf - t0) / float(nt)
 
@@ -445,8 +445,8 @@ program plasmaExpansion
   diag_high = 0.0_dp
   b         = 0.0_dp
   diag           = -2.0_dp / dr**2
-  diag_low       =  1.0_dp / dr**2! * r(1:nr-1) / r(2:nr)
-  diag_high      =  1.0_dp / dr**2! * r(2:nr)   / r(1:nr-1)
+  diag_low       =  1.0_dp / dr**2 - 1.0_dp / (r(2:nr)   * 2.0_dp * dr)
+  diag_high      =  1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
   diag(1)        =  1.0_dp ! Dirichlet
   diag_high(1)   =  0.0_dp ! Dirichlet
   ! diag_high(1)   =  2.0_dp / dr**2  ! Neumann
@@ -519,13 +519,13 @@ program plasmaExpansion
     do i = 1, nr
       ! Advect negative velocity
       if (i < nr) then
-        f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr*(f_i_old(i+1,1:j0-1) - &
-                                                             f_i_old(i  ,1:j0-1))
+        f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1) - &
+                                                                     r(i  )**2*f_i_old(i  ,1:j0-1))
       end if
       ! Advect positive velocity
       if (i >  1) then
-        f_i(i,j0:nv)  = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr*(f_i_old(i  , j0:nv) - &
-                                                             f_i_old(i-1, j0:nv))
+        f_i(i,j0:nv)  = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i  )**2*f_i_old(i  , j0:nv) - &
+                                                                     r(i-1)**2*f_i_old(i-1, j0:nv))
       end if
 
       n_i(i) = sum(f_i(i,:)*dv)
@@ -551,8 +551,8 @@ program plasmaExpansion
 
       ! Diagonal matrix for Newton integration scheme
       diag           = -2.0_dp / dr**2 - n_e
-      diag_low       =  1.0_dp / dr**2! * r(1:nr-1) / r(2:nr)
-      diag_high      =  1.0_dp / dr**2! * r(2:nr)   / r(1:nr-1)
+      diag_low       =  1.0_dp / dr**2 - 1.0_dp / (r(2:nr)   * 2.0_dp * dr)
+      diag_high      =  1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
       diag(1)        =  1.0_dp ! Dirichlet
       diag_high(1)   =  0.0_dp ! Dirichlet
       ! diag_high(1)   =  2.0_dp / dr**2 - n_e(1) ! Neumann
@@ -574,13 +574,14 @@ program plasmaExpansion
       ! Iterate system
       call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
       phi = phi_old + Res
+      ! phi0=phi(1)
 
       ! Calculate distribution of electrons
       n_e = Zave(1) * n_i(1) * exp((phi - phi0) / T_i(1))
 
       ! Check if the solution has converged
       phiConv = maxval(abs(Res),1)
-      if (phiConv < 1.0e-5_dp) then
+      if (phiConv < 1.0e-3_dp) then
         exit
 
       end if