!> Module with finite volume methods for fluxes module mod_finite_volume_mhd 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 use mod_usr, only: gravity_field use mod_radiative_cooling, only: rc_fl, radiative_cooling_add_source use mod_global_parameters, only:cmax_global 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) call addsource_usr(qdt, qt, wCT, wCTprim, wnew, x, .false.) else !--------------------------------- ! split sources !--------------------------------- wnew(psi_)=wnew(psi_)*dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dr)) 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) use mod_global_parameters, only: dt, ndir, smalldouble, courantpar 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 ! .. local .. real(dp) :: laplb_cd2 real(dp) :: Jdir1,Jdir2,Jdir3 integer :: idir real(dp) :: Te_c, rho_c, pth_c real(dp) :: gradT(3), Bmag2, Bmag, Bmag2_safe, bgradT real(dp) :: bhat(3), gradT_perp(3), gradTperp_mag real(dp) :: chi, sig_par, sig_perp real(dp) :: q_sat, f_sat_par, f_sat_perp, tau_par, tau_perp real(dp) :: cf2_tc, cmax2_tc real(dp) :: Qpar_proj, Qperp_proj_k integer :: k_tc if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- ! using the compact stencil formulation and adopting constant eta ! > eta*(laplacian of B)_idir added to B_idir ! > B_idir*[eta*(laplacian of B)_idir] added to e_ do idir = 1,ndim laplb_cd2 = & ( & wCTprim1(iw_mag(1)-1+idir,3) & - 2*wCTprim1(iw_mag(1)-1+idir,2) & + wCTprim1(iw_mag(1)-1+idir,1) & ) & / dx(1)**2 + & ( & wCTprim2(iw_mag(1)-1+idir,3) & - 2*wCTprim2(iw_mag(1)-1+idir,2) & + wCTprim2(iw_mag(1)-1+idir,1) & ) & / dx(2)**2 + & ( & wCTprim3(iw_mag(1)-1+idir,3) & - 2*wCTprim3(iw_mag(1)-1+idir,2) & + wCTprim3(iw_mag(1)-1+idir,1) & ) & / dx(3)**2 wnew(iw_mag(1)-1+idir) = wnew(iw_mag(1)-1+idir) + qdt*mhd_eta*laplb_cd2 wnew(iw_e) = wnew(iw_e) + qdt * mhd_eta * laplb_cd2 * wCTprim1(iw_mag(1)-1+idir,2) enddo ! > eta* J**2 added to e Jdir1 = (wCTprim2(iw_mag(3),3) - wCTprim2(iw_mag(3),1)) & / 2.0_dp/dx(2) & - (wCTprim3(iw_mag(2),3) - wCTprim3(iw_mag(2),1)) & / 2.0_dp/dx(3) Jdir2 = (wCTprim3(iw_mag(1),3) - wCTprim3(iw_mag(1),1)) & /2.0_dp/dx(3) & - (wCTprim1(iw_mag(3),3) - wCTprim1(iw_mag(3),1)) & / 2.0_dp/dx(1) Jdir3 = (wCTprim1(iw_mag(2),3) - wCTprim1(iw_mag(2),1)) & /2.0_dp/dx(1) & - (wCTprim2(iw_mag(1),3) - wCTprim2(iw_mag(1),1)) & / 2.0_dp/dx(2) wnew(iw_e) = wnew(iw_e) + qdt*mhd_eta*(Jdir1**2+Jdir2**2+Jdir3**2) ! ---- 1. Local thermodynamic state ---- Te_c = wCTprim1(iw_e,2) / wCTprim1(iw_rho,2) rho_c = wCTprim1(iw_rho,2) pth_c = wCTprim1(iw_e,2) ! ---- 2. Temperature gradient ---- gradT(1) = (wCTprim1(iw_e,3)/wCTprim1(iw_rho,3) - wCTprim1(iw_e,1)/wCTprim1(iw_rho,1)) / (2.d0*dx(1)) gradT(2) = (wCTprim2(iw_e,3)/wCTprim2(iw_rho,3) - wCTprim2(iw_e,1)/wCTprim2(iw_rho,1)) / (2.d0*dx(2)) gradT(3) = (wCTprim3(iw_e,3)/wCTprim3(iw_rho,3) - wCTprim3(iw_e,1)/wCTprim3(iw_rho,1)) / (2.d0*dx(3)) !--------------------------------- ! Anisotropic TC code here !--------------------------------- ! ---- 3. Field direction ---- Bmag2 = 0.0d0 do k_tc = 1, ndir Bmag2 = Bmag2 + wCTprim1(iw_mag(k_tc),2)**2 end do Bmag = sqrt(Bmag2) Bmag2_safe = max(Bmag2, smalldouble**2) do k_tc = 1, ndir bhat(k_tc) = wCTprim1(iw_mag(k_tc),2) * Bmag / Bmag2_safe end do ! ---- 4. Decompose gradT into parallel scalar + perpendicular vector ---- bgradT = bhat(1)*gradT(1) + bhat(2)*gradT(2) + bhat(3)*gradT(3) do k_tc = 1, ndir gradT_perp(k_tc) = gradT(k_tc) - bgradT * bhat(k_tc) end do gradTperp_mag = sqrt(gradT_perp(1)**2 + gradT_perp(2)**2 + gradT_perp(3)**2) ! ---- 5. Conductivities ---- if (tc_kappa_par > 0.0d0) then sig_par = tc_kappa_par else sig_par = tc_kappa0_par * Te_c**2.5d0 end if if (tc_kappa_perp > 0.0d0) then sig_perp = tc_kappa_perp else chi = htc_Cchi * Bmag * Te_c**1.5d0 / rho_c sig_perp = sig_par / (1.0d0 + chi**2) end if ! ---- 6. Saturation and relaxation times (per channel) ---- cf2_tc = (Bmag2 + mhd_gamma*pth_c) / rho_c cmax2_tc = 0.0d0 do k_tc = 1, ndim cmax2_tc = max(cmax2_tc, 0.5d0*(cf2_tc + sqrt(abs( & cf2_tc**2 - 4.0d0*mhd_gamma*pth_c*wCTprim1(iw_mag(k_tc),2)**2/rho_c**2)))) end do q_sat = 1.5d0 * rho_c * (pth_c / rho_c)**1.5d0 f_sat_par = 1.0d0 / (1.0d0 + abs(sig_par * bgradT ) / q_sat) f_sat_perp = 1.0d0 / (1.0d0 + abs(sig_perp * gradTperp_mag) / q_sat) tau_par = max(4.d0*dt, f_sat_par *sig_par *Te_c*courantpar**2*(mhd_gamma-1.0d0) / (pth_c*cmax2_tc)) tau_perp = max(4.d0*dt, f_sat_perp*sig_perp*Te_c*courantpar**2*(mhd_gamma-1.0d0) / (pth_c*cmax2_tc)) ! ---- 7. Decompose current q onto current field direction ---- Qpar_proj = wCTprim1(q_(1),2)*bhat(1) + wCTprim1(q_(2),2)*bhat(2) + wCTprim1(q_(3),2)*bhat(3) ! ---- 8. Relax parallel and perpendicular parts, each at its own rate ---- do k_tc = 1, ndir Qperp_proj_k = wCTprim1(q_(k_tc),2) - Qpar_proj * bhat(k_tc) wnew(q_(k_tc)) = wnew(q_(k_tc)) & - qdt*(f_sat_par *sig_par *bgradT + Qpar_proj)*bhat(k_tc)/tau_par & - qdt*(f_sat_perp*sig_perp*gradT_perp(k_tc) + Qperp_proj_k)/tau_perp end do !--------------------------------- ! Field-aligned TC code here !--------------------------------- else !--------------------------------- ! split sources !--------------------------------- ! Not yet implemented end if end subroutine addsource_compact pure subroutine to_primitive(u) real(dp), intent(inout) :: u(nw_phys) 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) u(iw_e)=gamma_1*(u(iw_e)-0.5_dp*& (u(iw_rho)*(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)+& u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**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)*(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)+& u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**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 !> Returns maximum local signal speed |v_n| + c_f (fast magnetosonic) 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) :: inv_rho, b2, cfast2 inv_rho=1.0_dp/u(iw_rho) wC=mhd_gamma*u(iw_e)*inv_rho cfast2=(u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**2)*inv_rho+wC b2=cfast2**2-4.0_dp*wC*u(iw_mag(flux_dim))**2*inv_rho if(b2<0.d0) b2=0.d0 wC=sqrt(0.5_dp*(cfast2+sqrt(b2)))+abs(u(iw_mom(flux_dim))) end function get_cmax subroutine get_flux(u, xC, flux_dim, flux) use mod_global_parameters, only: cmax_global use mod_global_parameters, only: smalldouble 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) :: ptotal ! Hyperbolic TC field geometry ! 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))-& u(iw_mag(flux_dim))*u(iw_mag(1)) flux(iw_mom(2))=u(iw_rho)*u(iw_mom(2))*u(iw_mom(flux_dim))-& u(iw_mag(flux_dim))*u(iw_mag(2)) flux(iw_mom(3))=u(iw_rho)*u(iw_mom(3))*u(iw_mom(flux_dim))-& u(iw_mag(flux_dim))*u(iw_mag(3)) ptotal=u(iw_e)+0.5_dp*(u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**2) flux(iw_mom(flux_dim))=flux(iw_mom(flux_dim))+ptotal ! Energy flux flux(iw_e)=u(iw_mom(flux_dim))*(u(iw_e)*inv_gamma_1+0.5_dp*& u(iw_rho)*(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)+& 2.0_dp*ptotal-u(iw_e))-u(iw_mag(flux_dim))*& (u(iw_mag(1))*u(iw_mom(1))+u(iw_mag(2))*u(iw_mom(2))+u(iw_mag(3))*u(iw_mom(3))) ! Magnetic flux flux(iw_mag(1))=u(iw_mom(flux_dim))*u(iw_mag(1))-u(iw_mag(flux_dim))*u(iw_mom(1)) flux(iw_mag(2))=u(iw_mom(flux_dim))*u(iw_mag(2))-u(iw_mag(flux_dim))*u(iw_mom(2)) flux(iw_mag(3))=u(iw_mom(flux_dim))*u(iw_mag(3))-u(iw_mag(flux_dim))*u(iw_mom(3)) ! GLM psi flux flux(iw_mag(flux_dim))=u(psi_) !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645 flux(psi_)=cmax_global**2*u(iw_mag(flux_dim)) ! Tracer flux. Note that tracers stay conservative. flux(tracer(1)) = u(tracer(1)) * u(iw_mom(flux_dim)) ! Hyperbolic TC fluxes flux(iw_e) = flux(iw_e) + u(iw_qvec(flux_dim)) flux(iw_qvec(1)) = 0.0_dp flux(iw_qvec(2)) = 0.0_dp flux(iw_qvec(3)) = 0.0_dp end subroutine get_flux !> Davis (1988) min/max wave speed estimates wL = min(v_n - c_f), wR = max(v_n + c_f) (fast magnetosonic) over left/right states; !> used in HLL flux estimation. 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) :: inv_rho, cs2, cA2, cAn2, sum2, discriminant, cfL, cfR ! Left state inv_rho = 1.0_dp / uL(iw_rho) cs2 = mhd_gamma * uL(iw_e) * inv_rho cA2 = (uL(iw_mag(1))**2 + uL(iw_mag(2))**2 + uL(iw_mag(3))**2) * inv_rho cAn2 = uL(iw_mag(flux_dim))**2 * inv_rho sum2 = cs2 + cA2 discriminant = max(sum2**2 - 4.0_dp*cs2*cAn2, 0.0_dp) cfL = sqrt(0.5_dp * (sum2 + sqrt(discriminant))) ! Right state inv_rho = 1.0_dp / uR(iw_rho) cs2 = mhd_gamma * uR(iw_e) * inv_rho cA2 = (uR(iw_mag(1))**2 + uR(iw_mag(2))**2 + uR(iw_mag(3))**2) * inv_rho cAn2 = uR(iw_mag(flux_dim))**2 * inv_rho sum2 = cs2 + cA2 discriminant = max(sum2**2 - 4.0_dp*cs2*cAn2, 0.0_dp) cfR = sqrt(0.5_dp * (sum2 + sqrt(discriminant))) wL = min(uL(iw_mom(flux_dim)) - cfL, uR(iw_mom(flux_dim)) - cfR) wR = max(uL(iw_mom(flux_dim)) + cfL, uR(iw_mom(flux_dim)) + cfR) 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) ! Add non-local compact source terms: wnew = bgb%w(ix1, ix2, ix3, 1:nw_phys, n) tmp1 = uprim(1:nw_phys, ix1-1:ix1+1, ix2, ix3) tmp2 = uprim(1:nw_phys, ix1, ix2-1:ix2+1, ix3) tmp3 = uprim(1:nw_phys, ix1, ix2, ix3-1:ix3+1) call addsource_compact(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), qtC, tmp1,tmp2,tmp3, & 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) ! Add non-local compact source terms: wnew = bgb%w(ix1, ix2, ix3, 1:nw_phys, n) tmp1 = uprim(1:nw_phys, ix1-1:ix1+1, ix2, ix3) tmp2 = uprim(1:nw_phys, ix1, ix2-1:ix2+1, ix3) tmp3 = uprim(1:nw_phys, ix1, ix2, ix3-1:ix3+1) call addsource_compact(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), qtC, tmp1,tmp2,tmp3, & 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) ! Add non-local compact source terms: wnew = bgb%w(ix1, ix2, ix3, 1:nw_phys, n) tmp1 = uprim(1:nw_phys, ix1-1:ix1+1, ix2, ix3) tmp2 = uprim(1:nw_phys, ix1, ix2-1:ix2+1, ix3) tmp3 = uprim(1:nw_phys, ix1, ix2, ix3-1:ix3+1) call addsource_compact(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),& dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), qtC, tmp1,tmp2,tmp3, & 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_mhd