From 422b9e84e2fbd40f14f285da7f139699322f7f73 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 3 Oct 2024 16:13:04 +0200 Subject: [PATCH] Able to read BC from file --- .gitignore | 1 + makefile | 14 +- moduleConstantParameters.f90 | 18 ++ moduleOutput.f90 | 233 +++++++++++++++++++++ moduleReferenceValues.f90 | 9 + moduleTableBC.f90 | 145 +++++++++++++ plasmaExpansion.f90 | 391 ++++++----------------------------- 7 files changed, 484 insertions(+), 327 deletions(-) create mode 100644 moduleConstantParameters.f90 create mode 100644 moduleOutput.f90 create mode 100644 moduleReferenceValues.f90 create mode 100644 moduleTableBC.f90 diff --git a/.gitignore b/.gitignore index eb6b4d9..265aa31 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ plasmaExpansion *.pyc *.csv *.mod +*.o diff --git a/makefile b/makefile index 868f7cc..5a0927f 100644 --- a/makefile +++ b/makefile @@ -1,3 +1,13 @@ -make all: - gfortran plasmaExpansion.f90 -lopenblas -Ofast -fopenmp -Wall -o plasmaExpansion +all: + gfortran moduleConstantParameters.f90 -c + gfortran moduleReferenceValues.f90 -c + gfortran moduleOutput.f90 -c + gfortran moduleTableBC.f90 -c + gfortran plasmaExpansion.f90 moduleConstantParameters.o moduleReferenceValues.o moduleOutput.o moduleTableBC.o -lopenblas -Ofast -fopenmp -Wall -o plasmaExpansion + +clean: + rm -f plasmaExpansion + rm -f *.mod + rm -f *.smod + rm -f *.o diff --git a/moduleConstantParameters.f90 b/moduleConstantParameters.f90 new file mode 100644 index 0000000..ee3a1b8 --- /dev/null +++ b/moduleConstantParameters.f90 @@ -0,0 +1,18 @@ +!Physical and mathematical constants +module constantParameters + implicit none + + public + + integer, parameter:: dp = kind(0.d0) ! Precision + real(dp), parameter:: PI = 4.0_dp*ATAN(1.0_dp) ! Number pi + real(dp), parameter:: qe = 1.60217662e-19_dp ! Elementary charge + real(dp), parameter:: kb = 1.38064852e-23_dp ! Boltzmann constants SI + real(dp), parameter:: eV2J = qe ! Electron volt to Joule conversion + real(dp), parameter:: eps_0 = 8.8542e-12_dp ! Epsilon_0 + + real(dp), parameter:: eV_to_K = 11604.5_dp ! Convert eV to K + real(dp), parameter:: cm3_to_m3 = 1.0e6_dp ! Convert cm^-3 to m^-3 + +end module constantParameters + diff --git a/moduleOutput.f90 b/moduleOutput.f90 new file mode 100644 index 0000000..aeac030 --- /dev/null +++ b/moduleOutput.f90 @@ -0,0 +1,233 @@ +module output + use constantParameters, only: dp + implicit none + + public + private:: dataRef_id, dataBC_id, dataF_id, dataPhi_id, dataCum_id + private:: formatFloat, formatSep, formatTime + + integer:: everyOutput, everyWrite + character(:), allocatable:: pathOutput + integer, parameter:: dataRef_id = 10 + integer, parameter:: dataBC_id = 20 + integer, parameter:: dataF_id = 30 + integer, parameter:: dataPhi_id = 40 + integer, parameter:: dataCum_id = 50 + character(len=7), parameter:: formatFloat = 'ES0.6e3' + character(len=3), parameter:: formatSep = '","' + character(len=7):: formatTime + + contains + subroutine setTimeFormat(nt) + integer, intent(in):: nt + integer:: l + + l=max(1,ceiling(log10(real(nt)))) + write(formatTime, '(A2,I0,".",I0,A1)') '(I',l,l,')' + + end subroutine setTimeFormat + + subroutine createPath() + character(8) :: date_now + character(10) :: time_now + + 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) // '/' + + call system('mkdir ' // pathOutput) + + end subroutine createPath + + subroutine writeOutputF(t, dt, nr, r, nv, v, f) + use referenceValues, only: L_ref, n_ref, u_ref, t_ref + + integer, intent(in):: t + integer, intent(in):: nr, nv + real(dp), intent(in):: dt + real(dp), intent(in):: r(1:nr) + real(dp), intent(in):: v(1:nv) + real(dp), intent(in):: f(1:nr,1:nv) + character(len=30) :: myfmt + character(:), allocatable:: filename + integer:: i + character(len=10):: timeString + + write (timeString, formatTime) t + filename = 'time_' // trim(timeString) // '_f_i.csv' + write (*, '(A, A)') 'Writing: ', filename + + open(unit=dataF_id, file=pathOutput // filename) + write(dataF_id, '(A)') "t (s)" + write(dataF_id, '('//formatFloat//')') t*dt*t_ref + write(myfmt, "(I0)") nr + myfmt = '(A,' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' + write(dataF_id, myfmt) "v (m/s) / r (m)", r*L_ref + write(myfmt, "(I0)") nr + myfmt = '(' // formatFloat // ',' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' + do i = 1, nv + write(dataF_id, myfmt) v(i)*u_ref, f(:,i)*n_ref/u_ref + + end do + + close(unit=dataF_id) + + end subroutine writeOutputF + + subroutine writeOutputPhi(t, dt, nr, r, phi, E, n_e) + use constantParameters, only: eV_to_K + use referenceValues, only: L_ref, phi_ref, t_ref, n_ref + + integer, intent(in):: t + integer, intent(in):: nr + real(dp), intent(in):: dt + real(dp), intent(in):: r(1:nr) + real(dp), intent(in):: phi(1:nr) + real(dp), intent(in):: E(1:nr) + real(dp), intent(in):: n_e(1:nr) + character(:), allocatable:: filename + integer:: i + character(len=10):: timeString + + write (timeString, formatTime) t + filename = 'time_' // trim(timeString)//'_phi.csv' + write (*, '(A, A)') 'Writing: ', filename + + open(unit=dataPhi_id, file=pathOutput//filename) + write(dataPhi_id, '(A)') "t (s)" + write(dataPhi_id, '('//formatFloat//')') t*dt*t_ref + write(dataPhi_id, '(A,3('//formatSep//',A))') "r (m)","phi (V)","E (V m^-1)","n_e (m^-3)" + do i = 1, nr + write(dataPhi_id, '('//formatFloat//',3('//formatSep //','//formatFloat//'))') & + r(i)*L_ref, & + phi(i)*phi_ref, & + E(i)*phi_ref/L_ref, & + n_e(i)*n_ref + + end do + + close(unit=dataPhi_id) + + end subroutine writeOutputPhi + + subroutine writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Z) + use constantParameters, only: eV_to_K + use referenceValues, only: L_ref, t_ref, n_ref, u_ref, Temp_ref + + integer, intent(in):: t + integer, intent(in):: nr + real(dp), intent(in):: dt + real(dp), intent(in):: r(1:nr) + real(dp), intent(in):: n_i(1:nr) + real(dp), intent(in):: u_i(1:nr) + real(dp), intent(in):: T_i(1:nr) + real(dp), intent(in):: Z(1:nr) + character(:), allocatable:: filename + integer:: i + character(len=10):: timeString + + write (timeString, formatTime) t + filename = 'time_' // trim(timeString)//'_mom_i.csv' + write (*, '(A, A)') 'Writing: ', filename + + open(unit=dataPhi_id, file=pathOutput//filename) + write(dataPhi_id, '(A)') "t (s)" + write(dataPhi_id, '('//formatFloat//')') t*dt*t_ref + write(dataPhi_id, '(A,4('//formatSep//',A))') "r (m)","n_i (m^-3)","u_i (m s^-1)", "T_i (eV)","Zave" + do i = 1, nr + write(dataPhi_id, '('//formatFloat//',4('//formatSep //','//formatFloat//'))') & + r(i)*L_ref, & + n_i(i)*n_ref, & + u_i(i)*u_ref, & + T_i(i)*Temp_ref/ev_to_K, & + Z(i) + + end do + + close(unit=dataPhi_id) + + end subroutine writeOutputMom + + subroutine writeOutputBoundary(t, dt, n, u, Temp, Z) + use constantParameters, only: eV_to_K + use referenceValues, only: t_ref, n_ref, u_ref, Temp_ref + + integer, intent(in):: t + real(dp), intent(in):: dt + real(dp), intent(in):: n, u, Temp, Z + character(len=6), parameter:: filename = 'bc.csv' + logical:: res + + inquire(file=pathOutput // filename, exist=res) + if (.not. res) then + write (*, '(A, A)') 'Writing: ', filename + open(unit=dataBC_id, file=pathOutput // filename, action='write', position='append') + write(dataBC_id, '(A,4(' // formatSep // ',A))') 't (s)', 'n (m^-3)', 'u (m s^-1)', 'T (eV)', 'Z' + close(dataBC_id) + + end if + + open(unit=dataBC_id, file=pathOutput // filename, action='write', position='append') + write(dataBC_id, '(' // formatFloat // ',4('// formatSep // ',' // formatFloat // '))') & + t*dt*t_ref, n*n_ref, u*u_ref, Temp*Temp_ref/eV_to_K, Z + + close(dataBC_id) + + end subroutine writeOutputBoundary + + subroutine writeOutputRef() + use referenceValues, only: t_ref, L_ref, n_ref, u_ref, Temp_ref, phi_ref + + character(len=7), parameter:: filename = 'ref.csv' + + write (*, '(A, A)') 'Writing: ', filename + open(unit=dataRef_id, file=pathOutput // filename) + write(dataRef_id, '(A,5(' // formatSep // ',A))') 't_ref (s)', 'L_ref (m)', 'n_ref (m^-3)', & + 'u_ref (m s^-1)', 'T_ref (K)', 'phi_ref (V)' + write(dataRef_id, '(' // formatFloat // ',5('// formatSep // ',' // formatFloat // '))') & + t_ref, L_ref, n_ref, & + u_ref, Temp_ref, phi_ref + + close(dataRef_id) + + end subroutine writeOutputRef + + subroutine writeOutputFCum(t, dt, r, nv, v, f) + use referenceValues, only: L_ref, n_ref, u_ref, t_ref + + integer, intent(in):: t + real(dp), intent(in):: dt + integer, intent(in):: nv + real(dp), intent(in):: r + real(dp), intent(in):: v(1:nv) + real(dp), intent(in):: f(1:nv) + character(len=30) :: myfmt + character(:), allocatable:: filename + integer:: i + character(len=10):: timeString + + write (timeString, formatTime) t + filename = 'time_' // trim(timeString) // '_fCum_i.csv' + write (*, '(A, A)') 'Writing: ', filename + + open(unit=dataCum_id, file=pathOutput // filename) + write(dataCum_id, '(A)') "t (s)" + write(dataCum_id, '('//formatFloat//')') t*dt*t_ref + write(myfmt, "(I0)") 1 + myfmt = '(A,' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' + write(dataCum_id, myfmt) "v (m/s) / r (m)", r*L_ref + write(myfmt, "(I0)") 1 + myfmt = '(' // formatFloat // ',' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' + do i = 1, nv + write(dataCum_id, myfmt) v(i)*u_ref, f(i)*n_ref/u_ref + + end do + + close(unit=dataCum_id) + + end subroutine writeOutputFCum + +end module output + diff --git a/moduleReferenceValues.f90 b/moduleReferenceValues.f90 new file mode 100644 index 0000000..c93679e --- /dev/null +++ b/moduleReferenceValues.f90 @@ -0,0 +1,9 @@ +module referenceValues + use constantParameters, only: dp + implicit none + + real(dp):: L_ref, t_ref, n_ref, u_ref, Temp_ref ! Reference values + real(dp):: phi_ref ! Reference values + +end module referenceValues + diff --git a/moduleTableBC.f90 b/moduleTableBC.f90 new file mode 100644 index 0000000..aa8a09d --- /dev/null +++ b/moduleTableBC.f90 @@ -0,0 +1,145 @@ +module moduleTableBC + use constantParameters, only: dp + + type:: tableBC + real(dp):: t_min, t_max + real(dp):: n_min, n_max + real(dp):: u_min, u_max + real(dp):: Temp_min, Temp_max + real(dp):: Zave_min, Zave_max + real(dp), allocatable, dimension(:):: t + real(dp), allocatable, dimension(:):: n, u, Temp, Zave + real(dp), allocatable, dimension(:):: n_k, u_k, Temp_k, Zave_k + contains + procedure, pass:: init => initTableBC + procedure, pass:: get => getValueTableBC + + end type tableBC + + contains + subroutine initTableBC(self, tableFile) + use constantParameters, only: eV_to_K + use referenceValues, only: t_ref, n_ref, u_ref, Temp_ref + implicit none + + class(tableBC), 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%t(1:amount)) + allocate( self%n(1:amount), self%u(1:amount), self%Temp(1:amount), self%Zave(1:amount)) + allocate(self%n_k(1:amount), self%u_k(1:amount), self%Temp_k(1:amount), self%Zave_k(1:amount)) + self%t = 0.0_dp + self%n = 0.0_dp + self%u = 0.0_dp + self%Temp = 0.0_dp + self%Zave = 0.0_dp + self%n_k = 0.0_dp + self%u_k = 0.0_dp + self%Temp_k = 0.0_dp + self%Zave_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%t(i), self%n(i), self%u(i), self%Temp(i), self%Zave(i) + + end do + + self%t = self%t / t_ref + self%n = self%n / n_ref + self%u = self%u / u_ref + self%Temp = self%Temp * eV_to_K / Temp_ref + + close(id) + + self%t_min = self%t(1) + self%t_max = self%t(amount) + self%n_min = self%n(1) + self%n_max = self%n(amount) + self%u_min = self%u(1) + self%u_max = self%u(amount) + self%Temp_min = self%Temp(1) + self%Temp_max = self%Temp(amount) + self%Zave_min = self%Zave(1) + self%Zave_max = self%Zave(amount) + + do i = 1, amount - 1 + self%n_k(i) = ( self%n(i+1) - self%n(i))/(self%t(i+1) - self%t(i)) + self%u_k(i) = ( self%u(i+1) - self%u(i))/(self%t(i+1) - self%t(i)) + self%Temp_k(i) = (self%Temp(i+1) - self%Temp(i))/(self%t(i+1) - self%t(i)) + self%Zave_k(i) = (self%Zave(i+1) - self%Zave(i))/(self%t(i+1) - self%t(i)) + + end do + + end subroutine initTableBC + + subroutine getValueTableBC(self, t, n, u, Temp, Zave) + implicit none + + class(tableBC), intent(in):: self + real(dp), intent(in):: t + real(dp), intent(out):: n, u, Temp, Zave + real(dp):: delta_t + integer:: i + + if (t <= self%t_min) THEN + n = self%n_min + u = self%u_min + Temp = self%Temp_min + Zave = self%Zave_min + + elseif (t >= self%t_max) THEN + n = self%n_max + u = self%u_max + Temp = self%Temp_max + Zave = self%Zave_max + + else + i = minloc(abs(t - self%t), 1) + delta_t = t - self%t(i) + if (delta_t < 0 ) THEN + i = i - 1 + delta_t = t - self%t(i) + + end if + + n = self%n(i) + self%n_k(i)*delta_t + u = self%u(i) + self%u_k(i)*delta_t + Temp = self%Temp(i) + self%Temp_k(i)*delta_t + Zave = self%Zave(i) + self%Zave_k(i)*delta_t + + end if + + end subroutine getValueTableBC + +end module moduleTableBC + diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index f8f9120..97e4309 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -1,321 +1,61 @@ -!Physical and mathematical constants -module constantParameters - implicit none - - public - - integer, parameter:: dp = kind(0.d0) ! Precision - real(dp), parameter:: PI = 4.0_dp*ATAN(1.0_dp) ! Number pi - real(dp), parameter:: qe = 1.60217662e-19_dp ! Elementary charge - real(dp), parameter:: kb = 1.38064852e-23_dp ! Boltzmann constants SI - real(dp), parameter:: eV2J = qe ! Electron volt to Joule conversion - real(dp), parameter:: eps_0 = 8.8542e-12_dp ! Epsilon_0 - - real(dp), parameter:: eV_to_K = 11604.5_dp ! Convert eV to K - real(dp), parameter:: cm3_to_m3 = 1.0e6_dp ! Convert cm^-3 to m^-3 - -end module constantParameters - -module referenceValues - use constantParameters, only: dp - implicit none - - real(dp):: L_ref, t_ref, n_ref, u_ref, Temp_ref ! Reference values - real(dp):: phi_ref ! Reference values - -end module referenceValues - -module output - use constantParameters, only: dp - implicit none - - public - private:: dataRef_id, dataBC_id, dataF_id, dataPhi_id, dataCum_id - private:: formatFloat, formatSep, formatTime - - integer:: everyOutput, everyWrite - character(:), allocatable:: pathOutput - integer, parameter:: dataRef_id = 10 - integer, parameter:: dataBC_id = 20 - integer, parameter:: dataF_id = 30 - integer, parameter:: dataPhi_id = 40 - integer, parameter:: dataCum_id = 50 - character(len=7), parameter:: formatFloat = 'ES0.6e3' - character(len=3), parameter:: formatSep = '","' - character(len=7):: formatTime - - contains - subroutine setTimeFormat(nt) - integer, intent(in):: nt - integer:: l - - l=max(1,ceiling(log10(real(nt)))) - write(formatTime, '(A2,I0,".",I0,A1)') '(I',l,l,')' - - end subroutine setTimeFormat - - subroutine createPath() - character(8) :: date_now - character(10) :: time_now - - 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) // '/' - - call system('mkdir ' // pathOutput) - - end subroutine createPath - - subroutine writeOutputF(t, dt, nr, r, nv, v, f) - use referenceValues, only: L_ref, n_ref, u_ref, t_ref - - integer, intent(in):: t - integer, intent(in):: nr, nv - real(dp), intent(in):: dt - real(dp), intent(in):: r(1:nr) - real(dp), intent(in):: v(1:nv) - real(dp), intent(in):: f(1:nr,1:nv) - character(len=30) :: myfmt - character(:), allocatable:: filename - integer:: i - character(len=10):: timeString - - write (timeString, formatTime) t - filename = 'time_' // trim(timeString) // '_f_i.csv' - write (*, '(A, A)') 'Writing: ', filename - - open(unit=dataF_id, file=pathOutput // filename) - write(dataF_id, '(A)') "t (s)" - write(dataF_id, '('//formatFloat//')') t*dt*t_ref - write(myfmt, "(I0)") nr - myfmt = '(A,' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' - write(dataF_id, myfmt) "v (m/s) / r (m)", r*L_ref - write(myfmt, "(I0)") nr - myfmt = '(' // formatFloat // ',' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' - do i = 1, nv - write(dataF_id, myfmt) v(i)*u_ref, f(:,i)*n_ref/u_ref - - end do - - close(unit=dataF_id) - - end subroutine writeOutputF - - subroutine writeOutputPhi(t, dt, nr, r, phi, E, n_e) - use constantParameters, only: eV_to_K - use referenceValues, only: L_ref, phi_ref, t_ref, n_ref - - integer, intent(in):: t - integer, intent(in):: nr - real(dp), intent(in):: dt - real(dp), intent(in):: r(1:nr) - real(dp), intent(in):: phi(1:nr) - real(dp), intent(in):: E(1:nr) - real(dp), intent(in):: n_e(1:nr) - character(:), allocatable:: filename - integer:: i - character(len=10):: timeString - - write (timeString, formatTime) t - filename = 'time_' // trim(timeString)//'_phi.csv' - write (*, '(A, A)') 'Writing: ', filename - - open(unit=dataPhi_id, file=pathOutput//filename) - write(dataPhi_id, '(A)') "t (s)" - write(dataPhi_id, '('//formatFloat//')') t*dt*t_ref - write(dataPhi_id, '(A,3('//formatSep//',A))') "r (m)","phi (V)","E (V m^-1)","n_e (m^-3)" - do i = 1, nr - write(dataPhi_id, '('//formatFloat//',3('//formatSep //','//formatFloat//'))') & - r(i)*L_ref, & - phi(i)*phi_ref, & - E(i)*phi_ref/L_ref, & - n_e(i)*n_ref - - end do - - close(unit=dataPhi_id) - - end subroutine writeOutputPhi - - subroutine writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Z) - use constantParameters, only: eV_to_K - use referenceValues, only: L_ref, t_ref, n_ref, u_ref, Temp_ref - - integer, intent(in):: t - integer, intent(in):: nr - real(dp), intent(in):: dt - real(dp), intent(in):: r(1:nr) - real(dp), intent(in):: n_i(1:nr) - real(dp), intent(in):: u_i(1:nr) - real(dp), intent(in):: T_i(1:nr) - real(dp), intent(in):: Z(1:nr) - character(:), allocatable:: filename - integer:: i - character(len=10):: timeString - - write (timeString, formatTime) t - filename = 'time_' // trim(timeString)//'_mom_i.csv' - write (*, '(A, A)') 'Writing: ', filename - - open(unit=dataPhi_id, file=pathOutput//filename) - write(dataPhi_id, '(A)') "t (s)" - write(dataPhi_id, '('//formatFloat//')') t*dt*t_ref - write(dataPhi_id, '(A,4('//formatSep//',A))') "r (m)","n_i (m^-3)","u_i (m s^-1)", "T_i (eV)","Zave" - do i = 1, nr - write(dataPhi_id, '('//formatFloat//',4('//formatSep //','//formatFloat//'))') & - r(i)*L_ref, & - n_i(i)*n_ref, & - u_i(i)*u_ref, & - T_i(i)*Temp_ref/ev_to_K, & - Z(i) - - end do - - close(unit=dataPhi_id) - - end subroutine writeOutputMom - - subroutine writeOutputBoundary(t, dt, n, u, Temp, Z) - use constantParameters, only: eV_to_K - use referenceValues, only: t_ref, n_ref, u_ref, Temp_ref - - integer, intent(in):: t - real(dp), intent(in):: dt - real(dp), intent(in):: n, u, Temp, Z - character(len=6), parameter:: filename = 'bc.csv' - logical:: res - - inquire(file=pathOutput // filename, exist=res) - if (.not. res) then - write (*, '(A, A)') 'Writing: ', filename - open(unit=dataBC_id, file=pathOutput // filename, action='write', position='append') - write(dataBC_id, '(A,4(' // formatSep // ',A))') 't (s)', 'n (m^-3)', 'u (m s^-1)', 'T (eV)', 'Z' - close(dataBC_id) - - end if - - open(unit=dataBC_id, file=pathOutput // filename, action='write', position='append') - write(dataBC_id, '(' // formatFloat // ',4('// formatSep // ',' // formatFloat // '))') & - t*dt*t_ref, n*n_ref, u*u_ref, Temp*Temp_ref/eV_to_K, Z - - close(dataBC_id) - - end subroutine writeOutputBoundary - - subroutine writeOutputRef() - use referenceValues, only: t_ref, L_ref, n_ref, u_ref, Temp_ref, phi_ref - - character(len=7), parameter:: filename = 'ref.csv' - - write (*, '(A, A)') 'Writing: ', filename - open(unit=dataRef_id, file=pathOutput // filename) - write(dataRef_id, '(A,5(' // formatSep // ',A))') 't_ref (s)', 'L_ref (m)', 'n_ref (m^-3)', & - 'u_ref (m s^-1)', 'T_ref (K)', 'phi_ref (V)' - write(dataRef_id, '(' // formatFloat // ',5('// formatSep // ',' // formatFloat // '))') & - t_ref, L_ref, n_ref, & - u_ref, Temp_ref, phi_ref - - close(dataRef_id) - - end subroutine writeOutputRef - - subroutine writeOutputFCum(t, dt, r, nv, v, f) - use referenceValues, only: L_ref, n_ref, u_ref, t_ref - - integer, intent(in):: t - real(dp), intent(in):: dt - integer, intent(in):: nv - real(dp), intent(in):: r - real(dp), intent(in):: v(1:nv) - real(dp), intent(in):: f(1:nv) - character(len=30) :: myfmt - character(:), allocatable:: filename - integer:: i - character(len=10):: timeString - - write (timeString, formatTime) t - filename = 'time_' // trim(timeString) // '_fCum_i.csv' - write (*, '(A, A)') 'Writing: ', filename - - open(unit=dataCum_id, file=pathOutput // filename) - write(dataCum_id, '(A)') "t (s)" - write(dataCum_id, '('//formatFloat//')') t*dt*t_ref - write(myfmt, "(I0)") 1 - myfmt = '(A,' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' - write(dataCum_id, myfmt) "v (m/s) / r (m)", r*L_ref - write(myfmt, "(I0)") 1 - myfmt = '(' // formatFloat // ',' // trim(myfmt) // '(' // formatSep // ',' // formatFloat // '))' - do i = 1, nv - write(dataCum_id, myfmt) v(i)*u_ref, f(i)*n_ref/u_ref - - end do - - close(unit=dataCum_id) - - end subroutine writeOutputFCum - -end module output - -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 = (Temp_ref * T / eV_to_K)**0.6 - ! Z = max(Z, 1.0_dp) - ! Z = min(Z, 22.0_dp) - Z = 12.0_dp - - end function T_to_Z - -end module eos +! 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 = (Temp_ref * T / eV_to_K)**0.6 +! ! Z = max(Z, 1.0_dp) +! ! Z = min(Z, 22.0_dp) +! Z = 12.0_dp +! +! end function T_to_Z +! +! end module eos program plasmaExpansion use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI use output use referenceValues - use eos, only: T_to_Z + ! use eos, only: T_to_Z + use moduleTableBC use omp_lib implicit none - real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin atom mass in kg - real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg - real(dp), parameter:: gam = 1.0_dp ! Adiabatic coefficient + 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 = 5.0_dp/3.0_dp ! Adiabatic coefficient for electrons real(dp):: r0, rf real(dp), allocatable, dimension(:):: r real(dp):: v0, vf real(dp), allocatable, dimension(:):: v real(dp):: t0, tf + real(dp):: time real(dp):: dr, dv, dt integer:: nr, nv, nt integer:: i, j, t integer:: j0 ! First integer of positive velocity real(dp):: Temp_bc ! Temperature - real(dp):: Temp0, TempF real(dp):: Zave_bc ! Average charge state - real(dp):: Zave0, ZaveF - real(dp):: n_ecr ! Electron critical density for the laser - real(dp):: c_s ! Ion sound speed real(dp):: u_bc ! Injection velocity - real(dp):: u_bc0, u_bcF real(dp):: n_bc ! Injection density - real(dp):: n_bc0, n_bcF - integer:: t_bc0, t_bcF + ! real(dp):: n_ecr ! Electron critical density for the laser + real(dp):: c_s ! Ion sound speed + type(tableBC):: boundaryConditions + character(:), allocatable:: bc_file real(dp), allocatable, dimension(:,:):: f_i, f_i_old real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0 @@ -333,7 +73,8 @@ program plasmaExpansion real(dp), allocatable, dimension(:):: phi, phi_old, E real(dp):: phiConv - real(dp):: phi0, phiF + real(dp):: phi0 + ! real(dp):: phiF integer:: k real(dp), allocatable, dimension(:):: fCum_i @@ -344,24 +85,26 @@ program plasmaExpansion call omp_set_num_threads(8) ! Set reference numbers (in SI units) - Temp_ref = 60.0_dp * eV_to_K - n_ref = 1.0e19_dp * cm3_to_m3 + 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) - Temp0 = 60.0_dp * eV_to_K / Temp_ref - TempF = 5.0_dp * eV_to_K / Temp_ref - Zave0 = 12.0_dp - ZaveF = 3.0_dp - n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref - c_s = sqrt(Zave0 * gam * Temp0) - u_bc0 = sqrt(Temp0) - u_bcF = sqrt(TempF) - n_bc0 = n_ecr / Zave0 - n_bcF = 1.0e-1_dp*n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) + ! Temp0 = 30.0_dp * eV_to_K / Temp_ref + ! TempF = 5.0_dp * eV_to_K / Temp_ref + ! Zave0 = 9.0_dp + ! ZaveF = 9.0_dp + ! n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref + c_s = sqrt(9.0_dp * gamma_i * 1.0_dp) + ! u_bc0 = sqrt(Temp0) + ! u_bcF = sqrt(TempF) + ! n_bc0 = n_ecr / Zave0 + ! n_bcF = n_bc0*1.0e-2_dp!n_ecr*1.0e-1 / T_to_Z(Temp0) + bc_file = 'bc.csv' + call boundaryConditions%init(bc_file) ! Set domain boundaries (non-dimensional units) r0 = 200.0e-6_dp / L_ref @@ -403,8 +146,8 @@ program plasmaExpansion nt = nint((tf - t0) / dt) dt = (tf - t0) / float(nt) - t_bc0 = nint(100.0e-9_dp / t_ref / dt) - t_bcF = nint(200.0e-9_dp / t_ref / dt) + ! t_bc0 = nint(100.0e-9_dp / t_ref / dt) + ! t_bcF = nint(150.0e-9_dp / t_ref / dt) everyOutput = nint(1.0e-9_dp/t_ref/dt) if (everyOutput == 0) then @@ -490,22 +233,20 @@ program plasmaExpansion call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave) ! Main loop - Temp_bc = Temp0 - Zave_bc = Zave0 - u_bc = u_bc0 - n_bc = n_bc0 do t = 1, nt - if (t > t_bc0 .and. t <= t_bcF) then - Temp_bc = (TempF - Temp0) / float(t_bcF - t_bc0)*(t - t_bc0) + Temp0 - Zave_bc = (ZaveF - Zave0) / float(t_bcF - t_bc0)*(t - t_bc0) + Zave0 - u_bc = (u_bcF - u_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + u_bc0 - n_bc = (n_bcF - n_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + n_bc0 - else if (t > t_bcF) then - Temp_bc = TempF - Zave_bc = ZaveF - u_bc = u_bcF - n_bc = n_bcF - end if + time = t * dt + t0 + ! if (t > t_bc0 .and. t <= t_bcF) then + ! Temp_bc = (TempF - Temp0) / float(t_bcF - t_bc0)*(t - t_bc0) + Temp0 + ! Zave_bc = (ZaveF - Zave0) / float(t_bcF - t_bc0)*(t - t_bc0) + Zave0 + ! u_bc = (u_bcF - u_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + u_bc0 + ! n_bc = (n_bcF - n_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + n_bc0 + ! else if (t > t_bcF) then + ! Temp_bc = TempF + ! Zave_bc = ZaveF + ! u_bc = u_bcF + ! n_bc = n_bcF + ! end if + call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc) call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc) f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0 = f0 * n_bc / (sum(f0)*dv)