!> Module with finite volume methods for fluxes module mod_finite_volume_srhd use mod_variables use mod_physics_vars use mod_global_parameters, only: ndim use mod_global_parameters, only: flux_ad_min, flux_ad_scale use mod_physicaldata implicit none private public :: finite_volume_local contains ! instantiate the templated functions here for inlining: subroutine addsource_local(qdt, dtfactor, qtC, wCT, wCTprim, qt, wnew, x, dr, & qsourcesplit) use mod_usr, only: addsource_usr real(dp), intent(in) :: qdt, dtfactor, qtC, qt real(dp), intent(in) :: wCT(nw_phys), wCTprim(nw_phys) real(dp), intent(in) :: x(1:ndim), dr(ndim) real(dp), intent(inout) :: wnew(nw_phys) logical, intent(in) :: qsourcesplit if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- call addsource_usr(qdt, qt, wCT, wCTprim, wnew, x, .false.) !!!else !--------------------------------- ! split sources !--------------------------------- end if end subroutine addsource_local subroutine addsource_nonlocal(qdt, dtfactor, qtC, wCTprim, qt, wnew, x, dx, idir, & qsourcesplit) real(dp), intent(in) :: qdt, dtfactor, qtC, qt real(dp), intent(in) :: wCTprim(nw_phys,5) real(dp), intent(in) :: x(1:ndim), dx(1:ndim) real(dp), intent(inout) :: wnew(nw_phys) integer, intent(in) :: idir logical, intent(in) :: qsourcesplit end subroutine addsource_nonlocal subroutine addsource_compact(qdt, dtfactor, qtC, wCTprim1, wCTprim2, wCTprim3, qt, wnew, x, dx, & qsourcesplit) real(dp), intent(in) :: qdt, dtfactor, qtC, qt real(dp), intent(in) :: wCTprim1(nw_phys,3),wCTprim2(nw_phys,3),wCTprim3(nw_phys,3) real(dp), intent(in) :: x(1:ndim), dx(1:ndim) real(dp), intent(inout) :: wnew(nw_phys) logical, intent(in) :: qsourcesplit end subroutine addsource_compact pure subroutine to_primitive(u) use mod_con2prim real(dp), intent(inout) :: u(nw_phys) real(dp) :: rho,rhoh,pth,E real(dp) :: ssqr !! begin: call srhd_get_auxiliary(u) ssqr=(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2) if(srhd_eos)then call con2prim_eos(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e)) else call con2prim(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e)) endif !! end: call srhd_get_auxiliary(u) rho=u(iw_rho)/u(lfac_) rhoh=u(xi_)/u(lfac_)**2 !! begin: call srhd_get_pressure_eos(rho,rhoh,pth,E) if(srhd_eos)then E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)*rho**2)) & /(srhd_gamma+1.0_dp) pth = 0.5d0*gamma_1* (E-rho*(rho/E)) else pth =(rhoh-rho)/gamma_to_gamma_1 endif !! end: call srhd_get_pressure_eos(rho,rhoh,pth,E) u(iw_rho)=rho u(iw_mom(1))=u(lfac_)*u(iw_mom(1))/u(xi_) u(iw_mom(2))=u(lfac_)*u(iw_mom(2))/u(xi_) u(iw_mom(3))=u(lfac_)*u(iw_mom(3))/u(xi_) u(iw_e)=pth end subroutine to_primitive pure subroutine to_conservative(u) real(dp), intent(inout) :: u(nw_phys) real(dp) :: rho,rhoh,pth real(dp) :: E,E_th rhoh=(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2) u(lfac_)=dsqrt(1.0d0+rhoh) rho=u(iw_rho) pth=u(iw_e) !! begin: call srhd_get_enthalpy_eos(rho,pth,rhoh) if(srhd_eos) then E_th = pth*inv_gamma_1 E = E_th+dsqrt(E_th**2+rho**2) ! writing rho/E on purpose, for numerics rhoh = 0.5_dp*((srhd_gamma+1.0_dp)*E & - gamma_1*rho*(rho/E)) else rhoh = rho+gamma_to_gamma_1*pth end if !! end: call srhd_get_enthalpy_eos(rho,pth,rhoh) rhoh=rhoh*u(lfac_) u(xi_)=u(lfac_)*rhoh u(iw_rho)=u(lfac_)*rho u(iw_mom(1))=rhoh*u(iw_mom(1)) u(iw_mom(2))=rhoh*u(iw_mom(2)) u(iw_mom(3))=rhoh*u(iw_mom(3)) u(iw_e)=u(xi_)-pth-u(iw_rho) end subroutine to_conservative !> Make a reconstructed primitive face state admissible before it is handed !> to the Riemann solver. !> !> Two things happen here, and only the second is optional. The MUSCL !> reconstruction limits every slot of the state independently, including the !> auxiliaries xi and lfac -- but those are not independent variables, they !> are functions of (rho, u^i, p). A reconstructed face state therefore does !> not in general satisfy lfac = sqrt(1+|u|^2) or xi = lfac^2*rho*h, while !> get_flux, get_cmax and estimate_speeds_minmax all read them. Left alone, !> v2 = 1-1/lfac^2 and v_n = u_n/lfac come from different sources, v_n can !> exceed v (even exceed 1), the radicand of the relativistic characteristic !> speeds goes negative and the HLL wave speeds collapse onto each other -- !> which removes the solver's dissipation exactly where it is needed. So the !> auxiliaries are always recomputed from the reconstructed primitives, by !> the same expressions to_conservative uses. !> !> The floors are the optional part. Both default to zero, in which case the !> max() calls only clamp states that were already unphysical; set !> srhd_small_pressure (and srhd_small_density) positive in &srhd_list to get !> an actual floor, which is the useful configuration -- a state floored to !> exactly zero still gives csound2 = gamma*p/(rho*h) a zero denominator. pure subroutine fix_prim_state(u) real(dp), intent(inout) :: u(nw_phys) real(dp) :: rho, pth, rhoh, E, E_th rho = max(u(iw_rho), srhd_small_density) pth = max(u(iw_e), srhd_small_pressure) u(iw_rho) = rho u(iw_e) = pth u(lfac_) = dsqrt(1.0d0 + u(iw_mom(1))**2 + u(iw_mom(2))**2 & + u(iw_mom(3))**2) !! begin: call srhd_get_enthalpy_eos(rho,pth,rhoh) if (srhd_eos) then E_th = pth*inv_gamma_1 E = E_th + dsqrt(E_th**2 + rho**2) ! writing rho/E on purpose, for numerics rhoh = 0.5_dp*((srhd_gamma + 1.0_dp)*E - gamma_1*rho*(rho/E)) else rhoh = rho + gamma_to_gamma_1*pth end if !! end: call srhd_get_enthalpy_eos(rho,pth,rhoh) u(xi_) = u(lfac_)**2 * rhoh end subroutine fix_prim_state !> Returns maximum local signal speed from primitive state u in direction flux_dim; !> used in LLF/TVDLF flux estimation. pure real(dp) function get_cmax(u, x, flux_dim) result(wC) real(dp), intent(in) :: u(nw_phys) real(dp), intent(in) :: x(1:ndim) integer, intent(in) :: flux_dim real(dp) :: tmp1,tmp2,v2,csound2,vidim,cmax,cmin real(dp) :: rho,rhoh,pth,E ! input u is in primitive form ! auxiliaries are filled here rho=u(iw_rho) rhoh=u(xi_)/u(lfac_)**2 pth=u(iw_e) !!tmp1=u(iw_rho) !!tmp2=u(xi_)/u(lfac_)**2.0d0 !!v2=u(iw_e) !! begin: call srhd_get_csound2_prim_eos(tmp1,tmp2,v2,csound2) if(srhd_eos) then E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)& *rho**2))/(srhd_gamma+1.0_dp) csound2=(pth*((srhd_gamma+1.0_dp)& +gamma_1*(rho/E)**2))& /(2.0d0*rhoh) else csound2=srhd_gamma*pth/rhoh end if !! end call srhd_get_csound2_prim_eos(tmp1,tmp2,v2,csound2) v2 = 1.0d0-1.0d0/u(lfac_)**2 vidim = u(iw_mom(flux_dim))/u(lfac_) tmp2 = vidim**2 tmp1 = 1.0d0-v2*csound2-tmp2*(1.0d0-csound2) ! the radicand is positive analytically, as in estimate_speeds_minmax; clip ! it so an inadmissible state (c_s^2 > 1 from a collapsed rhoh) yields a ! degenerate speed rather than a NaN. get_cmax feeds setdt, where a single ! NaN survives the global MPI_ALLREDUCE and takes dt down for every rank. tmp2 = dsqrt(max(csound2*(1.0_dp-v2)*tmp1, 0.0_dp)) tmp1 = vidim*(1.0_dp-csound2) cmax = (tmp1+tmp2)/(1.0_dp-v2*csound2) cmin = (tmp1-tmp2)/(1.0_dp-v2*csound2) ! Limit by speed of light cmin = max(cmin, - 1.0d0) cmin = min(cmin, 1.0d0) cmax = max(cmax, - 1.0d0) cmax = min(cmax, 1.0d0) ! now take extremal value only for dt limit wC = max(dabs(cmax),dabs(cmin)) end function get_cmax subroutine get_flux(u, xC, flux_dim, flux) use mod_global_parameters, only:cmax_global real(dp), intent(in) :: u(nw_phys) real(dp), intent(in) :: xC(1:ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: flux(nw_flux) real(dp) :: pth,vel(1:ndim) pth = u(iw_e) vel(1) = u(iw_mom(1)) / u(lfac_) vel(2) = u(iw_mom(2)) / u(lfac_) vel(3) = u(iw_mom(3)) / u(lfac_) ! Density flux flux(iw_rho) = u(iw_rho) * u(lfac_) * vel(flux_dim) ! Momentum flux with pressure term flux(iw_mom(1)) = u(xi_) * vel(1) * vel(flux_dim) flux(iw_mom(2)) = u(xi_) * vel(2) * vel(flux_dim) flux(iw_mom(3)) = u(xi_) * vel(3) * vel(flux_dim) flux(iw_mom(flux_dim)) = flux(iw_mom(flux_dim)) + pth ! Energy flux flux(iw_e) = vel(flux_dim) * ( u(xi_) - u(iw_rho) * u(lfac_) ) ! Tracer flux. Note that tracers stay conservative. flux(tracer(1)) = u(tracer(1)) * vel(flux_dim) end subroutine get_flux !> Davis (1988) min/max wave speed estimates for HLL, in their relativistic !> form: wL = min(lambda^-), wR = max(lambda^+) over the left and right states. !> !> The relativistic characteristic speeds in direction n are !> !> lambda^+- = ( v_n (1 - c_s^2) !> +- sqrt( c_s^2 (1 - v^2) (1 - v^2 c_s^2 - v_n^2 (1 - c_s^2)) ) ) !> / (1 - v^2 c_s^2) !> !> (Anile 1989; Mignone & Bodo 2005), which reduce to v_n +- c_s in the !> Newtonian limit. Unlike the HD and MHD versions of this routine, the !> velocity is already contained in lambda^+-, so the min/max is taken over the !> characteristic speeds themselves rather than over v_n +- c. !> !> This is the same algebra get_cmax uses; get_cmax collapses the two roots !> into max(|lambda^+|, |lambda^-|) for the time-step limit, whereas HLL needs !> them separately and signed. wprim(iw_mom(:)) holds the spatial four-velocity !> u^i = lfac*v^i, so v^i is recovered as u^i/lfac. subroutine estimate_speeds_minmax(uL, uR, xC, flux_dim, wL, wR) real(dp), intent(in) :: uL(nw_phys), uR(nw_phys) real(dp), intent(in) :: xC(ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: wL, wR real(dp) :: rho, rhoh, pth, E, csound2 real(dp) :: v2, vidim, root, denom real(dp) :: cminL, cmaxL, cminR, cmaxR ! Left state rho = uL(iw_rho) rhoh = uL(xi_) / uL(lfac_)**2 pth = uL(iw_e) if (srhd_eos) then E = (rhoh + dsqrt(rhoh**2 + (srhd_gamma**2 - 1.0_dp)*rho**2)) & /(srhd_gamma + 1.0_dp) csound2 = (pth*((srhd_gamma + 1.0_dp) + gamma_1*(rho/E)**2))/(2.0_dp*rhoh) else csound2 = srhd_gamma*pth/rhoh end if v2 = 1.0_dp - 1.0_dp/uL(lfac_)**2 vidim = uL(iw_mom(flux_dim))/uL(lfac_) denom = 1.0_dp - v2*csound2 ! the radicand is positive analytically; clip the roundoff-level negatives ! that would otherwise produce a NaN wave speed root = dsqrt(max(csound2*(1.0_dp - v2) & *(1.0_dp - v2*csound2 - vidim**2*(1.0_dp - csound2)), 0.0_dp)) cmaxL = (vidim*(1.0_dp - csound2) + root)/denom cminL = (vidim*(1.0_dp - csound2) - root)/denom ! Right state rho = uR(iw_rho) rhoh = uR(xi_) / uR(lfac_)**2 pth = uR(iw_e) if (srhd_eos) then E = (rhoh + dsqrt(rhoh**2 + (srhd_gamma**2 - 1.0_dp)*rho**2)) & /(srhd_gamma + 1.0_dp) csound2 = (pth*((srhd_gamma + 1.0_dp) + gamma_1*(rho/E)**2))/(2.0_dp*rhoh) else csound2 = srhd_gamma*pth/rhoh end if v2 = 1.0_dp - 1.0_dp/uR(lfac_)**2 vidim = uR(iw_mom(flux_dim))/uR(lfac_) denom = 1.0_dp - v2*csound2 root = dsqrt(max(csound2*(1.0_dp - v2) & *(1.0_dp - v2*csound2 - vidim**2*(1.0_dp - csound2)), 0.0_dp)) cmaxR = (vidim*(1.0_dp - csound2) + root)/denom cminR = (vidim*(1.0_dp - csound2) - root)/denom wL = min(cminL, cminR) wR = max(cmaxL, cmaxR) ! no signal can outrun light wL = min(max(wL, -1.0_dp), 1.0_dp) wR = min(max(wR, -1.0_dp), 1.0_dp) end subroutine estimate_speeds_minmax subroutine estimate_speeds_toro_pvrs(uL, uR, xC, flux_dim, sL, sR) #ifndef _CRAYFTN use mod_comm_lib, only: mpistop #endif real(dp), intent(in) :: uL(nw_phys), uR(nw_phys) real(dp), intent(in) :: xC(ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: sL, sR #ifndef _CRAYFTN call mpistop("estimate_speeds_toro_pvrs not implemented for this physics module") #endif sL = -1._dp sR = -1._dp end subroutine estimate_speeds_toro_pvrs ! flux scheme list : (scheme_tag, method_enum, faceflux_proc) subroutine finite_volume_local(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) use mod_global_parameters use mod_comm_lib, only: mpistop double precision, intent(in) :: qdt,& dtfactor, qtC, qt integer, intent(in) :: & ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,& ixOmin3,ixOmax1,ixOmax2,ixOmax3, idimsmin,idimsmax ! remember, old names map as: wCT => bga, wnew => bgb type(block_grid_t) :: bga type(block_grid_t) :: bgb ! .. local .. integer :: n, iigrid, ix1,ix2,ix3 double precision :: uprim(nw_phys, ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) real(dp) :: tmp(nw_phys,5) real(dp) :: f(nw_flux, 2) real(dp) :: inv_dr(ndim) real(dp) :: dr(ndim) integer :: typelim real(dp) :: xloc(ndim) real(dp) :: xlocC(ndim,2) real(dp) :: wprim(nw_phys), wCT(nw_phys), wnew(nw_phys) integer :: method !----------------------------------------------------------------------------- method = flux_method(1) ! TODO: implement per grid level select case (method) case (fs_tvdlf) call finite_volume_local_muscl_llf(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) case (fs_hll) call finite_volume_local_muscl_hll(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) case (fs_hllc) call finite_volume_local_muscl_hllc(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) case default call mpistop("finite_volume_local: unknown flux scheme") end select end subroutine finite_volume_local subroutine finite_volume_local_muscl_llf(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) use mod_global_parameters use mod_fix_conserve double precision, intent(in) :: qdt,& dtfactor, qtC, qt integer, intent(in) :: & ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,& ixOmin3,ixOmax1,ixOmax2,ixOmax3, idimsmin,idimsmax ! remember, old names map as: wCT => bga, wnew => bgb type(block_grid_t) :: bga type(block_grid_t) :: bgb ! .. local .. integer :: n, iigrid, ix1,ix2,ix3 integer :: nx1,nx2,nx3, nxCo1,nxCo2,nxCo3 double precision :: uprim(nw_phys, ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) real(dp) :: tmp(nw_phys,5) real(dp) :: tmp1(nw_phys,3),tmp2(nw_phys,3),tmp3(nw_phys,3) real(dp) :: f(nw_flux, 2) real(dp) :: inv_dr(ndim) real(dp) :: dr(ndim) integer :: typelim real(dp) :: xloc(ndim) real(dp) :: xlocC(ndim,2) real(dp) :: wprim(nw_phys), wCT(nw_phys), wnew(nw_phys) ! Block corner in the logical coordinate, for building face positions. real(dp) :: xlo(ndim) ! Physical cell size handed to the optional source terms. dr below is the ! *logical* spacing, which is not a length in a curvilinear system and is ! not even proportional to one under LOG_RADIUS. real(dp) :: dloc(ndim) real(dp) :: fC1(2,ixImin2:ixOmax2,ixImin3:ixOmax3,1:nwflux) real(dp) :: fC2(2,ixImin1:ixOmax1,ixImin3:ixOmax3,1:nwflux) real(dp) :: fC3(2,ixImin1:ixOmax1,ixImin2:ixOmax2,1:nwflux) integer, parameter :: max_batch=4096 integer :: nbatches, ibatch, igrid_beg, igrid_end integer :: neighbor_type_1m, neighbor_type_1p, neighbor_type_2m,& neighbor_type_2p, neighbor_type_3m, neighbor_type_3p !----------------------------------------------------------------------------- nx1=ixMhi1-ixMlo1+1 nx2=ixMhi2-ixMlo2+1 nx3=ixMhi3-ixMlo3+1 nxCo1=nx1/2 nxCo2=nx2/2 nxCo3=nx3/2 ! batch-launch the kernels (oom issue when many ~30000 blocks): nbatches = (igridstail_active + max_batch - 1) / max_batch ! ceiling division do ibatch = 1, nbatches igrid_beg = (ibatch-1) * max_batch + 1 igrid_end = min(ibatch * max_batch, igridstail_active) do iigrid = igrid_beg, igrid_end n = igrids_active(iigrid) dr = rnode(rpdx1_:rnodehi, n) inv_dr = 1/dr xlo = rnode(rpxmin1_:rpxmin1_+ndim-1, n) typelim = type_limiter(node(plevel_, n)) !JESSENEW These were called many times in the 3D loop below !but they technically only have an igrid dependence so I !thought it would be better to simply declare them via local !variable which will reduce the call overhead ~NX^3 neighbor_type_1m = neighbor_type(-1,0,0,n) neighbor_type_1p = neighbor_type( 1,0,0,n) neighbor_type_2m = neighbor_type(0,-1,0,n) neighbor_type_2p = neighbor_type(0, 1,0,n) neighbor_type_3m = neighbor_type(0,0,-1,n) neighbor_type_3p = neighbor_type(0,0, 1,n) do ix3=ixImin3,ixImax3 do ix2=ixImin2,ixImax2 do ix1=ixImin1,ixImax1 ! Convert to primitive uprim(1:nw_phys, ix1,ix2,ix3) = bga%w(ix1,ix2,ix3, 1:nw_phys, n) call to_primitive(uprim(1:nw_phys, ix1,ix2,ix3)) end do end do end do do ix3=ixOmin3,ixOmax3 do ix2=ixOmin2,ixOmax2 do ix1=ixOmin1,ixOmax1 ! Compute fluxes in all dimensions ! This cell's stored position - its volume barycentre - read ! once and reused by all three face constructions and by every ! source term below. The component index of bgeo%x sits between ! the cell indices and the grid index, so each of these is a ! strided gather; fetching it once rather than per direction ! takes ten of them per cell down to one. xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, n) tmp = uprim(1:nw_phys, ix1-2:ix1+2, ix2, ix3) ! transverse components stay on the barycentre; only the normal ! one is moved to the interface xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) ! Face positions, built from the block corner and the cell ! index. This is the one form that is right in every geometry, ! so it is used unconditionally. It cannot be shortened to ! x -/+ dr/2: bgeo%x is the cell's volume *barycentre* in a ! curvilinear build, and under LOG_RADIUS dr(1) is a ratio in ! ln(r) rather than a length, so that form would not merely be ! imprecise but wrong. Where x is the midpoint - Cartesian in ! every direction, and phi and cylindrical z besides - the two ! agree to round-off. ! ! The lower-corner form suffices because this loop covers the ! mesh interior only; fill_geometry_device's two-sided rule ! exists to make neighbouring blocks' *ghost* cells agree ! bit-for-bit, which nothing here needs. xlocC(1,1) = xlo(1) + dble(ix1-nghostcells-1)*dr(1) xlocC(1,2) = xlocC(1,1) + dr(1) call reconflux_muscl_llf_prim(tmp, xlocC, 1, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(1) ! Store fluxes for flux fixing in direction 1. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_1m) case (neighbor_fine) if (ix1.eq.ixOmin1) pflux(1,1)%flux(1,ix2-nghostcells,ix3-nghostcells,1:nw_flux,n) & = qdt * inv_dr(1) * f(:,1) case (neighbor_coarse) if (ix1.eq.ixOmin1) fC1(1,ix2,ix3,1:nw_flux) = - qdt * inv_dr(1) * f(:,1) end select select case (neighbor_type_1p) case (neighbor_fine) if (ix1.eq.ixOmax1) pflux(2,1)%flux(1,ix2-nghostcells,ix3-nghostcells,1:nw_flux,n) & = - qdt * inv_dr(1) * f(:,2) case (neighbor_coarse) if (ix1.eq.ixOmax1) fC1(2,ix2,ix3,1:nw_flux) = qdt * inv_dr(1) * f(:,2) end select tmp = uprim(1:nw_phys, ix1, ix2-2:ix2+2, ix3) xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) xlocC(2,1) = xlo(2) + dble(ix2-nghostcells-1)*dr(2) xlocC(2,2) = xlocC(2,1) + dr(2) call reconflux_muscl_llf_prim(tmp, xlocC, 2, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(2) ! Store fluxes for flux fixing in direction 2. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_2m) case (neighbor_fine) if (ix2.eq.ixOmin2) pflux(1,2)%flux(ix1-nghostcells,1,ix3-nghostcells,1:nw_flux,n) & = qdt * inv_dr(2) * f(:,1) case (neighbor_coarse) if (ix2.eq.ixOmin2) fC2(1,ix1,ix3,1:nw_flux) = - qdt * inv_dr(2) * f(:,1) end select select case (neighbor_type_2p) case (neighbor_fine) if (ix2.eq.ixOmax2) pflux(2,2)%flux(ix1-nghostcells,1,ix3-nghostcells,1:nw_flux,n) & = - qdt * inv_dr(2) * f(:,2) case (neighbor_coarse) if (ix2.eq.ixOmax2) fC2(2,ix1,ix3,1:nw_flux) = qdt * inv_dr(2) * f(:,2) end select tmp = uprim(1:nw_phys, ix1, ix2, ix3-2:ix3+2) xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) xlocC(3,1) = xlo(3) + dble(ix3-nghostcells-1)*dr(3) xlocC(3,2) = xlocC(3,1) + dr(3) call reconflux_muscl_llf_prim(tmp, xlocC, 3, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(3) ! Store fluxes for flux fixing in direction 3. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_3m) case (neighbor_fine) if (ix3.eq.ixOmin3) pflux(1,3)%flux(ix1-nghostcells,ix2-nghostcells,1,1:nw_flux,n) & = qdt * inv_dr(3) * f(:,1) case (neighbor_coarse) if (ix3.eq.ixOmin3) fC3(1,ix1,ix2,1:nw_flux) = - qdt * inv_dr(3) * f(:,1) end select select case (neighbor_type_3p) case (neighbor_fine) if (ix3.eq.ixOmax3) pflux(2,3)%flux(ix1-nghostcells,ix2-nghostcells,1,1:nw_flux,n) & = - qdt * inv_dr(3) * f(:,2) case (neighbor_coarse) if (ix3.eq.ixOmax3) fC3(2,ix1,ix2,1:nw_flux) = qdt * inv_dr(3) * f(:,2) end select ! The optional source terms want a physical cell size. In a ! curvilinear build that is bgeo%ds, not the logical spacing. dloc = dr ! Add local source terms: wprim = uprim(1:nw_phys, ix1, ix2, ix3) wCT = bga%w(ix1, ix2, ix3, 1:nw_phys, n) wnew = bgb%w(ix1, ix2, ix3, 1:nw_phys, n) call addsource_local(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), qtC, wCT,& wprim, qt, wnew, xloc, dloc, .false. ) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = wnew(1:nw_flux) end do end do end do ! Reduce fluxes to be stored for the coarse neighbor case ! Direction 1 select case (neighbor_type_1m) case (neighbor_coarse) do ix3=1,nxCo3 do ix2=1,nxCo2 pflux(1,1)%flux(1,ix2,ix3,1:nw_flux,n) = & fC1(1,1+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,2+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,1+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,2+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_1p) case (neighbor_coarse) do ix3=1,nxCo3 do ix2=1,nxCo2 pflux(2,1)%flux(1,ix2,ix3,1:nw_flux,n) = & fC1(2,1+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,2+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,1+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,2+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select ! Direction 2 select case (neighbor_type_2m) case (neighbor_coarse) do ix3=1,nxCo3 do ix1=1,nxCo1 pflux(1,2)%flux(ix1,1,ix3,1:nw_flux,n) = & fC2(1,1+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,2+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,1+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,2+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_2p) case (neighbor_coarse) do ix3=1,nxCo3 do ix1=1,nxCo1 pflux(2,2)%flux(ix1,1,ix3,1:nw_flux,n) = & fC2(2,1+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,2+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,1+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,2+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select ! Direction 3 select case (neighbor_type_3m) case (neighbor_coarse) do ix2=1,nxCo2 do ix1=1,nxCo1 pflux(1,3)%flux(ix1,ix2,1,1:nw_flux,n) = & fC3(1,1+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,2+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,1+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,2+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_3p) case (neighbor_coarse) do ix2=1,nxCo2 do ix1=1,nxCo1 pflux(2,3)%flux(ix1,ix2,1,1:nw_flux,n) = & fC3(2,1+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,2+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,1+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,2+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) end do end do end select end do end do end subroutine finite_volume_local_muscl_llf subroutine finite_volume_local_muscl_hll(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) use mod_global_parameters use mod_fix_conserve double precision, intent(in) :: qdt,& dtfactor, qtC, qt integer, intent(in) :: & ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,& ixOmin3,ixOmax1,ixOmax2,ixOmax3, idimsmin,idimsmax ! remember, old names map as: wCT => bga, wnew => bgb type(block_grid_t) :: bga type(block_grid_t) :: bgb ! .. local .. integer :: n, iigrid, ix1,ix2,ix3 integer :: nx1,nx2,nx3, nxCo1,nxCo2,nxCo3 double precision :: uprim(nw_phys, ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) real(dp) :: tmp(nw_phys,5) real(dp) :: tmp1(nw_phys,3),tmp2(nw_phys,3),tmp3(nw_phys,3) real(dp) :: f(nw_flux, 2) real(dp) :: inv_dr(ndim) real(dp) :: dr(ndim) integer :: typelim real(dp) :: xloc(ndim) real(dp) :: xlocC(ndim,2) real(dp) :: wprim(nw_phys), wCT(nw_phys), wnew(nw_phys) ! Block corner in the logical coordinate, for building face positions. real(dp) :: xlo(ndim) ! Physical cell size handed to the optional source terms. dr below is the ! *logical* spacing, which is not a length in a curvilinear system and is ! not even proportional to one under LOG_RADIUS. real(dp) :: dloc(ndim) real(dp) :: fC1(2,ixImin2:ixOmax2,ixImin3:ixOmax3,1:nwflux) real(dp) :: fC2(2,ixImin1:ixOmax1,ixImin3:ixOmax3,1:nwflux) real(dp) :: fC3(2,ixImin1:ixOmax1,ixImin2:ixOmax2,1:nwflux) integer, parameter :: max_batch=4096 integer :: nbatches, ibatch, igrid_beg, igrid_end integer :: neighbor_type_1m, neighbor_type_1p, neighbor_type_2m,& neighbor_type_2p, neighbor_type_3m, neighbor_type_3p !----------------------------------------------------------------------------- nx1=ixMhi1-ixMlo1+1 nx2=ixMhi2-ixMlo2+1 nx3=ixMhi3-ixMlo3+1 nxCo1=nx1/2 nxCo2=nx2/2 nxCo3=nx3/2 ! batch-launch the kernels (oom issue when many ~30000 blocks): nbatches = (igridstail_active + max_batch - 1) / max_batch ! ceiling division do ibatch = 1, nbatches igrid_beg = (ibatch-1) * max_batch + 1 igrid_end = min(ibatch * max_batch, igridstail_active) do iigrid = igrid_beg, igrid_end n = igrids_active(iigrid) dr = rnode(rpdx1_:rnodehi, n) inv_dr = 1/dr xlo = rnode(rpxmin1_:rpxmin1_+ndim-1, n) typelim = type_limiter(node(plevel_, n)) !JESSENEW These were called many times in the 3D loop below !but they technically only have an igrid dependence so I !thought it would be better to simply declare them via local !variable which will reduce the call overhead ~NX^3 neighbor_type_1m = neighbor_type(-1,0,0,n) neighbor_type_1p = neighbor_type( 1,0,0,n) neighbor_type_2m = neighbor_type(0,-1,0,n) neighbor_type_2p = neighbor_type(0, 1,0,n) neighbor_type_3m = neighbor_type(0,0,-1,n) neighbor_type_3p = neighbor_type(0,0, 1,n) do ix3=ixImin3,ixImax3 do ix2=ixImin2,ixImax2 do ix1=ixImin1,ixImax1 ! Convert to primitive uprim(1:nw_phys, ix1,ix2,ix3) = bga%w(ix1,ix2,ix3, 1:nw_phys, n) call to_primitive(uprim(1:nw_phys, ix1,ix2,ix3)) end do end do end do do ix3=ixOmin3,ixOmax3 do ix2=ixOmin2,ixOmax2 do ix1=ixOmin1,ixOmax1 ! Compute fluxes in all dimensions ! This cell's stored position - its volume barycentre - read ! once and reused by all three face constructions and by every ! source term below. The component index of bgeo%x sits between ! the cell indices and the grid index, so each of these is a ! strided gather; fetching it once rather than per direction ! takes ten of them per cell down to one. xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, n) tmp = uprim(1:nw_phys, ix1-2:ix1+2, ix2, ix3) ! transverse components stay on the barycentre; only the normal ! one is moved to the interface xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) ! Face positions, built from the block corner and the cell ! index. This is the one form that is right in every geometry, ! so it is used unconditionally. It cannot be shortened to ! x -/+ dr/2: bgeo%x is the cell's volume *barycentre* in a ! curvilinear build, and under LOG_RADIUS dr(1) is a ratio in ! ln(r) rather than a length, so that form would not merely be ! imprecise but wrong. Where x is the midpoint - Cartesian in ! every direction, and phi and cylindrical z besides - the two ! agree to round-off. ! ! The lower-corner form suffices because this loop covers the ! mesh interior only; fill_geometry_device's two-sided rule ! exists to make neighbouring blocks' *ghost* cells agree ! bit-for-bit, which nothing here needs. xlocC(1,1) = xlo(1) + dble(ix1-nghostcells-1)*dr(1) xlocC(1,2) = xlocC(1,1) + dr(1) call reconflux_muscl_hll_prim(tmp, xlocC, 1, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(1) ! Store fluxes for flux fixing in direction 1. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_1m) case (neighbor_fine) if (ix1.eq.ixOmin1) pflux(1,1)%flux(1,ix2-nghostcells,ix3-nghostcells,1:nw_flux,n) & = qdt * inv_dr(1) * f(:,1) case (neighbor_coarse) if (ix1.eq.ixOmin1) fC1(1,ix2,ix3,1:nw_flux) = - qdt * inv_dr(1) * f(:,1) end select select case (neighbor_type_1p) case (neighbor_fine) if (ix1.eq.ixOmax1) pflux(2,1)%flux(1,ix2-nghostcells,ix3-nghostcells,1:nw_flux,n) & = - qdt * inv_dr(1) * f(:,2) case (neighbor_coarse) if (ix1.eq.ixOmax1) fC1(2,ix2,ix3,1:nw_flux) = qdt * inv_dr(1) * f(:,2) end select tmp = uprim(1:nw_phys, ix1, ix2-2:ix2+2, ix3) xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) xlocC(2,1) = xlo(2) + dble(ix2-nghostcells-1)*dr(2) xlocC(2,2) = xlocC(2,1) + dr(2) call reconflux_muscl_hll_prim(tmp, xlocC, 2, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(2) ! Store fluxes for flux fixing in direction 2. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_2m) case (neighbor_fine) if (ix2.eq.ixOmin2) pflux(1,2)%flux(ix1-nghostcells,1,ix3-nghostcells,1:nw_flux,n) & = qdt * inv_dr(2) * f(:,1) case (neighbor_coarse) if (ix2.eq.ixOmin2) fC2(1,ix1,ix3,1:nw_flux) = - qdt * inv_dr(2) * f(:,1) end select select case (neighbor_type_2p) case (neighbor_fine) if (ix2.eq.ixOmax2) pflux(2,2)%flux(ix1-nghostcells,1,ix3-nghostcells,1:nw_flux,n) & = - qdt * inv_dr(2) * f(:,2) case (neighbor_coarse) if (ix2.eq.ixOmax2) fC2(2,ix1,ix3,1:nw_flux) = qdt * inv_dr(2) * f(:,2) end select tmp = uprim(1:nw_phys, ix1, ix2, ix3-2:ix3+2) xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) xlocC(3,1) = xlo(3) + dble(ix3-nghostcells-1)*dr(3) xlocC(3,2) = xlocC(3,1) + dr(3) call reconflux_muscl_hll_prim(tmp, xlocC, 3, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(3) ! Store fluxes for flux fixing in direction 3. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_3m) case (neighbor_fine) if (ix3.eq.ixOmin3) pflux(1,3)%flux(ix1-nghostcells,ix2-nghostcells,1,1:nw_flux,n) & = qdt * inv_dr(3) * f(:,1) case (neighbor_coarse) if (ix3.eq.ixOmin3) fC3(1,ix1,ix2,1:nw_flux) = - qdt * inv_dr(3) * f(:,1) end select select case (neighbor_type_3p) case (neighbor_fine) if (ix3.eq.ixOmax3) pflux(2,3)%flux(ix1-nghostcells,ix2-nghostcells,1,1:nw_flux,n) & = - qdt * inv_dr(3) * f(:,2) case (neighbor_coarse) if (ix3.eq.ixOmax3) fC3(2,ix1,ix2,1:nw_flux) = qdt * inv_dr(3) * f(:,2) end select ! The optional source terms want a physical cell size. In a ! curvilinear build that is bgeo%ds, not the logical spacing. dloc = dr ! Add local source terms: wprim = uprim(1:nw_phys, ix1, ix2, ix3) wCT = bga%w(ix1, ix2, ix3, 1:nw_phys, n) wnew = bgb%w(ix1, ix2, ix3, 1:nw_phys, n) call addsource_local(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), qtC, wCT,& wprim, qt, wnew, xloc, dloc, .false. ) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = wnew(1:nw_flux) end do end do end do ! Reduce fluxes to be stored for the coarse neighbor case ! Direction 1 select case (neighbor_type_1m) case (neighbor_coarse) do ix3=1,nxCo3 do ix2=1,nxCo2 pflux(1,1)%flux(1,ix2,ix3,1:nw_flux,n) = & fC1(1,1+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,2+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,1+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,2+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_1p) case (neighbor_coarse) do ix3=1,nxCo3 do ix2=1,nxCo2 pflux(2,1)%flux(1,ix2,ix3,1:nw_flux,n) = & fC1(2,1+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,2+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,1+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,2+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select ! Direction 2 select case (neighbor_type_2m) case (neighbor_coarse) do ix3=1,nxCo3 do ix1=1,nxCo1 pflux(1,2)%flux(ix1,1,ix3,1:nw_flux,n) = & fC2(1,1+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,2+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,1+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,2+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_2p) case (neighbor_coarse) do ix3=1,nxCo3 do ix1=1,nxCo1 pflux(2,2)%flux(ix1,1,ix3,1:nw_flux,n) = & fC2(2,1+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,2+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,1+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,2+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select ! Direction 3 select case (neighbor_type_3m) case (neighbor_coarse) do ix2=1,nxCo2 do ix1=1,nxCo1 pflux(1,3)%flux(ix1,ix2,1,1:nw_flux,n) = & fC3(1,1+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,2+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,1+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,2+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_3p) case (neighbor_coarse) do ix2=1,nxCo2 do ix1=1,nxCo1 pflux(2,3)%flux(ix1,ix2,1,1:nw_flux,n) = & fC3(2,1+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,2+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,1+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,2+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) end do end do end select end do end do end subroutine finite_volume_local_muscl_hll subroutine finite_volume_local_muscl_hllc(qdt, dtfactor, ixImin1,ixImin2,& ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,& ixOmax3, idimsmin,idimsmax, qtC, bga, qt, bgb) use mod_global_parameters use mod_fix_conserve double precision, intent(in) :: qdt,& dtfactor, qtC, qt integer, intent(in) :: & ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3, ixOmin1,ixOmin2,& ixOmin3,ixOmax1,ixOmax2,ixOmax3, idimsmin,idimsmax ! remember, old names map as: wCT => bga, wnew => bgb type(block_grid_t) :: bga type(block_grid_t) :: bgb ! .. local .. integer :: n, iigrid, ix1,ix2,ix3 integer :: nx1,nx2,nx3, nxCo1,nxCo2,nxCo3 double precision :: uprim(nw_phys, ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) real(dp) :: tmp(nw_phys,5) real(dp) :: tmp1(nw_phys,3),tmp2(nw_phys,3),tmp3(nw_phys,3) real(dp) :: f(nw_flux, 2) real(dp) :: inv_dr(ndim) real(dp) :: dr(ndim) integer :: typelim real(dp) :: xloc(ndim) real(dp) :: xlocC(ndim,2) real(dp) :: wprim(nw_phys), wCT(nw_phys), wnew(nw_phys) ! Block corner in the logical coordinate, for building face positions. real(dp) :: xlo(ndim) ! Physical cell size handed to the optional source terms. dr below is the ! *logical* spacing, which is not a length in a curvilinear system and is ! not even proportional to one under LOG_RADIUS. real(dp) :: dloc(ndim) real(dp) :: fC1(2,ixImin2:ixOmax2,ixImin3:ixOmax3,1:nwflux) real(dp) :: fC2(2,ixImin1:ixOmax1,ixImin3:ixOmax3,1:nwflux) real(dp) :: fC3(2,ixImin1:ixOmax1,ixImin2:ixOmax2,1:nwflux) integer, parameter :: max_batch=4096 integer :: nbatches, ibatch, igrid_beg, igrid_end integer :: neighbor_type_1m, neighbor_type_1p, neighbor_type_2m,& neighbor_type_2p, neighbor_type_3m, neighbor_type_3p !----------------------------------------------------------------------------- nx1=ixMhi1-ixMlo1+1 nx2=ixMhi2-ixMlo2+1 nx3=ixMhi3-ixMlo3+1 nxCo1=nx1/2 nxCo2=nx2/2 nxCo3=nx3/2 ! batch-launch the kernels (oom issue when many ~30000 blocks): nbatches = (igridstail_active + max_batch - 1) / max_batch ! ceiling division do ibatch = 1, nbatches igrid_beg = (ibatch-1) * max_batch + 1 igrid_end = min(ibatch * max_batch, igridstail_active) do iigrid = igrid_beg, igrid_end n = igrids_active(iigrid) dr = rnode(rpdx1_:rnodehi, n) inv_dr = 1/dr xlo = rnode(rpxmin1_:rpxmin1_+ndim-1, n) typelim = type_limiter(node(plevel_, n)) !JESSENEW These were called many times in the 3D loop below !but they technically only have an igrid dependence so I !thought it would be better to simply declare them via local !variable which will reduce the call overhead ~NX^3 neighbor_type_1m = neighbor_type(-1,0,0,n) neighbor_type_1p = neighbor_type( 1,0,0,n) neighbor_type_2m = neighbor_type(0,-1,0,n) neighbor_type_2p = neighbor_type(0, 1,0,n) neighbor_type_3m = neighbor_type(0,0,-1,n) neighbor_type_3p = neighbor_type(0,0, 1,n) do ix3=ixImin3,ixImax3 do ix2=ixImin2,ixImax2 do ix1=ixImin1,ixImax1 ! Convert to primitive uprim(1:nw_phys, ix1,ix2,ix3) = bga%w(ix1,ix2,ix3, 1:nw_phys, n) call to_primitive(uprim(1:nw_phys, ix1,ix2,ix3)) end do end do end do do ix3=ixOmin3,ixOmax3 do ix2=ixOmin2,ixOmax2 do ix1=ixOmin1,ixOmax1 ! Compute fluxes in all dimensions ! This cell's stored position - its volume barycentre - read ! once and reused by all three face constructions and by every ! source term below. The component index of bgeo%x sits between ! the cell indices and the grid index, so each of these is a ! strided gather; fetching it once rather than per direction ! takes ten of them per cell down to one. xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, n) tmp = uprim(1:nw_phys, ix1-2:ix1+2, ix2, ix3) ! transverse components stay on the barycentre; only the normal ! one is moved to the interface xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) ! Face positions, built from the block corner and the cell ! index. This is the one form that is right in every geometry, ! so it is used unconditionally. It cannot be shortened to ! x -/+ dr/2: bgeo%x is the cell's volume *barycentre* in a ! curvilinear build, and under LOG_RADIUS dr(1) is a ratio in ! ln(r) rather than a length, so that form would not merely be ! imprecise but wrong. Where x is the midpoint - Cartesian in ! every direction, and phi and cylindrical z besides - the two ! agree to round-off. ! ! The lower-corner form suffices because this loop covers the ! mesh interior only; fill_geometry_device's two-sided rule ! exists to make neighbouring blocks' *ghost* cells agree ! bit-for-bit, which nothing here needs. xlocC(1,1) = xlo(1) + dble(ix1-nghostcells-1)*dr(1) xlocC(1,2) = xlocC(1,1) + dr(1) call reconflux_muscl_hllc_prim(tmp, xlocC, 1, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(1) ! Store fluxes for flux fixing in direction 1. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_1m) case (neighbor_fine) if (ix1.eq.ixOmin1) pflux(1,1)%flux(1,ix2-nghostcells,ix3-nghostcells,1:nw_flux,n) & = qdt * inv_dr(1) * f(:,1) case (neighbor_coarse) if (ix1.eq.ixOmin1) fC1(1,ix2,ix3,1:nw_flux) = - qdt * inv_dr(1) * f(:,1) end select select case (neighbor_type_1p) case (neighbor_fine) if (ix1.eq.ixOmax1) pflux(2,1)%flux(1,ix2-nghostcells,ix3-nghostcells,1:nw_flux,n) & = - qdt * inv_dr(1) * f(:,2) case (neighbor_coarse) if (ix1.eq.ixOmax1) fC1(2,ix2,ix3,1:nw_flux) = qdt * inv_dr(1) * f(:,2) end select tmp = uprim(1:nw_phys, ix1, ix2-2:ix2+2, ix3) xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) xlocC(2,1) = xlo(2) + dble(ix2-nghostcells-1)*dr(2) xlocC(2,2) = xlocC(2,1) + dr(2) call reconflux_muscl_hllc_prim(tmp, xlocC, 2, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(2) ! Store fluxes for flux fixing in direction 2. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_2m) case (neighbor_fine) if (ix2.eq.ixOmin2) pflux(1,2)%flux(ix1-nghostcells,1,ix3-nghostcells,1:nw_flux,n) & = qdt * inv_dr(2) * f(:,1) case (neighbor_coarse) if (ix2.eq.ixOmin2) fC2(1,ix1,ix3,1:nw_flux) = - qdt * inv_dr(2) * f(:,1) end select select case (neighbor_type_2p) case (neighbor_fine) if (ix2.eq.ixOmax2) pflux(2,2)%flux(ix1-nghostcells,1,ix3-nghostcells,1:nw_flux,n) & = - qdt * inv_dr(2) * f(:,2) case (neighbor_coarse) if (ix2.eq.ixOmax2) fC2(2,ix1,ix3,1:nw_flux) = qdt * inv_dr(2) * f(:,2) end select tmp = uprim(1:nw_phys, ix1, ix2, ix3-2:ix3+2) xlocC(1:ndim,1) = xloc(1:ndim) xlocC(1:ndim,2) = xloc(1:ndim) xlocC(3,1) = xlo(3) + dble(ix3-nghostcells-1)*dr(3) xlocC(3,2) = xlocC(3,1) + dr(3) call reconflux_muscl_hllc_prim(tmp, xlocC, 3, f, typelim) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = bgb%w(ix1, ix2, ix3, 1:nw_flux,& n) + qdt * (f(:, 1) - f(:, 2)) * inv_dr(3) ! Store fluxes for flux fixing in direction 3. What is ! stored differs by geometry: qdt*f/dx, the face's own ! contribution to the cell update, on a Cartesian mesh, and the ! extensive qdt*f*A on a curvilinear one. The extensive form is ! what lets fix_conserve divide by the *coarse* cell's volume ! and sum the four fine faces without a volume ratio, since the ! fine face areas add up to the coarse one exactly. select case (neighbor_type_3m) case (neighbor_fine) if (ix3.eq.ixOmin3) pflux(1,3)%flux(ix1-nghostcells,ix2-nghostcells,1,1:nw_flux,n) & = qdt * inv_dr(3) * f(:,1) case (neighbor_coarse) if (ix3.eq.ixOmin3) fC3(1,ix1,ix2,1:nw_flux) = - qdt * inv_dr(3) * f(:,1) end select select case (neighbor_type_3p) case (neighbor_fine) if (ix3.eq.ixOmax3) pflux(2,3)%flux(ix1-nghostcells,ix2-nghostcells,1,1:nw_flux,n) & = - qdt * inv_dr(3) * f(:,2) case (neighbor_coarse) if (ix3.eq.ixOmax3) fC3(2,ix1,ix2,1:nw_flux) = qdt * inv_dr(3) * f(:,2) end select ! The optional source terms want a physical cell size. In a ! curvilinear build that is bgeo%ds, not the logical spacing. dloc = dr ! Add local source terms: wprim = uprim(1:nw_phys, ix1, ix2, ix3) wCT = bga%w(ix1, ix2, ix3, 1:nw_phys, n) wnew = bgb%w(ix1, ix2, ix3, 1:nw_phys, n) call addsource_local(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), qtC, wCT,& wprim, qt, wnew, xloc, dloc, .false. ) bgb%w(ix1, ix2, ix3, 1:nw_flux, n) = wnew(1:nw_flux) end do end do end do ! Reduce fluxes to be stored for the coarse neighbor case ! Direction 1 select case (neighbor_type_1m) case (neighbor_coarse) do ix3=1,nxCo3 do ix2=1,nxCo2 pflux(1,1)%flux(1,ix2,ix3,1:nw_flux,n) = & fC1(1,1+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,2+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,1+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(1,2+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_1p) case (neighbor_coarse) do ix3=1,nxCo3 do ix2=1,nxCo2 pflux(2,1)%flux(1,ix2,ix3,1:nw_flux,n) = & fC1(2,1+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,2+2*(ix2-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,1+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC1(2,2+2*(ix2-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select ! Direction 2 select case (neighbor_type_2m) case (neighbor_coarse) do ix3=1,nxCo3 do ix1=1,nxCo1 pflux(1,2)%flux(ix1,1,ix3,1:nw_flux,n) = & fC2(1,1+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,2+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,1+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(1,2+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_2p) case (neighbor_coarse) do ix3=1,nxCo3 do ix1=1,nxCo1 pflux(2,2)%flux(ix1,1,ix3,1:nw_flux,n) = & fC2(2,1+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,2+2*(ix1-1)+nghostcells,1+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,1+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) & +fC2(2,2+2*(ix1-1)+nghostcells,2+2*(ix3-1)+nghostcells,1:nw_flux) end do end do end select ! Direction 3 select case (neighbor_type_3m) case (neighbor_coarse) do ix2=1,nxCo2 do ix1=1,nxCo1 pflux(1,3)%flux(ix1,ix2,1,1:nw_flux,n) = & fC3(1,1+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,2+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,1+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(1,2+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) end do end do end select select case (neighbor_type_3p) case (neighbor_coarse) do ix2=1,nxCo2 do ix1=1,nxCo1 pflux(2,3)%flux(ix1,ix2,1,1:nw_flux,n) = & fC3(2,1+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,2+2*(ix1-1)+nghostcells,1+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,1+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) & +fC3(2,2+2*(ix1-1)+nghostcells,2+2*(ix2-1)+nghostcells,1:nw_flux) end do end do end select end do end do end subroutine finite_volume_local_muscl_hllc !> MUSCL reconstruction in primitive variables for two faces using a 5-point stencil. !> Returns uL(:,iface), uR(:,iface) for iface=1 (between cells 2-3) and iface=2 (between 3-4). pure subroutine muscl_reconstruct_prim(u, typelim, uL, uR) use mod_limiter, only: limiter_minmod, limiter_vanleer, limiter_mcbeta, limiter_koren real(dp), intent(in) :: u(nw_phys,5) integer, intent(in) :: typelim real(dp), intent(out) :: uL(nw_phys,2), uR(nw_phys,2) real(dp) :: sig(nw_phys,3) ! slopes at cells 2,3,4 integer :: iw select case (typelim) case (limiter_minmod) do iw=1,nw_phys sig(iw,1) = minmod (u(iw,2)-u(iw,1), u(iw,3)-u(iw,2)) ! cell 2 sig(iw,2) = minmod (u(iw,3)-u(iw,2), u(iw,4)-u(iw,3)) ! cell 3 sig(iw,3) = minmod (u(iw,4)-u(iw,3), u(iw,5)-u(iw,4)) ! cell 4 end do case (limiter_vanleer) do iw=1,nw_phys sig(iw,1) = vanleer(u(iw,2)-u(iw,1), u(iw,3)-u(iw,2)) sig(iw,2) = vanleer(u(iw,3)-u(iw,2), u(iw,4)-u(iw,3)) sig(iw,3) = vanleer(u(iw,4)-u(iw,3), u(iw,5)-u(iw,4)) end do case (limiter_mcbeta) do iw=1,nw_phys sig(iw,1) = mcbeta(u(iw,2)-u(iw,1), u(iw,3)-u(iw,2)) sig(iw,2) = mcbeta(u(iw,3)-u(iw,2), u(iw,4)-u(iw,3)) sig(iw,3) = mcbeta(u(iw,4)-u(iw,3), u(iw,5)-u(iw,4)) end do case (limiter_koren) do iw=1,nw_phys sig(iw,1) = koren(u(iw,2)-u(iw,1), u(iw,3)-u(iw,2)) sig(iw,2) = koren(u(iw,3)-u(iw,2), u(iw,4)-u(iw,3)) sig(iw,3) = koren(u(iw,4)-u(iw,3), u(iw,5)-u(iw,4)) end do case default ! Fallback: Godunov (piecewise constant) do iw=1,nw_phys sig(iw,1)=0._dp sig(iw,2)=0._dp sig(iw,3)=0._dp end do end select ! Face 1: between cells 2 and 3 uL(:,1) = u(:,2) + 0.5_dp*sig(:,1) uR(:,1) = u(:,3) - 0.5_dp*sig(:,2) ! Face 2: between cells 3 and 4 uL(:,2) = u(:,3) + 0.5_dp*sig(:,2) uR(:,2) = u(:,4) - 0.5_dp*sig(:,3) ! The slopes above were limited slot by slot, so a state can come out ! inadmissible or -- where the physics carries auxiliary variables -- ! internally inconsistent. Repair it here, once, rather than in each of the ! three Riemann solvers: this is the single producer of face states. call fix_prim_state(uL(:,1)) call fix_prim_state(uR(:,1)) call fix_prim_state(uL(:,2)) call fix_prim_state(uR(:,2)) end subroutine muscl_reconstruct_prim !> One-face LLF/Rusanov numerical flux from primitive L/R states. !> phi(nw_flux) is the adaptive-diffusion reduction factor; only present when FLUX_AD is defined. subroutine riemann_llf_prim(uL, uR, xC, flux_dim, F, phi) real(dp), intent(inout) :: uL(nw_phys), uR(nw_phys) real(dp), intent(in) :: xC(ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: F(nw_flux) real(dp), intent(in) :: phi(nw_flux) real(dp) :: flux_l(nw_flux), flux_r(nw_flux) real(dp) :: wL, wR, wmax call get_flux(uL, xC, flux_dim, flux_l) call get_flux(uR, xC, flux_dim, flux_r) wL = get_cmax(uL, xC, flux_dim) wR = get_cmax(uR, xC, flux_dim) wmax = max(abs(wL), abs(wR)) call to_conservative(uL) call to_conservative(uR) F = 0.5_dp * ((flux_l + flux_r) - wmax * (uR(1:nw_flux) - uL(1:nw_flux)) * phi) end subroutine riemann_llf_prim !> One-face HLL numerical flux from primitive L/R states. !> takes flux_type(flux_dim,nw_flux) into account (fallback to LLF for selected variables) subroutine riemann_hll_prim(uL, uR, xC, flux_dim, F) real(dp), intent(inout) :: uL(nw_phys), uR(nw_phys) real(dp), intent(in) :: xC(ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: F(nw_flux) real(dp) :: flux_l(nw_flux), flux_r(nw_flux) real(dp) :: sL, sR, ds, wmax real(dp), parameter :: eps = 1e-14_dp integer :: iw call get_flux(uL, xC, flux_dim, flux_l) call get_flux(uR, xC, flux_dim, flux_r) call estimate_speeds_minmax(uL, uR, xC, flux_dim, sL, sR) wmax = max(abs(sL), abs(sR)) ! for fallback call to_conservative(uL) call to_conservative(uR) if (sL .ge. 0._dp) then F = flux_l else if (sR .le. 0._dp) then F = flux_r else ds = sR - sL if (ds > eps * (abs(sL) + abs(sR))) then F = (sR*flux_l - sL*flux_r + sL*sR*(uR(1:nw_flux) - uL(1:nw_flux))) / ds else ! Fall back to LLF F = 0.5_dp * ((flux_l + flux_r) - wmax * (uR(1:nw_flux) - uL(1:nw_flux))) end if end if ! overwrite with tvdlf flux if requested (faster this way) if ( any(flux_type(flux_dim,:) == flux_tvdlf) ) then do iw = 1, nw_flux if (flux_type(flux_dim,iw) == flux_tvdlf) then F(iw) = 0.5_dp * ((flux_l(iw) + flux_r(iw)) - wmax * (uR(iw) - uL(iw))) end if end do end if end subroutine riemann_hll_prim !> One-face HLLC numerical flux from primitive L/R states (HD only). !> Uses scalar pressure assumption; not valid for MHD (use HLLD instead). !> Reference: Toro (2010), chapter 10 (Variant 2) !> does not yet take flux_type(flux_dim,nw_flux) into account subroutine riemann_hllc_prim(uL, uR, xC, flux_dim, F) real(dp), intent(inout) :: uL(nw_phys), uR(nw_phys) real(dp), intent(in) :: xC(ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: F(nw_flux) real(dp) :: flux_l(nw_flux), flux_r(nw_flux) real(dp) :: sL, sR, s_star real(dp) :: dL, dR, denom real(dp) :: rhoL, rhoR, unL, unR, pL, pR real(dp) :: PLR real(dp) :: Dstar real(dp) :: wmax real(dp), parameter :: eps = 1e-14_dp integer :: i call get_flux(uL, xC, flux_dim, flux_l) call get_flux(uR, xC, flux_dim, flux_r) call estimate_speeds_toro_pvrs(uL, uR, xC, flux_dim, sL, sR) ! cache primitive scalars before uL/uR are overwritten to conservative rhoL = uL(iw_rho) rhoR = uR(iw_rho) unL = uL(iw_mom(flux_dim)) unR = uR(iw_mom(flux_dim)) pL = uL(iw_e) pR = uR(iw_e) ! we re-use these later (Eq. 10.70, Eq. 10.76) dL = rhoL*(sL - unL) dR = rhoR*(sR - unR) denom = dL - dR ! needed because (Eq. 10.75) uses U_K call to_conservative(uL) call to_conservative(uR) ! fall back to LLF/Rusanov if denom is ~0 if (abs(denom) <= eps*(abs(dL) + abs(dR) + 1._dp)) then wmax = max(abs(sL), abs(sR)) F = 0.5_dp * ((flux_l + flux_r) - wmax * (uR(1:nw_flux) - uL(1:nw_flux))) return end if ! contact speed S* (Eq. 10.70) s_star = ((pR - pL) + unL*dL - unR*dR) / denom ! pressure-like term (Eq. 10.76) PLR = 0.5_dp * (pL + pR + dL*(s_star - unL) + dR*(s_star - unR)) ! HLLC flux (Eq. 10.71) + variant 2 (Eq. 10.75) if (0._dp .le. sL) then F = flux_l return elseif (0._dp .ge. sR) then F = flux_r return end if if (0._dp .le. s_star) then ! F_*L do i = 1, nw_flux Dstar = 0._dp if (i .eq. iw_mom(flux_dim)) Dstar = 1._dp if (i .eq. iw_e) Dstar = s_star F(i) = (s_star*(sL*uL(i) - flux_l(i)) + sL*PLR*Dstar) / (sL - s_star) end do else ! F_*R do i = 1, nw_flux Dstar = 0._dp if (i .eq. iw_mom(flux_dim)) Dstar = 1._dp if (i .eq. iw_e) Dstar = s_star F(i) = (s_star*(sR*uR(i) - flux_r(i)) + sR*PLR*Dstar) / (sR - s_star) end do end if end subroutine riemann_hllc_prim !> MUSCL (primitive-variable) reconstruction with slope limiter; HLL two-wave approximate Riemann flux at faces. !> Uses estimated left/right signal speeds (Davis (1988)) for less diffusion than LLF, no contact resolution. subroutine reconflux_muscl_hll_prim(u, xlocC, flux_dim, flux, typelim) real(dp), intent(in) :: u(nw_phys, 5) real(dp), intent(in) :: xlocC(1:ndim, 2) integer, intent(in) :: flux_dim, typelim real(dp), intent(out) :: flux(nw_flux, 2) real(dp) :: uL(nw_phys,2), uR(nw_phys,2) integer :: iface call muscl_reconstruct_prim(u, typelim, uL, uR) do iface=1,2 call riemann_hll_prim(uL(:,iface), uR(:,iface), xlocC(:,iface), flux_dim, flux(:,iface)) end do end subroutine reconflux_muscl_hll_prim !> MUSCL (primitive-variable) reconstruction with slope limiter; LLF/Rusanov numerical flux at faces. !> Robust and diffusive; uses local max wavespeed for upwinding. !> Adaptive diffusion (Rempel et al. 2009) is compiled in only when FLUX_AD is defined. subroutine reconflux_muscl_llf_prim(u, xlocC, flux_dim, flux, typelim) real(dp), intent(in) :: u(nw_phys, 5) real(dp), intent(in) :: xlocC(1:ndim, 2) integer, intent(in) :: flux_dim, typelim real(dp), intent(out) :: flux(nw_flux, 2) real(dp) :: uL(nw_phys,2), uR(nw_phys,2) integer :: iface real(dp) :: phi(nw_flux,2) real(dp) :: delta_rc, delta_ct, phi_v integer :: iw call muscl_reconstruct_prim(u, typelim, uL, uR) ! Rempel et al. (2009): reduce dissipation where reconstruction and cell-centre agree in sign phi = 1.0_dp do iface = 1, 2 do iw = 1, nw_flux delta_rc = uR(iw,iface) - uL(iw,iface) delta_ct = u(iw, iface+2) - u(iw, iface+1) if (delta_rc * delta_ct > 1.0e-18_dp) then phi_v = flux_ad_scale * delta_rc**2 / (delta_ct**2 + 1.0e-18_dp) phi(iw,iface) = max(flux_ad_min, min(phi_v, 1.0_dp)) end if end do end do do iface = 1, 2 call riemann_llf_prim(uL(:,iface), uR(:,iface), xlocC(:,iface), flux_dim, flux(:,iface), phi(:,iface)) end do end subroutine reconflux_muscl_llf_prim !> MUSCL (primitive-variable) reconstruction with slope limiter; HLLC approximate Riemann flux at faces. !> Restores the contact wave (and shear in Euler/HD), typically sharper than HLL for similar cost. subroutine reconflux_muscl_hllc_prim(u, xlocC, flux_dim, flux, typelim) real(dp), intent(in) :: u(nw_phys, 5) real(dp), intent(in) :: xlocC(1:ndim, 2) integer, intent(in) :: flux_dim, typelim real(dp), intent(out) :: flux(nw_flux, 2) real(dp) :: uL(nw_phys,2), uR(nw_phys,2) integer :: iface call muscl_reconstruct_prim(u, typelim, uL, uR) do iface=1,2 call riemann_hllc_prim(uL(:,iface), uR(:,iface), xlocC(:,iface), flux_dim, flux(:,iface)) end do end subroutine reconflux_muscl_hllc_prim pure real(dp) function vanleer(a, b) result(phi) real(dp), intent(in) :: a, b real(dp) :: ab ab = a * b if (ab > 0) then phi = 2.0_dp * ab / (a + b) else phi = 0.0_dp end if end function vanleer pure real(dp) function minmod(a, b) real(dp), intent(in) :: a, b if (a * b <= 0) then minmod = 0.0_dp else if (abs(a) < abs(b)) then minmod = a else minmod = b end if end function minmod ! scavenged from afivo, thanks Jannis :-) ! !> Modified implementation of Koren limiter, to avoid division and the min/max !> functions, which can be problematic / expensive. In most literature, you !> have r = a / b (ratio of gradients). Then the limiter phi(r) is multiplied !> with b. With this implementation, you get phi(r) * b pure real(dp) function koren(a, b) real(dp), intent(in) :: a !< Density gradient (numerator) real(dp), intent(in) :: b !< Density gradient (denominator) real(dp), parameter :: third = 1/3.0_dp real(dp) :: aa, ab aa = a * a ab = a * b if (ab <= 0) then ! a and b have different sign or one of them is zero, so r is either 0, ! inf or negative (special case a == b == 0 is ignored) koren = 0 else if (aa <= 0.25_dp * ab) then ! 0 < a/b <= 1/4, limiter has value 2*a/b koren = 2*a else if (aa <= 2.5_dp * ab) then ! 1/4 < a/b <= 2.5, limiter has value (1+2*a/b)/3 koren = third * (b + 2*a) else ! (1+2*a/b)/6 >= 1, limiter has value 2 koren = 2*b end if end function koren !> Monotonised central-difference limiter with tunable beta (AMRVAC's !> mcbeta; beta=2 recovers the classic MC limiter, van Leer 1979) pure real(dp) function mcbeta(a, b) result(phi) real(dp), intent(in) :: a, b real(dp), parameter :: c_mcbeta = 1.4_dp real(dp) :: ab ab = a * b if (ab > 0) then phi = sign(1.0_dp, a + b) * min(c_mcbeta*abs(a), c_mcbeta*abs(b), abs(a+b)*0.5_dp) else phi = 0.0_dp end if end function mcbeta end module mod_finite_volume_srhd