/*@@ @file GRHydro_HLLE.F90 @date Sat Jan 26 01:40:14 2002 @author Ian Hawke, Pedro Montero, Toni Font @desc The HLLE solver. Called from the wrapper function, so works in all directions. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" #include "GRHydro_Macros.h" #include "SpaceMask.h" /*@@ @routine GRHydro_HLLEGeneral @date Sat Jan 26 01:41:02 2002 @author Ian Hawke, Pedro Montero, Toni Font @desc The HLLE solver. Sufficiently simple that its just one big routine. Rewritten for the new EOS interface. @enddesc @calls @calledby @history Altered from Cactus 3 routines originally written by Toni Font. @endhistory @@*/ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) USE GRHydro_Eigenproblem implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer :: i, j, k, m CCTK_REAL, dimension(5) :: cons_p,cons_m,fplus,fminus,lamplus CCTK_REAL, dimension(5) :: f1,lamminus CCTK_REAL, dimension(5) :: qdiff CCTK_REAL, dimension(6) :: prim_p, prim_m CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & cs2_p, cs2_m, dpdeps_p, dpdeps_m CCTK_INT :: type_bits, trivial, not_trivial integer tadmor if(CCTK_EQUALS(HLLE_type,"Tadmor")) then tadmor = 1 else tadmor = 0 endif if (flux_direction == 1) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & &"trivial") else if (flux_direction == 2) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", & &"trivial") else if (flux_direction == 3) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else call CCTK_WARN(0, "Flux direction not x,y,z") end if do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil f1 = 0.d0 lamminus = 0.d0 lamplus = 0.d0 cons_p = 0.d0 cons_m = 0.d0 fplus = 0.d0 fminus = 0.d0 qdiff = 0.d0 !!$ Set the left (p for plus) and right (m_i for minus, i+1) states cons_p(1) = densplus(i,j,k) cons_p(2) = sxplus(i,j,k) cons_p(3) = syplus(i,j,k) cons_p(4) = szplus(i,j,k) cons_p(5) = tauplus(i,j,k) cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset) cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset) cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset) cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) prim_p(1) = rhoplus(i,j,k) prim_p(2) = velxplus(i,j,k) prim_p(3) = velyplus(i,j,k) prim_p(4) = velzplus(i,j,k) prim_p(5) = epsplus(i,j,k) prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset) prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset) prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) prim_p(6) = pressplus(i,j,k) cs2_p = eos_cs2_p(i,j,k) dpdeps_p = eos_dpdeps_p(i,j,k) prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset) dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset) !!$ Calculate various metric terms here. !!$ Note also need the average of the lapse at the !!$ left and right points. if (shift_state .ne. 0) then if (flux_direction == 1) then avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + & betax(i,j,k)) else if (flux_direction == 2) then avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + & betay(i,j,k)) else if (flux_direction == 3) then avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + & betaz(i,j,k)) else call CCTK_WARN(0, "Flux direction not x,y,z") end if else avg_beta = 0.d0 end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k)) gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k)) gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k)) gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k)) gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) !!$ If the Riemann problem is trivial, just calculate the fluxes from the !!$ left state and skip to the next cell if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then if (flux_direction == 1) then call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& f1(1),f1(2),f1(3),f1(4),f1(5),& prim_m(2),prim_m(6), & avg_det,avg_alp,avg_beta) else if (flux_direction == 2) then call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& f1(1),f1(3),f1(4),f1(2),f1(5),& prim_m(3),prim_m(6), & avg_det,avg_alp,avg_beta) else if (flux_direction == 3) then call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& f1(1),f1(4),f1(2),f1(3),f1(5), & prim_m(4),prim_m(6), & avg_det,avg_alp,avg_beta) else call CCTK_WARN(0, "Flux direction not x,y,z") end if else !!! The end of this branch is right at the bottom of the routine call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & avg_det,gxxh,gxyh,gxzh, & gyyh,gyzh, gzzh) if (flux_direction == 1) then usendh = uxxh else if (flux_direction == 2) then usendh = uyyh else if (flux_direction == 3) then usendh = uzzh else call CCTK_WARN(0, "Flux direction not x,y,z") end if !!$ Calculate the jumps in the conserved variables qdiff(1) = cons_m(1) - cons_p(1) qdiff(2) = cons_m(2) - cons_p(2) qdiff(3) = cons_m(3) - cons_p(3) qdiff(4) = cons_m(4) - cons_p(4) qdiff(5) = cons_m(5) - cons_p(5) !!$ Eigenvalues and fluxes either side of the cell interface if (flux_direction == 1) then call eigenvalues_general(& prim_m(2),prim_m(3),prim_m(4),cs2_m, & lamminus,& gxxh,gxyh,gxzh,& gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(2),prim_p(3),prim_p(4),cs2_p, & lamplus,& gxxh,gxyh,gxzh,& gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),& prim_p(2),prim_p(6), & avg_det,avg_alp,avg_beta) call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),& prim_m(2),prim_m(6), & avg_det,avg_alp,avg_beta) else if (flux_direction == 2) then call eigenvalues_general(& prim_m(3),prim_m(4),prim_m(2),cs2_m, & lamminus,& gyyh,gyzh,gxyh,& gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(3),prim_p(4),prim_p(2),cs2_p, & lamplus,& gyyh,gyzh,gxyh,& gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),& prim_p(3),prim_p(6), & avg_det,avg_alp,avg_beta) call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),& prim_m(3),prim_m(6), & avg_det,avg_alp,avg_beta) else if (flux_direction == 3) then call eigenvalues_general(& prim_m(4),prim_m(2),prim_m(3),cs2_m, & lamminus,& gzzh,gxzh,gyzh,& gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(4),prim_p(2),prim_p(3),cs2_p, & lamplus,& gzzh,gxzh,gyzh,& gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& fplus(1),fplus(4),fplus(2),fplus(3),fplus(5), & prim_p(4),prim_p(6), & avg_det,avg_alp,avg_beta) call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), & prim_m(4),prim_m(6), & avg_det,avg_alp,avg_beta) else call CCTK_WARN(0, "Flux direction not x,y,z") end if if(tadmor.eq.0) then !!$ Find minimum and maximum wavespeeds charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& lamminus(4),lamminus(5)) charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& lamminus(4),lamminus(5)) charpm = charmax - charmin !!$ Calculate flux by standard formula do m = 1,5 qdiff(m) = cons_m(m) - cons_p(m) f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & charmax * charmin * qdiff(m)) / charpm end do else ! Tadmor's semi-discrete scheme: JcP 160, 241 (2000) charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), & abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), & abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5))) do m = 1,5 qdiff(m) = cons_m(m) - cons_p(m) f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* & qdiff(m) end do end if end if !!! The end of the SpaceMask check for a trivial RP. densflux(i, j, k) = f1(1) sxflux(i, j, k) = f1(2) syflux(i, j, k) = f1(3) szflux(i, j, k) = f1(4) tauflux(i, j, k) = f1(5) if(evolve_Y_e.ne.0) then if (densflux(i, j, k) > 0.d0) then Y_e_con_flux(i, j, k) = & Y_e_plus(i, j, k) * & densflux(i, j, k) else Y_e_con_flux(i, j, k) = & Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * & densflux(i, j, k) endif endif end do end do end do end subroutine GRHydro_HLLEGeneral subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) USE GRHydro_Eigenproblem implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer :: i, j, k, m CCTK_REAL, dimension(number_of_tracers) :: cons_p,cons_m,fplus,fminus,f1 CCTK_REAL, dimension(5) :: lamminus,lamplus CCTK_REAL, dimension(number_of_tracers) :: qdiff CCTK_REAL, dimension(6) :: prim_p, prim_m CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & cs2_p, cs2_m, dpdeps_p, dpdeps_m integer tadmor if(CCTK_EQUALS(HLLE_type,"Tadmor")) then tadmor = 1 else tadmor = 0 endif do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil f1 = 0.d0 lamminus = 0.d0 lamplus = 0.d0 cons_p = 0.d0 cons_m = 0.d0 fplus = 0.d0 fminus = 0.d0 qdiff = 0.d0 !!$ Set the left (p for plus) and right (m_i for minus, i+1) states cons_p(:) = cons_tracerplus(i,j,k,:) cons_m(:) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) prim_p(1) = rhoplus(i,j,k) prim_p(2) = velxplus(i,j,k) prim_p(3) = velyplus(i,j,k) prim_p(4) = velzplus(i,j,k) prim_p(5) = epsplus(i,j,k) prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset) prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset) prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) prim_p(6) = pressplus(i,j,k) cs2_p = eos_cs2_p(i,j,k) dpdeps_p = eos_dpdeps_p(i,j,k) prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset) dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset) !!$ Calculate various metric terms here. !!$ Note also need the average of the lapse at the !!$ left and right points. if (shift_state .ne. 0) then if (flux_direction == 1) then avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + & betax(i,j,k)) else if (flux_direction == 2) then avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + & betay(i,j,k)) else if (flux_direction == 3) then avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + & betaz(i,j,k)) else call CCTK_WARN(0, "Flux direction not x,y,z") end if else avg_beta = 0.d0 end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k)) gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k)) gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k)) gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k)) gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) !!$ If the Riemann problem is trivial, just calculate the fluxes from the !!$ left state and skip to the next cell call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & avg_det,gxxh, gxyh, gxzh, & gyyh, gyzh, gzzh) if (flux_direction == 1) then usendh = uxxh else if (flux_direction == 2) then usendh = uyyh else if (flux_direction == 3) then usendh = uzzh else call CCTK_WARN(0, "Flux direction not x,y,z") end if !!$ Eigenvalues and fluxes either side of the cell interface if (flux_direction == 1) then call eigenvalues_general(& prim_m(2),prim_m(3),prim_m(4),cs2_m, & lamminus,& gxxh,gxyh,gxzh,& gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(2),prim_p(3),prim_p(4),cs2_p, & lamplus,& gxxh,gxyh,gxzh,& gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * & cons_tracerplus(i,j,k,:) fminus(:) = (velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else if (flux_direction == 2) then call eigenvalues_general(& prim_m(3),prim_m(4),prim_m(2),cs2_m, & lamminus,& gyyh,gyzh,gxyh,& gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(3),prim_p(4),prim_p(2),cs2_p, & lamplus,& gyyh,gyzh,gxyh,& gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * & cons_tracerplus(i,j,k,:) fminus(:) = (velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else if (flux_direction == 3) then call eigenvalues_general(& prim_m(4),prim_m(2),prim_m(3),cs2_m, & lamminus,& gzzh,gxzh,gyzh,& gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(4),prim_p(2),prim_p(3),cs2_p, & lamplus,& gzzh,gxzh,gyzh,& gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * & cons_tracerplus(i,j,k,:) fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else call CCTK_WARN(0, "Flux direction not x,y,z") end if if(tadmor.eq.0) then !!$ Find minimum and maximum wavespeeds charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& lamminus(4),lamminus(5)) charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& lamminus(4),lamminus(5)) charpm = charmax - charmin !!$ Calculate flux by standard formula do m = 1,number_of_tracers qdiff(m) = cons_m(m) - cons_p(m) f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & charmax * charmin * qdiff(m)) / charpm end do else ! Tadmor's semi-descrite scheme: JcP 160, 241 (2000) charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), & abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), & abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5))) do m = 1,number_of_tracers qdiff(m) = cons_m(m) - cons_p(m) f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* & qdiff(m) end do end if cons_tracerflux(i,j,k,:) = f1(:) end do end do end do end subroutine GRHydro_HLLE_TracerGeneral