module heat_conduction_module use boundary_module implicit none integer,private,parameter :: xm=(xmar(1)-1)*udim(1),ym=(ymar(1)-1)*udim(2),& zm=(zmar(1)-1)*udim(3) integer,private,parameter :: nx=xbox+2*xm,ny=ybox+2*ym,nz=zbox+2*zm integer,private,parameter :: nxnynz=nx*ny*nz real(DP),private,dimension(nxnynz) :: Tmp_o,Tmp_n,Tmp_n1 real(DP),private :: dt contains subroutine heat_conduction_RKL2(U0,dU,t0,dt0) implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN) :: U0 real(DP),dimension(xlen,ylen,zlen,8),intent(OUT) :: dU real(DP),dimension(xlen,ylen,zlen) :: Fc real(DP),dimension(xlen,ylen,zlen) :: divFc,divFc_0 real(DP),intent(IN) :: t0,dt0 real(DP),allocatable,dimension(:,:,:,:) :: dYs real(DP) :: dtH_dtP,t_step real(DP),allocatable :: cb(:),cmu(:),cnu(:),cmu_(:),cgm_(:) integer :: nsteps,cnt,cnt1,axis dtH_dtP=dt0/CFL_parabolic(U0) nsteps=ceiling(.25_DP*(1._DP+sqrt(1._DP+.125_DP*(1._DP+2._DP*dtH_dtP)))) nsteps=max(nsteps,3) nsteps=nsteps+nsteps-1 dU=0._DP; divFc_0=div_conduction_flux(U0) rkstep=-1 dU(:,:,:,Sax(2))=divFc_0*4._DP*dt0/3._DP/dble(nsteps**2+nsteps-2) allocate(dYs(xlen,ylen,zlen,2)) allocate(cb(nsteps+1)); allocate(cmu(nsteps)); allocate(cnu(nsteps)) allocate(cmu_(nsteps)); allocate(cgm_(nsteps)) !coefficients prepared ******************************************* do cnt=1,nsteps+1 cb(cnt)=max(.5_DP*(1._DP-2._DP/dble(cnt-1)/(dble(cnt-1)+1._DP)),& 1._DP/3._DP) end do do cnt=2,nsteps cmu(cnt)=(2._DP-1._DP/dble(cnt+1))*cb(cnt+1)/cb(cnt) cnu(cnt)=(1._DP/dble(cnt+1)-1._DP)*cb(cnt+1)/cb(cnt-1) cmu_(cnt)=(8._DP*dble(cnt)-4._DP)*cb(cnt+1)/cb(cnt)& /dble(cnt*(nsteps*nsteps+nsteps-2)) cgm_(cnt)=(cb(cnt)-1._DP)*cmu_(cnt) end do !***************************************************************** dYs=0._DP; dYs(:,:,:,2)=dU(:,:,:,Sax(2)) !$omp parallel do cnt=2,nsteps rkstep=-cnt t_step=t0+dt0*dble(cnt*cnt+cnt-2)/dble(nsteps*nsteps+nsteps-2) divFc=div_conduction_flux(bc(U0,dU,t_step)) !$omp workshare dU(:,:,:,Sax(2))=cmu(cnt)*dYs(:,:,:,2)+cnu(cnt)*dYs(:,:,:,1)& +cmu_(cnt)*dt0*divFc+cgm_(cnt)*dt0*divFc_0 dYs(:,:,:,1)=dYs(:,:,:,2); dYs(:,:,:,2)=dU(:,:,:,Sax(2)) !$omp end workshare end do !$omp end parallel end subroutine heat_conduction_RKL2 function div_conduction_flux(U0) implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN) :: U0 real(DP),dimension(xlen,ylen,zlen) :: Fc,div_conduction_flux real(DP) :: hc(3),pfR(3),hfR(3),sfR,pfL(3),hfL(3),sfL,dcell_f integer :: axis,kdlt(3),cnt,cnt1,szi,szj,szk,s(3) szi=xr(2)-xr(1)+1; szj=yr(2)-yr(1)+1; szk=zr(2)-zr(1)+1 do cnt1=1,sum(udim) axis=udir(cnt1); kdlt=0**abs(axis-(/1,2,3/)) call conduction_flux(U0,Fc,axis) !$omp parallel do private(s,hc,pfR,hfR,sfR,pfL,hfL,sfL,dcell_f) do cnt=1,szi*szj*szk call loopcounter(cnt,(/szk,szj,szi/),(/zr(1),yr(1),xr(1)/),s) hc=scale_factors(s(3),s(2),s(1)) call cell_face(s(3),s(2),s(1),axis, 1,pfR,hfR,sfR) call cell_face(s(3),s(2),s(1),axis,-1,pfL,hfL,sfL) dcell_f=pfR(axis)-pfL(axis) div_conduction_flux(s(3),s(2),s(1))& =div_conduction_flux(s(3),s(2),s(1))*ge0(dble(axis-2))& +(sfR*Fc(s(3),s(2),s(1))& -sfL*Fc(s(3)-kdlt(1),s(2)-kdlt(2),s(1)-kdlt(3)))& /dcell_f/hc(1)/hc(2)/hc(3) end do !$omp end parallel do end do contains subroutine loopcounter(count,sizes,org,subs) implicit none integer,intent(IN) :: count,sizes(:),org(:) integer,intent(OUT) :: subs(:) integer :: dim,s,quo(size(subs)) dim=size(sizes) quo=count do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1 end do do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s) end do end subroutine loopcounter end function div_conduction_flux subroutine conduction_flux(U0,Fc,axis) implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN) :: U0 real(DP),dimension(xlen,ylen,zlen),intent(OUT) :: Fc integer,intent(IN) :: axis integer :: axis2,axis3,axes(3),axes2(3),axes3(3) integer :: ir(2),jr(2),kr(2) integer :: cnt,szi,szj,szk,s(3),sR(3) real(DP) :: Uc8(8),UR8(8) real(DP) :: eth_c,eth_R,temp_c,temp_R,eth_c2(2),temp_c2(2) real(DP) :: eth_R2(2),temp_R2(2),eth_c3(2),temp_c3(2),eth_R3(2),temp_R3(2) real(DP) :: temp_f,eth_f,mag_f(3),kappa_f(3,3) real(DP) :: dTc(3),Fc_dTc(3),maxFc(3),case real(DP) :: pc(3),hc(3),pL(3),pR(3) real(DP) :: dcell_e2,dcell_e3,dcell_c,dcell_f real(DP) :: denomi real(DP) :: pf(3),hf(3),sf Fc=0._DP if (udim(axis).eq.0) return axes=0**abs(axis-(/1,2,3/)) axis2=mod(axis,3)+1; axis3=mod(axis2,3)+1 axes=0**abs((/1,2,3/)-axis); axes2=0**abs((/1,2,3/)-axis2) axes3=(/1,1,1/)-axes-axes2 ir=xr+((/-1,1/)+(/-1,0/)*axes(1))*udim(1) jr=yr+((/-1,1/)+(/-1,0/)*axes(2))*udim(2) kr=zr+((/-1,1/)+(/-1,0/)*axes(3))*udim(3) szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1 !$omp parallel do private(s,pc,hc,pL,pR,dcell_e2,dcell_e3,dcell_c,denomi,& !$omp pf,hf,sf,sR,Uc8,UR8,eth_c,eth_R,temp_c,temp_R,eth_c2,temp_c2,eth_R2,& !$omp temp_R2,eth_c3,temp_c3,eth_R3,temp_R3,temp_f,eth_f,mag_f,kappa_f,& !$omp dTc,Fc_dTc,maxFc,case) do cnt=1,szi*szj*szk call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s) call cell_center(s(3),s(2),s(1),pc,hc) call cell_center(s(3)-axes2(1),s(2)-axes2(2),s(1)-axes2(3),pL,hc) call cell_center(s(3)+axes2(1),s(2)+axes2(2),s(1)+axes2(3),pR,hc) dcell_e2=(pR(axis2)-pL(axis2))*.5_DP call cell_center(s(3)-axes3(1),s(2)-axes3(2),s(1)-axes3(3),pL,hc) call cell_center(s(3)+axes3(1),s(2)+axes3(2),s(1)+axes3(3),pR,hc) dcell_e3=(pR(axis3)-pL(axis3))*.5_DP call cell_center(s(3)+axes(1),s(2)+axes(2),s(1)+axes(3),pR,hc) dcell_c=pR(axis)-pc(axis) call cell_face(s(3),s(2),s(1),axis, 1,pf,hf,sf) !calculation kappa_f +++++++++++++++++++++++++++++++++++++++++++++++++++ sR=s+(/axes(3),axes(2),axes(1)/) Uc8=CT_preprocess(U0,s(3),s(2),s(1)) UR8=CT_preprocess(U0,sR(3),sR(2),sR(1)) eth_c=max(U2eth(Uc8),pressure_inf*gamma1i) eth_R=max(U2eth(UR8),pressure_inf*gamma1i) temp_c=eq_state(eth_c/gamma1i,Uc8(Sax(1)),-1._DP) temp_R=eq_state(eth_R/gamma1i,UR8(Sax(1)),-1._DP) eth_c2(1)=U2eth(U0(max(s(3)-axes2(1),1),max(s(2)-axes2(2),1),& max(s(1)-axes2(3),1),:)) eth_c2(2)=U2eth(U0(min(s(3)+axes2(1),xlen),min(s(2)+axes2(2),ylen),& min(s(1)+axes2(3),zlen),:)) temp_c2(1)=eq_state(eth_c2(1)/gamma1i,& U0(max(s(3)-axes2(1),1),max(s(2)-axes2(2),1),& max(s(1)-axes2(3),1),Sax(1)),-1._DP) temp_c2(2)=eq_state(eth_c2(1)/gamma1i,& U0(min(s(3)+axes2(1),xlen),min(s(2)+axes2(2),ylen),& min(s(1)+axes2(3),zlen),Sax(1)),-1._DP) eth_R2(1)=U2eth(U0(max(sR(3)-axes2(1),1),max(sR(2)-axes2(2),1),& max(sR(1)-axes2(3),1),:)) eth_R2(2)=U2eth(U0(min(sR(3)+axes2(1),xlen),min(sR(2)+axes2(2),ylen),& min(sR(1)+axes2(3),zlen),:)) temp_R2(1)=eq_state(eth_R2(1)/gamma1i,& U0(max(sR(3)-axes2(1),1),max(sR(2)-axes2(2),1),& max(sR(1)-axes2(3),1),Sax(1)),-1._DP) temp_R2(2)=eq_state(eth_R2(2)/gamma1i,& U0(min(sR(3)+axes2(1),xlen),min(sR(2)+axes2(2),ylen),& min(sR(1)+axes2(3),zlen),Sax(1)),-1._DP) eth_c3(1)=U2eth(U0(max(s(3)-axes3(1),1),max(s(2)-axes3(2),1),& max(s(1)-axes3(3),1),:)) eth_c3(2)=U2eth(U0(min(s(3)+axes3(1),xlen),min(s(2)+axes3(2),ylen),& min(s(1)+axes3(3),zlen),:)) temp_c3(1)=eq_state(eth_c3(1)/gamma1i,& U0(max(s(3)-axes3(1),1),max(s(2)-axes3(2),1),& max(s(1)-axes3(3),1),Sax(1)),-1._DP) temp_c3(2)=eq_state(eth_c3(1)/gamma1i,& U0(min(s(3)+axes3(1),xlen),min(s(2)+axes3(2),ylen),& min(s(1)+axes3(3),zlen),Sax(1)),-1._DP) eth_R3(1)=U2eth(U0(max(sR(3)-axes3(1),1),max(sR(2)-axes3(2),1),& max(sR(1)-axes3(3),1),:)) eth_R3(2)=U2eth(U0(min(sR(3)+axes3(1),xlen),min(sR(2)+axes3(2),ylen),& min(sR(1)+axes3(3),zlen),:)) temp_R3(1)=eq_state(eth_R3(1)/gamma1i,& U0(max(sR(3)-axes3(1),1),max(sR(2)-axes3(2),1),& max(sR(1)-axes3(3),1),Sax(1)),-1._DP) temp_R3(2)=eq_state(eth_R3(2)/gamma1i,& U0(min(sR(3)+axes3(1),xlen),min(sR(2)+axes3(2),ylen),& min(sR(1)+axes3(3),zlen),Sax(1)),-1._DP) temp_f=(temp_c+temp_R)*.5_DP eth_f=(eth_c+eth_R)*.5_DP mag_f=(UR8(Bax)+UC8(Bax))*.5_DP mag_f(axis)=U0(s(3),s(2),s(1),Bax(axis)) kappa_f=conductivity_mtrx(mag_f,temp_f) !calculation kappa_fR,kappa_fL dTc(axis)=temp_R-temp_c dTc(axis2)=(temp_c2(2)-temp_c2(1)+temp_R2(2)-temp_R2(1))*.25_DP dTc(axis3)=(temp_c3(2)-temp_c3(1)+temp_R3(2)-temp_R3(1))*.25_DP denomi=hf(axis)*dcell_c; denomi=ne0(denomi)/(denomi+eq0(denomi)) Fc_dTc(axis)=kappa_f(axis,axis)*denomi denomi=hf(axis2)*dcell_e2; denomi=ne0(denomi)/(denomi+eq0(denomi)) Fc_dTc(axis2)=kappa_f(axis,axis2)*denomi denomi=hf(axis3)*dcell_e3; denomi=ne0(denomi)/(denomi+eq0(denomi)) Fc_dTc(axis3)=kappa_f(axis,axis3)*denomi !Conductive flux limitter +++++++++++++++++++++++++++++++++++ maxFc=saturated_flux(temp_f,eth_f,dTc,mag_f) case=sqrt(sum((Fc_dTc*dTc)**2)/sum(maxFc**2)) case=max(0._DP,min(1._DP,1._DP-.5_DP*case)) !Conductive flux limitter Fc(s(3),s(2),s(1))=sum(Fc_dTc*dTc*case+maxFc*(1._DP-case)) end do !$omp end parallel do contains subroutine loopcounter(count,sizes,org,subs) implicit none integer,intent(IN) :: count,sizes(:),org(:) integer,intent(OUT) :: subs(:) integer :: dim,s,quo(size(subs)) dim=size(sizes) quo=count do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1 end do do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s) end do end subroutine loopcounter end subroutine conduction_flux function CFL_parabolic(U0) implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0 real(DP) :: CFL_parabolic real(DP) :: V8(8),kappa_mtrx(3,3) real(DP) :: pc(3),hc(3),dcell_f(sum(udim)) real(DP) :: pfR(3),pfL(3),hfR(3),hfL(3),sfR,sfL real(DP) :: eps=p_yocto*p_yocto real(DP) :: dt_a(sum(udim)) integer :: ir(2),jr(2),kr(2),szi,szj,szk,s(3),kdlt(3),cnt,cnt1 integer :: axis CFL_parabolic=p_giga ir=xr+(/-1,1/)*udim(1); jr=yr+(/-1,1/)*udim(2) kr=zr+(/-1,1/)*udim(3) szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1 !$omp parallel do private(s,pc,hc,V8,kappa_mtrx,cnt1,axis,pfR,hfR,sfR,pfL,hfL,sfL,dcell_f,kdlt,dt_a),reduction(min:CFL_parabolic) do cnt=1,szi*szj*szk call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s) call cell_center(s(3),s(2),s(1),pc,hc) V8=cnsrvtv2prmtv(CT_preprocess(U0,s(3),s(2),s(1)),gamma1i) kappa_mtrx=conductivity_mtrx(V8(Bax),& eq_state(V8(Sax(2)),V8(Sax(1)),-1._DP)) do cnt1=0,8 kappa_mtrx(mod(cnt1,3)+1,cnt1/3+1)& =max(abs(kappa_mtrx(mod(cnt1,3)+1,cnt1/3+1)),p_yocto) end do do cnt1=1,sum(udim) axis=udir(cnt1) call cell_face(s(3),s(2),s(1),axis, 1,pfR,hfR,sfR) call cell_face(s(3),s(2),s(1),axis,-1,pfL,hfL,sfL) dcell_f(cnt1)=hc(axis)*abs(pfR(axis)-pfL(axis)) end do do cnt1=1,sum(udim) axis=udir(cnt1) dt_a(cnt1)=1._DP/sum(2._DP*kappa_mtrx(axis,udir)/dcell_f/dcell_f) end do CFL_parabolic=min(CFL_parabolic,1._DP/sum(1._DP/dt_a)) end do !$omp end parallel do contains subroutine loopcounter(count,sizes,org,subs) implicit none integer,intent(IN) :: count,sizes(:),org(:) integer,intent(OUT) :: subs(:) integer :: dim,s,quo(size(subs)) dim=size(sizes) quo=count do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1 end do do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s) end do end subroutine loopcounter end function CFL_parabolic subroutine heat_conduction_sor(U0,dU,dt0) !$ use omp_lib implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0 real(DP),dimension(xlen,ylen,zlen,8),intent(OUT)::dU real(DP),intent(IN) :: dt0 real(DP) :: psum,pre1,rho1,tmp1 integer :: cnt,cnt1,s(3) !cnt=i+nx*(j-1)+nx*ny*(k-1) integer :: iteration_number real(DP) :: iteration_thr=p_mili real(DP),parameter :: omega=1.5_DP dt=dt0 dU=0._DP !$omp parallel do private(s,pre1,rho1) do cnt=1,nxnynz call loopcounter(cnt,(/nz,ny,nx/),(/1,1,1/),s) pre1=U2p(CT_preprocess(U0,s(3)-xm+xmar(1),s(2)-ym+ymar(1),& s(1)-zm+zmar(1)),gamma1i) pre1=max(pre1,pressure_inf) rho1=U0(s(3)-xm+xmar(1),s(2)-ym+ymar(1),s(1)-zm+zmar(1),Sax(1)) Tmp_o(s(3)+nx*(s(2)-1)+nx*ny*(s(1)-1))=eq_state(pre1,rho1,-1._DP) end do !$omp end parallel do Tmp_n=Tmp_o iteration_number=0 do psum=0._DP !$omp parallel !$omp do reduction(+:psum) do cnt=1,nxnynz psum=psum+sor_coeff_mtrx(1,cnt,U0,dt0)*Tmp_n(cnt) end do !$omp end do Tmp_n1(1)=Tmp_n(1)+omega*(Tmp_o(1)-psum)/sor_coeff_mtrx(1,1,U0,dt0) !$omp do private(psum,cnt1) do cnt=2,nxnynz psum=0._DP do cnt1=cnt,nxnynz psum=psum+sor_coeff_mtrx(cnt,cnt1,U0,dt0)*Tmp_n(cnt1) end do Tmp_n1(cnt)=0._DP do cnt1=1,cnt-1 Tmp_n1(cnt)=Tmp_n1(cnt)& +sor_coeff_mtrx(cnt,cnt1,U0,dt0)*Tmp_n1(cnt1) end do Tmp_n1(cnt)=Tmp_n(cnt)& +omega*(Tmp_o(cnt)-psum-Tmp_n1(cnt))& /sor_coeff_mtrx(cnt,cnt,U0,dt0) end do !$omp end do !$omp end parallel if (maxval(abs(Tmp_n1/Tmp_n-1._DP)) .le. iteration_thr) exit if (iteration_number .ge. 10) then print*,'% Iteration over! @ heat_conduction_module.f90',n,& maxval(abs(Tmp_n1/Tmp_n-1._DP)) exit end if Tmp_n=Tmp_n1 iteration_number=iteration_number+1 end do !$omp parallel do private(s) do cnt=1,xbox*ybox*zbox call loopcounter(cnt,(/zbox,ybox,xbox/),(/1,1,1/),s) dU(s(3)+xmar(1),s(2)+ymar(1),s(1)+zmar(1),Sax(2))=& (eq_state(-1._DP,U0(s(3)+xmar(1),s(2)+ymar(1),s(1)+zmar(1),Sax(1)),& Tmp_n1(s(3)+xm+nx*(s(2)+ym-1)+nx*ny*(s(1)+zm-1)))-& eq_state(-1._DP,U0(s(3)+xmar(1),s(2)+ymar(1),s(1)+zmar(1),Sax(1)),& Tmp_o(s(3)+xm+nx*(s(2)+ym-1)+nx*ny*(s(1)+zm-1))))*gamma1i end do !$omp end parallel do contains subroutine loopcounter(count,sizes,org,subs) implicit none integer,intent(IN) :: count,sizes(:),org(:) integer,intent(OUT) :: subs(:) integer :: dim,s,quo(size(subs)) dim=size(sizes) quo=count do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1 end do do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s) end do end subroutine loopcounter end subroutine heat_conduction_sor function conductivity_mtrx(mag,tmp) implicit none real(DP),intent(IN) :: mag(3),tmp real(DP) :: denomi,nonmag real(DP) :: conductivity_mtrx(3,3) integer s1,s2 denomi=sum(mag*mag); nonmag=le0(denomi-p_yocto*p_yocto) denomi=denomi*(1._DP-nonmag)+nonmag do s1=1,3; do s2=1,3 conductivity_mtrx(s1,s2)=0._DP**(abs(s1-s2))*nonmag& +(1._DP-nonmag)*mag(s1)*mag(s2)/denomi end do; end do conductivity_mtrx=heat_conductivity(tmp)*conductivity_mtrx end function conductivity_mtrx function saturated_flux(tmp,Eth,gradT,mag) implicit none real(DP),intent(IN) :: tmp,Eth,gradT(3),mag(3) real(DP) :: denomi,nonmag,umag(3) real(DP) :: saturated_flux(3) integer s1,s2 denomi=sum(mag*mag); nonmag=le0(denomi-p_yocto*p_yocto) denomi=denomi*(1._DP-nonmag)+nonmag umag=(mag+nonmag)/denomi saturated_flux=sign(1._DP,sum(umag*gradT))*umag saturated_flux=saturated_flux*Eth*sqrt(gamma*tmp*scTmp*c_kb/c_me)/scV end function saturated_flux function sor_coeff_mtrx(mu,nu,U0,dt0) implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0 real(DP),intent(IN) :: dt0 integer,intent(IN) :: mu,nu integer :: axis,axis2,axis3,axes(3),axes2(3),axes3(3),s,s1 integer :: mu_i,mu_j,mu_k,mu1i,mu1j,mu1k,nu_ijk(3) real(DP) :: UR8(8),UC8(8),UL8(8) real(DP) :: pc(3),pR(3),pL(3),pfR(3),pfL(3),hc(3),hfR(3),hfL(3),sfR,sfL real(DP) :: eth_c(3),temp_c(3),dTR(3),dTL(3),eth_f(2),temp_f(2) real(DP) :: mag_fR(3),mag_fL(3),kappa_fR(3,3),kappa_fL(3,3) real(DP) :: temp_c2(3,2),temp_c3(3,2) real(DP) :: FR_dT(3),FL_dT(3),hhFR_dT(3),hhFL_dT(3) real(DP) :: maxFR(3),maxFL(3),maxFR_dT(3),maxFL_dT(3),F_Fsat,case real(DP) :: alpha(3) real(DP) :: dcell_cR,dcell_cL,dcell_e2,dcell_e3,dcell_f real(DP) :: sor_coeff_mtrx mu_i=mod(mu-1,nx)+1; mu_j=mod((mu-1)/nx,ny)+1; mu_k=(mu-1)/nx/ny+1 nu_ijk=(/mod(nu-1,nx)+1-mu_i,mod((nu-1)/nx,ny)+1-mu_j,(nu-1)/nx/ny+1-mu_k/) mu1i=mu_i-xm+xmar(1); mu1j=mu_j-ym+ymar(1); mu1k=mu_k-zm+zmar(1) if (sum(abs(nu_ijk)) .ge. 3 .or. maxval(abs(nu_ijk)) .ge. 2) then sor_coeff_mtrx=0._DP; return end if alpha=0._DP call cell_center(mu1i,mu1j,mu1k,pc,hc) do s1=1,sum(udim) axis=udir(s1); axis2=mod(axis,3)+1; axis3=mod(axis2,3)+1 axes=0**abs((/1,2,3/)-axis); axes2=0**abs((/1,2,3/)-axis2) axes3=(/1,1,1/)-axes-axes2 call cell_center(mu1i-axes2(1),mu1j-axes2(2),mu1k-axes2(3),pL,hc) call cell_center(mu1i+axes2(1),mu1j+axes2(2),mu1k+axes2(3),pR,hc) dcell_e2=(pR(axis2)-pL(axis2))*.5_DP call cell_center(mu1i-axes3(1),mu1j-axes3(2),mu1k-axes3(3),pL,hc) call cell_center(mu1i+axes3(1),mu1j+axes3(2),mu1k+axes3(3),pR,hc) dcell_e3=(pR(axis3)-pL(axis3))*.5_DP call cell_center(mu1i-axes(1),mu1j-axes(2),mu1k-axes(3),pL,hc) call cell_center(mu1i+axes(1),mu1j+axes(2),mu1k+axes(3),pR,hc) dcell_cR=pR(axis)-pc(axis); dcell_cL=pc(axis)-pL(axis) call cell_face(mu1i,mu1j,mu1k,axis, 1,pfR,hfR,sfR) call cell_face(mu1i,mu1j,mu1k,axis,-1,pfL,hfL,sfL) dcell_f=pfR(axis)-pfL(axis) !calculation kappa_fR, kappa_fL ++++++++++++++++++++++++++++++++++++++++ do s=-1,1 eth_c(2+s)=U2eth(CT_preprocess(U0,mu1i+s*axes(1),mu1j+s*axes(2),& mu1k+s*axes(3))) temp_c(2+s)=temperature_array(mu1i+s*axes(1),mu1j+s*axes(2),& mu1k+s*axes(3),Tmp_o) temp_c2(2+s,1)=temperature_array(mu1i+s*axes(1)-axes2(1),& mu1j+s*axes(2)-axes2(2),& mu1k+s*axes(3)-axes2(3),Tmp_o) temp_c2(2+s,2)=temperature_array(mu1i+s*axes(1)+axes2(1),& mu1j+s*axes(2)+axes2(2),& mu1k+s*axes(3)+axes2(3),Tmp_o) temp_c3(2+s,1)=temperature_array(mu1i+s*axes(1)-axes3(1),& mu1j+s*axes(2)-axes3(2),& mu1k+s*axes(3)-axes3(3),Tmp_o) temp_c3(2+s,2)=temperature_array(mu1i+s*axes(1)+axes3(1),& mu1j+s*axes(2)+axes3(2),& mu1k+s*axes(3)+axes3(3),Tmp_o) end do temp_f(1)=Lag_interpol2(pL(axis),pc(axis),pR(axis),& temp_c(1),temp_c(2),temp_c(3),pfL(axis)) temp_f(2)=Lag_interpol2(pL(axis),pc(axis),pR(axis),& temp_c(1),temp_c(2),temp_c(3),pfR(axis)) eth_f(1)=Lag_interpol2(pL(axis),pc(axis),pR(axis),& eth_c(1),eth_c(2),eth_c(3),pfL(axis)) eth_f(2)=Lag_interpol2(pL(axis),pc(axis),pR(axis),& eth_c(1),eth_c(2),eth_c(3),pfR(axis)) do s=1,2 temp_f(s)=max(min(temp_f(s),maxval(temp_c)),minval(temp_c)) eth_f(s)=max(min(eth_f(s),maxval(eth_c)),minval(eth_c)) end do UL8=CT_preprocess(U0,mu1i-axes(1),mu1j-axes(2),mu1k-axes(3)) UC8=CT_preprocess(U0,mu1i,mu1j,mu1k) UR8=CT_preprocess(U0,mu1i+axes(1),mu1j+axes(2),mu1k+axes(3)) mag_fL=(UL8(Bax)+UC8(Bax))*.5_DP mag_fL(axis)=U0(mu1i-axes(1),mu1j-axes(2),mu1k-axes(3),Bax(axis)) mag_fR=(UR8(Bax)+UC8(Bax))*.5_DP mag_fR(axis)=U0(mu1i,mu1j,mu1k,Bax(axis)) kappa_fL=conductivity_mtrx(mag_fL,temp_f(1)) kappa_fR=conductivity_mtrx(mag_fR,temp_f(2)) !calculation kappa_fR,kappa_fL dTL(axis)=temp_c(2)-temp_c(1) dTL(axis2)=(temp_c2(1,2)-temp_c2(1,1)+temp_c2(2,2)-temp_c2(2,1))*.25_DP dTL(axis3)=(temp_c3(1,2)-temp_c3(1,1)+temp_c3(2,2)-temp_c3(2,1))*.25_DP dTR(axis)=temp_c(3)-temp_c(2) dTR(axis2)=(temp_c2(2,2)-temp_c2(2,1)+temp_c2(3,2)-temp_c2(3,1))*.25_DP dTR(axis3)=(temp_c3(2,2)-temp_c3(2,1)+temp_c3(3,2)-temp_c3(3,1))*.25_DP FL_dT(axis)=kappa_fL(axis,axis)/(hfL(axis)*dcell_cL) FL_dT(axis2)=kappa_fL(axis,axis2)/(hfL(axis2)*dcell_e2) FL_dT(axis3)=kappa_fL(axis,axis3)/(hfL(axis3)*dcell_e3) FR_dT(axis)=kappa_fR(axis,axis)/(hfR(axis)*dcell_cR) FR_dT(axis2)=kappa_fR(axis,axis2)/(hfR(axis2)*dcell_e2) FR_dT(axis3)=kappa_fR(axis,axis3)/(hfR(axis3)*dcell_e3) !Conductive flux limiter +++++++++++++++++++++++++++++++++++ maxFL=saturated_flux(temp_f(1),eth_f(1),dTL,mag_fL) maxFL_dT=maxFL*sign(1._DP,dTL)/(/max(abs(dTL(1)),p_yocto),& max(abs(dTL(2)),p_yocto),& max(abs(dTL(3)),p_yocto)/) case=sqrt(sum(FL_dT**2)/sum(maxFL_dT**2)) case=max(0._DP,min(1._DP,1._DP-.5_DP*case)) FL_dT=FL_dT*case+maxFL_dT*(1._DP-case) maxFR=saturated_flux(temp_f(2),eth_f(2),dTR,mag_fR) maxFR_dT=maxFR*sign(1._DP,dTR)/(/max(abs(dTR(1)),p_yocto),& max(abs(dTR(2)),p_yocto),& max(abs(dTR(3)),p_yocto)/) case=sqrt(sum(FR_dT**2)/sum(maxFR_dT**2)) case=max(0._DP,min(1._DP,1._DP-.5_DP*case)) FR_dT=FR_dT*case+maxFR_dT*(1._DP-case) !Conductive flux limiter hhFR_dT=sfR*FR_dT; hhFL_dT=sfL*FL_dT alpha(axis)=0._DP alpha(axis)=alpha(axis)+(hhFR_dT(axis)+hhFL_dT(axis))& *(0._DP**sum(abs(nu_ijk))) ! 1. 000 alpha(axis)=alpha(axis)-hhFR_dT(axis)& *(0._DP**sum(abs(nu_ijk-axes))) ! 2. +00 alpha(axis)=alpha(axis)-(hhFR_dT(axis2)-hhFL_dT(axis2))& *(0._DP**sum(abs(nu_ijk-axes2))) ! 3. 0+0 alpha(axis)=alpha(axis)-(hhFR_dT(axis3)-hhFL_dT(axis3))& *(0._DP**sum(abs(nu_ijk-axes3))) ! 4. 00+ alpha(axis)=alpha(axis)-hhFL_dT(axis)& *(0._DP**sum(abs(nu_ijk+axes))) ! 5. -00 alpha(axis)=alpha(axis)+(hhFR_dT(axis2)-hhFL_dT(axis2))& *(0._DP**sum(abs(nu_ijk+axes2))) ! 6. 0-0 alpha(axis)=alpha(axis)+(hhFR_dT(axis3)-hhFL_dT(axis3))& *(0._DP**sum(abs(nu_ijk+axes3))) ! 7. 00- alpha(axis)=alpha(axis)-0._DP& *(0._DP**sum(abs(nu_ijk-axes2-axes3))) ! 8. 0++ alpha(axis)=alpha(axis)-hhFR_dT(axis3)& *(0._DP**sum(abs(nu_ijk-axes3-axes))) ! 9. +0+ alpha(axis)=alpha(axis)-hhFR_dT(axis2)& *(0._DP**sum(abs(nu_ijk-axes-axes2))) !10. ++0 alpha(axis)=alpha(axis)-0._DP& *(0._DP**sum(abs(nu_ijk+axes2+axes3))) !11. 0-- alpha(axis)=alpha(axis)-hhFL_dT(axis3)& *(0._DP**sum(abs(nu_ijk+axes3+axes))) !12. -0- alpha(axis)=alpha(axis)-hhFL_dT(axis2)& *(0._DP**sum(abs(nu_ijk+axes+axes2))) !13. --0 alpha(axis)=alpha(axis)-0._DP& *(0._DP**sum(abs(nu_ijk+axes2-axes3))) !14. 0-+ alpha(axis)=alpha(axis)+hhFR_dT(axis3)& *(0._DP**sum(abs(nu_ijk+axes3-axes))) !15. +0- alpha(axis)=alpha(axis)+hhFL_dT(axis2)& *(0._DP**sum(abs(nu_ijk+axes-axes2))) !16. -+0 alpha(axis)=alpha(axis)-0._DP& *(0._DP**sum(abs(nu_ijk-axes2+axes3))) !17. 0+- alpha(axis)=alpha(axis)+hhFL_dT(axis3)& *(0._DP**sum(abs(nu_ijk-axes3+axes))) !18. -0+ alpha(axis)=alpha(axis)+hhFR_dT(axis2)& *(0._DP**sum(abs(nu_ijk-axes+axes2))) !19. +-0 alpha(axis)=alpha(axis)/dcell_f end do sor_coeff_mtrx=0._DP**abs(mu-nu)& +dt0*sum(alpha)/(hc(1)*hc(2)*hc(3)*gamma1i*U0(mu1i,mu1j,mu1k,Sax(1))) end function sor_coeff_mtrx function temperature_array(xi,yj,zk,Tarr) implicit none real(DP),intent(IN) :: Tarr(nxnynz) integer,intent(IN) :: xi,yj,zk real(DP) :: temperature_array temperature_array=Tarr((xi-xmar(1)+xm)+nx*(yj-ymar(1)+ym-1)& +nx*ny*(zk-zmar(1)+zm-1)) end function temperature_array end module heat_conduction_module