subroutine NumericalFlux1 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(i, Pleftc1, Pleftc2, Plefte,Prigtc1, Prigtc2, Prigte,dsv,dsvp,dsvm,etot,ptl,css,cts,leftco,rigtco,nflux) do k=ks,ke do j=js,je do i=is,ie+1 Pleftc1(:) = svc(:,i-2,j,k) Pleftc2(:) = svc(:,i-1,j,k) Prigtc1(:) = svc(:,i ,j,k) Prigtc2(:) = svc(:,i+1,j,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) ! rho leftco(muvu)=Plefte(nve1)*Plefte(nden) ! rho v_x leftco(muvv)=Plefte(nve2)*Plefte(nden) ! rho v_y leftco(muvw)=Plefte(nve3)*Plefte(nden) ! rho v_z etot =Plefte(nene)*Plefte(nden) & ! e_i & +0.5d0*Plefte(nden)*( & & +Plefte(nve1)**2 & & +Plefte(nve2)**2 & & +Plefte(nve3)**2) & ! + rho v^2/2 & +0.5d0* ( & & +Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2) ! + B^2/2 leftco(muet) = etot leftco(mubu)=Plefte(nbm1) ! b_x leftco(mubv)=Plefte(nbm2) ! b_y leftco(mubw)=Plefte(nbm3) ! b_z leftco(mubp)=Plefte(nbps) ! psi leftco(must:mued)=Plefte(nst:ned)*Plefte(nden) ! rho X ! Flux ptl = Plefte(npre) + ( Plefte(nbm1)**2 & & +Plefte(nbm2)**2 & & +Plefte(nbm3)**2)/2.0d0 leftco(mfdn)=Plefte(nden) *Plefte(nve1) leftco(mfvu)=Plefte(nden)*Plefte(nve1)*Plefte(nve1) & & +ptl-Plefte(nbm1)*Plefte(nbm1) leftco(mfvv)=Plefte(nden)*Plefte(nve2)*Plefte(nve1) & & -Plefte(nbm2)*Plefte(nbm1) leftco(mfvw)=Plefte(nden)*Plefte(nve3)*Plefte(nve1) & & -Plefte(nbm3)*Plefte(nbm1) leftco(mfet)= (etot +ptl )*Plefte(nve1) & & -( Plefte(nbm1)*Plefte(nve1) & & +Plefte(nbm2)*Plefte(nve2) & & +Plefte(nbm3)*Plefte(nve3))*Plefte(nbm1) leftco(mfbu) = 0.0d0 leftco(mfbv) = Plefte(nbm2)*Plefte(nve1) & & -Plefte(nve2)*Plefte(nbm1) leftco(mfbw) = Plefte(nbm3)*Plefte(nve1) & & -Plefte(nve3)*Plefte(nbm1) leftco(mfbp) = 0.0d0 ! psi leftco(mfst:mfed)=Plefte(nden)*Plefte(nst:ned)*Plefte(nve1) ! 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(nbm1)**2 & & /Plefte(nden)) & & )/2.0d0) leftco(mvel)= Plefte(nve1) 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(:) ! Consvative variables rigtco(mudn)=Prigte(nden) ! rho rigtco(muvu)=Prigte(nve1)*Prigte(nden) ! rho v_x rigtco(muvv)=Prigte(nve2)*Prigte(nden) ! rho v_y rigtco(muvw)=Prigte(nve3)*Prigte(nden) ! rho v_z etot =Prigte(nene)*Prigte(nden) & ! e_i & +0.5d0*Prigte(nden)*( & & +Prigte(nve1)**2 & & +Prigte(nve2)**2 & & +Prigte(nve3)**2) & ! + rho v^2/2 & +0.5d0* ( & & +Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2) ! + B^2/2 rigtco(muet) = etot rigtco(mubu)=Prigte(nbm1) ! b_x rigtco(mubv)=Prigte(nbm2) ! b_y rigtco(mubw)=Prigte(nbm3) ! b_z rigtco(mubp)=Prigte(nbps) ! psi rigtco(must:mued)=Prigte(nden)*Prigte(nst:ned) ! rho X ! Flux ptl = Prigte(npre) + ( Prigte(nbm1)**2 & & +Prigte(nbm2)**2 & & +Prigte(nbm3)**2)/2.0d0 rigtco(mfdn)=Prigte(nden) *Prigte(nve1) rigtco(mfvu)=Prigte(nden)*Prigte(nve1)*Prigte(nve1) & & +ptl-Prigte(nbm1)*Prigte(nbm1) rigtco(mfvv)=Prigte(nden)*Prigte(nve2)*Prigte(nve1) & & -Prigte(nbm2)*Prigte(nbm1) rigtco(mfvw)=Prigte(nden)*Prigte(nve3)*Prigte(nve1) & & -Prigte(nbm3)*Prigte(nbm1) rigtco(mfet)=(etot+ptl)*Prigte(nve1) & & -( Prigte(nbm1)*Prigte(nve1) & & +Prigte(nbm2)*Prigte(nve2) & & +Prigte(nbm3)*Prigte(nve3))*Prigte(nbm1) rigtco(mfbu) = 0.0d0 rigtco(mfbv) = Prigte(nbm2)*Prigte(nve1) & & -Prigte(nve2)*Prigte(nbm1) rigtco(mfbw) = Prigte(nbm3)*Prigte(nve1) & & -Prigte(nve3)*Prigte(nbm1) rigtco(mfbp) = 0.0d0 ! b_z rigtco(mfst:mfed)=Prigte(nden)*Prigte(nst:ned)*Prigte(nve1) ! 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(nbm1)**2 & & /Prigte(nden)) & & )/2.0d0) rigtco(mvel)= Prigte(nve1) rigtco(mpre)= ptl ! call HLLE(leftco,rigtco,nflux) ! call HLLC(leftco,rigtco,nflux) call HLLD(leftco,rigtco,nflux) nflux1(mden,i,j,k)=nflux(mden) nflux1(mrv1,i,j,k)=nflux(mrvu) nflux1(mrv2,i,j,k)=nflux(mrvv) nflux1(mrv3,i,j,k)=nflux(mrvw) nflux1(meto,i,j,k)=nflux(meto) nflux1(mbm1,i,j,k)=nflux(mbmu) nflux1(mbm2,i,j,k)=nflux(mbmv) nflux1(mbm3,i,j,k)=nflux(mbmw) nflux1(mst:med,i,j,k)=nflux(mst:med) nflux1(mbm1,i,j,k) = 0.5d0*(leftco(mubp)+rigtco(mubp)) & & -0.5d0*chg*(rigtco(mubu)-leftco(mubu)) ! finite volume nflux1(mbps,i,j,k) = (0.5d0*(leftco(mubu)+rigtco(mubu)) & & -0.5d0/chg*(rigtco(mubp)-leftco(mubp)))*chg**2 ! finite volume ! print *, "bpf1",nflux1(mbps,i,j,k) enddo enddo enddo !$acc end kernels return end subroutine Numericalflux1