From 50b4258c8fc74a5253cb2825dedd78325fbc88c6 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 10 Apr 2025 11:56:05 +0200 Subject: [PATCH] Input file This should not have been a single commit. Now we have input files. Also, I've restructured the code and renamed modules. --- .gitignore | 2 + .../boundary/bc_80ns_T10.csv | 0 .../boundary/bc_80ns_T30.csv | 0 .../boundary/bc_80ns_T6.csv | 0 .../boundary/bc_80ns_T60.csv | 0 bc_fa_T30.csv => data/boundary/bc_fa_T30.csv | 0 input/input_80ns_T10.nml | 33 ++++ input/input_80ns_T30.nml | 33 ++++ input/input_80ns_T6.nml | 33 ++++ input/input_80ns_T60.nml | 33 ++++ input/input_fa_T30.nml | 33 ++++ makefile | 44 +++-- moduleTableTtoZ.f90 | 116 ------------ src/makefile | 16 ++ .../modules/constantParameters.f90 | 0 src/modules/input.f90 | 142 ++++++++++++++ src/modules/makefile | 15 ++ moduleOutput.f90 => src/modules/output.f90 | 0 .../modules/referenceValues.f90 | 1 + .../modules/tableBoundary.f90 | 6 +- .../modules/tableTNZ.f90 | 12 +- vlaplex.f90 => src/vlaplex.f90 | 178 +++++++----------- 22 files changed, 446 insertions(+), 251 deletions(-) rename bc_80ns_T10.csv => data/boundary/bc_80ns_T10.csv (100%) rename bc_80ns_T30.csv => data/boundary/bc_80ns_T30.csv (100%) rename bc_80ns_T6.csv => data/boundary/bc_80ns_T6.csv (100%) rename bc_80ns_T60.csv => data/boundary/bc_80ns_T60.csv (100%) rename bc_fa_T30.csv => data/boundary/bc_fa_T30.csv (100%) create mode 100644 input/input_80ns_T10.nml create mode 100644 input/input_80ns_T30.nml create mode 100644 input/input_80ns_T6.nml create mode 100644 input/input_80ns_T60.nml create mode 100644 input/input_fa_T30.nml delete mode 100644 moduleTableTtoZ.f90 create mode 100644 src/makefile rename moduleConstantParameters.f90 => src/modules/constantParameters.f90 (100%) create mode 100644 src/modules/input.f90 create mode 100644 src/modules/makefile rename moduleOutput.f90 => src/modules/output.f90 (100%) rename moduleReferenceValues.f90 => src/modules/referenceValues.f90 (77%) rename moduleTableBC.f90 => src/modules/tableBoundary.f90 (97%) rename moduleTableTtoZne.f90 => src/modules/tableTNZ.f90 (95%) rename vlaplex.f90 => src/vlaplex.f90 (73%) diff --git a/.gitignore b/.gitignore index 3320179..55fa2bd 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ vlaplex *.mod *.o *.tar.gz +obj/ +mod/ diff --git a/bc_80ns_T10.csv b/data/boundary/bc_80ns_T10.csv similarity index 100% rename from bc_80ns_T10.csv rename to data/boundary/bc_80ns_T10.csv diff --git a/bc_80ns_T30.csv b/data/boundary/bc_80ns_T30.csv similarity index 100% rename from bc_80ns_T30.csv rename to data/boundary/bc_80ns_T30.csv diff --git a/bc_80ns_T6.csv b/data/boundary/bc_80ns_T6.csv similarity index 100% rename from bc_80ns_T6.csv rename to data/boundary/bc_80ns_T6.csv diff --git a/bc_80ns_T60.csv b/data/boundary/bc_80ns_T60.csv similarity index 100% rename from bc_80ns_T60.csv rename to data/boundary/bc_80ns_T60.csv diff --git a/bc_fa_T30.csv b/data/boundary/bc_fa_T30.csv similarity index 100% rename from bc_fa_T30.csv rename to data/boundary/bc_fa_T30.csv diff --git a/input/input_80ns_T10.nml b/input/input_80ns_T10.nml new file mode 100644 index 0000000..c1373b5 --- /dev/null +++ b/input/input_80ns_T10.nml @@ -0,0 +1,33 @@ +&reference + m_ref = 1.9712258e-25 ! Ion mass [kg] + Temp_ref = 30.0 ! Plasma temperature [eV] + n_ref = 1.0e24 ! Plasma density [m^-3] +/ + +&grid + r0 = 10.0e-6 ! Initial position [m] + rf = 2.0e-3 ! Final position [m] + dr = 1.0e-6 ! Spatial step [m] + v0 = -80e3 ! Lower limit for velocity space [m s^-1] + vf = 240e3 ! Upper limit for velocity space [m s^-1] + nv = 301 ! Number of grid points in the velocity space +/ + +&time + t0 = 0.0 ! Initial time [s] + tf = 3.0e-7 ! Final time [s] + cfl = 0.5 ! CFL condition (dt = CFL * dr / vf) +/ + +&detector + rCum = 1.0e-3 ! Position of the detector [m] +/ + +&boundary + filename = 'data/boundary/bc_80ns_T10.csv' ! File for boundary value +/ + +&Zbins + 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 +/ diff --git a/input/input_80ns_T30.nml b/input/input_80ns_T30.nml new file mode 100644 index 0000000..6a9d7fb --- /dev/null +++ b/input/input_80ns_T30.nml @@ -0,0 +1,33 @@ +&reference + m_ref = 1.9712258e-25 ! Ion mass [kg] + Temp_ref = 30.0 ! Plasma temperature [eV] + n_ref = 1.0e24 ! Plasma density [m^-3] +/ + +&grid + r0 = 10.0e-6 ! Initial position [m] + rf = 2.0e-3 ! Final position [m] + dr = 1.0e-6 ! Spatial step [m] + v0 = -80e3 ! Lower limit for velocity space [m s^-1] + vf = 240e3 ! Upper limit for velocity space [m s^-1] + nv = 301 ! Number of grid points in the velocity space +/ + +&time + t0 = 0.0 ! Initial time [s] + tf = 3.0e-7 ! Final time [s] + cfl = 0.5 ! CFL condition (dt = CFL * dr / vf) +/ + +&detector + rCum = 1.0e-3 ! Position of the detector [m] +/ + +&boundary + filename = 'data/boundary/bc_80ns_T30.csv' ! File for boundary value +/ + +&Zbins + 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 +/ diff --git a/input/input_80ns_T6.nml b/input/input_80ns_T6.nml new file mode 100644 index 0000000..30f01f5 --- /dev/null +++ b/input/input_80ns_T6.nml @@ -0,0 +1,33 @@ +&reference + m_ref = 1.9712258e-25 ! Ion mass [kg] + Temp_ref = 30.0 ! Plasma temperature [eV] + n_ref = 1.0e24 ! Plasma density [m^-3] +/ + +&grid + r0 = 10.0e-6 ! Initial position [m] + rf = 2.0e-3 ! Final position [m] + dr = 1.0e-6 ! Spatial step [m] + v0 = -80e3 ! Lower limit for velocity space [m s^-1] + vf = 240e3 ! Upper limit for velocity space [m s^-1] + nv = 301 ! Number of grid points in the velocity space +/ + +&time + t0 = 0.0 ! Initial time [s] + tf = 3.0e-7 ! Final time [s] + cfl = 0.5 ! CFL condition (dt = CFL * dr / vf) +/ + +&detector + rCum = 1.0e-3 ! Position of the detector [m] +/ + +&boundary + filename = 'data/boundary/bc_80ns_T6.csv' ! File for boundary value +/ + +&Zbins + 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 +/ diff --git a/input/input_80ns_T60.nml b/input/input_80ns_T60.nml new file mode 100644 index 0000000..b9ef995 --- /dev/null +++ b/input/input_80ns_T60.nml @@ -0,0 +1,33 @@ +&reference + m_ref = 1.9712258e-25 ! Ion mass [kg] + Temp_ref = 30.0 ! Plasma temperature [eV] + n_ref = 1.0e24 ! Plasma density [m^-3] +/ + +&grid + r0 = 10.0e-6 ! Initial position [m] + rf = 2.0e-3 ! Final position [m] + dr = 1.0e-6 ! Spatial step [m] + v0 = -80e3 ! Lower limit for velocity space [m s^-1] + vf = 240e3 ! Upper limit for velocity space [m s^-1] + nv = 301 ! Number of grid points in the velocity space +/ + +&time + t0 = 0.0 ! Initial time [s] + tf = 3.0e-7 ! Final time [s] + cfl = 0.5 ! CFL condition (dt = CFL * dr / vf) +/ + +&detector + rCum = 1.0e-3 ! Position of the detector [m] +/ + +&boundary + filename = 'data/boundary/bc_80ns_T60.csv' ! File for boundary value +/ + +&Zbins + 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 +/ diff --git a/input/input_fa_T30.nml b/input/input_fa_T30.nml new file mode 100644 index 0000000..5f548e6 --- /dev/null +++ b/input/input_fa_T30.nml @@ -0,0 +1,33 @@ +&reference + m_ref = 1.9712258e-25 ! Ion mass [kg] + Temp_ref = 30.0 ! Plasma temperature [eV] + n_ref = 1.0e24 ! Plasma density [m^-3] +/ + +&grid + r0 = 10.0e-6 ! Initial position [m] + rf = 2.0e-3 ! Final position [m] + dr = 1.0e-6 ! Spatial step [m] + v0 = -80e3 ! Lower limit for velocity space [m s^-1] + vf = 240e3 ! Upper limit for velocity space [m s^-1] + nv = 301 ! Number of grid points in the velocity space +/ + +&time + t0 = 0.0 ! Initial time [s] + tf = 3.0e-7 ! Final time [s] + cfl = 0.5 ! CFL condition (dt = CFL * dr / vf) +/ + +&detector + rCum = 1.0e-3 ! Position of the detector [m] +/ + +&boundary + filename = 'data/boundary/bc_fa_T30.csv' ! File for boundary value +/ + +&Zbins + 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 +/ diff --git a/makefile b/makefile index e210cc5..385af6c 100644 --- a/makefile +++ b/makefile @@ -1,15 +1,33 @@ -all: - gfortran moduleConstantParameters.f90 -c -lopenblas -Ofast -fopenmp -Wall - gfortran moduleReferenceValues.f90 -c -lopenblas -Ofast -fopenmp -Wall - gfortran moduleOutput.f90 -c -lopenblas -Ofast -fopenmp -Wall - gfortran moduleTableBC.f90 -c -lopenblas -Ofast -fopenmp -Wall - gfortran moduleTableTtoZ.f90 -c -lopenblas -Ofast -fopenmp -Wall - gfortran moduleTableTtoZne.f90 -c -lopenblas -Ofast -fopenmp -Wall - gfortran vlaplex.f90 moduleConstantParameters.o moduleReferenceValues.o moduleOutput.o moduleTableBC.o moduleTableTtoZ.o moduleTableTtoZne.o -lopenblas -Ofast -fopenmp -Wall -o vlaplex +# set folders +TOPDIR = $(PWD)# top directory +MODDIR := $(TOPDIR)/mod# module folder +OBJDIR := $(TOPDIR)/obj# object folder +SRCDIR := $(TOPDIR)/src# source folder + +# compiler +# gfortran: +FC := gfortran + +# compiler flags +# gfortran: +FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -Wall -march=native -g + +#Output file +OUTPUT = vlaplex + +export + +all: $(OUTPUT) + +$(OUTPUT): src.o + +src.o: + @mkdir -p $(MODDIR) + @mkdir -p $(OBJDIR) + $(MAKE) -C src $(OUTPUT) clean: - rm -f vlaplex - rm -f *.mod - rm -f *.smod - rm -f *.o - + rm -f $(OUTPUT) + rm -f $(MODDIR)/*.mod + rm -f $(MODDIR)/*.smod + rm -f $(OBJDIR)/*.o diff --git a/moduleTableTtoZ.f90 b/moduleTableTtoZ.f90 deleted file mode 100644 index 61106de..0000000 --- a/moduleTableTtoZ.f90 +++ /dev/null @@ -1,116 +0,0 @@ -module moduleTableTtoZ - use constantParameters, only: dp - - type:: tableTtoZ - real(dp):: Z_min, Z_max - real(dp):: Temp_min, Temp_max - real(dp), allocatable, dimension(:):: Temp - real(dp), allocatable, dimension(:):: Z - real(dp), allocatable, dimension(:):: Z_k - contains - procedure, pass:: init => initTableTtoZ - procedure, pass:: get => getValueTableTtoZ - - end type tableTtoZ - - contains - subroutine initTableTtoZ(self, tableFile) - use constantParameters, only: eV_to_K - use referenceValues, only: Temp_ref - implicit none - - class(tableTtoZ), intent(inout):: self - character(:), allocatable, intent(in):: tableFile - character(100):: dummy - integer:: amount - integer:: i - integer:: stat - integer:: id = 20 - - open(id, file = tableFile) - amount = -1 ! Remove header - do - read(id, '(A)', iostat = stat) dummy - !If EOF or error, exit file - if (stat /= 0) EXIT - ! !Skip comment - ! if (index(dummy,'#') /= 0) CYCLE - !Add row - amount = amount + 1 - - end do - - !Go bback to initial point - rewind(id) - - !Allocate table arrays - allocate(self%Temp(1:amount)) - allocate(self%Z(1:amount)) - allocate(self%Z_k(1:amount)) - self%Temp = 0.0_dp - self%Z = 0.0_dp - self%Z_k = 0.0_dp - - i = 0 - read(id, *) ! Skip header - do - read(id, '(A)', iostat = stat) dummy - !TOdo: Make this a function - if (stat /= 0) EXIT - !Add data - !TODO: substitute with extracting information from dummy - backspace(id) - i = i + 1 - read(id, *) self%Temp(i), self%Z(i) - - end do - self%Temp = self%Temp * eV_to_K / Temp_ref - - close(id) - - self%Temp_min = self%Temp(1) - self%Temp_max = self%Temp(amount) - self%Z_min = self%Z(1) - self%Z_max = self%Z(amount) - - do i = 1, amount - 1 - self%Z_k(i) = ( self%Z(i+1) - self%Z(i))/(self%Temp(i+1) - self%Temp(i)) - - end do - - end subroutine initTableTtoZ - - subroutine getValueTableTtoZ(self, Temp, Z) - implicit none - - class(tableTtoZ), intent(in):: self - real(dp), intent(in):: Temp - real(dp), intent(out):: Z - real(dp):: delta_Temp - integer:: i - - if (Temp <= self%Temp_min) THEN - Z = self%Z_min - - elseif (Temp >= self%Temp_max) THEN - Z = self%Z_max - - else - i = minloc(abs(Temp - self%Temp), 1) - delta_Temp = Temp - self%Temp(i) - if (delta_Temp < 0 ) THEN - i = i - 1 - delta_Temp = Temp - self%Temp(i) - - end if - - Z = self%Z(i) + self%Z_k(i)*delta_Temp - - end if - - end subroutine getValueTableTtoZ - - - -end module moduleTableTtoZ - diff --git a/src/makefile b/src/makefile new file mode 100644 index 0000000..8baded5 --- /dev/null +++ b/src/makefile @@ -0,0 +1,16 @@ +OBJECTS = $(OBJDIR)/constantParameters.o \ + $(OBJDIR)/input.o \ + $(OBJDIR)/output.o \ + $(OBJDIR)/referenceValues.o \ + $(OBJDIR)/tableBoundary.o \ + $(OBJDIR)/tableTNZ.o + + +all: $(OUTPUT) + +$(OUTPUT): modules.o $(OUTPUT).f90 + $(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f90 + $(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(OBJDIR)/$(OUTPUT).o -lopenblas + +modules.o: + $(MAKE) -C modules all diff --git a/moduleConstantParameters.f90 b/src/modules/constantParameters.f90 similarity index 100% rename from moduleConstantParameters.f90 rename to src/modules/constantParameters.f90 diff --git a/src/modules/input.f90 b/src/modules/input.f90 new file mode 100644 index 0000000..371f1e8 --- /dev/null +++ b/src/modules/input.f90 @@ -0,0 +1,142 @@ +module input + use constantParameters, only: dp + implicit none + + character(:), allocatable:: inputFile + integer:: inputFile_id, inputFile_io + + + contains + subroutine openInputFile() + + inquire(file=inputFile, iostat=inputfile_io) + + if (inputFile_io /= 0) then + write (*, '("Error: input file ", a, " does not exist")') inputFile + return + + end if + + open(action='read', file=inputFile, iostat=inputFile_io, newunit=inputFile_id) + + end subroutine openInputFile + + subroutine readReference() + use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI + use referenceValues + + namelist /reference/ m_ref, Temp_ref, n_ref + + read(nml=reference, unit=inputFile_id, iostat=inputFile_io) + + ! Set reference numbers (in SI units) + Temp_ref = Temp_ref * eV_to_K + n_ref = n_ref + t_ref = sqrt(eps_0 * m_ref / (n_ref * 1.0_dp * qe**2)) ! 1.0_dp represents Z = 1 for reference values + u_ref = sqrt(kb * Temp_ref / m_ref) + L_ref = u_ref * t_ref + phi_ref = kb * Temp_ref / qe + + end subroutine readReference + + subroutine readGrid(r0, rf, dr, v0, vf, nv) + use referenceValues, only: L_ref, u_ref + + real(dp), intent(out):: r0, rf, dr + real(dp), intent(out):: v0, vf + integer, intent(out):: nv + namelist /grid/ r0, rf, dr, v0, vf, nv + + read(nml=grid, unit=inputFile_id, iostat=inputFile_io) + + r0 = r0/L_ref + rf = rf/L_ref + dr = dr/L_ref + + v0 = v0/u_ref + vf = vf/u_ref + + end subroutine readGrid + + subroutine readTime(t0, tf, CFL) + use referenceValues, only: t_ref + + real(dp), intent(out):: t0, tf, CFL + + namelist /time/ t0, tf, CFL + + read(nml=time, unit=inputFile_id, iostat=inputFile_io) + + t0 = t0/t_ref + tf = tf/t_ref + + end subroutine readTime + + subroutine readDetector(rCum) + use referenceValues, only: L_ref + + real(dp), intent(out):: rCum + namelist /detector/ rCum + + read(nml=detector, unit=inputFile_id, iostat=inputFile_io) + + ! Set position to calculate cumulative sum of f (non-dimensional units) + rCum = rCum/L_ref + + end subroutine readDetector + + subroutine readBoundary(bc) + use tableBoundary + + type(tableBC), intent(out):: bc + character(len=128):: filename + character(:), allocatable:: filename_dynamic ! Needed to pass as an argument, might delete + namelist /boundary/ filename + + read(nml=boundary, unit=inputFile_id, iostat=inputFile_io) + + filename_dynamic = trim(filename) + + call bc%init(filename_dynamic) + + end subroutine readBoundary + subroutine readZ(Zlist, nz, Tene_to_Z) + use tableTNZ + + real(dp), allocatable, intent(out):: Zlist(:) + integer, intent(out):: nz + type(tableTn_to_Z), intent(out):: Tene_to_Z + real(dp), allocatable:: ZList_dummy(:) + character(len=128):: filename + character(:), allocatable:: filename_dynamic ! Needed to pass as an argument, might delete + namelist /Zbins/ ZList, filename + + ! List of dummy size + allocate(ZList(1:99)) + ZList = -1.0_dp + + ! Read namelist + read(nml=Zbins, unit=inputFile_id, iostat=inputFile_io) + + ! Get the real Z bins + nz = count(ZList > 0.0_dp) + allocate(ZList_dummy(1:nz)) + ZList_dummy = ZList(1:nz) + + ! Copy to final list + deallocate(ZList) + ZList = ZList_dummy + deallocate(ZList_dummy) + + filename_dynamic = trim(filename) + call Tene_to_Z%init(filename_dynamic) + + end subroutine readZ + + subroutine closeInputFile() + + close(inputFile_id) + + end subroutine closeInputFile + +end module input diff --git a/src/modules/makefile b/src/modules/makefile new file mode 100644 index 0000000..30e0954 --- /dev/null +++ b/src/modules/makefile @@ -0,0 +1,15 @@ +OBJS = constantParameters.o \ + input.o \ + output.o \ + referenceValues.o \ + tableBoundary.o \ + tableTNZ.o + +all: $(OBJS) + +input.o: referenceValues.o tableBoundary.o tableTNZ.o input.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +%.o: %.f90 + $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ + diff --git a/moduleOutput.f90 b/src/modules/output.f90 similarity index 100% rename from moduleOutput.f90 rename to src/modules/output.f90 diff --git a/moduleReferenceValues.f90 b/src/modules/referenceValues.f90 similarity index 77% rename from moduleReferenceValues.f90 rename to src/modules/referenceValues.f90 index c93679e..5a7e6d9 100644 --- a/moduleReferenceValues.f90 +++ b/src/modules/referenceValues.f90 @@ -2,6 +2,7 @@ module referenceValues use constantParameters, only: dp implicit none + real(dp):: m_ref ! Reference values real(dp):: L_ref, t_ref, n_ref, u_ref, Temp_ref ! Reference values real(dp):: phi_ref ! Reference values diff --git a/moduleTableBC.f90 b/src/modules/tableBoundary.f90 similarity index 97% rename from moduleTableBC.f90 rename to src/modules/tableBoundary.f90 index a79337c..dfcd081 100644 --- a/moduleTableBC.f90 +++ b/src/modules/tableBoundary.f90 @@ -1,4 +1,4 @@ -module moduleTableBC +module tableBoundary use constantParameters, only: dp type:: tableBC @@ -61,7 +61,7 @@ module moduleTableBC read(id, *) ! Skip header do read(id, '(A)', iostat = stat) dummy - !TOdo: Make this a function + !TODO: Make this a function if (stat /= 0) EXIT !Add data !TODO: substitute with extracting information from dummy @@ -132,5 +132,5 @@ module moduleTableBC end subroutine getValueTableBC -end module moduleTableBC +end module tableBoundary diff --git a/moduleTableTtoZne.f90 b/src/modules/tableTNZ.f90 similarity index 95% rename from moduleTableTtoZne.f90 rename to src/modules/tableTNZ.f90 index be9abe3..fc45963 100644 --- a/moduleTableTtoZne.f90 +++ b/src/modules/tableTNZ.f90 @@ -1,7 +1,7 @@ -module moduleTableTtoZne +module tableTNZ use constantParameters, only: dp - type:: tableTtoZne + type:: tableTn_to_Z real(dp) :: Z_min, Z_max real(dp) :: Temp_min, Temp_max real(dp), allocatable, dimension(:):: Temp @@ -12,7 +12,7 @@ module moduleTableTtoZne procedure, pass:: init => initTableTtoZne procedure, pass:: get => getValueTableTtoZne - end type tableTtoZne + end type tableTn_to_Z contains subroutine initTableTtoZne(self, tableFile) @@ -20,7 +20,7 @@ module moduleTableTtoZne use referenceValues, only: Temp_ref, n_ref implicit none - class(tableTtoZne), intent(inout):: self + class(tableTn_to_Z), intent(inout):: self character(:), allocatable, intent(in):: tableFile character(100):: dummy character(len=512) :: line @@ -110,7 +110,7 @@ module moduleTableTtoZne subroutine getValueTableTtoZne(self, Temp, ne, Z) implicit none - class(tableTtoZne), intent(in):: self + class(tableTn_to_Z), intent(in):: self real(dp), intent(in):: Temp, ne real(dp), intent(out):: Z real(dp):: delta_Temp @@ -144,5 +144,5 @@ module moduleTableTtoZne -end module moduleTableTtoZne +end module tableTNZ diff --git a/vlaplex.f90 b/src/vlaplex.f90 similarity index 73% rename from vlaplex.f90 rename to src/vlaplex.f90 index 4682862..5ce22f8 100644 --- a/vlaplex.f90 +++ b/src/vlaplex.f90 @@ -1,37 +1,13 @@ - module eos - use constantParameters, only: dp - implicit none - - private - public:: T_to_Z - - contains - pure function T_to_Z(T) result(Z) - use constantParameters, only: eV_to_K - use referenceValues, only: Temp_ref - implicit none - - real(dp), intent(in):: T - real(dp):: Z - - Z = 22.5 * (Temp_ref * T / eV_to_K / 100.0)**0.6 - - end function T_to_Z - -end module eos - program VlaPlEx use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI + use input use output use referenceValues - use eos, only: T_to_Z - use moduleTableBC - use moduleTableTtoZ - use moduleTableTtoZne + use tableBoundary + use tableTNZ use omp_lib implicit none - real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin atom mass in kg real(dp), parameter:: gamma_i = 1.0_dp ! Adiabatic coefficient for ions real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg real(dp), parameter:: gamma_e = 4.0_dp / 3.0_dp ! Adiabatic coefficient for electrons @@ -44,10 +20,10 @@ program VlaPlEx real(dp):: v0, vf real(dp), allocatable, dimension(:):: v real(dp):: t0, tf + real(dp):: CFL real(dp):: time real(dp):: dr, dv, dt integer:: nr, nv, nt, nz - integer:: nzMin, nzMax integer:: i, iz, j, t, z_inj integer:: j0 ! First integer of positive velocity @@ -55,13 +31,8 @@ program VlaPlEx real(dp):: Zave_bc, Zave_bc_old ! Average charge state real(dp):: u_bc ! Injection velocity real(dp):: n_bc ! Injection density - real(dp):: c_s ! Ion sound speed type(tableBC):: boundaryConditions - character(:), allocatable:: bc_file - type(tableTtoZ):: TtoZ - character(:), allocatable:: TtoZ_file - type(tableTtoZne):: TtoZne - character(:), allocatable:: TtoZne_file + type(tableTn_to_Z):: Tene_to_Z real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0 @@ -72,7 +43,7 @@ program VlaPlEx real(dp), allocatable, dimension(:,:):: T_i real(dp), allocatable, dimension(:):: n_e real(dp), allocatable, dimension(:):: Zave - real(dp), allocatable, dimension(:):: Z_list + real(dp), allocatable, dimension(:):: Zlist real(dp), allocatable, dimension(:):: diag, diag_low, diag_high real(dp), allocatable, dimension(:,:):: A real(dp), allocatable, dimension(:):: Res @@ -86,54 +57,48 @@ program VlaPlEx ! real(dp):: phiF integer:: k + integer:: nThreads ! number of threads for OpenMP + real(dp), allocatable, dimension(:,:):: fCum_i real(dp):: rCum integer:: rCum_index + character(len=128) arg + + CALL get_command_argument(1, arg) + if (arg == '') then + write (*, '(a)') "No input file provided" + return + + end if + inputFile = trim(arg) + call openInputFile() + call readReference() + call readGrid(r0, rf, dr, v0, vf, nv) + call readTime(t0, tf, CFL) + call readDetector(rCum) + call readBoundary(boundaryConditions) + nThreads = 16 + call readZ(Zlist, nz, Tene_to_Z) + call closeInputFile() + ! Set number of threads - call omp_set_num_threads(16) + call omp_set_num_threads(nThreads) - ! Set reference numbers (in SI units) - Temp_ref = 30.0_dp * eV_to_K - n_ref = 1.0e20_dp * cm3_to_m3 - t_ref = sqrt(eps_0 * m_i / (n_ref * 1.0_dp * qe**2)) ! 1.0_dp represents Z = 1 for reference values - u_ref = sqrt(kb * Temp_ref / m_i) - L_ref = u_ref * t_ref - phi_ref = kb * Temp_ref / qe - - ! Set input parameters (remember these have to be in non-dimensional units) - c_s = sqrt(11.0_dp * gamma_i * 1.0_dp) - bc_file = 'bc.csv' - call boundaryConditions%init(bc_file) - TtoZ_file = 'average_TtoZ_curve.csv' - call TtoZ%init(TtoZ_file) - TtoZne_file = 'Zave_values_per_Ne.csv' - call TtoZne%init(TtoZne_file) - - ! Set domain boundaries (non-dimensional units) - r0 = 10.0e-6_dp / L_ref - rf = 2.0e-3_dp / L_ref - dr = 1.0e-6_dp / L_ref + ! Grid in the position space nr = nint((rf - r0) / dr) + 1 dr = (rf - r0) / float(nr-1) - allocate(r(1:nr)) do i = 1, nr r(i) = dr * float(i-1) + r0 end do - ! Set position to calculate cumulative sum of f (non-dimensional units) - rCum = 1.0e-3 / L_ref - ! Index for cumulative sum rCum_index = minloc(abs(r - rCum), 1) - v0 =-0.5e1_dp*c_s - vf = 1.5e1_dp*c_s - dv = 2.0e-1_dp - nv = nint((vf - v0) / dv) + 1 - dv = (vf - v0) / float(nv-1) - + ! Grid in the velocity space + print*, v0, vf + dv = (vf - v0) / float(nv-1) allocate(v(1:nv)) do j = 1, nv v(j) = dv * float(j-1) + v0 @@ -145,9 +110,7 @@ program VlaPlEx j0 = j0 + 1 end if - t0 = 0.0_dp - tf = 3.0e-7_dp / t_ref - dt = 5.0e-2_dp*dr/c_s + dt = CFL*dr/vf nt = nint((tf - t0) / dt) dt = (tf - t0) / float(nt) @@ -163,41 +126,30 @@ program VlaPlEx end if - write(*, '(A,ES0.4e3)') 'CFL: ', dt*vf/dr - - nzMin = 3 - nzMax = 13 - nz = nzMax - nzMin + 1 - nz = nz + 1 ! Add bin for low Z plasma ! Allocate vectors allocate(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv)) allocate(n_i(1:nz,1:nr)) allocate(sum_ni(1:nr)) allocate(u_i(1:nz,1:nr), E_i(1:nr), T_i(1:nz,1:nr)) allocate(Zave(1:nr)) - allocate(Z_list(1:nz)) allocate(n_e(1:nr)) allocate(phi(1:nr), phi_old(1:nr), E(1:nr)) allocate(fCum_i(1:nz,1:nv)) - f_i = 0.0_dp - f_i_old = 0.0_dp - n_i = 0.0_dp - sum_ni = 0.0_dp - u_i = 0.0_dp - E_i = 0.0_dp - T_i = 0.0_dp - n_e = 0.0_dp - T_e = 0.0_dp - Zave = 0.0_dp - Z_list(1) = 0.1_dp ! Low Z bin - do iz = nzMin, nzMax - Z_list(iz-nzMin+1+1) = float(iz) - end do + f_i = 0.0_dp + f_i_old = 0.0_dp + n_i = 0.0_dp + sum_ni = 0.0_dp + u_i = 0.0_dp + E_i = 0.0_dp + T_i = 0.0_dp + n_e = 0.0_dp + T_e = 0.0_dp + Zave = 0.0_dp Zave_bc_old = 0.0_dp - phi = 0.0_dp - phi_old = 0.0_dp - E = 0.0_dp - fCum_i = 0.0_dp + phi = 0.0_dp + phi_old = 0.0_dp + E = 0.0_dp + fCum_i = 0.0_dp ! Allocate matrix for Poisson equation allocate(diag(1:nr), diag_low(1:nr-1), diag_high(1:nr-1)) @@ -246,10 +198,10 @@ program VlaPlEx t = 0 call writeOutputRef() ! call writeOutputF(t, dt, nr, r, nv, v, f_i_old) - call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Z_list) + call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Zlist) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) - call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Z_list) - call writeOutputZList(nz, Z_list) + call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Zlist) + call writeOutputZList(nz, Zlist) ! Main loop do t = 1, nt @@ -258,14 +210,14 @@ program VlaPlEx ! Get boundary conditions for specific time call boundaryConditions%get(time, n_bc, u_bc, Temp_bc) ! Find new \bar{Z}_i based on T_e = Temp_bc and n_e = n_bc - call TtoZne%get(Temp_bc, n_bc, Zave_bc) + call Tene_to_Z%get(Temp_bc, n_bc, Zave_bc) ! Assign Z(T,n) to bin - z_inj = minloc(abs(Z_list - Zave_bc),1) + z_inj = minloc(abs(Zlist - Zave_bc),1) ! Calculate inject (sonic) speed - u_bc = sqrt(Z_list(z_inj)* Temp_bc) + u_bc = sqrt(Zlist(z_inj)* Temp_bc) ! Calculate ion density to inject - n_bc = n_bc / Z_list(z_inj) - call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc, Z_list(z_inj)) + n_bc = n_bc / Zlist(z_inj) + call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc, Zlist(z_inj)) ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) @@ -316,7 +268,7 @@ program VlaPlEx end do !$omp end parallel do - sum_ni = sum_ni + Z_list(iz) * n_i(iz,:) + sum_ni = sum_ni + Zlist(iz) * n_i(iz,:) end do ! Assume quasi-neutrality to start iterating @@ -402,18 +354,18 @@ program VlaPlEx ! i = 1, v<0 i = 1 if (E(i) >= 0.0_dp) then - f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,2:j0-2) - f_i_old(iz,i,1:j0-3)) + f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zlist(iz)*E(i)*dt/dv*(f_i_old(iz,i,2:j0-2) - f_i_old(iz,i,1:j0-3)) else - f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:j0-1) - f_i_old(iz,i,2:j0-2)) + f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zlist(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:j0-1) - f_i_old(iz,i,2:j0-2)) end if ! i = 2, nr-1; all v !$omp parallel do do i = 2, nr-1 if (E(i) >= 0.0_dp) then - f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,2:nv-1) - f_i_old(iz,i,1:nv-2)) + f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zlist(iz)*E(i)*dt/dv*(f_i_old(iz,i,2:nv-1) - f_i_old(iz,i,1:nv-2)) else - f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1)) + f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zlist(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1)) end if @@ -422,9 +374,9 @@ program VlaPlEx ! i = nr, v>=0 i = nr if (E(i) >= 0.0_dp) then - f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,j0+1:nv-1) - f_i_old(iz,i,j0:nv-2)) + f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zlist(iz)*E(i)*dt/dv*(f_i_old(iz,i,j0+1:nv-1) - f_i_old(iz,i,j0:nv-2)) else - f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,j0+2:nv) - f_i_old(iz,i,j0+1:nv-1)) + f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zlist(iz)*E(i)*dt/dv*(f_i_old(iz,i,j0+2:nv) - f_i_old(iz,i,j0+1:nv-1)) end if end do @@ -439,10 +391,10 @@ program VlaPlEx end do ! Write output if (mod(t,everyOutput) == 0 .or. t == nt) then - ! call writeOutputF(t, dt, nz, nr, r, nv, v, f_i_old, Z_list) + ! call writeOutputF(t, dt, nz, nr, r, nv, v, f_i_old, Zlist) call writeOutputPhi(t, dt, nr, r, phi, E, n_e) - call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Z_list) - call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Z_list) + call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Zlist) + call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Zlist) end if