From f25abb3213c16ec13218c93da2de31144ea3fcba Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 17 Apr 2025 17:07:07 +0200 Subject: [PATCH] New input options and fCUm now calculates the number of ions passing by (to plot dN/dE easily --- README.md | 12 +++++++ input/input_80ns_T10.nml | 12 ++++++- input/input_80ns_T30.nml | 12 ++++++- input/input_80ns_T6.nml | 12 ++++++- input/input_80ns_T60.nml | 12 ++++++- input/input_fa_T30.nml | 12 ++++++- makefile | 2 +- scripts_python/plotBC.py | 4 +-- scripts_python/plotCumF.py | 15 +++++---- scripts_python/plotMom.py | 6 ++-- scripts_python/readF.py | 2 +- src/modules/input.f90 | 65 ++++++++++++++++++++++++++++++++++++++ src/modules/output.f90 | 10 +++--- src/modules/tableTNZ.f90 | 7 ++++ src/vlaplex.f90 | 13 +++++--- 15 files changed, 169 insertions(+), 27 deletions(-) diff --git a/README.md b/README.md index f7cbda4..e30de96 100644 --- a/README.md +++ b/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 diff --git a/input/input_80ns_T10.nml b/input/input_80ns_T10.nml index c1373b5..d2dc557 100644 --- a/input/input_80ns_T10.nml +++ b/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 / + +¶llel + nThreads = 16 ! Number of threads for OpenMP +/ + diff --git a/input/input_80ns_T30.nml b/input/input_80ns_T30.nml index 6a9d7fb..787acd3 100644 --- a/input/input_80ns_T30.nml +++ b/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 / + +¶llel + nThreads = 16 ! Number of threads for OpenMP +/ + diff --git a/input/input_80ns_T6.nml b/input/input_80ns_T6.nml index 30f01f5..7c9c101 100644 --- a/input/input_80ns_T6.nml +++ b/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 / + +¶llel + nThreads = 16 ! Number of threads for OpenMP +/ + diff --git a/input/input_80ns_T60.nml b/input/input_80ns_T60.nml index b9ef995..d8106c1 100644 --- a/input/input_80ns_T60.nml +++ b/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 / + +¶llel + nThreads = 16 ! Number of threads for OpenMP +/ + diff --git a/input/input_fa_T30.nml b/input/input_fa_T30.nml index 5f548e6..ddbc986 100644 --- a/input/input_fa_T30.nml +++ b/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 / + +¶llel + nThreads = 16 ! Number of threads for OpenMP +/ + diff --git a/makefile b/makefile index 385af6c..de0165b 100644 --- a/makefile +++ b/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 diff --git a/scripts_python/plotBC.py b/scripts_python/plotBC.py index 207869e..9e756d6 100644 --- a/scripts_python/plotBC.py +++ b/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() diff --git a/scripts_python/plotCumF.py b/scripts_python/plotCumF.py index 01a8d31..5c5b500 100644 --- a/scripts_python/plotCumF.py +++ b/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.ylabel('Sum f(e) / sqrt(e) (m^-3 eV^-1)') +plt.ylim([1e8,5e11]) +plt.ylabel('dN / dE (eV^-1)') plt.xscale('log') plt.xlim([1e0,1e4]) plt.xlabel('e (eV)') diff --git a/scripts_python/plotMom.py b/scripts_python/plotMom.py index 2f1dacb..0da99b8 100644 --- a/scripts_python/plotMom.py +++ b/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() diff --git a/scripts_python/readF.py b/scripts_python/readF.py index 0c8cdf6..37ac329 100644 --- a/scripts_python/readF.py +++ b/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 = [] diff --git a/src/modules/input.f90 b/src/modules/input.f90 index 371f1e8..5b2e832 100644 --- a/src/modules/input.f90 +++ b/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) diff --git a/src/modules/output.f90 b/src/modules/output.f90 index 9a3a882..0d3351a 100644 --- a/src/modules/output.f90 +++ b/src/modules/output.f90 @@ -37,15 +37,17 @@ module output end function setZFormat - subroutine createPath() - character(8) :: date_now - character(10) :: time_now + subroutine createPath(folder) + character(8) :: date_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) diff --git a/src/modules/tableTNZ.f90 b/src/modules/tableTNZ.f90 index fc45963..12cb341 100644 --- a/src/modules/tableTNZ.f90 +++ b/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 diff --git a/src/vlaplex.f90 b/src/vlaplex.f90 index 49cbe41..7fc9ebe 100644 --- a/src/vlaplex.f90 +++ b/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