aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorcott <cott@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45>2010-12-06 21:13:20 +0000
committercott <cott@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45>2010-12-06 21:13:20 +0000
commite6322a50c9ceb945b7c2ec9966b33da831662882 (patch)
tree07f7e189e781b22ce0448aa6c99c8db7a0c7c15c
parent9b147a8c397b5bb7a0c1e8c1537f870474c71a74 (diff)
* 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
Diffstat
-rw-r--r--src/GRHydro_HLLE.F90669
-rw-r--r--src/GRHydro_HLLEM.F90383
-rw-r--r--src/GRHydro_HLLEPoly.F90621
-rw-r--r--src/GRHydro_HLLE_EOSG.F90576
-rw-r--r--src/GRHydro_HLLE_EOSGM.F90 (renamed from src/GRHydro_HLLEPolyM.F90)383
-rw-r--r--src/make.code.defn6
6 files changed, 1319 insertions, 1319 deletions
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
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90
index 8fe2ef0..aec2899 100644
--- a/src/GRHydro_HLLEM.F90
+++ b/src/GRHydro_HLLEM.F90
@@ -1,5 +1,5 @@
/*@@
- @file GRHydro_HLLEM.F90
+ @file GRHydro_HLLEPolyM.F90
@date Aug 30, 2010
@author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke, Pedro Montero, Toni Font
@desc
@@ -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_HLLEGeneralM
+ @routine GRHydro_HLLEM
@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
@@ -31,18 +31,17 @@
@@*/
-subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
+subroutine GRHydro_HLLEM(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(7) :: 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, &
@@ -58,15 +57,6 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
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", &
@@ -126,28 +116,22 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
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_p(7) = w_lorentzplus(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)
-
+ prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
+ prim_m(7) = w_lorentzminus(i+xoffset,j+yoffset,k+zoffset)
+
if(clean_divergence.ne.0) then
psidcp = psidcplus(i,j,k)
psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset)
endif
-
+
!!$ Calculate various metric terms here.
!!$ Note also need the average of the lapse at the
!!$ left and right points.
@@ -186,7 +170,7 @@ subroutine GRHydro_HLLEGeneralM(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)
vxtp = prim_p(2)-avg_betax/avg_alp
vytp = prim_p(3)-avg_betay/avg_alp
@@ -225,30 +209,30 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
!!$ 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),&
+ 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),&
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, &
+ vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
psidcf = cons_p(6)
endif
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),&
+ 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),&
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, &
+ vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
psidcf = cons_p(7)
endif
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),&
+ 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),&
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, &
+ vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
psidcf = cons_p(8)
@@ -292,12 +276,14 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
!!$ 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, &
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), &
+ cons_m(6),cons_m(7),cons_m(8),&
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, &
+ call eigenvaluesM(GRHydro_eos_handle, &
+ prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), &
+ cons_p(6),cons_p(7),cons_p(8),&
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),&
@@ -318,12 +304,14 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
endif
else if (flux_direction == 2) then
- call eigenvalues_generalM(&
- prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, &
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), &
+ cons_m(7),cons_m(8),cons_m(6),&
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, &
+ call eigenvaluesM(GRHydro_eos_handle, &
+ prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), &
+ cons_p(7),cons_p(8),cons_p(6),&
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),&
@@ -343,12 +331,14 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
endif
else if (flux_direction == 3) then
- call eigenvalues_generalM(&
- prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, &
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), &
+ cons_m(8),cons_m(6),cons_m(7),&
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, &
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), &
+ cons_p(8),cons_p(6),cons_p(7),&
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),&
@@ -370,66 +360,38 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
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,8
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
-
- end do
-
- if(clean_divergence.ne.0) then
- psidcdiff = psidcm - psidcp
- psidcf = (charmax * psidcfp - charmin * psidcfm + &
- charmax * charmin * psidcdiff) / charpm
- endif
-
- 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)
- 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
- if(clean_divergence.ne.0) then
- psidcdiff = psidcm - psidcp
- psidcf = 0.5d0 * (psidcfp + psidcfm) - 0.5d0*charmax* &
- psidcdiff
- endif
-
- end if
-
+ end do
+ if(clean_divergence.ne.0) then
+ psidcdiff = psidcm - psidcp
+ psidcf = (charmax * psidcfp - charmin * psidcfm + &
+ charmax * charmin * psidcdiff) / charpm
+ endif
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)
@@ -439,13 +401,13 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
Bvecyflux(i, j, k) = f1(7)
Bveczflux(i, j, k) = f1(8)
psidcflux(i,j,k) = psidcf
-
+
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
+ else
Y_e_con_flux(i, j, k) = &
Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * &
densflux(i, j, k)
@@ -456,10 +418,24 @@ subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
end do
end do
-end subroutine GRHydro_HLLEGeneralM
+end subroutine GRHydro_HLLEM
+
+ /*@@
+ @routine GRHydro_HLLE_TracerM
+ @date Aug 30, 2010
+ @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
+ @desc
+ HLLE just for the tracer.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+@@*/
-subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
+subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS)
USE GRHydro_EigenproblemM
@@ -467,28 +443,36 @@ subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
- integer :: i, j, k, m
+ 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(7) :: 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
+
+ CCTK_INT :: type_bits, trivial, not_trivial
- integer tadmor
-
- 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
@@ -523,21 +507,17 @@ subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
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_p(7) = w_lorentzplus(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)
+ prim_m(7) = w_lorentzminus(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.
@@ -568,84 +548,77 @@ subroutine GRHydro_HLLE_TracerGeneralM(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
!!$ 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)) + &
+ b2p=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_p(1),mag_p(2),mag_p(3))/prim_p(7)**2 + &
(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)) + &
+ b2m=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3))/prim_m(7)**2 + &
(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)
-
+
+!!$ Calculate the jumps in the conserved variables
+
+ qdiff = cons_m - cons_p
+
!!$ 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,&
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), &
+ cons_m(6),cons_m(7),cons_m(8),&
+ 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,&
+ call eigenvaluesM(GRHydro_eos_handle, &
+ prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), &
+ cons_p(6),cons_p(7),cons_p(8),&
+ 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,&
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), &
+ cons_m(7),cons_m(8),cons_m(6),&
+ 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,&
+ call eigenvaluesM(GRHydro_eos_handle, &
+ prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), &
+ cons_p(7),cons_p(8),cons_p(6),&
+ 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,&
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), &
+ cons_m(8),cons_m(6),cons_m(7),&
+ 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,&
+ call eigenvaluesM(GRHydro_eos_handle,&
+ prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), &
+ cons_p(8),cons_p(6),cons_p(7),&
+ 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,:)
@@ -654,57 +627,43 @@ subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
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
- charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
+ qdiff(m) = cons_m(m) - cons_p(m)
- 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
+ f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
+ charmax * charmin * qdiff(m)) / charpm
+
+ end do
+
+ 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), cons_m(1), cons_p(1)
+!!$ end if
+
+ end do
+ end do
+end do
+
+
+end subroutine GRHydro_HLLE_TracerM
diff --git a/src/GRHydro_HLLEPoly.F90 b/src/GRHydro_HLLEPoly.F90
deleted file mode 100644
index bd976e1..0000000
--- a/src/GRHydro_HLLEPoly.F90
+++ /dev/null
@@ -1,621 +0,0 @@
- /*@@
- @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_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.
- @enddesc
- @calls
- @calledby
- @history
- Altered from Cactus 3 routines originally written by Toni Font.
- @endhistory
-
-@@*/
-
-subroutine GRHydro_HLLE(CCTK_ARGUMENTS)
- USE GRHydro_Eigenproblem
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
-
- integer :: i, j, k, m
- CCTK_REAL, dimension(5) :: consp,consm_i,fplus,fminus,lamplus
- CCTK_REAL, dimension(5) :: f1,lamminus
- CCTK_REAL, dimension(5) :: qdiff
- 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
- CCTK_REAL :: xtemp
-
- CCTK_INT :: type_bits, trivial, not_trivial
-
- 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
- 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
-
- 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)
-
- 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
-!!$ 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(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(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(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")
- 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) = 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
- 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(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
- 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(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
- 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
-
-!!$ 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) = consm_i(m) - consp(m)
-
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
-
- 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)
- 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_HLLE
-
- /*@@
- @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_Tracer(CCTK_ARGUMENTS)
-
- USE GRHydro_Eigenproblem
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
-
- 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 :: 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
-
- CCTK_INT :: type_bits, trivial, not_trivial
-
- 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
- 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
-
- 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.
-
- 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)
-
- 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 = 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))
-
- 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) = 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(:)
-!!$
-!!$ 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_Tracer
-
diff --git a/src/GRHydro_HLLE_EOSG.F90 b/src/GRHydro_HLLE_EOSG.F90
new file mode 100644
index 0000000..7a3f5f8
--- /dev/null
+++ b/src/GRHydro_HLLE_EOSG.F90
@@ -0,0 +1,576 @@
+ /*@@
+ @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
+
+
+
+
diff --git a/src/GRHydro_HLLEPolyM.F90 b/src/GRHydro_HLLE_EOSGM.F90
index aec2899..8fe2ef0 100644
--- a/src/GRHydro_HLLEPolyM.F90
+++ b/src/GRHydro_HLLE_EOSGM.F90
@@ -1,5 +1,5 @@
/*@@
- @file GRHydro_HLLEPolyM.F90
+ @file GRHydro_HLLEM.F90
@date Aug 30, 2010
@author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke, Pedro Montero, Toni Font
@desc
@@ -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_HLLEM
+ @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
@@ -31,17 +31,18 @@
@@*/
-subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
+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(7) :: prim_p, prim_m
+ 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, &
@@ -57,6 +58,15 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
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", &
@@ -116,22 +126,28 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
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_p(7) = w_lorentzplus(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)
- prim_m(7) = w_lorentzminus(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)
+
if(clean_divergence.ne.0) then
psidcp = psidcplus(i,j,k)
psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset)
endif
-
+
!!$ Calculate various metric terms here.
!!$ Note also need the average of the lapse at the
!!$ left and right points.
@@ -170,7 +186,7 @@ subroutine GRHydro_HLLEM(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)
vxtp = prim_p(2)-avg_betax/avg_alp
vytp = prim_p(3)-avg_betay/avg_alp
@@ -209,30 +225,30 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
!!$ alp, and beta in the flux dir
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),&
+ 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),&
- vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
+ vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
psidcf = cons_p(6)
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),&
+ 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),&
- vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
+ vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
psidcf = cons_p(7)
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),&
+ 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), &
- vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
+ vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
psidcf = cons_p(8)
@@ -276,14 +292,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
!!$ Eigenvalues and fluxes either side of the cell interface
if (flux_direction == 1) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), &
- cons_m(6),cons_m(7),cons_m(8),&
+ 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 eigenvaluesM(GRHydro_eos_handle, &
- prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), &
- cons_p(6),cons_p(7),cons_p(8),&
+ 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),&
@@ -304,14 +318,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
endif
else if (flux_direction == 2) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), &
- cons_m(7),cons_m(8),cons_m(6),&
+ 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 eigenvaluesM(GRHydro_eos_handle, &
- prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), &
- cons_p(7),cons_p(8),cons_p(6),&
+ 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),&
@@ -331,14 +343,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
endif
else if (flux_direction == 3) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), &
- cons_m(8),cons_m(6),cons_m(7),&
+ 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 eigenvaluesM(GRHydro_eos_handle,&
- prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), &
- cons_p(8),cons_p(6),cons_p(7),&
+ 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),&
@@ -360,38 +370,66 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
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,8
- qdiff(m) = cons_m(m) - cons_p(m)
+ qdiff(m) = cons_m(m) - cons_p(m)
+
+ f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
+ charmax * charmin * qdiff(m)) / charpm
+
+ end do
+
+ if(clean_divergence.ne.0) then
+ psidcdiff = psidcm - psidcp
+ psidcf = (charmax * psidcfp - charmin * psidcfm + &
+ charmax * charmin * psidcdiff) / charpm
+ endif
+
+ 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
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
+ qdiff(m) = cons_m(m) - cons_p(m)
+
+ f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
+ qdiff(m)
+
+ end do
- end do
+ if(clean_divergence.ne.0) then
+ psidcdiff = psidcm - psidcp
+ psidcf = 0.5d0 * (psidcfp + psidcfm) - 0.5d0*charmax* &
+ psidcdiff
+ endif
+
+ end if
+
- if(clean_divergence.ne.0) then
- psidcdiff = psidcm - psidcp
- psidcf = (charmax * psidcfp - charmin * psidcfm + &
- charmax * charmin * psidcdiff) / charpm
- endif
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)
@@ -401,13 +439,13 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
Bvecyflux(i, j, k) = f1(7)
Bveczflux(i, j, k) = f1(8)
psidcflux(i,j,k) = psidcf
-
+
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
+ else
Y_e_con_flux(i, j, k) = &
Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * &
densflux(i, j, k)
@@ -418,24 +456,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
end do
end do
-end subroutine GRHydro_HLLEM
-
- /*@@
- @routine GRHydro_HLLE_TracerM
- @date Aug 30, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
- @desc
- HLLE just for the tracer.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
+end subroutine GRHydro_HLLEGeneralM
-@@*/
-subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS)
+subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
USE GRHydro_EigenproblemM
@@ -443,36 +467,28 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
- integer :: i, j, k, m
+ 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(7) :: prim_p, prim_m
+ 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
-
- CCTK_INT :: type_bits, trivial, not_trivial
- 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")
+ integer tadmor
+
+ if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
+ tadmor = 1
else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
+ tadmor = 0
+ endif
+
do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
@@ -507,17 +523,21 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS)
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_p(7) = w_lorentzplus(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)
- prim_m(7) = w_lorentzminus(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.
@@ -548,77 +568,84 @@ subroutine GRHydro_HLLE_TracerM(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
!!$ b^2 = (1-v^2)B^2+(B dot v)^2
- b2p=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_p(1),mag_p(2),mag_p(3))/prim_p(7)**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=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3))/prim_m(7)**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)
-
-!!$ Calculate the jumps in the conserved variables
-
- qdiff = cons_m - cons_p
-
+
!!$ Eigenvalues and fluxes either side of the cell interface
-
+
if (flux_direction == 1) then
- call eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), &
- cons_m(6),cons_m(7),cons_m(8),&
- lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
+ 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 eigenvaluesM(GRHydro_eos_handle, &
- prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), &
- cons_p(6),cons_p(7),cons_p(8),&
- lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
+ 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 eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), &
- cons_m(7),cons_m(8),cons_m(6),&
- lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
+ 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 eigenvaluesM(GRHydro_eos_handle, &
- prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), &
- cons_p(7),cons_p(8),cons_p(6),&
- lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
+ 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 eigenvaluesM(GRHydro_eos_handle,&
- prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), &
- cons_m(8),cons_m(6),cons_m(7),&
- lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
+ 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 eigenvaluesM(GRHydro_eos_handle,&
- prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), &
- cons_p(8),cons_p(6),cons_p(7),&
- lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
+ 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,:)
@@ -627,43 +654,57 @@ subroutine GRHydro_HLLE_TracerM(CCTK_ARGUMENTS)
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))
-
- 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
+
+ if(tadmor.eq.0) then
- qdiff(m) = cons_m(m) - cons_p(m)
+!!$ Find minimum and maximum wavespeeds
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
+ charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
+ lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
+ lamminus(4),lamminus(5))
- end do
-
- 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), cons_m(1), cons_p(1)
-!!$ end if
-
- end do
- end do
-end do
-
-
-end subroutine GRHydro_HLLE_TracerM
+ 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
diff --git a/src/make.code.defn b/src/make.code.defn
index 53fe17c..2a00a1b 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -14,8 +14,8 @@ SRCS = Utils.F90 \
GRHydro_EOS.c \
GRHydro_Flux.F90 \
GRHydro_FluxSplit.F90 \
+ GRHydro_HLLE_EOSG.F90 \
GRHydro_HLLE.F90 \
- GRHydro_HLLEPoly.F90 \
GRHydro_Loop.F90 \
GRHydro_Minima.F90 \
GRHydro_Minima.cc \
@@ -52,9 +52,9 @@ SRCS = Utils.F90 \
GRHydro_Con2PrimM_pt.c \
GRHydro_EigenproblemM.F90 \
GRHydro_FluxM.F90 \
+ GRHydro_HLLE_EOSGM.F90 \
GRHydro_HLLEM.F90 \
- GRHydro_HLLEPolyM.F90 \
- GRHydro_PPMM.F90 \
+ GRHydro_PPMM.F90 \
GRHydro_ParamCheckM.F90 \
GRHydro_Prim2ConM.F90 \
GRHydro_RiemannSolveM.F90 \
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 12:02:43 +0000