!!****if* source/physics/Hydro/HydroMain/unsplit/hy_uhd_getFaceFlux !! !! NAME !! !! hy_uhd_getFaceFlux !! !! SYNOPSIS !! !! hy_uhd_getFaceFlux( integer(IN) :: blockID, !! integer(IN) :: blkLimits(2,MDIM) !! integer(IN) :: blkLimitsGC(2,MDIM), !! integer(IN) :: datasize(MDIM), !! real(IN) :: del(MDIM), !! real(IN) :: lambda(:,:,:,:,:), !! real(OUT) :: xflux(:,:,:,:), !! real(OUT) :: yflux(:,:,:,:), !! real(OUT) :: zflux(:,:,:,:), !! real(OUT) :: vint(:,:,:,:), !! logical,optional(IN) :: lastCall ) !! !! ARGUMENTS !! !! blockID - a current block ID !! blkLimits - an array that holds the lower and upper indices of the section !! of block without the guard cells !! blkLimitsGC - an array that holds the lower and upper indices of the section !! of block with the guard cells !! datasize - data size for boundary extrapolated data, boundary data !! del - grid deltas !! lambda - eigenvalues !! xflux,yflux,zflux - face fluxes at each {x,y,z} direction !! vint - interface velocity !! lastCall - if true then store flux data in scratch array !! !! DESCRIPTION !! !! This routine computes high-order Godunov fluxes at cell interface centers !! for each spatial direction using a choice of Riemann solvers. !! Choices of Riemann solvers are Roe-type, HLL(E), HLLC, HLLD, !! Hybrid (HLLC+HLL, Hydro solver only), and local Lax-Friedrichs solvers. !! !!*** !!REORDER(4):U, scrch_Ctr, scrch_[XYZ], [xyz]flux #include "Flash.h" subroutine hy_uhd_getFaceFlux ( blockID,blkLimits,blkLimitsGC,datasize,del,lambda,& xflux,yflux,zflux,vint,lastCall) use Hydro_data, ONLY : hy_nref,hy_kref, & hy_shockInstabilityFix, & hy_RiemannSolver, & hy_useDiffuse, & hy_useViscosity, & hy_useConductivity, & hy_EOSforRiemann, & hy_use_avisc, & hy_cvisc use hy_uhd_interface, ONLY : hy_uhd_addViscousFluxes, & hy_uhd_addThermalFluxes, & hy_uhd_Roe, & hy_uhd_LLF, & hy_uhd_HLL, & hy_uhd_HLLC,& hy_uhd_Marquina,& hy_uhd_shockInstabilityFix use hy_uhd_slopeLimiters, ONLY : signum use Grid_interface, ONLY : Grid_getBlkPtr, Grid_releaseBlkPtr use Eos_interface, ONLY : Eos use Conductivity_interface, ONLY : Conductivity use Viscosity_interface, ONLY : Viscosity use Timers_interface, ONLY : Timers_start, Timers_stop #ifdef FLASH_USM_MHD use Hydro_data, ONLY : hy_mref, hy_useResistivity, hy_useMagneticResistivity use hy_uhd_interface, ONLY : hy_uhd_addResistiveFluxes, hy_uhd_HLLD use MagneticResistivity_interface, ONLY : MagneticResistivity #endif implicit none #include "constants.h" #include "Eos.h" #include "UHD.h" !! Arguments type declaration ------------------------------ integer, intent(IN) :: blockID integer, dimension(LOW:HIGH,MDIM),intent(IN) :: blkLimits, blkLimitsGC integer, dimension(MDIM), intent(IN) :: datasize real, dimension(MDIM), intent(IN) :: del #ifdef FIXEDBLOCKSIZE real, intent(IN) :: lambda(NDIM,HY_WAVENUM,& GRID_IHI_GC,GRID_JHI_GC,GRID_KHI_GC) real, intent(out) :: xflux (NFLUXES,& GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) real, intent(out) :: yflux (NFLUXES,& GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) real, intent(out) :: zflux (NFLUXES,& GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) real, intent(out) :: vint (NDIM,& GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) #else real, intent(IN) :: lambda(NDIM,HY_WAVENUM,& blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) real, intent(OUT) :: xflux (NFLUXES,& blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) real, intent(OUT) :: yflux (NFLUXES,& blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) real, intent(OUT) :: zflux (NFLUXES,& blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) real, intent(OUT) :: vint (NDIM,& blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) #endif logical, optional, intent(IN) :: lastCall !! --------------------------------------------------------- integer :: i0,imax,j0,jmax,k0,kmax,jbeg,jend,kbeg,kend, i,j,k real, pointer, dimension(:,:,:,:) :: scrch_Ctr,scrch_X,scrch_Y,scrch_Z,U !(DENS,VELX,VELY,VELZ,(MAGX,MAGY,MAGZ),PRES + GAMC,GAME,EINT,TEMP) real, dimension(HY_VARINUM4) :: VL,VR real, dimension(NSPECIES) :: speciesArr #ifdef FIXEDBLOCKSIZE real, dimension(GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) & :: viscDynamic,viscKinematic,cond,dcff #ifdef FLASH_USM_MHD real, dimension(GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) & :: magResist #endif real, dimension(NDIM,2,& GRID_ILO_GC:GRID_IHI_GC, & GRID_JLO_GC:GRID_JHI_GC, & GRID_KLO_GC:GRID_KHI_GC) :: LambdaFix #else real, dimension(blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) & :: viscDynamic,viscKinematic,cond,dcff #ifdef FLASH_USM_MHD real, dimension(blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) & :: magResist #endif real, dimension(NDIM,2, & blkLimitsGC(LOW,IAXIS):blkLimitsGC(HIGH,IAXIS), & blkLimitsGC(LOW,JAXIS):blkLimitsGC(HIGH,JAXIS), & blkLimitsGC(LOW,KAXIS):blkLimitsGC(HIGH,KAXIS)) & :: LambdaFix #endif real,dimension(EOS_NUM)::eosData logical,dimension(EOS_VARS+1:EOS_NUM)::eosMask=.false. real, dimension(NSPECIES) :: eosMf integer :: vecLen integer :: interp_eosMode = MODE_DENS_PRES real :: sig,absig real :: cvisc !! initialize vecLen=1 eosMf =1. VL = 0. VR = 0. !! Call scratch for storing fluxes call Grid_getBlkPtr(blockID,scrch_Ctr,SCRATCH_CTR) call Grid_getBlkPtr(blockID,scrch_X,SCRATCH_CTR) if (NDIM > 1) call Grid_getBlkPtr(blockID,scrch_Y,SCRATCH_CTR) if (NDIM > 2) call Grid_getBlkPtr(blockID,scrch_Z,SCRATCH_CTR) call Grid_getBlkPtr(blockID,U,CENTER) i0 = blkLimits(LOW, IAXIS) imax = blkLimits(HIGH,IAXIS) j0 = blkLimits(LOW, JAXIS) jmax = blkLimits(HIGH,JAXIS) k0 = blkLimits(LOW, KAXIS) kmax = blkLimits(HIGH,KAXIS) if (NDIM == 1) then jbeg= 3 jend=-1 kbeg= 3 kend=-1 elseif (NDIM == 2) then jbeg= j0 jend= jmax kbeg= 3 kend=-1 else jbeg= j0 jend= jmax kbeg= k0 kend= kmax endif if (hy_useDiffuse) then ! call Timers_start("get diffusion") !! Initialize #ifdef FLASH_USM_MHD magResist = 0.0 #endif viscDynamic = 0.0 viscKinematic = 0.0 cond = 0.0 dcff = 0.0 do k=kbeg-2,kend+2 do j=jbeg-2,jend+2 do i=i0-2,imax+2 !! copy species to a temporal array speciesArr(:) = U(SPECIES_BEGIN:SPECIES_END,i,j,k) if (hy_useViscosity) then !! Get viscosity call Viscosity(U(TEMP_VAR,i,j,k),U(DENS_VAR,i,j,k),& speciesArr,viscDynamic(i,j,k),viscKinematic(i,j,k)) endif if (hy_useConductivity) then !! Get heat conductivity call Conductivity(U(TEMP_VAR,i,j,k),U(DENS_VAR,i,j,k),& speciesArr,cond(i,j,k),dcff(i,j,k),2) endif #ifdef FLASH_USM_MHD if (hy_useMagneticResistivity) then !! Get magnetic viscosity call MagneticResistivity(U(TEMP_VAR,i,j,k),U(DENS_VAR,i,j,k),& speciesArr,magResist(i,j,k)) endif #endif enddo enddo enddo !! For non-ideal fluxes #ifdef FLASH_USM_MHD magResist = magResist/hy_mref #endif viscDynamic= viscDynamic/hy_nref endif !! Initialize flux arrays xflux = 0. yflux = 0. zflux = 0. !! Fix for shock-instabilities such as odd-even, carbuncle instabilities if (hy_shockInstabilityFix) then call hy_uhd_shockInstabilityFix(blockID,blkLimits,blkLimitsGC,lambda,LambdaFix) endif !! Compute intercell fluxes using the updated left & right states !! Calculate x-flux first do k=kbeg-2,kend+2 do j=jbeg-2,jend+2 do i=i0-1,imax+2 VL(HY_DENS:HY_GAME)=scrch_X(XP01_SCRATCH_CENTER_VAR:XP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i-1,j,k) VR(HY_DENS:HY_GAME)=scrch_X(XN01_SCRATCH_CENTER_VAR:XN01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i, j,k) if (hy_EOSforRiemann) then !! Left states eosData(EOS_DENS)=VL(HY_DENS) eosData(EOS_PRES)=VL(HY_PRES) eosData(EOS_TEMP)=U(TEMP_VAR,i-1,j,k) #ifdef YE_MSCALAR eosData(EOS_ABAR)=1./U(SUMY_MSCALAR,i-1,j,k) eosData(EOS_ZBAR)= U(YE_MSCALAR, i-1,j,k)/U(SUMY_MSCALAR,i-1,j,k) #endif eosMf=U(SPECIES_BEGIN:SPECIES_END,i-1,j,k) call Eos(interp_eosMode,vecLen,eosData,eosMf,eosMask) VL(HY_GAMC) =eosData(EOS_GAMC) VL(HY_EINT) =eosData(EOS_EINT) VL(HY_GAME) =1.+eosData(EOS_PRES)/(eosData(EOS_DENS)*eosData(EOS_EINT)) ! game VL(HY_TEMP) =eosData(EOS_TEMP) !! Right states eosData(EOS_DENS)=VR(HY_DENS) eosData(EOS_PRES)=VR(HY_PRES) eosData(EOS_TEMP)=U(TEMP_VAR,i,j,k) #ifdef YE_MSCALAR eosData(EOS_ABAR)=1./U(SUMY_MSCALAR,i,j,k) eosData(EOS_ZBAR)= U(YE_MSCALAR, i,j,k)/U(SUMY_MSCALAR,i,j,k) #endif eosMf=U(SPECIES_BEGIN:SPECIES_END,i,j,k) call Eos(interp_eosMode,vecLen,eosData,eosMf,eosMask) VR(HY_GAMC) =eosData(EOS_GAMC) VR(HY_EINT) =eosData(EOS_EINT) VR(HY_GAME) =1.+eosData(EOS_PRES)/(eosData(EOS_DENS)*eosData(EOS_EINT)) ! game VR(HY_TEMP) =eosData(EOS_TEMP) endif !end of EOS if (hy_RiemannSolver == ROE) then call hy_uhd_Roe(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k),& i,j,k,blkLimitsGC,LambdaFix) elseif (hy_RiemannSolver == HLL) then call hy_uhd_HLL(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) elseif (hy_RiemannSolver == HLLC) then call hy_uhd_HLLC(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) #ifdef FLASH_USM_MHD elseif (hy_RiemannSolver == HLLD) then call hy_uhd_HLLD(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) #endif elseif (hy_RiemannSolver == MARQ) then call hy_uhd_Marquina(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) elseif (hy_RiemannSolver == LLF) then call hy_uhd_LLF(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) #ifndef FLASH_USM_MHD elseif (hy_RiemannSolver == HYBR) then if (scrch_Ctr(FLAG_SCRATCH_CENTER_VAR,i-1,j,k) + scrch_Ctr(FLAG_SCRATCH_CENTER_VAR,i,j,k) > 0.) then !! use diffusive HLL solver for a local strong shock/rarefaction region call hy_uhd_HLL(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) else call hy_uhd_HLLC(DIR_X,VL,VR,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_X,i,j,k)) endif #endif endif !! Flux for internal energy density !! Reference: "Simple Method to Track Pressure Accurately", S. Li, Astronum Proceeding, 2007 !! Note that there is a typo in Li's paper related on the left and right states. sig=signum(xflux(HY_DENS_FLUX,i,j,k)) absig=abs(sig) xflux(HY_EINT_FLUX,i,j,k) = & 0.5*absig*( xflux(HY_DENS_FLUX,i,j,k)*VR(HY_EINT)*(1.-sig)& +xflux(HY_DENS_FLUX,i,j,k)*VL(HY_EINT)*(1.+sig)) xflux(HY_PRES_FLUX,i,j,k) = & 0.5*absig*( xflux(HY_DENS_FLUX,i,j,k)/VR(HY_DENS)*(1.-sig)& +xflux(HY_DENS_FLUX,i,j,k)/VL(HY_DENS)*(1.+sig)) !! Artificial viscosity as in PPM, Colella and Woodward, 1984. if (hy_use_avisc) then if (NDIM == 1) then cvisc = hy_cvisc*max( -(U(VELX_VAR,i, j, k )-U(VELX_VAR,i-1,j, k )),0.) elseif (NDIM == 2) then cvisc = hy_cvisc*max(-( U(VELX_VAR,i, j, k )-U(VELX_VAR,i-1,j, k ) + & 0.25*(U(VELY_VAR,i, j+1,k )-U(VELY_VAR,i, j-1,k ) + & U(VELY_VAR,i-1,j+1,k )-U(VELY_VAR,i-1,j-1,k ))*del(DIR_X)/del(DIR_Y)),& 0.) else cvisc = hy_cvisc*max(-( U(VELX_VAR,i, j, k )-U(VELX_VAR,i-1,j, k ) + & 0.25*( U(VELY_VAR,i, j+1,k )-U(VELY_VAR,i, j-1,k ) + & U(VELY_VAR,i-1,j+1,k )-U(VELY_VAR,i-1,j-1,k ))*del(DIR_X)/del(DIR_Y) + & 0.25*( U(VELZ_VAR,i, j, k+1)-U(VELZ_VAR,i, j, k-1) + & U(VELZ_VAR,i-1,j, k+1)-U(VELZ_VAR,i-1,j, k-1))*del(DIR_X)/del(DIR_Z)), & 0.) endif xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k) = & xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k) & +cvisc*(/U(DENS_VAR,i-1,j,k) -U(DENS_VAR,i,j,k),& U(DENS_VAR,i-1,j,k)*U(VELX_VAR,i-1,j,k)-U(DENS_VAR,i,j,k)*U(VELX_VAR,i,j,k),& U(DENS_VAR,i-1,j,k)*U(VELY_VAR,i-1,j,k)-U(DENS_VAR,i,j,k)*U(VELY_VAR,i,j,k),& U(DENS_VAR,i-1,j,k)*U(VELZ_VAR,i-1,j,k)-U(DENS_VAR,i,j,k)*U(VELZ_VAR,i,j,k),& U(DENS_VAR,i-1,j,k)*U(ENER_VAR,i-1,j,k)-U(DENS_VAR,i,j,k)*U(ENER_VAR,i,j,k) & #ifdef FLASH_USM_MHD ,& U(MAGX_VAR,i-1,j,k) -U(MAGX_VAR,i,j,k),& U(MAGY_VAR,i-1,j,k) -U(MAGY_VAR,i,j,k),& U(MAGZ_VAR,i-1,j,k) -U(MAGZ_VAR,i,j,k) & #endif /) endif if (hy_useDiffuse) then if (hy_useViscosity) then call hy_uhd_addViscousFluxes& (blockID,blkLimitsGC,i,j,k,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),viscDynamic,DIR_X) endif if (hy_useConductivity) then call hy_uhd_addThermalFluxes& (blockID,blkLimitsGC,i,j,k,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),cond,DIR_X) endif #ifdef FLASH_USM_MHD if (hy_useMagneticResistivity) then call hy_uhd_addResistiveFluxes& (blockID,blkLimitsGC,i,j,k,xflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),magResist,DIR_X) endif #endif endif enddo enddo enddo #ifdef FLASH_USM_MHD if (present(lastCall)) then #endif !! Make another copy of pressure for internal energy update in hy_uhd_unsplitUpdate scrch_X(XN04_SCRATCH_CENTER_VAR,:,:,:) = scrch_X(XP05_SCRATCH_CENTER_VAR,:,:,:) #ifndef FLASH_GRID_UG !! initialize with zero scrch_X(XP01_SCRATCH_CENTER_VAR:XP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,:,:,:) = 0. !! re-use XP array for storing fluxes do k=kbeg-2,kend+2 do j=jbeg-2,jend+2 do i=i0-1,imax+2 scrch_X(XP01_SCRATCH_CENTER_VAR:XP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j,k) = & xflux(HY_DENS_FLUX:HY_PRES_FLUX,i,j,k) enddo enddo enddo #endif #ifdef FLASH_USM_MHD endif !end of if-lastCall #endif #if NDIM >= 2 !! Calculate y-flux do k=kbeg-2,kend+2 do j=jbeg-1,jend+2 do i=i0-2,imax+2 VL(HY_DENS:HY_GAME)=scrch_Y(YP01_SCRATCH_CENTER_VAR:YP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j-1,k) VR(HY_DENS:HY_GAME)=scrch_Y(YN01_SCRATCH_CENTER_VAR:YN01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j, k) if (hy_EOSforRiemann) then !! Left states eosData(EOS_DENS)=VL(HY_DENS) eosData(EOS_PRES)=VL(HY_PRES) eosData(EOS_TEMP)=U(TEMP_VAR,i,j-1,k) #ifdef YE_MSCALAR eosData(EOS_ABAR)=1./U(SUMY_MSCALAR,i,j-1,k) eosData(EOS_ZBAR)= U(YE_MSCALAR, i,j-1,k)/U(SUMY_MSCALAR,i,j-1,k) #endif eosMf=U(SPECIES_BEGIN:SPECIES_END,i,j-1,k) call Eos(interp_eosMode,vecLen,eosData,eosMf,eosMask) VL(HY_GAMC) =eosData(EOS_GAMC) VL(HY_EINT) =eosData(EOS_EINT) VL(HY_GAME) =1.+eosData(EOS_PRES)/(eosData(EOS_DENS)*eosData(EOS_EINT)) ! game VL(HY_TEMP) =eosData(EOS_TEMP) !! Right states eosData(EOS_DENS)=VR(HY_DENS) eosData(EOS_PRES)=VR(HY_PRES) eosData(EOS_TEMP)=U(TEMP_VAR,i,j,k) #ifdef YE_MSCALAR eosData(EOS_ABAR)=1./U(SUMY_MSCALAR,i,j,k) eosData(EOS_ZBAR)= U(YE_MSCALAR, i,j,k)/U(SUMY_MSCALAR,i,j,k) #endif eosMf=U(SPECIES_BEGIN:SPECIES_END,i,j,k) call Eos(interp_eosMode,vecLen,eosData,eosMf,eosMask) VR(HY_GAMC) =eosData(EOS_GAMC) VR(HY_EINT) =eosData(EOS_EINT) VR(HY_GAME) =1.+eosData(EOS_PRES)/(eosData(EOS_DENS)*eosData(EOS_EINT)) ! game VR(HY_TEMP) =eosData(EOS_TEMP) endif !end of EOS if (hy_RiemannSolver == ROE) then call hy_uhd_Roe(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k),& i,j,k,blkLimitsGC,LambdaFix) elseif (hy_RiemannSolver == HLL) then call hy_uhd_HLL(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) elseif (hy_RiemannSolver == HLLC) then call hy_uhd_HLLC(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) #ifdef FLASH_USM_MHD elseif (hy_RiemannSolver == HLLD) then call hy_uhd_HLLD(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) #endif elseif (hy_RiemannSolver == MARQ) then call hy_uhd_Marquina(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) elseif (hy_RiemannSolver == LLF) then call hy_uhd_LLF(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) #ifndef FLASH_USM_MHD elseif (hy_RiemannSolver == HYBR) then if (scrch_Ctr(FLAG_SCRATCH_CENTER_VAR,i,j-1,k) + scrch_Ctr(FLAG_SCRATCH_CENTER_VAR,i,j,k) > 0.) then !! use diffusive HLL solver for a local strong shock/rarefaction region call hy_uhd_HLL(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) else call hy_uhd_HLLC(DIR_Y,VL,VR,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Y,i,j,k)) endif #endif endif !! Flux for internal energy density !! Reference: "Simple Method to Track Pressure Accurately", S. Li, Astronum Proceeding, 2007 sig=signum(yflux(HY_DENS_FLUX,i,j,k)) absig=abs(sig) yflux(HY_EINT_FLUX,i,j,k) = & 0.5*absig*( yflux(HY_DENS_FLUX,i,j,k)*VR(HY_EINT)*(1.-sig)& +yflux(HY_DENS_FLUX,i,j,k)*VL(HY_EINT)*(1.+sig)) yflux(HY_PRES_FLUX,i,j,k) = & 0.5*absig*( yflux(HY_DENS_FLUX,i,j,k)/VR(HY_DENS)*(1.-sig)& +yflux(HY_DENS_FLUX,i,j,k)/VL(HY_DENS)*(1.+sig)) !! Artificial viscosity as in PPM, Colella and Woodward, 1984. if (hy_use_avisc) then if (NDIM == 2) then cvisc = hy_cvisc*max(-(0.25*( U(VELX_VAR,i+1,j, k )-U(VELX_VAR,i-1,j, k ) + & U(VELX_VAR,i+1,j-1,k )-U(VELX_VAR,i-1,j-1,k ))*del(DIR_Y)/del(DIR_X) + & U(VELY_VAR,i, j, k )-U(VELY_VAR,i, j-1,k )), & 0.) else cvisc = hy_cvisc*max(-(0.25*( U(VELX_VAR,i+1,j, k )-U(VELX_VAR,i-1,j, k ) + & U(VELX_VAR,i+1,j-1,k )-U(VELX_VAR,i-1,j-1,k ))*del(DIR_Y)/del(DIR_X) + & U(VELY_VAR,i, j, k )-U(VELY_VAR,i, j-1,k ) + & 0.25*( U(VELZ_VAR,i, j, k+1)-U(VELZ_VAR,i-1,j, k-1) + & U(VELZ_VAR,i, j-1,k+1)-U(VELZ_VAR,i, j-1,k-1))*del(DIR_Y)/del(DIR_Z)), & 0.) endif yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k) = & yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k) & +cvisc*(/U(DENS_VAR,i,j-1,k) -U(DENS_VAR,i,j,k),& U(DENS_VAR,i,j-1,k)*U(VELX_VAR,i,j-1,k)-U(DENS_VAR,i,j,k)*U(VELX_VAR,i,j,k),& U(DENS_VAR,i,j-1,k)*U(VELY_VAR,i,j-1,k)-U(DENS_VAR,i,j,k)*U(VELY_VAR,i,j,k),& U(DENS_VAR,i,j-1,k)*U(VELZ_VAR,i,j-1,k)-U(DENS_VAR,i,j,k)*U(VELZ_VAR,i,j,k),& U(DENS_VAR,i,j-1,k)*U(ENER_VAR,i,j-1,k)-U(DENS_VAR,i,j,k)*U(ENER_VAR,i,j,k) & #ifdef FLASH_USM_MHD ,& U(MAGX_VAR,i,j-1,k) -U(MAGX_VAR,i,j,k),& U(MAGY_VAR,i,j-1,k) -U(MAGY_VAR,i,j,k),& U(MAGZ_VAR,i,j-1,k) -U(MAGZ_VAR,i,j,k) & #endif /) endif if (hy_useDiffuse) then if (hy_useViscosity) then call hy_uhd_addViscousFluxes& (blockID,blkLimitsGC,i,j,k,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),viscDynamic,DIR_Y) endif if (hy_useConductivity) then call hy_uhd_addThermalFluxes& (blockID,blkLimitsGC,i,j,k,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),cond,DIR_Y) endif #ifdef FLASH_USM_MHD if (hy_useMagneticResistivity) then call hy_uhd_addResistiveFluxes& (blockID,blkLimitsGC,i,j,k,yflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),magResist,DIR_Y) endif #endif endif enddo enddo enddo #ifdef FLASH_USM_MHD if (present(lastCall)) then #endif !! Make another copy of pressure for internal energy update in hy_uhd_unsplitUpdate scrch_Y(YN04_SCRATCH_CENTER_VAR,:,:,:) = scrch_Y(YP05_SCRATCH_CENTER_VAR,:,:,:) #ifndef FLASH_GRID_UG !! initialize with zero scrch_Y(YP01_SCRATCH_CENTER_VAR:YP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,:,:,:) = 0. !! re-use YP array for storing fluxes do k=kbeg-2,kend+2 do j=jbeg-1,jend+2 do i=i0-2,imax+2 scrch_Y(YP01_SCRATCH_CENTER_VAR:YP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j,k) = & yflux(HY_DENS_FLUX:HY_PRES_FLUX,i,j,k) enddo enddo enddo #endif #ifdef FLASH_USM_MHD endif !end of if-lastCall #endif #if NDIM == 3 !! Calculate z-flux do k=kbeg-1,kend+2 do j=jbeg-2,jend+2 do i=i0-2,imax+2 VL(HY_DENS:HY_GAME)=scrch_Z(ZP01_SCRATCH_CENTER_VAR:ZP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j,k-1) VR(HY_DENS:HY_GAME)=scrch_Z(ZN01_SCRATCH_CENTER_VAR:ZN01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j,k) if (hy_EOSforRiemann) then !! Left states eosData(EOS_DENS)=VL(HY_DENS) eosData(EOS_PRES)=VL(HY_PRES) eosData(EOS_TEMP)=U(TEMP_VAR,i,j,k-1) #ifdef YE_MSCALAR eosData(EOS_ABAR)=(1./U(SUMY_MSCALAR,i,j,k-1) eosData(EOS_ZBAR)= U(YE_MSCALAR, i,j,k-1)/U(SUMY_MSCALAR,i,j,k-1) #endif eosMf=U(SPECIES_BEGIN:SPECIES_END,i,j,k-1) call Eos(interp_eosMode,vecLen,eosData,eosMf,eosMask) VL(HY_GAMC) =eosData(EOS_GAMC) VL(HY_EINT) =eosData(EOS_EINT) VL(HY_GAME) =1.+eosData(EOS_PRES)/(eosData(EOS_DENS)*eosData(EOS_EINT)) ! game VL(HY_TEMP) =eosData(EOS_TEMP) !! Right states eosData(EOS_DENS)=VR(HY_DENS) eosData(EOS_PRES)=VR(HY_PRES) eosData(EOS_TEMP)=U(TEMP_VAR,i,j,k) #ifdef YE_MSCALAR eosData(EOS_ABAR)=1./U(SUMY_MSCALAR,i,j,k) eosData(EOS_ZBAR)= U(YE_MSCALAR, i,j,k)/U(SUMY_MSCALAR,i,j,k) #endif eosMf=U(SPECIES_BEGIN:SPECIES_END,i,j,k) call Eos(interp_eosMode,vecLen,eosData,eosMf,eosMask) VR(HY_GAMC) =eosData(EOS_GAMC) VR(HY_EINT) =eosData(EOS_EINT) VR(HY_GAME) =1.+eosData(EOS_PRES)/(eosData(EOS_DENS)*eosData(EOS_EINT)) ! game VR(HY_TEMP) =eosData(EOS_TEMP) endif !end of EOS if (hy_RiemannSolver == ROE) then call hy_uhd_Roe(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k),& i,j,k,blkLimitsGC,LambdaFix) elseif (hy_RiemannSolver == HLL) then call hy_uhd_HLL(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) elseif (hy_RiemannSolver == HLLC) then call hy_uhd_HLLC(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) #ifdef FLASH_USM_MHD elseif (hy_RiemannSolver == HLLD) then call hy_uhd_HLLD(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) #endif elseif (hy_RiemannSolver == MARQ) then call hy_uhd_Marquina(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) elseif (hy_RiemannSolver == LLF) then call hy_uhd_LLF(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) #ifndef FLASH_USM_MHD elseif (hy_RiemannSolver == HYBR) then if (scrch_Ctr(FLAG_SCRATCH_CENTER_VAR,i,j,k-1) + scrch_Ctr(FLAG_SCRATCH_CENTER_VAR,i,j,k) > 0.) then !! use diffusive HLL solver for a local strong shock/rarefaction region call hy_uhd_HLL(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) else call hy_uhd_HLLC(DIR_Z,VL,VR,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),vint(DIR_Z,i,j,k)) endif #endif endif !! Flux for internal energy density !! Reference: "Simple Method to Track Pressure Accurately", S. Li, Astronum Proceeding, 2007 sig=signum(zflux(HY_DENS_FLUX,i,j,k)) absig=abs(sig) zflux(HY_EINT_FLUX,i,j,k) = & 0.5*absig*( zflux(HY_DENS_FLUX,i,j,k)*VR(HY_EINT)*(1.-sig)& +zflux(HY_DENS_FLUX,i,j,k)*VL(HY_EINT)*(1.+sig)) zflux(HY_PRES_FLUX,i,j,k) = & 0.5*absig*( zflux(HY_DENS_FLUX,i,j,k)/VR(HY_DENS)*(1.-sig)& +zflux(HY_DENS_FLUX,i,j,k)/VL(HY_DENS)*(1.+sig)) !! Artificial viscosity as in PPM, Colella and Woodward, 1984. if (hy_use_avisc) then cvisc = hy_cvisc*max(-(0.25*( U(VELX_VAR,i+1,j, k )-U(VELX_VAR,i-1,j, k ) + & U(VELX_VAR,i+1,j, k-1)-U(VELX_VAR,i-1,j, k-1))*del(DIR_Z)/del(DIR_X) + & 0.25*( U(VELY_VAR,i, j+1,k )-U(VELY_VAR,i, j-1,k ) + & U(VELY_VAR,i, j+1,k-1)-U(VELY_VAR,i, j-1,k-1))*del(DIR_Z)/del(DIR_Y) + & U(VELZ_VAR,i, j, k )-U(VELZ_VAR,i, j, k-1)),& 0.) zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k) = & zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k) & +cvisc*(/U(DENS_VAR,i,j,k-1) -U(DENS_VAR,i,j,k),& U(DENS_VAR,i,j,k-1)*U(VELX_VAR,i,j,k-1)-U(DENS_VAR,i,j,k)*U(VELX_VAR,i,j,k),& U(DENS_VAR,i,j,k-1)*U(VELY_VAR,i,j,k-1)-U(DENS_VAR,i,j,k)*U(VELY_VAR,i,j,k),& U(DENS_VAR,i,j,k-1)*U(VELZ_VAR,i,j,k-1)-U(DENS_VAR,i,j,k)*U(VELZ_VAR,i,j,k),& U(DENS_VAR,i,j,k-1)*U(ENER_VAR,i,j,k-1)-U(DENS_VAR,i,j,k)*U(ENER_VAR,i,j,k) & #ifdef FLASH_USM_MHD ,& U(MAGX_VAR,i,j,k-1) -U(MAGX_VAR,i,j,k),& U(MAGY_VAR,i,j,k-1) -U(MAGY_VAR,i,j,k),& U(MAGZ_VAR,i,j,k-1) -U(MAGZ_VAR,i,j,k) & #endif /) endif if (hy_useDiffuse) then if (hy_useViscosity) then call hy_uhd_addViscousFluxes& (blockID,blkLimitsGC,i,j,k,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),viscDynamic,DIR_Z) endif if (hy_useConductivity) then call hy_uhd_addThermalFluxes& (blockID,blkLimitsGC,i,j,k,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),cond,DIR_Z) endif #ifdef FLASH_USM_MHD if (hy_useMagneticResistivity) then call hy_uhd_addResistiveFluxes& (blockID,blkLimitsGC,i,j,k,zflux(F01DENS_FLUX:F01DENS_FLUX+HY_VARINUM-1,i,j,k),magResist,DIR_Z) endif #endif endif enddo enddo enddo #ifdef FLASH_USM_MHD if (present(lastCall)) then #endif !! Make another copy of pressure for internal energy update in hy_uhd_unsplitUpdate scrch_Z(ZN04_SCRATCH_CENTER_VAR,:,:,:) = scrch_Z(ZP05_SCRATCH_CENTER_VAR,:,:,:) #ifndef FLASH_GRID_UG !! initialize with zero scrch_Z(ZP01_SCRATCH_CENTER_VAR:ZP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,:,:,:) = 0. !! re-use ZP array for storing fluxes do k=kbeg-1,kend+2 do j=jbeg-2,jend+2 do i=i0-2,imax+2 scrch_Z(ZP01_SCRATCH_CENTER_VAR:ZP01_SCRATCH_CENTER_VAR+HY_SCRATCH_NUM-1,i,j,k) = & zflux(HY_DENS_FLUX:HY_PRES_FLUX,i,j,k) enddo enddo enddo #endif #ifdef FLASH_USM_MHD endif !end of if-lastCall #endif #endif !end of #if NDIM==3 #endif !end of #if NDIM >= 2 !! Release pointer call Grid_releaseBlkPtr(blockID,scrch_Ctr,SCRATCH_CTR) call Grid_releaseBlkPtr(blockID,scrch_X,SCRATCH_CTR) if (NDIM > 1) call Grid_releaseBlkPtr(blockID,scrch_Y,SCRATCH_CTR) if (NDIM > 2) call Grid_releaseBlkPtr(blockID,scrch_Z,SCRATCH_CTR) call Grid_releaseBlkPtr(blockID,U,CENTER) End Subroutine hy_uhd_getFaceFlux