From d95f7bc4e19ff9d991e17417b63318ea63d18491 Mon Sep 17 00:00:00 2001 From: bmundim Date: Wed, 29 Sep 2010 21:47:21 +0000 Subject: Current RIT GRMHD code contributions: Add the magnetized counterparts for several GRHydro routines. Adjust interface.ccl, param.ccl and schedule.ccl appropriately. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@158 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_HLLEM.F90 | 642 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 642 insertions(+) create mode 100644 src/GRHydro_HLLEM.F90 (limited to 'src/GRHydro_HLLEM.F90') diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 new file mode 100644 index 0000000..98c4874 --- /dev/null +++ b/src/GRHydro_HLLEM.F90 @@ -0,0 +1,642 @@ + /*@@ + @file GRHydro_HLLEM.F90 + @date Aug 30, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, 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_HLLEGeneralM + @date Aug 30, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, 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_HLLEGeneralM(CCTK_ARGUMENTS) + USE GRHydro_EigenproblemM + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: i, j, k, m + CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff + CCTK_REAL, dimension(6) :: prim_p, prim_m + CCTK_REAL, dimension(5) :: lamminus,lamplus + 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_REAL :: rhoenth_p, rhoenth_m, 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,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm + + 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_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_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) + rhoenth_p = prim_p(1)*(1.0d0+prim_p(5))+prim_p(6) + + 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) + rhoenth_m = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6) + +!!$ Calculate various metric terms here. +!!$ Note also need the average of the lapse at the +!!$ left and right points. +!!$ +!!$ In MHD, we need all three shift components regardless of the flux direction + + 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 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 + +!!$ we now pass in the B-field as conserved and flux, the vtilde's instead of v's, +!!$ pressstar instead of P, b_i, alp b^0, w, metric determinant, +!!$ alp, and beta in the flux dir + + if (flux_direction == 1) then + 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),& + f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),& + vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + 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),& + f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),& + vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + 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),& + f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), & + vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & + 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) + qdiff(6) = cons_m(6) - cons_p(6) + qdiff(7) = cons_m(7) - cons_p(7) + qdiff(8) = cons_m(8) - cons_p(8) + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + call eigenvalues_generalM(& + prim_m(2),prim_m(3),prim_m(4),cs2_m,vA2m, & + lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues_generalM(& + prim_p(2),prim_p(3),prim_p(4),cs2_p,vA2p, & + lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + 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),& + fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),& + vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + 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),& + fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),& + fminus(6),fminus(7),fminus(8),& + vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + call eigenvalues_generalM(& + prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, & + lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues_generalM(& + prim_p(3),prim_p(4),prim_p(2),cs2_p,vA2p, & + lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + 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),& + fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),& + vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + 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),& + fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),& + fminus(7),fminus(8),fminus(6),& + vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + call eigenvalues_generalM(& + prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, & + lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues_generalM(& + prim_p(4),prim_p(2),prim_p(3),cs2_p,vA2p, & + lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + 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),& + fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), & + vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + 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),& + fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), & + fminus(8),fminus(6),fminus(7), & + vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & + 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,8 + + 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,8 + + 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) + Bvecxflux(i, j, k) = f1(6) + Bvecyflux(i, j, k) = f1(7) + Bveczflux(i, j, k) = f1(8) + + end do + end do + end do + +end subroutine GRHydro_HLLEGeneralM + + +subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS) + + USE GRHydro_EigenproblemM + + 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, dimension(3) :: mag_p, mag_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_REAL :: b2p,b2m,vA2m,vA2p + + 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 + mag_p = 0.d0 + mag_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,:) + + mag_p(1) = Bvecxplus(i,j,k) + mag_p(2) = Bvecyplus(i,j,k) + mag_p(3) = Bveczplus(i,j,k) + + mag_m(1) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) + mag_m(2) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) + mag_m(3) = 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_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 + +!!$ b^2 = (1-v^2)B^2+(B dot v)^2 + b2p=(1.d0-DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4)))* & + DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_p(1),mag_p(2),mag_p(3)) + & + (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4),mag_p(1),mag_p(2),mag_p(3)))**2 + b2m=(1.d0-DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4)))* & + DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3)) + & + (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3)))**2 + + vA2p = b2p/(prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)+b2p) + vA2m = b2m/(prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)+b2m) + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + call eigenvalues_generalM(& + prim_m(2),prim_m(3),prim_m(4),cs2_m,vA2m, & + lamminus,& + gxxh,gxyh,gxzh,& + gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues_generalM(& + prim_p(2),prim_p(3),prim_p(4),cs2_p,vA2p, & + 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_generalM(& + prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, & + lamminus,& + gyyh,gyzh,gxyh,& + gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues_generalM(& + prim_p(3),prim_p(4),prim_p(2),cs2_p,vA2p, & + 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_generalM(& + prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, & + lamminus,& + gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues_generalM(& + prim_p(4),prim_p(2),prim_p(3),cs2_p,vA2p, & + 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_TracerGeneralM + -- cgit v1.2.3