New input options and fCUm now calculates the number of ions passing by (to plot dN/dE easily

README.md

@@ -5,3 +5,15 @@ It solves the Vlasov equation for ions in a 1D spherical domain.
 Electrons are treated as a fluid either isothermal (Boltzmann) or polytropic.
 
 Self-consistent electric field is solved using a Newton-Raphson method to solve the Poisson equation.
+
+A 2D table linking T_e, n_e and the average Z must be provided. Examples can be found in 'data/TNZ'.
+
+Input files for different cases are in the 'input' folder.
+
+Files with the boundary conditions are provided in 'data/boundary'.
+
+Input files must specify which TNZ table and boundary file are using.
+
+To run the code execute
+
+./vlaplex <path-to-input-file>

input/input_80ns_T10.nml

@@ -13,9 +13,14 @@
   nv = 301     ! Number of grid points in the velocity space
 /
 
+&output
+  folder = 'polytropic_80ns_T10' ! Folder name
+  outputStep = 1e-9              ! Time step for file write [s]
+/
+
 &time
   t0  = 0.0    ! Initial time [s]
-  tf  = 3.0e-7 ! Final time [s]
+  tf  = 2.0e-7 ! Final time [s]
   cfl = 0.5    ! CFL condition (dt = CFL * dr / vf)
 /
 
@@ -31,3 +36,8 @@
   filename = 'data/TNZ/Sn.csv' ! File with table to get Z from Te and ne
   ZList = 0.1, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0 ! Bins of average Z
 /
+
+&parallel
+  nThreads = 16 ! Number of threads for OpenMP
+/
+

input/input_80ns_T30.nml

@@ -13,9 +13,14 @@
   nv = 301     ! Number of grid points in the velocity space
 /
 
+&output
+  folder = 'polytropic_80ns_T30' ! Folder name
+  outputStep = 1e-9              ! Time step for file write [s]
+/
+
 &time
   t0  = 0.0    ! Initial time [s]
-  tf  = 3.0e-7 ! Final time [s]
+  tf  = 2.0e-7 ! Final time [s]
   cfl = 0.5    ! CFL condition (dt = CFL * dr / vf)
 /
 
@@ -31,3 +36,8 @@
   filename = 'data/TNZ/Sn.csv' ! File with table to get Z from Te and ne
   ZList = 0.1, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0 ! Bins of average Z
 /
+
+&parallel
+  nThreads = 16 ! Number of threads for OpenMP
+/
+

input/input_80ns_T6.nml

@@ -13,9 +13,14 @@
   nv = 301     ! Number of grid points in the velocity space
 /
 
+&output
+  folder = 'polytropic_80ns_T6'  ! Folder name
+  outputStep = 1e-9              ! Time step for file write [s]
+/
+
 &time
   t0  = 0.0    ! Initial time [s]
-  tf  = 3.0e-7 ! Final time [s]
+  tf  = 2.0e-7 ! Final time [s]
   cfl = 0.5    ! CFL condition (dt = CFL * dr / vf)
 /
 
@@ -31,3 +36,8 @@
   filename = 'data/TNZ/Sn.csv' ! File with table to get Z from Te and ne
   ZList = 0.1, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0 ! Bins of average Z
 /
+
+&parallel
+  nThreads = 16 ! Number of threads for OpenMP
+/
+

input/input_80ns_T60.nml

@@ -13,9 +13,14 @@
   nv = 301     ! Number of grid points in the velocity space
 /
 
+&output
+  folder = 'polytropic_80ns_T60' ! Folder name
+  outputStep = 1e-9              ! Time step for file write [s]
+/
+
 &time
   t0  = 0.0    ! Initial time [s]
-  tf  = 3.0e-7 ! Final time [s]
+  tf  = 2.0e-7 ! Final time [s]
   cfl = 0.5    ! CFL condition (dt = CFL * dr / vf)
 /
 
@@ -31,3 +36,8 @@
   filename = 'data/TNZ/Sn.csv' ! File with table to get Z from Te and ne
   ZList = 0.1, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0 ! Bins of average Z
 /
+
+&parallel
+  nThreads = 16 ! Number of threads for OpenMP
+/
+

input/input_fa_T30.nml

@@ -13,9 +13,14 @@
   nv = 301     ! Number of grid points in the velocity space
 /
 
+&output
+  folder = 'polytropic_fa_T30' ! Folder name
+  outputStep = 1e-9              ! Time step for file write [s]
+/
+
 &time
   t0  = 0.0    ! Initial time [s]
-  tf  = 3.0e-7 ! Final time [s]
+  tf  = 2.0e-7 ! Final time [s]
   cfl = 0.5    ! CFL condition (dt = CFL * dr / vf)
 /
 
@@ -31,3 +36,8 @@
   filename = 'data/TNZ/Sn.csv' ! File with table to get Z from Te and ne
   ZList = 0.1, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0 ! Bins of average Z
 /
+
+&parallel
+  nThreads = 16 ! Number of threads for OpenMP
+/
+

makefile

@@ -10,7 +10,7 @@ FC := gfortran
 
 # compiler flags
 # gfortran:
-FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -Wall -march=native -g
+FCFLAGS := -fopenmp -Ofast -g -J $(MODDIR) -Wall -march=native
 
 #Output file
 OUTPUT = vlaplex

scripts_python/plotBC.py

@@ -4,7 +4,7 @@ import numpy as np
 import readBC
 
 
-paths = ['../2025-04-10_08.40.09/']
+paths = ['../2025-04-15_13.33.28/']
 time, n_i, u_i, T_i, Zinj, Z_Tne = readBC.read(paths[0] + 'bc.csv')
 
 fig, ax = plt.subplots()
@@ -14,7 +14,7 @@ n_e = n_i * Zinj
 plt.plot(time, n_e / n_e[0], label = f"$n_e$ ($\\times {n_e[0] * 1e-6:.0e} \\; cm^{{-3}})$")
 plt.plot(time, T_i / T_i[0], label = f"$T  \\; (\\times{T_i[0]:.1f} \\; eV)$")
 plt.plot(time, Zinj, label = "Zinj")
-plt.plot(time, Z_Tne, label = "Z_Tne")
+#  plt.plot(time, Z_Tne, label = "Z_Tne")
 plt.semilogy()
 plt.legend()
 plt.show()

scripts_python/plotCumF.py

@@ -16,27 +16,30 @@ m_i = 1.9712e-25
 #  paths = ['../2024-12-10_18.45.17/', '../Poisson_50ns_T30Z11/']
 #  paths = ['../2024-12-11_12.38.27/', '../Poisson_polytropic_fa_T30Z11/', '../Poisson_fa_T30Z11/']
 #  paths = ['../Poisson_partialAblation/','../Poisson_partialAblation_lowerT/','../Poisson_partialAblation_lowT/','../Poisson_partialAblation_highT/']
-paths = ['../polytropic_80ns_T30/']
+paths = ['../2025-04-17_14.39.34/']
 labels = [path[3:-1] for path in paths]
 
+m2s2_to_eV = m_i*0.5/e
+
 for path, label in zip(paths, labels):
     Zlist      = readZlist.read(path+'ZList.csv')
     filesCum_i = sorted(glob.glob(path+'time_*_fCum_i.csv'))
     _, _, v, _ = readF.read(filesCum_i[-1])
     sumF = np.zeros(len(v))
+    E = v**2*m2s2_to_eV
 
     for Z in Zlist:
         filename='time_*_Z_{:.1f}_fCum_i.csv'.format(Z)
         filesCum_i = sorted(glob.glob(path+filename))
-        time, x, v, f_i = readF.read(filesCum_i[-1])
+        time, rCum, v, f_i = readF.read(filesCum_i[-1])
         sumF += f_i[0]
-        plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label=Z)
+        plt.plot(E, 4.0*np.pi*rCum[0]**2*f_i[0]/E, label=Z)
+    plt.plot(E, 4.0*np.pi*rCum[0]**2*sumF/E, label='sum', color='k', linestyle='dashed')
     print(time)
