/*@@ @file GRHydro_RiemannSolve.F90 @date Sat Jan 26 02:20:25 2002 @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 RiemannSolve @date Sat Jan 26 02:20:48 2002 @author Pedro Montero, Ian Hawke @desc A wrapper routine to switch between the different Riemann solvers. @enddesc @calls @calledby @history @endhistory @@*/ subroutine RiemannSolve(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_HLLE(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then 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 end if end subroutine RiemannSolve /*@@ @routine RiemannSolvePolytype @date Tue Mar 19 11:40:20 2002 @author 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 RiemannSolvePolytype(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_HLLE(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then 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 end if end subroutine RiemannSolvePolytype /*@@ @routine RiemannSolveGeneral @date Tue Mar 19 11:40:20 2002 @author 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 RiemannSolveGeneral(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(5) :: tmp_flux, cons_p, cons_m CCTK_REAL, dimension(6) :: prim_p, prim_m CCTK_REAL :: avg_det, avg_alp, avg_beta CCTK_REAL :: psi4h, gxxh, gyyh, gzzh, gxyh, gxzh, gyzh densflux = 0.d0 sxflux = 0.d0 syflux = 0.d0 szflux = 0.d0 tauflux = 0.d0 !!$ 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_HLLEGeneral(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then 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_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_p(6) = pressplus(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_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) !!$ Set metric terms at interface 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 (flux_direction == 1) then call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), & prim_p(2),prim_p(6), & 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) call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), & prim_m(2),prim_m(6), & 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) else if (flux_direction == 2) then call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5), & prim_p(3),prim_p(6), & 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) call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5), & prim_m(3),prim_m(6), & 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) else if (flux_direction == 3) then call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), & prim_p(4),prim_p(6), & 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) call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), & prim_m(4),prim_m(6), & 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) 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 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 GRHydro_MarquinaGeneral(CCTK_PASS_FTOF) !!$ Tracers are built directly in to the Marquina solver end if end if end subroutine RiemannSolveGeneral