module ct_module use common_module implicit none contains subroutine CT(U0,Fxyz,SLRs) !$use omp_lib implicit none real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0 real(DP),dimension(xlen,ylen,zlen,8,3),intent(INOUT)::Fxyz real(DP),dimension(xlen,ylen,zlen,3,2),intent(IN)::SLRs real(DP),dimension(xlen,ylen,zlen,3)::Ee real(DP) :: Ec(4),Ef(8),Mc(4,2),Mf(12),slr(4,2),denomi ! e...cell edge; c...cell center; f...cell face integer axis,axis2,axis3,kdlt(3),kdlt2(3),kdlt3(3) integer axR(3),axL(3),axb(3),axv(3),f,dir integer :: ir(2),jr(2),kr(2),cnt,szi,szj,szk,s(3) real(DP) :: UR(8),UL(8) ! real(DP) :: vnR2,BnR2,vnL2,BnL2,pR,pL,CsR,CsL,CaR,CaL,CfR,CfL,SR,SL real(DP) :: Ohm_e,Ohm_f2,Ohm_f3,B_f2,B_f3 !------------------------! ! f8 f7 ! ! f9 c4 f3 c3 f6 ! axis2 ! f4 e f2 ! A ! f10 c1 f1 c2 f5 ! | ! f11 f12 ! --> axis !------------------------! Ee=0._DP !$omp parallel shared(U0,Fxyz,SLRs,Ee),default(firstprivate) do axis=1,3 axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1 if (max(udim(axis),udim(axis2)).eq.0) cycle kdlt=0**abs(axis-(/1,2,3/)) ; kdlt2=0**abs(axis2-(/1,2,3/)) kdlt3=(/1,1,1/)-kdlt-kdlt2 axv=(/Vax(axis),Vax(axis2),Vax(axis3)/) axb=(/Bax(axis),Bax(axis2),Bax(axis3)/) ir=xr-(/1-kdlt3(1),0/) jr=yr-(/1-kdlt3(2),0/) kr=zr-(/1-kdlt3(3),0/) ir=xr-(/1-kdlt3(1),0/)*udim(1) jr=yr-(/1-kdlt3(2),0/)*udim(2) kr=zr-(/1-kdlt3(3),0/)*udim(3) szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1 !$omp barrier !$omp do do cnt=1,szi*szj*szk call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s) if (0**(s(3)-xmar(1))+0**(s(2)-ymar(1))+0**(s(1)-zmar(1)) .ge. 2) & cycle Ef(1)=-Fxyz(s(3),s(2),s(1),axb(2),axis) Ef(2)=Fxyz(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axb(1),axis2) Ef(3)=-Fxyz(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),axb(2),axis) Ef(4)=Fxyz(s(3),s(2),s(1),axb(1),axis2) Ee(s(3),s(2),s(1),axis3)=0.25_DP*(Ef(1)+Ef(2)+Ef(3)+Ef(4)) if (miyoshi_kusano_2011 .ne. 1) cycle !------------------------! ! f8 f7 ! ! f9 c4 f3 c3 f6 ! axis2 ! f4 e f2 ! A ! f10 c1 f1 c2 f5 ! | ! f11 f12 ! --> axis !------------------------! Ef(5)=-Fxyz(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axb(2),axis) Ef(6)=-Fxyz(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),axb(2),axis) Ef(7)=Fxyz(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),axb(1),axis2) Ef(8)=Fxyz(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),axb(1),axis2) Mf(1)=U0(s(3),s(2),s(1),axb(1)) Mf(2)=U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axb(2)) Mf(3)=U0(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),axb(1)) Mf(4)=U0(s(3),s(2),s(1),axb(2)) Mf(5)=U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axb(1)) Mf(6)=U0(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),axb(1)) Mf(7)=U0(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),axb(2)) Mf(8)=U0(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),axb(2)) Mf(9)=U0(min(max(s(3)-kdlt(1)+kdlt2(1),1),xlen),& min(max(s(2)-kdlt(2)+kdlt2(2),1),ylen),& min(max(s(1)-kdlt(3)+kdlt2(3),1),zlen),axb(1)) Mf(10)=U0(max(s(3)-kdlt(1),1),max(s(2)-kdlt(2),1),& max(s(1)-kdlt(3),1),axb(1)) Mf(11)=U0(max(s(3)-kdlt2(1),1),max(s(2)-kdlt2(2),1),& max(s(1)-kdlt2(3),1),axb(2)) Mf(12)=U0(min(max(s(3)+kdlt(1)-kdlt2(1),1),xlen),& min(max(s(2)+kdlt(2)-kdlt2(2),1),ylen),& min(max(s(1)+kdlt(3)-kdlt2(3),1),zlen),axb(2)) Mc(1,:)=(/Mf(1)+Mf(10),Mf(4)+Mf(11)/)*0.5_DP Mc(2,:)=(/Mf(1)+Mf(5),Mf(2)+Mf(12)/)*0.5_DP Mc(3,:)=(/Mf(3)+Mf(6),Mf(2)+Mf(7)/)*0.5_DP Mc(4,:)=(/Mf(3)+Mf(9),Mf(4)+Mf(8)/)*0.5_DP Ec(1)=-(U0(s(3),s(2),s(1),axv(1))*Mc(1,2)& -U0(s(3),s(2),s(1),axv(2))*Mc(1,1))/U0(s(3),s(2),s(1),Sax(1)) Ec(2)=-(U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axv(1))*Mc(2,2)& -U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axv(2))*Mc(2,1))& /U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),Sax(1)) Ec(3)=-(U0(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),axv(1))*Mc(3,2)& -U0(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),axv(2))*Mc(3,1))& /U0(min(s(3)+1-kdlt3(1),xlen),min(s(2)+1-kdlt3(2),ylen),& min(s(1)+1-kdlt3(3),zlen),Sax(1)) Ec(4)=-(U0(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),axv(1))*Mc(4,2)& -U0(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),axv(2))*Mc(4,1))& /U0(min(s(3)+kdlt2(1),xlen),min(s(2)+kdlt2(2),ylen),& min(s(1)+kdlt2(3),zlen),Sax(1)) do f=1,4 !------------------------! ! f8 f7 ! ! f9 c4 f3 c3 f6 ! axis2 ! f4 e f2 ! A ! f10 c1 f1 c2 f5 ! | ! f11 f12 ! --> axis !------------------------! !f=1...R=c2, L=c1, dir=1 !f=2...R=c3, L=c2, dir=2 !f=3...R=c3, L=c4, dir=1 !f=4...R=c4, L=c1, dir=2 axL=kdlt*(0**abs(f-2))+kdlt2*(0**abs(f-3)) dir=axis*0**abs(1-mod(f,2))+axis2*0**abs(mod(f,2)) slr(f,:)=SLRs(min(s(3)+axL(1),xlen),min(s(2)+axL(2),ylen),& min(s(1)+axL(3),zlen),dir,:) end do !------------------------! ! f8 f7 ! ! f9 c4 f3 c3 f6 ! axis2 ! f4 e f2 ! A ! f10 c1 f1 c2 f5 ! | ! f11 f12 ! --> axis !------------------------! !f=1...R=c2, L=c1, dir=1 !f=2...R=c3, L=c2, dir=2 !f=3...R=c3, L=c4, dir=1 !f=4...R=c4, L=c1, dir=2 denomi=slr(4,2)-slr(4,1) denomi=.25_DP*ne0(denomi)/(denomi+eq0(denomi)) Ee(s(3),s(2),s(1),axis3)=Ee(s(3),s(2),s(1),axis3)+& (slr(4,2)*(Ef(1)-Ec(1))-slr(4,1)*(Ef(3)-Ec(4))& +slr(4,1)*slr(4,2)*(Mf(3)-Mf(1)-Mc(4,1)+Mc(1,1)))*denomi denomi=slr(2,2)-slr(2,1) denomi=.25_DP*ne0(denomi)/(denomi+eq0(denomi)) Ee(s(3),s(2),s(1),axis3)=Ee(s(3),s(2),s(1),axis3)-& (slr(2,2)*(Ef(5)-Ec(2))-slr(2,1)*(Ef(6)-Ec(3))& +slr(2,1)*slr(2,2)*(Mf(6)-Mf(5)-Mc(3,1)+Mc(2,1)))*denomi denomi=slr(1,2)-slr(1,1) denomi=.25_DP*ne0(denomi)/(denomi+eq0(denomi)) Ee(s(3),s(2),s(1),axis3)=Ee(s(3),s(2),s(1),axis3)+& (slr(1,2)*(Ef(4)-Ec(1))-slr(1,1)*(Ef(2)-Ec(2))& +slr(1,1)*slr(1,2)*(-Mf(2)+Mf(4)+Mc(2,2)-Mc(1,2)))*denomi denomi=slr(3,2)-slr(3,1) denomi=.25_DP*ne0(denomi)/(denomi+eq0(denomi)) Ee(s(3),s(2),s(1),axis3)=Ee(s(3),s(2),s(1),axis3)-& (slr(3,2)*(Ef(8)-Ec(4))-slr(3,1)*(Ef(7)-Ec(3))& +slr(3,1)*slr(3,2)*(-Mf(7)+Mf(8)+Mc(3,2)-Mc(4,2)))*denomi end do !$omp end do end do !Ohm's law ++++++++++++++++++++++++++++++++++++++ do axis=1,3 axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1 kdlt=0**abs(axis-(/1,2,3/)) ; kdlt2=0**abs(axis2-(/1,2,3/)) kdlt3=(/1,1,1/)-kdlt-kdlt2 axv=(/Vax(axis),Vax(axis2),Vax(axis3)/) axb=(/Bax(axis),Bax(axis2),Bax(axis3)/) ir=xr+(/-1,0/)*udim(1); jr=yr+(/-1,0/)*udim(2); kr=zr+(/-1,0/)*udim(3) szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1 !$omp do do cnt=1,szi*szj*szk call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s) Ohm_e=Ohm(s(3),s(2),s(1),axis3,U0) Ee(s(3),s(2),s(1),axis3)=Ee(s(3),s(2),s(1),axis3)+Ohm_e Ohm_f2=(Ohm(s(3),s(2),s(1),axis2,U0)+& Ohm(max(s(3)-kdlt3(1),1),max(s(2)-kdlt3(2),1),& max(s(1)-kdlt3(3),1),axis2,U0))*.5_DP Ohm_f3=(Ohm(s(3),s(2),s(1),axis3,U0)+& Ohm(max(s(3)-kdlt2(1),1),max(s(2)-kdlt2(2),1),& max(s(1)-kdlt2(3),1),axis3,U0))*.5_DP B_f2=(U0(s(3),s(2),s(1),axb(2))& +U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axb(2))& +U0(min(max(s(3)+kdlt(1)-kdlt2(1),1),xlen),& min(max(s(2)+kdlt(2)-kdlt2(2),1),ylen),& min(max(s(1)+kdlt(3)-kdlt2(3),1),zlen),axb(2))& +U0(max(s(3)-kdlt2(1),1),max(s(2)-kdlt2(2),1),& max(s(1)-kdlt2(3),1),axb(2)))*.25_DP B_f3=(U0(s(3),s(2),s(1),axb(3))& +U0(min(s(3)+kdlt(1),xlen),min(s(2)+kdlt(2),ylen),& min(s(1)+kdlt(3),zlen),axb(3))& +U0(min(max(s(3)+kdlt(1)-kdlt3(1),1),xlen),& min(max(s(2)+kdlt(2)-kdlt3(2),1),ylen),& min(max(s(1)+kdlt(3)-kdlt3(3),1),zlen),axb(3))& +U0(max(s(3)-kdlt3(1),1),max(s(2)-kdlt3(2),1),& max(s(1)-kdlt3(3),1),axb(3)))*.25_DP Fxyz(s(3),s(2),s(1),Sax(2),axis)=Fxyz(s(3),s(2),s(1),Sax(2),axis)& +Ohm_f2*B_f3-Ohm_f3*B_f2 end do !$omp end do end do !------------------------------------------------ do axis=1,3 if (udim(axis).eq.0) cycle axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1 axb=(/Bax(axis),Bax(axis2),Bax(axis3)/) kdlt=0**abs(axis-(/1,2,3/)) ir=xr+(/-kdlt(1),0/) jr=yr+(/-kdlt(2),0/) kr=zr+(/-kdlt(3),0/) ir=xr+(/-kdlt(1),0/)*udim(1) jr=yr+(/-kdlt(2),0/)*udim(2) kr=zr+(/-kdlt(3),0/)*udim(3) szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1 !$omp do do cnt=1,szi*szj*szk call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s) Fxyz(s(3),s(2),s(1),axb(2),axis)=-Ee(s(3),s(2),s(1),axis3) Fxyz(s(3),s(2),s(1),axb(3),axis)=Ee(s(3),s(2),s(1),axis2) end do !$omp end do end do !$omp end parallel 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 CT end module CT_MODULE