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