!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