subroutine NumericalFlux2 use basicmod, only: is,ie,in,js,je,jn,ks,ke,kn use fluxmod implicit none ! | Pleftc1 | Pleftc2 | Prigtc1 | Prigtc2 | ! You are here integer::i,j,k real(8),dimension(nhyd):: dsvp,dsvm,dsvc,dsv real(8),dimension(nhyd):: Pleftc1, Pleftc2, Plefte real(8),dimension(nhyd):: Prigtc1, Prigtc2, Prigte real(8),dimension(2*mflx+madd):: leftco,rigtco real(8),dimension(mflx):: nflux real(8):: ptl,css,cts,etot !$acc kernels !$acc loop collapse(2) independent private(j, Pleftc1, Pleftc2, Plefte,Prigtc1, Prigtc2, Prigte,dsv,dsvp,dsvm,etot,ptl,css,cts,leftco,rigtco,nflux) do k=ks,ke do i=is,ie do j=js,je+1 Pleftc1(:) = svc(:,i,j-2,k) Pleftc2(:) = svc(:,i,j-1,k) Prigtc1(:) = svc(:,i,j ,k) Prigtc2(:) = svc(:,i,j+1,k) ! | Pleftc1 | Pleftc2 |<= Prigtc1 | Prigtc2 | ! You are here !==================== ! Left !==================== dsvp(:) = Prigtc1(:) - Pleftc2(:) dsvm(:) = Pleftc2(:) - Pleftc1(:) call vanLeer(dsvp,dsvm,dsv) ! call minmod(dsvp,dsvm,dsv) Plefte(:) = Pleftc2(:) + 0.5d0*dsv(:) ! Consvative variables leftco(mudn)=Plefte(nden) leftco(muvw)=Plefte(nve1)*Plefte(nden) leftco(muvu)=Plefte(nve2)*Plefte(nden) ! rho v leftco(muvv)=Plefte(nve3)*Plefte(nden) etot =Plefte(nene)*Plefte(nden) & ! internal & +0.5d0*Plefte(nden)*( & & +Plefte(nve1)**2 & & +Plefte(nve2)**2 & & +Plefte(nve3)**2) & ! kinetic & +0.5d0* ( & & +Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2) ! magnetic leftco(muet)=etot leftco(mubw)=Plefte(nbm1) ! b_x leftco(mubu)=Plefte(nbm2) ! b_y leftco(mubv)=Plefte(nbm3) ! b_z leftco(mubp)=Plefte(nbps) ! psi leftco(must:mued)=Plefte(nst:ned)*Plefte(nden) ! rho X ptl = Plefte(npre) + ( Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2)/2.0d0 leftco(mfdn)=Plefte(nden) *Plefte(nve2) ! rho v leftco(mfvw)=Plefte(nden)*Plefte(nve1)*Plefte(nve2) & & -Plefte(nbm1)*Plefte(nbm2) leftco(mfvu)=Plefte(nden)*Plefte(nve2)*Plefte(nve2) & & +ptl-Plefte(nbm2)*Plefte(nbm2) leftco(mfvv)=Plefte(nden)*Plefte(nve3)*Plefte(nve2) & & -Plefte(nbm3)*Plefte(nbm2) leftco(mfet)=(etot+ptl)*Plefte(nve2) & & -( Plefte(nbm1)*Plefte(nve1) & & +Plefte(nbm2)*Plefte(nve2) & & +Plefte(nbm3)*Plefte(nve3))*Plefte(nbm2) leftco(mfbw) = Plefte(nbm1)*Plefte(nve2) & & -Plefte(nve1)*Plefte(nbm2) leftco(mfbu) = 0.0d0 leftco(mfbv) = Plefte(nbm3)*Plefte(nve2) & & -Plefte(nve3)*Plefte(nbm2) leftco(mfbp) = 0.0d0 ! psi leftco(mfst:mfed)=Plefte(nst:ned)*Plefte(nden)*Plefte(nve2) ! rho X v css = Plefte(ncsp)**2 cts = css & !c_s^2*c_a^2 & +( Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2)/Plefte(nden) leftco(mcsp)= sqrt((cts +sqrt(cts**2 & & -4.0d0*css*Plefte(nbm2)**2 & & /Plefte(nden)) & & )/2.0d0) leftco(mvel)= Plefte(nve2) leftco(mpre)= ptl ! | Pleftc1 | Pleftc2 =>| Prigtc1 | Prigtc2 | ! You are here !==================== ! Right !==================== dsvp(:) = Prigtc2(:) - Prigtc1(:) dsvm(:) = Prigtc1(:) - Pleftc2(:) call vanLeer(dsvp,dsvm,dsv) ! call minmod(dsvp,dsvm,dsv) Prigte(:) = Prigtc1(:) - 0.5d0*dsv(:) rigtco(mudn)=Prigte(nden) rigtco(muvw)=Prigte(nve1)*Prigte(nden) rigtco(muvu)=Prigte(nve2)*Prigte(nden) ! rho v rigtco(muvv)=Prigte(nve3)*Prigte(nden) etot =Prigte(nene)*Prigte(nden) & ! internal & +0.5d0*Prigte(nden)*( & & +Prigte(nve1)**2 & & +Prigte(nve2)**2 & & +Prigte(nve3)**2) & ! kinetic & +0.5d0* ( & & +Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2) ! magnetic rigtco(muet) = etot rigtco(mubw)=Prigte(nbm1) ! b_x rigtco(mubu)=Prigte(nbm2) ! b_y rigtco(mubv)=Prigte(nbm3) ! b_z rigtco(mubp)=Prigte(nbps) ! psi rigtco(must:mued)=Prigte(nden)*Prigte(nst:ned) ! rho X v ptl = Prigte(npre) + ( Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2)/2.0d0 rigtco(mfdn)=Prigte(nden) *Prigte(nve2) ! rho v rigtco(mfvw)=Prigte(nden)*Prigte(nve1)*Prigte(nve2) & & -Prigte(nbm1)*Prigte(nbm2) rigtco(mfvu)=Prigte(nden)*Prigte(nve2)*Prigte(nve2) & & +ptl-Prigte(nbm2)*Prigte(nbm2) rigtco(mfvv)=Prigte(nden)*Prigte(nve3)*Prigte(nve2) & & -Prigte(nbm3)*Prigte(nbm2) rigtco(mfet)=(etot+ptl)*Prigte(nve2) & & -( Prigte(nbm1)*Prigte(nve1) & & +Prigte(nbm2)*Prigte(nve2) & & +Prigte(nbm3)*Prigte(nve3))*Prigte(nbm2) rigtco(mfbw) = Prigte(nbm1)*Prigte(nve2) & & -Prigte(nve1)*Prigte(nbm2) rigtco(mfbu) = 0.0d0 rigtco(mfbv) = Prigte(nbm3)*Prigte(nve2) & & -Prigte(nve3)*Prigte(nbm2) rigtco(mfbp) = 0.0d0 ! psi rigtco(mfst:mfed)=Prigte(nden)*Prigte(nst:ned)*Prigte(nve2) ! rho X v css = Prigte(ncsp)**2 cts = css & !c_s^2*c_a^2 & +( Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2)/Prigte(nden) rigtco(mcsp)= sqrt((cts +sqrt(cts**2 & & -4.0d0*css*Prigte(nbm2)**2 & & /Prigte(nden)) & & )/2.0d0) rigtco(mvel)= Prigte(nve2) rigtco(mpre)= ptl ! call HLLE(leftco,rigtco,nflux) ! call HLLC(leftco,rigtco,nflux) call HLLD(leftco,rigtco,nflux) nflux2(mden,i,j,k)=nflux(mden) nflux2(mrv1,i,j,k)=nflux(mrvw) nflux2(mrv2,i,j,k)=nflux(mrvu) ! mrv2=3, mrvu=2 nflux2(mrv3,i,j,k)=nflux(mrvv) nflux2(meto,i,j,k)=nflux(meto) nflux2(mbm1,i,j,k)=nflux(mbmw) nflux2(mbm2,i,j,k)=nflux(mbmu) nflux2(mbm3,i,j,k)=nflux(mbmv) nflux2(mst:med,i,j,k)=nflux(mst:med) nflux2(mbm2,i,j,k) = 0.5d0*(leftco(mubp)+rigtco(mubp)) & & -0.5d0*chg*(rigtco(mubu)-leftco(mubu)) ! finite volume nflux2(mbps,i,j,k) = (0.5d0*(leftco(mubu)+rigtco(mubu)) & & -0.5d0/chg*(rigtco(mubp)-leftco(mubp)))*chg**2 ! finite volume ! print *, "bpf2",nflux2(mbps,i,j,k) enddo enddo enddo !$acc end kernels return end subroutine Numericalflux2