JGonzalez/vlaplex
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Able to read BC from file

Browse source
This commit is contained in:
Jorge Gonzalez 2024-10-03 16:13:04 +02:00
commit 422b9e84e2
7 changed files with 484 additions and 327 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -2,3 +2,4 @@ plasmaExpansion
*.pyc
*.csv
*.mod
*.o

14
makefile
View file

@ -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

18
moduleConstantParameters.f90 Normal file
View file

@ -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

233
moduleOutput.f90 Normal file
View file

@ -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

9
moduleReferenceValues.f90 Normal file
View file

@ -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

145
moduleTableBC.f90 Normal file
View file

@ -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

391
plasmaExpansion.f90
View file

@ -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)