!2017.6.9
module common_module
use sktk_fortran
use setup_template,only:vax,bax,sax,eq_state,eforce,eheat,heat_conductivity,ohm,iflux_damp
implicit none
character(20),parameter :: outdir="/nwork/sakaue/test1/"
integer n,rkstep
real(DP) dx,dy,dz
integer,parameter :: xbox=256
integer,parameter :: xmar(2)=(/ 10,10 /)
integer,parameter :: ybox=1
integer,parameter :: ymar(2)=(/ 10,10 /)*0
integer,parameter :: zbox=1
integer,parameter :: zmar(2)=(/ 10,10 /)*0
integer,parameter :: tlen=500
integer,parameter :: xlen=xbox+xmar(1)+xmar(2)
integer,parameter :: ylen=ybox+ymar(1)+ymar(2)
integer,parameter :: zlen=zbox+zmar(1)+zmar(2)
integer,parameter :: xr(2)=(/xmar(1)+1,xmar(1)+xbox/)
integer,parameter :: yr(2)=(/ymar(1)+1,ymar(1)+ybox/)
integer,parameter :: zr(2)=(/zmar(1)+1,zmar(1)+zbox/)
integer,parameter :: udim(3)=1-int(0._DP**(/xbox-1,ybox-1,zbox-1/))
integer :: udir(sum(udim))
real(DP),parameter :: gamma=5._DP/3._DP,gamma1i=1.0_DP/(gamma-1.0_DP)
real(DP),parameter :: pressure_inf=p_yocto
!Scale of physical quantities
real(DP),parameter :: scV=1._DP,scL=1._DP,scRho=1._DP
! [cm/s], [cm], [g/cm3]
real(DP),parameter :: scTmp=scV*scV/c_Rg ![K]
real(DP),parameter :: scPre=scRho*scV*scV
real(DP),parameter :: scMag=sqrt(4._DP*pi*scRho)*scV
real(DP),dimension(xlen,ylen,zlen,8)::U
integer,parameter :: taxlen=10
real(DP) :: tax(taxlen+1)=0._DP
!-------------------------------------------------------------!
! Coordinate System ------------------------------------------!
real(DP) :: xax(xlen),yax(ylen),zax(zlen)
real(DP),allocatable,dimension(:,:,:) :: cax1,cax2,cax3
real(DP),allocatable,dimension(:,:,:) ::hx,hy,hz
!-------------------------------------------------------------!
real(DP),parameter :: cflf=0.6_DP
integer :: HLLDerror=0,INTERPerror=0,pressure_blurred=-1
integer,parameter :: heat_conduction=0
integer,parameter :: miyoshi_kusano_2011=1
integer,parameter :: powell_et_al_1999=1
integer,parameter :: initial_flux_cancel=0
real(DP),allocatable,dimension(:,:,:,:) :: initial_flux
integer,parameter :: torder=3
!-------------------------------!
! torder !
! 1...Euler !
! 2...TVD Runge-Kutta 2nd !
! 3...TVD Runge-Kutta 3rd !
! 4...TVD Runge-Kutta 4th !
!-------------------------------!
contains
function outputtime(tout)
implicit none
real(DP) outputtime,dt_out,tout
integer,parameter :: nskip=3!.. xbox=100 10!...xbox=2048
dt_out=0.01_DP
! dt_out=0.0001_DP
outputtime=tout+dt_out*ge0(tax(mod(n,taxlen)+1)-tout)
end function outputtime
function U2eth(U8)
implicit none
real(DP) :: U8(8),U2eth
U2eth=U8(Sax(2))-.5_DP*(sum(U8(Vax)*U8(Vax))/U8(Sax(1))&
+sum(U8(Bax)*U8(Bax)))
end function U2eth
function U2p(U8,gamma1i)
implicit none
real(DP),intent(IN) :: gamma1i,U8(8)
real(DP) :: U2p
U2p=U2eth(U8)/gamma1i
end function U2p
function CT_preprocess(U0,xi,yj,zk)
implicit none
real(DP),dimension(:,:,:,:),intent(IN)::U0
integer,intent(IN) :: xi,yj,zk
real(DP),dimension(8) :: CT_preprocess
integer s,kdlt(3)
CT_preprocess=U0(xi,yj,zk,:)
do s=1,3
kdlt=0**abs(s-(/1,2,3/))
CT_preprocess(Bax(s))=.5_DP*(CT_preprocess(Bax(s))&
+U0(max(xi-kdlt(1),1),max(yj-kdlt(2),1),max(zk-kdlt(3),1),Bax(s)))
end do
end function CT_preprocess
function cnsrvtv2prmtv(conservative,gamma1i)
implicit none
real(DP),intent(IN) :: gamma1i,conservative(8)
real(DP) :: cnsrvtv2prmtv(8),Eth,Ekin,Emag
Ekin=.5_DP*sum(conservative(Vax)*conservative(Vax))/conservative(Sax(1))
Emag=.5_DP*sum(conservative(Bax)*conservative(Bax))
Eth=conservative(Sax(2))-Ekin-Emag
cnsrvtv2prmtv(Sax(1))=conservative(Sax(1))
cnsrvtv2prmtv(Sax(2))=max(Eth/gamma1i,pressure_inf)
cnsrvtv2prmtv(Vax)=conservative(Vax)/conservative(Sax(1))
cnsrvtv2prmtv(Bax)=conservative(Bax)
end function cnsrvtv2prmtv
function prmtv2cnsrvtv(primitive,gamma1i)
implicit none
real(DP),intent(IN) :: gamma1i,primitive(8)
real(DP) :: prmtv2cnsrvtv(8),Eth,Ekin,Emag
Eth=max(primitive(Sax(2)),pressure_inf)*gamma1i
Ekin=.5_DP*sum(primitive(Vax)*primitive(Vax))*primitive(Sax(1))
Emag=.5_DP*sum(primitive(Bax)*primitive(Bax))
prmtv2cnsrvtv(Sax(1))=primitive(Sax(1))
prmtv2cnsrvtv(Sax(2))=Eth+Ekin+Emag
prmtv2cnsrvtv(Vax)=primitive(Vax)*primitive(Sax(1))
prmtv2cnsrvtv(Bax)=primitive(Bax)
end function prmtv2cnsrvtv
function scale_factors(xi,yj,zk)
implicit none
real(DP) :: scale_factors(3)
integer :: xi,yj,zk,shh(3)
shh=shape(hx)
scale_factors&
=(/hx(min(max(xi,1),shh(1)),min(max(yj,1),shh(2)),&
min(max(zk,1),shh(3))),&
hy(min(max(xi,1),shh(1)),min(max(yj,1),shh(2)),&
min(max(zk,1),shh(3))),&
hz(min(max(xi,1),shh(1)),min(max(yj,1),shh(2)),&
min(max(zk,1),shh(3)))/)
end function scale_factors
subroutine cell_center(xi,yj,zk,pc,hc)
implicit none
integer,intent(IN) :: xi,yj,zk
real(DP),intent(OUT) :: pc(3),hc(3)
hc=scale_factors(xi,yj,zk)
pc=(/xax(min(max(xi,1),xlen)),yax(min(max(yj,1),ylen)),&
zax(min(max(zk,1),zlen))/)
end subroutine cell_center
subroutine cell_face(xi,yj,zk,axis,dir,pf,hf,sf)
implicit none
integer,intent(IN) :: xi,yj,zk,axis,dir
real(DP),intent(OUT) :: pf(3),hf(3),sf
integer :: kdlt(3),kdlt2(3),kdlt3(3),axis2,axis3,cnt
real(DP) :: hc(3),hR(3),hL(3),pc(3),pR(3),pL(3),sfR,sfL
axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
kdlt=0**abs(axis-(/1,2,3/))
call cell_center(xi,yj,zk,pc,hc)
if (udim(axis).eq.0) then
pf=pc; hf=hc
sf=2._DP*hc(axis2)*hc(axis3)*hf(axis2)*hf(axis3)&
/(hc(axis2)*hc(axis3)+hf(axis2)*hf(axis3))
return
end if
call cell_center(xi+kdlt(1),yj+kdlt(2),zk+kdlt(3),pR,hR)
call cell_center(xi-kdlt(1),yj-kdlt(2),zk-kdlt(3),pL,hL)
do cnt=1,3
pf(cnt)=Lag_interpol2(-1._DP,0._DP,1._DP,pL(cnt),pc(cnt),pR(cnt),&
.5_DP*dir)
end do
do cnt=1,3
hf(cnt)=Lag_interpol2(pL(axis),pc(axis),pR(axis),hL(cnt),hc(cnt),&
hR(cnt),pf(axis))
end do
sfR=2._DP*hc(axis2)*hc(axis3)*hR(axis2)*hR(axis3)&
/(hc(axis2)*hc(axis3)+hR(axis2)*hR(axis3))
sfL=2._DP*hc(axis2)*hc(axis3)*hL(axis2)*hL(axis3)&
/(hc(axis2)*hc(axis3)+hL(axis2)*hL(axis3))
sf=sfR*eq0(dir-1._DP)+sfL*eq0(dir+1._DP)
end subroutine cell_face
subroutine cell_edge(xi,yj,zk,axis1,axis2,pe,he)
implicit none
integer,intent(IN) :: xi,yj,zk,axis1,axis2
real(DP),intent(OUT) :: pe(3),he(3)
integer :: kdlt1(3),kdlt2(3),kdlt3(3),axis3
real(DP) :: hc(3),hp1(3),hp2(3),hp3(3)
kdlt1=0**abs(axis1-(/1,2,3/))
kdlt2=0**abs(axis2-(/1,2,3/))
kdlt3=1-kdlt1-kdlt2
pe(1)=(xax(min(max(xi+1-kdlt3(1),1),xlen))+xax(min(max(xi,1),xlen)))*.5_DP
pe(2)=(yax(min(max(yj+1-kdlt3(2),1),ylen))+yax(min(max(yj,1),ylen)))*.5_DP
pe(3)=(zax(min(max(zk+1-kdlt3(3),1),zlen))+zax(min(max(zk,1),zlen)))*.5_DP
hc=scale_factors(xi,yj,zk)
hp1=scale_factors(xi+kdlt1(1),yj+kdlt1(2),zk+kdlt1(3))
hp2=scale_factors(xi+kdlt2(1),yj+kdlt2(2),zk+kdlt2(3))
hp3=scale_factors(xi+1-kdlt3(1),yj+1-kdlt3(2),zk+1-kdlt3(3))
he(1)=(hc(1)+hp1(1)+hp2(1)+hp3(1))*.25_DP
he(2)=(hc(2)+hp1(2)+hp2(2)+hp3(2))*.25_DP
he(3)=(hc(3)+hp1(3)+hp2(3)+hp3(3))*.25_DP
end subroutine cell_edge
function divBf(vec,xi,yj,zk)
! vec is defined at the cell face
implicit none
real(DP) :: divBf
real(DP),dimension(xlen,ylen,zlen,3),intent(IN):: vec
integer,intent(IN) :: xi,yj,zk
integer :: s,axis,kdlt(3)
real(DP) :: pfR(3),pfL(3),hfR(3),hfL(3),hc(3)
real(DP) :: sfR,sfL,dBxdx,dBydy,dBzdz
hc=scale_factors(xi,yj,zk)
divBf=0._DP
do s=1,sum(udim)
axis=udir(s); kdlt=0**abs(axis-(/1,2,3/))
call cell_face(xi,yj,zk,axis, 1,pfR,hfR,sfR)
call cell_face(xi,yj,zk,axis,-1,pfL,hfL,sfL)
divBf=divBf+(sfR*vec(xi,yj,zk,axis)&
-sfL*vec(max(xi-kdlt(1),1),max(yj-kdlt(2),1),&
max(zk-kdlt(3),1),axis))&
/(pfR(axis)-pfL(axis))
end do
divBf=divBf/hc(1)/hc(2)/hc(3)
end function divBf
function rotAe(vec,xi,yj,zk,axis)
!vec is defined at the cell edge (right-up)
!rotAe is defined at the cell face
implicit none
real(DP) :: rotAe
real(DP),dimension(xlen,ylen,zlen,3),intent(IN):: vec
integer,intent(IN) :: xi,yj,zk,axis
real(DP) :: pfR(3),hfR(3),sfR,denomi
real(DP),dimension(3) :: peR,heR,peL,heL
integer :: axis2,axis3,kdlt2(3),kdlt3(3)
axis2=mod(axis,3)+1; axis3=mod(axis2,3)+1
kdlt2=0**abs(axis2-(/1,2,3/)); kdlt3=0**abs(axis3-(/1,2,3/))
call cell_face(xi,yj,zk,axis,1,pfR,hfR,sfR)
call cell_edge(xi,yj,zk,axis,axis2,peR,heR)
call cell_edge(xi-kdlt2(1),yj-kdlt2(2),zk-kdlt2(3),axis,axis2,peL,heL)
denomi=peR(axis2)-peL(axis2); denomi=ne0(denomi)/(denomi+eq0(denomi))
rotAe=(heR(axis3)*vec(xi,yj,zk,axis3)&
-heL(axis3)*vec(xi-kdlt2(1),yj-kdlt2(2),zk-kdlt2(3),axis3))*denomi
call cell_edge(xi,yj,zk,axis,axis3,peR,heR)
call cell_edge(xi-kdlt3(1),yj-kdlt3(2),zk-kdlt3(3),axis,axis3,peL,heL)
denomi=peR(axis3)-peL(axis3); denomi=ne0(denomi)/(denomi+eq0(denomi))
rotAe=rotAe-(heR(axis2)*vec(xi,yj,zk,axis2)&
-heL(axis2)*vec(xi-kdlt3(1),yj-kdlt3(2),zk-kdlt3(3),axis2))*denomi
rotAe=rotAe/sfR
end function rotAe
function rotAf(vec,xi,yj,zk,axis)
!vec is defined at the cell edge (right)
!rotAf is defined at the cell center
implicit none
real(DP) :: rotAf
real(DP),dimension(xlen,ylen,zlen,3),intent(IN):: vec
integer,intent(IN) :: xi,yj,zk,axis
real(DP) :: pc(3),hc(3),pfR(3),hfR(3),sfR,pfL(3),hfL(3),denomi
integer :: axis2,axis3,kdlt2(3),kdlt3(3)
axis2=mod(axis,3)+1; axis3=mod(axis2,3)+1
kdlt2=0**abs(axis2-(/1,2,3/)); kdlt3=0**abs(axis3-(/1,2,3/))
call cell_center(xi,yj,zk,pc,hc)
call cell_face(xi,yj,zk,axis2, 1,pfR,hfR,sfR)
call cell_face(xi,yj,zk,axis2,-1,pfL,hfL,sfR)
denomi=pfR(axis2)-pfL(axis2); denomi=ne0(denomi)/(denomi+eq0(denomi))
rotAf=(hfR(axis3)*vec(xi,yj,zk,axis3)&
-hfL(axis3)*vec(xi-kdlt2(1),yj-kdlt2(2),zk-kdlt2(3),axis3))*denomi
call cell_face(xi,yj,zk,axis3, 1,pfR,hfR,sfR)
call cell_face(xi,yj,zk,axis3,-1,pfL,hfL,sfR)
denomi=pfR(axis3)-pfL(axis3); denomi=ne0(denomi)/(denomi+eq0(denomi))
rotAf=rotAf-(hfR(axis2)*vec(xi,yj,zk,axis2)&
-hfL(axis2)*vec(xi-kdlt3(1),yj-kdlt3(2),zk-kdlt3(3),axis2))*denomi
rotAf=rotAf/(hc(axis2)*hc(axis3))
end function rotAf
end module common_module
module initial_module
use common_module
contains
subroutine ic
implicit none
real(DP),dimension(xlen,ylen,zlen)::rho,pre,ene
real(DP),dimension(xlen,ylen,zlen,3)::vel,mag
real(DP):: b0=0._DP
integer :: i,j,k
integer :: Brio_Wu_1988=1
integer :: Ryu_Jones_1995
j=1
do i=1,3
if (udim(i).eq.1) then
udir(j)=i; j=j+1
end if
end do
allocate(hx(1,1,1)) !allocate(hx(xlen,ylen,zlen))
allocate(hy(1,1,1)) !allocate(hy(xlen,ylen,zlen))
allocate(hz(1,1,1)) !allocate(hz(xlen,ylen,zlen))
allocate(cax1(1,1,1)) !allocate(cax1(xlen,ylen,zlen))
allocate(cax2(1,1,1)) !allocate(cax2(xlen,ylen,zlen))
allocate(cax3(1,1,1)) !allocate(cax3(xlen,ylen,zlen))
if (initial_flux_cancel .ge. 1) then
allocate(initial_flux(xlen,ylen,zlen,8))
else
allocate(initial_flux(1,1,1,1))
end if
Ryu_Jones_1995=1-Brio_Wu_1988
!Orthogonal Curvilinear Coordinate System ++++++++++++++++++++
hx=1._DP
hy=1._DP
hz=1._DP
cax1=1._DP
cax2=1._DP
cax3=1._DP
!-------------------------------------------------------------
!axis coordinate +++++++++++++++++++++++++++++++++++++++++++++
do i=1,xlen
xax(i)=1._DP/dble(xbox)*dble(i-xr(1))-.5_DP
end do
! non-uniform grid !---------------------------------------
! do i=3,xlen
! xax(i)=xax(i-1)+(xax(i-1)-xax(i-2))*1.01_DP
! end do
!----------------------------------------------------------
dx=minval(abs(xax(2:xlen)-xax(1:xlen-1)))
do j=1,ylen; yax(j)=1._DP/dble(ybox)*dble(j-yr(1))-.5_DP
enddo
dy=minval(abs(yax(2:ylen)-yax(1:ylen-1)))
do k=1,zlen; zax(k)=1._DP/dble(xbox)*dble(k-zr(1))
enddo
!-------------------------------------------------------------
if (Brio_Wu_1988 .eq. 1) then
do i=1,xlen
if (xax(i) .ge. 0) then
rho(i,:,:)=0.125_DP
pre(i,:,:)=0.1_DP
vel(i,:,:,:)=0._DP
mag(i,:,:,1)=.75_DP*b0
mag(i,:,:,2)=-1._DP*b0
else
rho(i,:,:)=1.0_DP
pre(i,:,:)=1.0_DP
vel(i,:,:,:)=0._DP
mag(i,:,:,1)=.75_DP*b0
mag(i,:,:,2)=1._DP*b0
end if
mag(i,:,:,3)=0._DP
enddo
end if
if (Ryu_Jones_1995 .eq. 1) then
do i=1,xlen
if (xax(i) .le. 0) then
rho(i,:,:)=1.08_DP
vel(i,:,:,1)=1.2_DP
vel(i,:,:,2)=0.01_DP
vel(i,:,:,3)=0.5_DP
mag(i,:,:,1)=2._DP
mag(i,:,:,2)=3.6_DP
mag(i,:,:,3)=2._DP
pre(i,:,:)=0.95_DP/gamma1i
if (pre(i,1,1) .le. 0._DP) then
print*,i,pre(i,1,1)
stop
end if
else
rho(i,:,:)=1.0_DP
vel(i,:,:,:)=0._DP
mag(i,:,:,1)=2._DP
mag(i,:,:,2)=4._DP
mag(i,:,:,3)=2._DP
pre(i,:,:)=1._DP/gamma1i
end if
enddo
end if
do i=1,xlen; do j=1,ylen; do k=1,zlen
U(i,j,k,Bax)=mag(i,j,k,:)
U(i,j,k,Vax)=vel(i,j,k,:)*rho(i,j,k)
U(i,j,k,Sax(1))=rho(i,j,k)
end do; end do; end do
do i=1,xlen; do j=1,ylen; do k=1,zlen
U(i,j,k,Sax(2))=pre(i,j,k)*gamma1i&
+sum(U(i,j,k,Vax)**2)/U(i,j,k,Sax(1))/2._DP&
+(((U(i,j,k,Bax(1))+U(max(i-1,1),j,k,Bax(1)))/2._DP)**2&
+((U(i,j,k,Bax(2))+U(i,max(j-1,1),k,Bax(2)))/2._DP)**2&
+((U(i,j,k,Bax(3))+U(i,j,max(k-1,1),Bax(3)))/2._DP)**2)/2._DP
end do; end do; end do
!-------------------------------!
! bc_select !
! 1...mirror boundary !
! 2...anti-mirror boundary !
! 3...periodic boundary !
! 4...derichret free boundary !
! 5...No consideration !
!-------------------------------!
end subroutine ic
end module initial_module