-    plt.plot(v**2*m_i*0.5/e, sumF*e/m_i/v, label='sum', color='k', linestyle='dashed')
 
 plt.yscale('log')
-plt.ylim([1e16,5e27])
+plt.ylim([1e8,5e11])
-plt.ylabel('Sum f(e) / sqrt(e) (m^-3 eV^-1)')
+plt.ylabel('dN / dE (eV^-1)')
 plt.xscale('log')
 plt.xlim([1e0,1e4])
 plt.xlabel('e (eV)')

scripts_python/plotMom.py

@@ -12,7 +12,8 @@ from scipy.constants import e, k
 #  paths = ['../quasiNeutral_partialAblation/','../Poisson_partialAblation/']
 #  paths = ['../2024-10-02_14.30.44/']
 #  paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/']
-paths = ['../2025-04-09_16.45.52/']
+#  paths = ['../2025-04-10_11.59.02/']
+paths = ['../polytropic_80ns_T30/']
 labels = [path[3:-1] for path in paths]
 
 for path, label in zip(paths, labels):
@@ -21,7 +22,7 @@ for path, label in zip(paths, labels):
     start = 80
     end   = 85#len(filesPhi)
     every = 1
-    fig, ax = plt.subplots(4, sharex='all')
+    fig, ax = plt.subplots(3, sharex='all')
     ax[1].set_yscale('log')
     ax[1].set_ylim(bottom=1e10, top=1e24)
     _, r, _, _, _         = readPhi.read(filesPhi[0])
