!> Module with finite volume methods for fluxes module mod_finite_volume_hd 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: gravity_field use mod_radiative_cooling, only: rc_fl, radiative_cooling_add_source 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 ! .. local .. integer :: idim real(dp) :: field if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- do idim = 1, ndim field = gravity_field(wCT, x, idim) wnew(iw_mom(idim)) = wnew(iw_mom(idim)) + qdt * field * wCT(iw_rho) wnew(iw_e) = wnew(iw_e) + qdt * field * wCT(iw_mom(idim)) end do call radiative_cooling_add_source(qdt,wCT,wCTprim,wnew,x) else !--------------------------------- ! split sources !--------------------------------- ! Not yet implemented 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) real(dp), intent(inout) :: u(nw_phys) ! Compute velocity from momentum u(iw_mom(1)) = u(iw_mom(1))/u(iw_rho) u(iw_mom(2)) = u(iw_mom(2))/u(iw_rho) u(iw_mom(3)) = u(iw_mom(3))/u(iw_rho) ! Compute pressure from energy u(iw_e) = gamma_1 * (u(iw_e) - 0.5_dp * u(iw_rho) * & sum(u(iw_mom(1:ndim))**2) ) end subroutine to_primitive pure subroutine to_conservative(u) real(dp), intent(inout) :: u(nw_phys) ! Compute energy from pressure and kinetic energy u(iw_e) = u(iw_e) * inv_gamma_1 + 0.5_dp * u(iw_rho) * & sum(u(iw_mom(1:ndim))**2) ! Compute momentum from density and velocity components u(iw_mom(1)) = u(iw_rho) * u(iw_mom(1)) u(iw_mom(2)) = u(iw_rho) * u(iw_mom(2)) u(iw_mom(3)) = u(iw_rho) * u(iw_mom(3)) end subroutine to_conservative !> Make a reconstructed primitive face state admissible before it is handed to !> the Riemann solver. Genuinely empty by default: a physics module without !> auxiliary variables and without a floor of its own needs nothing here, and !> the call is inlined away. pure subroutine fix_prim_state(u) real(dp), intent(inout) :: u(nw_phys) end subroutine fix_prim_state 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 wC = sqrt(hd_gamma * u(iw_e) / u(iw_rho)) + abs(u(iw_mom(flux_dim))) end function get_cmax subroutine get_flux(u, xC, flux_dim, flux) real(dp), intent(in) :: u(nw_phys) real(dp), intent(in) :: xC(ndim) integer, intent(in) :: flux_dim real(dp), intent(out) :: flux(nw_flux) ! Density flux flux(iw_rho) = u(iw_rho) * u(iw_mom(flux_dim)) ! Momentum flux with pressure term flux(iw_mom(1)) = u(iw_rho) * u(iw_mom(1)) * u(iw_mom(flux_dim)) flux(iw_mom(2)) = u(iw_rho) * u(iw_mom(2)) * u(iw_mom(flux_dim)) flux(iw_mom(3)) = u(iw_rho) * u(iw_mom(3)) * u(iw_mom(flux_dim)) flux(iw_mom(flux_dim)) = flux(iw_mom(flux_dim)) + u(iw_e) ! Energy flux flux(iw_e) = u(iw_mom(flux_dim)) * (u(iw_e) * inv_gamma_1 + 0.5_dp * & u(iw_rho) * sum(u(iw_mom(1:ndim))**2) + u(iw_e)) ! Tracer flux. Note that tracers stay conservative. flux(tracer(1)) = u(tracer(1)) * u(iw_mom(flux_dim)) end subroutine get_flux !> Wave speed estimates: min/max acoustic bounds (Davis 1988) !> Reference: Toro 2010, Chapter 10. 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) :: cL, cR cL = sqrt(hd_gamma * uL(iw_e) / uL(iw_rho)) cR = sqrt(hd_gamma * uR(iw_e) / uR(iw_rho)) wL = min(uL(iw_mom(flux_dim)) - cL, uR(iw_mom(flux_dim)) - cR) wR = max(uL(iw_mom(flux_dim)) + cL, uR(iw_mom(flux_dim)) + cR) end subroutine estimate_speeds_minmax !> Wave speed estimates for HLL/HLLC using Toro (2010) PVRS pressure estimate !> Implements Eq. (10.67)-(10.69) subroutine estimate_speeds_toro_pvrs(uL, uR, xC, flux_dim, sL, sR) 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 real(dp) :: rhoL, rhoR, pL, pR, unL, unR real(dp) :: aL, aR, abar, rhobar real(dp) :: ppvrs, pstar real(dp) :: qL, qR, pratio real(dp), parameter :: tiny = 1e-30_dp rhoL = uL(iw_rho); rhoR = uR(iw_rho) pL = uL(iw_e); pR = uR(iw_e) unL = uL(iw_mom(flux_dim)) unR = uR(iw_mom(flux_dim)) ! sound speeds aL = sqrt(hd_gamma * pL / max(rhoL, tiny)) aR = sqrt(hd_gamma * pR / max(rhoR, tiny)) ! PVRS pressure estimate (Eq. 10.67) rhobar = 0.5_dp*(rhoL + rhoR) abar = 0.5_dp*(aL + aR) ppvrs = 0.5_dp*(pL + pR) - 0.5_dp*(unR - unL)*rhobar*abar pstar = max(0._dp, ppvrs) ! (Eq. 10.69) if (pstar <= pL) then qL = 1._dp else pratio = pstar / max(pL, tiny) qL = sqrt(1._dp + 0.5_dp*(hd_gamma + 1._dp)/hd_gamma * (pratio - 1._dp)) end if if (pstar <= pR) then qR = 1._dp else pratio = pstar / max(pR, tiny) qR = sqrt(1._dp + 0.5_dp*(hd_gamma + 1._dp)/hd_gamma * (pratio - 1._dp)) end if sL = unL - aL*qL sR = unR + aR*qR 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_hd