Initialize the field.
subroutine init()
!---------------------------------------------------------------------------------------------------------------------------------
!< Initialize the field.
!---------------------------------------------------------------------------------------------------------------------------------
real(R_P), parameter :: pi=4._R_P * atan(1._R_P) !< Pi greek.
integer(I_P) :: i !< Space counter.
integer(I_P) :: s !< Species counter.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
allocate(xcenter(1:Ni))
allocate(xnode(0:Ni))
Dx=1._R_P/Ni
do i=1, Ni
xcenter(i) = Dx * i - 0.5_R_P * Dx
enddo
do i=0, Ni
xnode(i) = Dx * i
enddo
if ((av_Ni>0).and.(av_Ni/=Ni)) then
allocate(av_xnode(0:av_Ni))
do i=0, av_Ni
av_xnode(i) = 1._R_P/av_Ni * i
enddo
endif
select case(trim(adjustl(problem)))
case('sod')
print "(A)", 'Solving "'//trim(adjustl(problem))//'" problem'
t_final = 0.2_R_P
CFL = 0.7_R_P
Ns = 1
Nc = Ns + 2
Np = Ns + 4
allocate(initial_state(1:Np, 1:Ni))
allocate(cp0(1:Ns))
allocate(cv0(1:Ns))
variables = 'VARIABLES="x"'
do s=1, Ns
variables = variables//' "rho('//trim(str(s,.true.))//')"'
enddo
variables = variables//' "u" "p" "rho" "gamma"'
BC_L = 'TRA'
BC_R = 'TRA'
cp0(1) = 1040._R_P
cv0(1) = 743._R_P
do i=1, Ni/2
initial_state(:, i) = [1._R_P, & ! rho(s)
0._R_P, & ! u
1._R_P, & ! p
1._R_P, & ! sum(rho(s))
cp0/cv0] ! gamma = cp/cv
enddo
do i=Ni/2 + 1, Ni
initial_state(:, i) = [0.125_R_P, & ! rho(s)
0._R_P, & ! u
0.1_R_P, & ! p
0.125_R_P, & ! sum(rho(s))
cp0/cv0] ! gamma = cp/cv
enddo
case('smooth')
print "(A)", 'Solving "'//trim(adjustl(problem))//'" problem'
t_final = 0.1_R_P
CFL = 0.7_R_P
Ns = 1
Nc = Ns + 2
Np = Ns + 4
allocate(initial_state(1:Np, 1:Ni))
allocate(cp0(1:Ns))
allocate(cv0(1:Ns))
variables = 'VARIABLES="x"'
do s=1, Ns
variables = variables//' "rho('//trim(str(s,.true.))//')"'
enddo
variables = variables//' "u" "p" "rho" "gamma"'
BC_L = 'TRA'
BC_R = 'TRA'
cp0(1) = 1040._R_P
cv0(1) = 743._R_P
do i=1, Ni
initial_state(:, i) = [1._R_P + 4._R_P / 5._R_P * sin( pi * xcenter(i) * 0.5_R_P) + &
1._R_P / 10._R_P * sin(5._R_P * pi * xcenter(i) * 0.5_R_P), & ! rho(s)
0.5_R_P * (xcenter(i) - 0.5_R_P)**4, & ! u
10._R_P + 2._R_P * xcenter(i)**4, & ! p
1._R_P + 4._R_P / 5._R_P * sin( pi * xcenter(i) * 0.5_R_P) + &
1._R_P / 10._R_P * sin(5._R_P * pi * xcenter(i) * 0.5_R_P), & ! sum(rho(s))
cp0/cv0] ! gamma = cp/cv
enddo
case default
print "(A)", 'Error: unknown problem "'//trim(adjustl(problem))//'"'
print "(A)", 'Valid problem names are:'
print "(A)", ' + sod'
print "(A)", ' + smooth'
stop
endselect
if (trim(adjustl(output_cli))/='unset') then
output = trim(adjustl(output_cli))//'-'//trim(adjustl(solver))
else
output = 'euler_1D_integration-'//trim(adjustl(solver))
endif
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine init