From e6322a50c9ceb945b7c2ec9966b33da831662882 Mon Sep 17 00:00:00 2001 From: cott Date: Mon, 6 Dec 2010 21:13:20 +0000 Subject: * rename HLLE routine source code files to reflect true content git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@189 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_HLLE.F90 | 669 +++++++++++++++++++++++++++------------------------ 1 file changed, 357 insertions(+), 312 deletions(-) (limited to 'src/GRHydro_HLLE.F90') diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90 index 7a3f5f8..bd976e1 100644 --- a/src/GRHydro_HLLE.F90 +++ b/src/GRHydro_HLLE.F90 @@ -12,16 +12,16 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" + #include "GRHydro_Macros.h" #include "SpaceMask.h" /*@@ - @routine GRHydro_HLLEGeneral + @routine GRHydro_HLLE @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 @@ -31,36 +31,25 @@ @@*/ -subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) - +subroutine GRHydro_HLLE(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) :: consp,consm_i,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 + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r + CCTK_REAL :: xtemp 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", & @@ -84,45 +73,25 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) f1 = 0.d0 lamminus = 0.d0 lamplus = 0.d0 - cons_p = 0.d0 - cons_m = 0.d0 + consp = 0.d0 + consm_i = 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) + consp(1) = densplus(i,j,k) + consp(2) = sxplus(i,j,k) + consp(3) = syplus(i,j,k) + consp(4) = szplus(i,j,k) + consp(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) + consm_i(1) = densminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(3) = syminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(4) = szminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) !!$ Calculate various metric terms here. !!$ Note also need the average of the lapse at the @@ -154,27 +123,31 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) 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) + 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), & + call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),& + f1(1),f1(2),f1(3),& + f1(4),f1(5),& + velxplus(i,j,k),pressplus(i,j,k),& 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), & + call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),& + f1(1),f1(3),f1(4),& + f1(2),f1(5),& + velyplus(i,j,k),pressplus(i,j,k),& 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), & + call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),& + f1(1),f1(4),f1(2),& + f1(3),f1(5),& + velzplus(i,j,k),pressplus(i,j,k),& avg_det,avg_alp,avg_beta) else call CCTK_WARN(0, "Flux direction not x,y,z") @@ -183,8 +156,8 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) 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) + avg_det,gxxh, gxyh, gxzh, & + gyyh, gyzh, gzzh) if (flux_direction == 1) then usendh = uxxh @@ -198,127 +171,201 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) !!$ 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(1) = consm_i(1) - consp(1) + qdiff(2) = consm_i(2) - consp(2) + qdiff(3) = consm_i(3) - consp(3) + qdiff(4) = consm_i(4) - consp(4) + qdiff(5) = consm_i(5) - consp(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), & + if(evolve_temper.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,i,j,k,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + xtemp,& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,i,j,k,& + rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k), & + xtemp,y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + endif + call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),& + fplus(1),fplus(2),fplus(3),fplus(4),& + fplus(5),velxplus(i,j,k),pressplus(i,j,k),& 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), & + call num_x_flux(consm_i(1),consm_i(2),consm_i(3),& + consm_i(4),consm_i(5),fminus(1),fminus(2),fminus(3),& + fminus(4), fminus(5),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + pressminus(i+xoffset,j+yoffset,k+zoffset),& 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), & + if(evolve_temper.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,i,j,k,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + xtemp,& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,i,j,k,& + rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),& + xtemp,y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + endif + call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),& + fplus(1),fplus(3),fplus(4),fplus(2),& + fplus(5),velyplus(i,j,k),pressplus(i,j,k),& 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), & + call num_x_flux(consm_i(1),consm_i(3),consm_i(4),& + consm_i(2),consm_i(5),fminus(1),fminus(3),fminus(4),& + fminus(2), fminus(5),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + pressminus(i+xoffset,j+yoffset,k+zoffset),& 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), & + if(evolve_temper.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,i,j,k,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,i,j,k,& + rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),& + xtemp,y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + endif + call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),& + fplus(1),fplus(4),fplus(2),fplus(3),& + fplus(5),velzplus(i,j,k),pressplus(i,j,k),avg_det,& + avg_alp,avg_beta) + call num_x_flux(consm_i(1),consm_i(4),consm_i(2),& + consm_i(3),consm_i(5),fminus(1),fminus(4),fminus(2),& + fminus(3), fminus(5),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + pressminus(i+xoffset,j+yoffset,k+zoffset),& 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)) + 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)) + 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 + 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 + do m = 1,5 - 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) = consm_i(m) - consp(m) - qdiff(m) = cons_m(m) - cons_p(m) - - f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* & - qdiff(m) - - end do + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm - end if - + end do 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) @@ -330,22 +377,36 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) 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 + 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 do + +end subroutine GRHydro_HLLE -end subroutine GRHydro_HLLEGeneral + /*@@ + @routine GRHydro_HLLE_Tracer + @date Mon Mar 8 13:47:13 2004 + @author Ian Hawke + @desc + HLLE just for the tracer. + @enddesc + @calls + @calledby + @history + + @endhistory +@@*/ -subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) +subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) USE GRHydro_Eigenproblem @@ -353,26 +414,32 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) 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 + integer :: i, j, k, m + CCTK_REAL, dimension(number_of_tracers) :: consp,consm_i,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 + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r + + CCTK_INT :: type_bits, trivial, not_trivial - if(CCTK_EQUALS(HLLE_type,"Tadmor")) then - tadmor = 1 + 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 - tadmor = 0 - endif - + 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 @@ -381,37 +448,18 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) f1 = 0.d0 lamminus = 0.d0 lamplus = 0.d0 - cons_p = 0.d0 - cons_m = 0.d0 + consp = 0.d0 + consm_i = 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) - + do m=1,number_of_tracers + consp(m) = cons_tracerplus(i,j,k,m) + consm_i(m) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,m) + enddo !!$ Calculate various metric terms here. !!$ Note also need the average of the lapse at the !!$ left and right points. @@ -442,135 +490,132 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) 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) + 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 + usendh = uxxh else if (flux_direction == 2) then - usendh = uyyh + usendh = uyyh else if (flux_direction == 3) then - usendh = uzzh + usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_WARN(0, "Flux direction not x,y,z") end if -!!$ Eigenvalues and fluxes either side of the cell interface +!!$ Calculate the jumps in the conserved variables - 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 + qdiff = consm_i - consp + +!!$ Eigenvalues and fluxes either side of the cell interface + if (flux_direction == 1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + 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(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + 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(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + 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 + !!$ 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)) + 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 + 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 + do m = 1,number_of_tracers - - end if - + qdiff(m) = consm_i(m) - consp(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + + end do - cons_tracerflux(i,j,k,:) = f1(:) - + cons_tracerflux(i, j, k,:) = f1(:) +!!$ +!!$ if ( ((flux_direction.eq.3).and.(i.eq.4).and.(j.eq.4)).or.& +!!$ ((flux_direction.eq.2).and.(i.eq.4).and.(k.eq.4)).or.& +!!$ ((flux_direction.eq.1).and.(j.eq.4).and.(k.eq.4))& +!!$ ) then +!!$ write(*,*) flux_direction, i, j, k, f1(1), consm_i(1), consp(1) +!!$ end if + end do end do end do -end subroutine GRHydro_HLLE_TracerGeneral - - - + +end subroutine GRHydro_HLLE_Tracer -- cgit v1.2.3