/*@@ @file GRHydro_RiemannSolveM.F90 @date Sep 1, 2010 @author @desc A wrapper routine to call the correct Riemann solver @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" #include "GRHydro_Macros.h" /*@@ @routine RiemannSolveM @date Sep 1, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc A wrapper routine to switch between the different Riemann solvers. @enddesc @calls @calledby @history @endhistory @@*/ subroutine RiemannSolveM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLEM(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then !!$ There are no special calls for tracers, which care not one whit about B-fields! !!$ Just call the standard version... call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) end if !!$ else if (CCTK_EQUALS(riemann_solver,"Roe")) then !!$ !!$ call GRHydro_RoeSolveM(CCTK_PASS_FTOF) !!$ !!$ if (evolve_tracer .ne. 0) then !!$ !!$ call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) !!$ !!$ end if !!$ !!$ else if (CCTK_EQUALS(riemann_solver,"Marquina")) then !!$ !!$ call GRHydro_MarquinaM(CCTK_PASS_FTOF) !!$ Tracers are built directly in to the Marquina solver else call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!") end if end subroutine RiemannSolveM /*@@ @routine RiemannSolvePolytypeM @date Sep 1, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc The same as above, just specializing to polytropic type EOS. Currently there is no point to this routine right now. @enddesc @calls @calledby @history @endhistory @@*/ subroutine RiemannSolvePolytypeM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLEM(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then !!$ Call the non-MHD version - see above call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) end if !!$ else if (CCTK_EQUALS(riemann_solver,"Roe")) then !!$ !!$ call GRHydro_RoeSolve(CCTK_PASS_FTOF) !!$ !!$ if (evolve_tracer .ne. 0) then !!$ !!$ call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) !!$ !!$ end if !!$ !!$ else if (CCTK_EQUALS(riemann_solver,"Marquina")) then !!$ !!$ call GRHydro_Marquina(CCTK_PASS_FTOF) !!$ Tracers are built directly in to the Marquina solver else call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!") end if end subroutine RiemannSolvePolytypeM /*@@ @routine RiemannSolveGeneralM @date Tue Mar 19 11:40:20 2002 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc The Riemann solvers for the new general EOS routines. This sets the fluxes from the left and right reconstructed states, so that after this routine they are effectively scratch space. @enddesc @calls @calledby @history @endhistory @@*/ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS) USE GRHydro_Scalars implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k, ierr CCTK_REAL, dimension(8) :: tmp_flux, cons_p, cons_m CCTK_REAL, dimension(6) :: prim_p, prim_m CCTK_REAL :: avg_det, avg_alp, avg_beta CCTK_REAL :: gxxh, gyyh, gzzh, gxyh, gxzh, gyzh CCTK_REAL :: avg_betax, avg_betay, avg_betaz CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm CCTK_REAL :: pressstarp,pressstarm,rhoenth_p,rhoenth_m CCTK_REAL :: velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm densflux = 0.d0 sxflux = 0.d0 syflux = 0.d0 szflux = 0.d0 tauflux = 0.d0 Bvecxflux = 0.d0 Bvecyflux = 0.d0 Bveczflux = 0.d0 if(clean_divergence.ne.0) then psidcflux = 0.d0 endif !!$ Do the EOS call to set the pressure, derivative and cs2 ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallPlus) ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallMinus) if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLEGeneralM(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then !!$ No b-field component for tracers! call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF) end if else 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 !!$ Set the left (p for plus) and right (m 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_p(6) = Bvecxplus(i,j,k) cons_p(7) = Bvecyplus(i,j,k) cons_p(8) = Bveczplus(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) cons_m(6) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) cons_m(7) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) cons_m(8) = Bveczminus(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_p(6) = pressplus(i,j,k) rhoenth_p = prim_p(1)*(1.0d0+prim_p(5))+prim_p(6) 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_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) rhoenth_m = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6) if(clean_divergence.ne.0) then psidcp = psidcplus(i,j,k) psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset) endif !!$ Set metric terms at interface if (shift_state .ne. 0) then avg_betax = 0.5d0*(betax(i+xoffset,j+yoffset,k+zoffset)+betax(i,j,k)) avg_betay = 0.5d0*(betay(i+xoffset,j+yoffset,k+zoffset)+betay(i,j,k)) avg_betaz = 0.5d0*(betaz(i+xoffset,j+yoffset,k+zoffset)+betaz(i,j,k)) if (flux_direction == 1) then avg_beta = avg_betax else if (flux_direction == 2) then avg_beta = avg_betay else if (flux_direction == 3) then avg_beta = avg_betaz else call CCTK_WARN(0, "Flux direction not x,y,z") end if else avg_beta = 0.d0 avg_betax = 0.d0 avg_betay = 0.d0 avg_betaz = 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) vxtp = prim_p(2)-avg_betax/avg_alp vytp = prim_p(3)-avg_betay/avg_alp vztp = prim_p(4)-avg_betaz/avg_alp vxtm = prim_m(2)-avg_betax/avg_alp vytm = prim_m(3)-avg_betay/avg_alp vztm = prim_m(4)-avg_betaz/avg_alp call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & prim_p(2),prim_p(3),prim_p(4),cons_p(6),cons_p(7),cons_p(8), & velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, & bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp) call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), & velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, & bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm) ab0p = wp*bdotvp ab0m = wm*bdotvm vA2p = b2p/(rhoenth_p+b2p) vA2m = b2m/(rhoenth_m+b2m) !!$ p^* = p+b^2/2 in Anton et al. pressstarp = prim_p(6)+0.5d0*b2p pressstarm = prim_m(6)+0.5d0*b2m if (flux_direction == 1) then call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& cons_p(6),cons_p(7),cons_p(8),& tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), & tmp_flux(6),tmp_flux(7),tmp_flux(8), & vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) densflux(i,j,k) = 0.5d0 * tmp_flux(1) sxflux(i,j,k) = 0.5d0 * tmp_flux(2) syflux(i,j,k) = 0.5d0 * tmp_flux(3) szflux(i,j,k) = 0.5d0 * tmp_flux(4) tauflux(i,j,k) = 0.5d0 * tmp_flux(5) Bvecxflux(i,j,k) = 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k) = 0.5d0 * tmp_flux(7) Bveczflux(i,j,k) = 0.5d0 * tmp_flux(8) if(clean_divergence.ne.0) then psidcf = cons_p(6) psidcflux(i,j,k) = 0.5d0 * psidcf endif call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& cons_m(6),cons_m(7),cons_m(8),& tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5),& tmp_flux(6),tmp_flux(7),tmp_flux(8),& vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1) sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2) syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3) szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4) tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5) Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7) Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8) if(clean_divergence.ne.0) then psidcf = cons_m(6) psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf endif else if (flux_direction == 2) then call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& cons_p(7),cons_p(8),cons_p(6),& tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),& tmp_flux(7),tmp_flux(8),tmp_flux(6),& vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) densflux(i,j,k) = 0.5d0 * tmp_flux(1) sxflux(i,j,k) = 0.5d0 * tmp_flux(2) syflux(i,j,k) = 0.5d0 * tmp_flux(3) szflux(i,j,k) = 0.5d0 * tmp_flux(4) tauflux(i,j,k) = 0.5d0 * tmp_flux(5) Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7) Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8) if(clean_divergence.ne.0) then psidcf = cons_p(7) psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf endif call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& cons_m(7),cons_m(8),cons_m(6),& tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),& tmp_flux(7),tmp_flux(8),tmp_flux(6),& vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1) sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2) syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3) szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4) tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5) Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7) Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8) if(clean_divergence.ne.0) then psidcf = cons_m(7) psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf endif else if (flux_direction == 3) then call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& cons_p(8),cons_p(6),cons_p(7),& tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5),& tmp_flux(8),tmp_flux(6),tmp_flux(7), & vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) densflux(i,j,k) = 0.5d0 * tmp_flux(1) sxflux(i,j,k) = 0.5d0 * tmp_flux(2) syflux(i,j,k) = 0.5d0 * tmp_flux(3) szflux(i,j,k) = 0.5d0 * tmp_flux(4) tauflux(i,j,k) = 0.5d0 * tmp_flux(5) Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7) Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8) if(clean_divergence.ne.0) then psidcf = cons_p(8) psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf endif call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& cons_m(8),cons_m(6),cons_m(7),& tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), & tmp_flux(8),tmp_flux(6),tmp_flux(7), & vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & avg_det,avg_alp,avg_beta) densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1) sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2) syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3) szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4) tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5) Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6) Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7) Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8) if(clean_divergence.ne.0) then psidcf = cons_m(8) psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf endif else call CCTK_WARN(0, "Flux direction not x,y,z") end if end do end do end do if (CCTK_EQUALS(riemann_solver,"Roe")) then call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!") !!$ !!$ call GRHydro_RoeSolveGeneral(CCTK_PASS_FTOF) !!$ !!$ if (evolve_tracer .ne. 0) then !!$ !!$ call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF) !!$ !!$ end if !!$ else if (CCTK_EQUALS(riemann_solver,"Marquina")) then call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!") !!$ !!$ call GRHydro_MarquinaGeneral(CCTK_PASS_FTOF) !!$ !!$ Tracers are built directly in to the Marquina solver end if end if end subroutine RiemannSolveGeneralM