@@ -50,7 +51,6 @@ for path, label in zip(paths, labels):
         ax[1].plot(r, sum_Zni)
         ax[1].plot(r,      n_e, color='k', linestyle='dashed')
         ax[2].plot(r, ave_ui)
-        ax[3].plot(r, ave_Ti)
 
     ax[0].set_title(label)
     ax[0].legend()

scripts_python/readF.py

@@ -8,7 +8,7 @@ def read(filename):
     df = pandas.read_csv(filename,skiprows=2,nrows=1)
     Z = df['Z'].to_numpy()[0]
 
-    df = pandas.read_csv(filename,skiprows=2,nrows=1,header=None)
+    df = pandas.read_csv(filename,skiprows=4,nrows=1,header=None)
     x = df.to_numpy()[0][1:]
     df = pandas.read_csv(filename,skiprows=5,header=None)
     f = []

src/modules/input.f90

@@ -21,6 +21,24 @@ module input
 
     end subroutine openInputFile
 
+    subroutine checkIO(subroutineName)
+
+      character(len=*), intent(in):: subroutineName
+
+      if     (inputFile_io < 0) then
+        write (*, '("End-of-file found when reading in subroutine: ",a)') subroutineName
+        error stop
+
+      elseif (inputFile_io > 0) then
+        write (*, '("Error condition ", i3," found when reading in subroutine: ",a)') inputFile_io, subroutineName
+        error stop
+
+      end if
+
+      inputFile_io = 0
+
+    end subroutine checkIO
+
     subroutine readReference()
       use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI
       use referenceValues
