subroutine NumericalFlux3 use basicmod, only: is,ie,in,js,je,jn,ks,ke,kn use fluxmod implicit none 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(k, Pleftc1, Pleftc2, Plefte,Prigtc1, Prigtc2, Prigte,dsv,dsvp,dsvm,etot,ptl,css,cts,leftco,rigtco,nflux) do j=js,je do i=is,ie do k=ks,ke+1 Pleftc1(:) = svc(:,i,j,k-2) Pleftc2(:) = svc(:,i,j,k-1) Prigtc1(:) = svc(:,i,j,k ) Prigtc2(:) = svc(:,i,j,k+1) ! | 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(muvv)=Plefte(nve1)*Plefte(nden) leftco(muvw)=Plefte(nve2)*Plefte(nden) ! rho v leftco(muvu)=Plefte(nve3)*Plefte(nden) etot =Plefte(nene)*Plefte(nden) & ! thermal & +0.5d0*Plefte(nden)*( & ! kinetic & +Plefte(nve1)**2 & & +Plefte(nve2)**2 & & +Plefte(nve3)**2) & & +0.5d0* ( & ! magnetic & +Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2) leftco(muet) = etot leftco(mubv)=Plefte(nbm1) ! b_x leftco(mubw)=Plefte(nbm2) ! b_y leftco(mubu)=Plefte(nbm3) ! b_z leftco(mubp)=Plefte(nbps) ! psi leftco(must:mued)=Plefte(nden)*Plefte(nst:ned) ptl = Plefte(npre) + ( Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2)/2.0d0 leftco(mfdn)=Plefte(nden) *Plefte(nve3) ! rho v leftco(mfvv)=Plefte(nden)*Plefte(nve1)*Plefte(nve3) & & -Plefte(nbm1)*Plefte(nbm3) leftco(mfvw)=Plefte(nden)*Plefte(nve2)*Plefte(nve3) & & -Plefte(nbm2)*Plefte(nbm3) leftco(mfvu)=Plefte(nden)*Plefte(nve3)*Plefte(nve3) & & +ptl -Plefte(nbm3)*Plefte(nbm3) leftco(mfet)=(etot+ptl )*Plefte(nve3) & & -( Plefte(nbm1)*Plefte(nve1) & & +Plefte(nbm2)*Plefte(nve2) & & +Plefte(nbm3)*Plefte(nve3))*Plefte(nbm3) leftco(mfbv) = Plefte(nbm1)*Plefte(nve3) & & -Plefte(nve1)*Plefte(nbm3) leftco(mfbw) = Plefte(nbm2)*Plefte(nve3) & & -Plefte(nve2)*Plefte(nbm3) leftco(mfbu) = 0.d0 leftco(mfbp) = 0.0d0 ! psi leftco(mfst:mfed)=Plefte(nden)*Plefte(nst:ned)*Plefte(nve3) 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(nbm3)**2 & & /Plefte(nden)) & & )/2.0d0) leftco(mvel)= Plefte(nve3) 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(muvv)=Prigte(nve1)*Prigte(nden) rigtco(muvw)=Prigte(nve2)*Prigte(nden) rigtco(muvu)=Prigte(nve3)*Prigte(nden) etot =Prigte(nene)*Prigte(nden) & ! internal & +0.5d0*Prigte(nden)*( & ! kinetic & +Prigte(nve1)**2 & & +Prigte(nve2)**2 & & +Prigte(nve3)**2) & & +0.5d0* ( & ! magnetic & +Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2) rigtco(muet) = etot rigtco(mubv)=Prigte(nbm1) ! b_x rigtco(mubw)=Prigte(nbm2) ! b_y rigtco(mubu)=Prigte(nbm3) ! b_z rigtco(mubp)=Prigte(nbps) ! psi rigtco(must:mued)=Prigte(nden)*Prigte(nst:ned) ptl = Prigte(npre) + ( Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2)/2.0d0 rigtco(mfdn)=Prigte(nden) *Prigte(nve3) rigtco(mfvv)=Prigte(nden)*Prigte(nve1)*Prigte(nve3) & & -Prigte(nbm1)*Prigte(nbm3) rigtco(mfvw)=Prigte(nden)*Prigte(nve2)*Prigte(nve3) & & -Prigte(nbm2)*Prigte(nbm3) rigtco(mfvu)=Prigte(nden)*Prigte(nve3)*Prigte(nve3) & & +ptl -Prigte(nbm3)*Prigte(nbm3) rigtco(mfet)=(etot+ptl )*Prigte(nve3) & & -( Prigte(nbm1)*Prigte(nve1) & & +Prigte(nbm2)*Prigte(nve2) & & +Prigte(nbm3)*Prigte(nve3))*Prigte(nbm3) rigtco(mfbv) = Prigte(nbm1)*Prigte(nve3) & & -Prigte(nve1)*Prigte(nbm3) rigtco(mfbw) = Prigte(nbm2)*Prigte(nve3) & & -Prigte(nve2)*Prigte(nbm3) rigtco(mfbu) = 0.0d0 rigtco(mfbp) = 0.0d0 ! psi rigtco(mfst:mfed)=Prigte(nden)*Prigte(nst:ned)*Prigte(nve3) 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(nbm3)**2 & & /Prigte(nden)) & & )/2.0d0) rigtco(mvel)= Prigte(nve3) rigtco(mpre)= ptl ! call HLLE(leftco,rigtco,nflux) ! call HLLC(leftco,rigtco,nflux) call HLLD(leftco,rigtco,nflux) nflux3(mden,i,j,k)=nflux(mden) nflux3(mrv1,i,j,k)=nflux(mrvv) nflux3(mrv2,i,j,k)=nflux(mrvw) nflux3(mrv3,i,j,k)=nflux(mrvu) nflux3(meto,i,j,k)=nflux(meto) nflux3(mbm1,i,j,k)=nflux(mbmv) nflux3(mbm2,i,j,k)=nflux(mbmw) nflux3(mbm3,i,j,k)=nflux(mbmu) nflux3(mst:med,i,j,k)=nflux(mst:med) nflux3(mbm3,i,j,k) = 0.5d0*(leftco(mubp)+rigtco(mubp)) & & -0.5d0*chg*(rigtco(mubu)-leftco(mubu)) ! finite volume nflux3(mbps,i,j,k) = (0.5d0*(leftco(mubu)+rigtco(mubu)) & & -0.5d0/chg*(rigtco(mubp)-leftco(mubp)))*chg**2 ! finite volume enddo enddo enddo !$acc end kernels return end subroutine Numericalflux3