@@ -28,6 +46,7 @@ module input
       namelist /reference/ m_ref, Temp_ref, n_ref
 
       read(nml=reference, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
 
       ! Set reference numbers (in SI units)
       Temp_ref = Temp_ref * eV_to_K
@@ -37,6 +56,8 @@ module input
       L_ref    = u_ref * t_ref
       phi_ref  = kb * Temp_ref / qe
 
+      call checkIO('readReference')
+
     end subroutine readReference
 
     subroutine readGrid(r0, rf, dr, v0, vf, nv)
@@ -48,6 +69,9 @@ module input
       namelist /grid/ r0, rf, dr, v0, vf, nv
 
       read(nml=grid, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readGrid')
 
       r0 = r0/L_ref
       rf = rf/L_ref
@@ -58,6 +82,22 @@ module input
 
     end subroutine readGrid
 
+    subroutine readOutput(folder_alloc, outputStep)
+
+      character(:), allocatable, intent(out):: folder_alloc
+      real(dp), intent(out):: outputStep
+      character(len=128):: folder
+      namelist /output/ folder, outputStep
+
+      read(nml=output, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readOutput')
+
+      folder_alloc = trim(folder)
+
+    end subroutine readOutput
+
     subroutine readTime(t0, tf, CFL)
       use referenceValues, only: t_ref
 
@@ -66,6 +106,9 @@ module input
       namelist /time/ t0, tf, CFL
 
       read(nml=time, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readTime')
 
       t0 = t0/t_ref
       tf = tf/t_ref
@@ -79,6 +122,9 @@ module input
       namelist /detector/ rCum
 
       read(nml=detector, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readDetector')
 
       ! Set position to calculate cumulative sum of f (non-dimensional units)
       rCum = rCum/L_ref
@@ -94,12 +140,16 @@ module input
       namelist /boundary/ filename
 
       read(nml=boundary, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readBoundary')
 
       filename_dynamic = trim(filename)
 
       call bc%init(filename_dynamic)
 
     end subroutine readBoundary
+
     subroutine readZ(Zlist, nz, Tene_to_Z)
       use tableTNZ
       
@@ -117,6 +167,9 @@ module input
 
       ! Read namelist
       read(nml=Zbins, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readZ')
 
       ! Get the real Z bins
       nz = count(ZList > 0.0_dp)
@@ -133,6 +186,18 @@ module input
 
     end subroutine readZ
 
+    subroutine readParallel(nThreads)
+
+      integer, intent(out):: nThreads
+      namelist /parallel/ nThreads
+
+      read(nml=parallel, unit=inputFile_id, iostat=inputFile_io)
+      rewind(unit=inputFile_id)
+
+      call checkIO('readParallel')
+
+    end subroutine readParallel
+
     subroutine closeInputFile()
 
       close(inputFile_id)

src/modules/output.f90

@@ -37,15 +37,17 @@ module output
 
     end function setZFormat
 
-    subroutine createPath()
+    subroutine createPath(folder)
-      character(8)  :: date_now
+      character(8) :: date_now
-      character(10) :: time_now
+      character(10):: time_now
+      character(:), allocatable:: folder
 
       call date_and_time(date_now, time_now)
 
       !Compose the folder name
       pathOutput = date_now(1:4) // '-' // date_now(5:6) // '-' // date_now(7:8) // '_' // &
-                   time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6) // '/'
+                   time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6) // '_' // &
+                   folder        // '/'
 
       call system('mkdir ' // pathOutput)
 

src/modules/tableTNZ.f90

@@ -31,6 +31,13 @@ module tableTNZ
       integer:: id = 20
       num_ne = 0
 
+      inquire(file=tableFile, iostat=stat)
+      if (stat /= 0) then
+          write (*, '("Error: TNZ table file ", a, " does not exist")') tableFile
+          return
+
+      end if
+
       open(id, file = tableFile)
       amount = -1 ! Remove header
       do

src/vlaplex.f90

@@ -62,8 +62,10 @@ program VlaPlEx
   real(dp), allocatable, dimension(:,:):: fCum_i
   real(dp):: rCum
   integer:: rCum_index
+  real(dp):: outputStep
+  character(:), allocatable:: folder
 
-  character(len=128) arg
+  character(len=128):: arg
 
   CALL get_command_argument(1, arg)
   if (arg == '') then
@@ -74,11 +76,12 @@ program VlaPlEx
   inputFile = trim(arg)
   call openInputFile()
   call readReference()
+  call readParallel(nThreads)
   call readGrid(r0, rf, dr, v0, vf, nv)
+  call readOutput(folder, outputStep)
   call readTime(t0, tf, CFL)
   call readDetector(rCum)
   call readBoundary(boundaryConditions)
-  nThreads = 16
   call readZ(Zlist, nz, Tene_to_Z)
   call closeInputFile()
 
@@ -113,7 +116,7 @@ program VlaPlEx
   nt = nint((tf - t0) / dt) 
   dt =      (tf - t0) / float(nt)
 
-  everyOutput = nint(1.0e-9_dp/t_ref/dt)
+  everyOutput = nint(outputStep/t_ref/dt)
   if (everyOutput  == 0) then
     everyOutput = 1
 
@@ -192,7 +195,7 @@ program VlaPlEx
   f0 = 0.0_dp
 
   ! Output initial values
-  call createPath()
+  call createPath(folder)
   call setTimeFormat(nt)
   t = 0
   call writeOutputRef()
@@ -386,7 +389,7 @@ program VlaPlEx
       if (all(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then
           cycle
       end if
-      fCum_i(iz,:)  = fCum_i(iz,:) + f_i_old(iz,rCum_index,:)
+      fCum_i(iz,:)  = fCum_i(iz,:) + f_i_old(iz,rCum_index,:)*dt*(v+dv/2.0)*dv
     end do
     ! Write output
     if (mod(t,everyOutput) == 0 .or. t == nt) then