aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorrhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45>2011-09-15 16:49:53 +0000
committerrhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45>2011-09-15 16:49:53 +0000
commite0dc2af4862d5ddb874328bd097f7f516231dd8c (patch)
treed2ee226c88cbe662530d9b847b2a4fe6d72ef32b
parent374f96eabf2c562985209711c05b68624cd866e1 (diff)
add Multipatch support to GRHydro
* not all features of GRHydro are supported yet, in particular only the HLLE solver supports Mulitpatch yet. Original commit by Christian Reisswig and Christian Ott git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@273 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat
-rw-r--r--interface.ccl43
-rw-r--r--param.ccl4
-rw-r--r--schedule.ccl57
-rw-r--r--src/GRHydro_Boundaries.F9024
-rw-r--r--src/GRHydro_Con2Prim.F90114
-rw-r--r--src/GRHydro_ENOReconstruct_drv.F9012
-rw-r--r--src/GRHydro_EoSChangeGamma.F9022
-rw-r--r--src/GRHydro_HLLE.F9048
-rw-r--r--src/GRHydro_Jacobian_state.c32
-rw-r--r--src/GRHydro_PPMMReconstruct_drv.F9024
-rw-r--r--src/GRHydro_PPMReconstruct_drv.F9018
-rw-r--r--src/GRHydro_Prim2Con.F90224
-rw-r--r--src/GRHydro_RegisterVars.cc8
-rw-r--r--src/GRHydro_Source.F90432
-rw-r--r--src/GRHydro_TVDReconstruct_drv.F9024
-rw-r--r--src/GRHydro_TransformTensorBasis.F90211
-rw-r--r--src/GRHydro_UpdateMask.F9076
-rw-r--r--src/make.code.defn13
18 files changed, 861 insertions, 525 deletions
diff --git a/interface.ccl b/interface.ccl
index fc74a19..3a347bd 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -7,7 +7,7 @@
####################################################################
implements: GRHydro
-inherits: ADMBase, Boundary, SpaceMask, ADMMacros, Tmunubase, HydroBase
+inherits: ADMBase, Boundary, SpaceMask, ADMMacros, Tmunubase, HydroBase, Coordinates
#INCLUDES SOURCE: GRHydro_EMTensor.inc in CalcTmunu.inc
#INCLUDES: GRHydro_EMTensor_temps.inc in CalcTmunu_temps.inc
@@ -325,15 +325,15 @@ CCTK_REAL GRHydro_minima type = SCALAR
# GRHydro_dens_min
} "Atmosphere values"
-real dens type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "generalized particle number"
+real dens type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "generalized particle number"
-real tau type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "internal energy"
+real tau type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "internal energy"
-real scon[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "generalized momenta"
+real scon[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "generalized momenta"
-real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "B-field conservative variable"
+real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "B-field conservative variable"
-real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "Conserved electron fraction"
+real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Conserved electron fraction"
#real bcom[3] type = GF Timelevels = 3 tags='Prolongation="none" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "comoving magnetic field components"
@@ -346,7 +346,7 @@ real GRHydro_tracers[number_of_tracers] type = GF Timelevels = 3 tags='Prolongat
#real w_lorentz type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "Lorentz factor"
-real psidc type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 interpolator="matter"' "Psi parameter for divergence cleaning"
+real psidc type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Psi parameter for divergence cleaning"
real densrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for dens"
real taurhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for tau"
@@ -358,6 +358,35 @@ real psidcrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"
real divB type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Magnetic field constraint"
##################################################
+### variables in the local tensor basis
+##################################################
+
+CCTK_REAL lvel[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" jacobian="jacobian" interpolator="matter"' "local velocity v^i"
+CCTK_REAL lBvec[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" jacobian="jacobian" tensorparity=-1 interpolator="matter"' "local Magnetic field components B^i"
+
+# We should only need one timelevel. However, InitialAtmosphereReset also wants to set past timelevels!
+# There must be a better way!
+CCTK_REAL local_metric type = GF Timelevels = 3 tags='Prolongation="None" checkpoint="no"'
+{
+ gaa, gab, gac
+ gbb, gbc
+ gcc
+} "local ADM metric g_ij"
+
+CCTK_REAL local_extrinsic_curvature type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"'
+{
+ kaa, kab, kac
+ kbb, kbc
+ kcc
+} "local extrinsic curvature K_ij"
+
+
+CCTK_REAL local_shift type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"'
+{
+ betaa, betab, betac
+} "local ADM shift \beta^i"
+
+##################################################
### These variables are only protected so that ###
### the tests in init_data work. Should fix. ###
##################################################
diff --git a/param.ccl b/param.ccl
index 8d5f459..7b5049a 100644
--- a/param.ccl
+++ b/param.ccl
@@ -81,8 +81,8 @@ CCTK_INT GRHydro_MaxNumEvolvedVars "The maximum number of evolved variables used
CCTK_INT GRHydro_MaxNumConstrainedVars "The maximum number of constrained variables used by GRHydro" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars
{
- 7:15 :: "A small range, depending on testing or not"
-} 9
+ 7:33 :: "A small range, depending on testing or not"
+} 27
CCTK_INT GRHydro_MaxNumSandRVars "The maximum number of save and restore variables used by GRHydro" ACCUMULATOR-BASE=MethodofLines::MoL_Num_SaveAndRestore_Vars
{
diff --git a/schedule.ccl b/schedule.ccl
index f1de45f..1c64d4f 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -181,6 +181,18 @@ if (!CCTK_Equals(initial_shift,"none"))
}
##############################################
+### Storage for local tensor quantities ###
+##############################################
+
+STORAGE: lvel[3]
+if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) {
+ STORAGE: lBvec[3]
+}
+STORAGE: local_metric[3]
+STORAGE: local_extrinsic_curvature
+STORAGE: local_shift
+
+##############################################
### Storage for the conformal state scalar ###
##############################################
@@ -214,6 +226,12 @@ schedule GRHydro_ParamCheck AT PARAMCHECK
LANG: Fortran
} "Check parameters"
+schedule GRHydro_check_Jacobian_state AT BASEGRID AFTER (TmunuBase_SetStressEnergyState Coordinates_SetGlobalCoords_Group)
+{
+ LANG: C
+ OPTIONS: GLOBAL
+} "Test state of Jacobians"
+
######################################
### Standard symmetry registration ###
######################################
@@ -385,6 +403,7 @@ if (EoS_Change)
SYNC: hydrobase::press
SYNC: hydrobase::eps
SYNC: hydrobase::vel
+ SYNC: lvel
SYNC: hydrobase::temperature
SYNC: hydrobase::entropy
SYNC: hydrobase::Y_e
@@ -481,6 +500,42 @@ schedule GRHydro_SetupDescriptors AT CCTK_Initial BEFORE HydroBase_Initial
### the calculation of the source terms and the flux update. ###
#####################################################################
+
+schedule GRHydroTransformPrimToLocalBasis AT INITIAL AFTER (HydroBase_Initial, ADMBase_PostInitial) BEFORE HydroBase_Prim2ConInitial
+{
+ LANG: FORTRAN
+} "Transform primitive vars to local tensor basis."
+
+
+schedule GRHydroTransformADMToLocalBasis AT INITIAL AFTER HydroBase_Initial BEFORE GRHydroTransformPrimToLocalBasis
+{
+ LANG: FORTRAN
+} "Transform ADM metric, extr. curv. and shift to local tensor basis."
+
+
+schedule GRHydroTransformADMToLocalBasis IN ADMBase_SetADMVars
+{
+ LANG: FORTRAN
+} "Transform metric and shift to local tensor basis."
+
+#schedule GRHydroTransformADMToLocalBasis IN CTG_Convert_to_ADM AFTER CTGBase_Convert_CTG_to_ADM
+#{
+# LANG: FORTRAN
+#} "Transform metric and shift to local tensor basis."
+
+
+#schedule GRHydroTransformADMToLocalBasis IN MoL_Step BEFORE MoL_CalcRHS
+#{
+# LANG: FORTRAN
+#} "Transform metric and shift to local tensor basis."
+
+
+schedule GRHydroTransformPrimToGlobalBasis IN HydroBase_PostStep AFTER HydroBase_Con2Prim
+{
+ LANG: FORTRAN
+} "Transform primitive vars to global tensor basis."
+
+
schedule group GRHydroRHS IN HydroBase_RHS
{
STORAGE:GRHydro_scalars
@@ -942,6 +997,8 @@ schedule GRHydro_Boundaries IN HydroBase_Select_Boundaries AS GRHydro_Bound
SYNC: hydrobase::entropy
SYNC: hydrobase::Y_e
SYNC: Y_e_con
+ SYNC: lvel
+ SYNC: lBvec
} "Select GRHydro boundary conditions"
############################################################
diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90
index 5bfd457..0a1b3d8 100644
--- a/src/GRHydro_Boundaries.F90
+++ b/src/GRHydro_Boundaries.F90
@@ -15,15 +15,15 @@
#include "util_Table.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
-#define Bvecx(i,j,k) Bvec(i,j,k,1)
-#define Bvecy(i,j,k) Bvec(i,j,k,2)
-#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bvecx(i,j,k) lBvec(i,j,k,1)
+#define Bvecy(i,j,k) lBvec(i,j,k,2)
+#define Bvecz(i,j,k) lBvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)
@@ -94,7 +94,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
sym(2) = 1
sym(3) = 1
- call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]")
+ call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[0]")
call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]")
if(evolve_mhd.ne.0) then
call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]")
@@ -105,7 +105,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
sym(2) = -1
sym(3) = 1
- call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]")
+ call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[1]")
call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]")
if(evolve_mhd.ne.0) then
call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]")
@@ -116,7 +116,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS)
sym(2) = 1
sym(3) = -1
- call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]")
+ call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[2]")
call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]")
if(evolve_mhd.ne.0) then
call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]")
@@ -182,10 +182,10 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
"HydroBase::eps", "Flat")
ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
- "HydroBase::vel", "Flat")
+ "GRHydro::lvel", "Flat")
if(evolve_mhd.ne.0) then
ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
- "HydroBase::Bvec", "Flat")
+ "GRHydro::lBvec", "Flat")
ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
"GRHydro::Bcons", "Flat")
if(clean_divergence.ne.0) then
@@ -237,7 +237,7 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)
"HydroBase::vel", "None")
if(evolve_mhd.ne.0) then
ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
- "HydroBase::Bvec", "None")
+ "GRHydro::lBvec", "None")
ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, &
"GRHydro::Bcons", "None")
if(clean_divergence.ne.0) then
diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index d3c9732..3c51302 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -101,11 +101,11 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
epsnegative = .false.
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
- gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
- gyz(i,j,k),gzz(i,j,k))
+ gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),&
+ gbc(i,j,k),gcc(i,j,k))
!!$ if (det < 0.e0) then
!!$ call CCTK_WARN(0, "nan produced (det negative)")
@@ -135,7 +135,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
rho(i,j,k) = GRHydro_rho_min
scon(i,j,k,:) = 0.d0
- vel(i,j,k,:) = 0.d0
+ lvel(i,j,k,:) = 0.d0
w_lorentz(i,j,k) = 1.d0
if(evolve_temper.ne.0) then
@@ -148,7 +148,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
press(i,j,k),keyerr,anyerr)
else
- ! w_lorentz=1, so the expression for tau reduces to:
keytemp = 0
call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
@@ -157,6 +156,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
endif
+ ! w_lorentz=1, so the expression for tau reduces to:
tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
cycle
@@ -165,15 +165,15 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
if(evolve_temper.eq.0) then
call Con2Prim_pt(GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2), &
- scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), &
- vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+ scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), &
+ lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
else
call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),&
- scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vel(i,j,k,1),&
- vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+ scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),lvel(i,j,k,1),&
+ lvel(i,j,k,2),lvel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
press(i,j,k),w_lorentz(i,j,k), &
uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
@@ -186,7 +186,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
local_perc_ptol = GRHydro_perc_ptol*100.0d0
call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),&
scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),&
- vel(i,j,k,1),vel(i,j,k,2), vel(i,j,k,3),eps(i,j,k),&
+ lvel(i,j,k,1),lvel(i,j,k,2), lvel(i,j,k,3),eps(i,j,k),&
temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, &
@@ -237,7 +237,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
rho(i,j,k) = GRHydro_rho_min
scon(i,j,k,:) = 0.d0
- vel(i,j,k,:) = 0.d0
+ lvel(i,j,k,:) = 0.d0
w_lorentz(i,j,k) = 1.d0
call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
@@ -259,8 +259,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS)
!$OMP END CRITICAL
call Con2Prim_ptPolytype(GRHydro_polytrope_handle, &
dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
- tau(i,j,k), rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), &
- vel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), &
+ tau(i,j,k), rho(i,j,k), lvel(i,j,k,1), lvel(i,j,k,2), &
+ lvel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), &
uxx, uxy, uxz, uyy, uyz, uzz, det, x(i,j,k), y(i,j,k), &
z(i,j,k), r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
#endif
@@ -1048,18 +1048,18 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS)
if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+ gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
epsnegative = .false.
@@ -1251,11 +1251,11 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
- gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
- gyz(i,j,k),gzz(i,j,k))
+ gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),&
+ gbc(i,j,k),gcc(i,j,k))
if (evolve_tracer .ne. 0) then
do itracer=1,number_of_tracers
@@ -1277,8 +1277,8 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS)
call Con2Prim_ptPolytype(GRHydro_eos_handle, dens(i,j,k),&
scon(i,j,k,1),scon(i,j,k,2), &
- scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), &
- vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
+ scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), &
+ lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), &
z(i,j,k),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k))
@@ -1600,18 +1600,18 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS)
if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+ gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
@@ -1686,18 +1686,18 @@ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS)
do j = GRHydro_stencil, ny - GRHydro_stencil + 1
do i = GRHydro_stencil, nx - GRHydro_stencil + 1
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+ gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
@@ -1878,7 +1878,7 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS)
write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k)
call CCTK_WARN(1,warnline)
write(warnline,'(a20,3g16.7)') 'velocities: ',&
- vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3)
+ lvel(i,j,k,1),lvel(i,j,k,2),lvel(i,j,k,3)
call CCTK_WARN(0,warnline)
!$OMP END CRITICAL
diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90
index 0fcf71f..247b398 100644
--- a/src/GRHydro_ENOReconstruct_drv.F90
+++ b/src/GRHydro_ENOReconstruct_drv.F90
@@ -14,15 +14,15 @@
#include "SpaceMask.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
-#define Bvecx(i,j,k) Bvec(i,j,k,1)
-#define Bvecy(i,j,k) Bvec(i,j,k,2)
-#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bvecx(i,j,k) lBvec(i,j,k,1)
+#define Bvecy(i,j,k) lBvec(i,j,k,2)
+#define Bvecz(i,j,k) lBvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)
diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90
index 2c5776b..d1c022c 100644
--- a/src/GRHydro_EoSChangeGamma.F90
+++ b/src/GRHydro_EoSChangeGamma.F90
@@ -14,9 +14,9 @@
#include "cctk_Arguments.h"
#include "GRHydro_Macros.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
@@ -78,10 +78,10 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
- call prim2conpolytype(GRHydro_polytrope_handle,gxx(i,j,k),gxy(i,j,k),&
- gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
+ call prim2conpolytype(GRHydro_polytrope_handle,gaa(i,j,k),gab(i,j,k),&
+ gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
@@ -242,10 +242,10 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
- call prim2con(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
- gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
+ call prim2con(GRHydro_eos_handle,gaa(i,j,k),gab(i,j,k),&
+ gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90
index 2074b60..ee6c0ad 100644
--- a/src/GRHydro_HLLE.F90
+++ b/src/GRHydro_HLLE.F90
@@ -104,26 +104,26 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS)
!!$ left and right points.
if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
+ avg_beta = 0.5d0 * (betaa(i+xoffset,j+yoffset,k+zoffset) + &
+ betaa(i,j,k))
else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
+ avg_beta = 0.5d0 * (betab(i+xoffset,j+yoffset,k+zoffset) + &
+ betab(i,j,k))
else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
+ avg_beta = 0.5d0 * (betac(i+xoffset,j+yoffset,k+zoffset) + &
+ betac(i,j,k))
else
call CCTK_WARN(0, "Flux direction not x,y,z")
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))
+ gxxh = 0.5d0 * (gaa(i+xoffset,j+yoffset,k+zoffset) + gaa(i,j,k))
+ gxyh = 0.5d0 * (gab(i+xoffset,j+yoffset,k+zoffset) + gab(i,j,k))
+ gxzh = 0.5d0 * (gac(i+xoffset,j+yoffset,k+zoffset) + gac(i,j,k))
+ gyyh = 0.5d0 * (gbb(i+xoffset,j+yoffset,k+zoffset) + gbb(i,j,k))
+ gyzh = 0.5d0 * (gbc(i+xoffset,j+yoffset,k+zoffset) + gbc(i,j,k))
+ gzzh = 0.5d0 * (gcc(i+xoffset,j+yoffset,k+zoffset) + gcc(i,j,k))
avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\
gyyh,gyzh,gzzh)
@@ -467,26 +467,26 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS)
!!$ left and right points.
if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
+ avg_beta = 0.5d0 * (betaa(i+xoffset,j+yoffset,k+zoffset) + &
+ betaa(i,j,k))
else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
+ avg_beta = 0.5d0 * (betab(i+xoffset,j+yoffset,k+zoffset) + &
+ betab(i,j,k))
else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
+ avg_beta = 0.5d0 * (betac(i+xoffset,j+yoffset,k+zoffset) + &
+ betac(i,j,k))
else
call CCTK_WARN(0, "Flux direction not x,y,z")
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))
+ gxxh = 0.5d0 * (gaa(i+xoffset,j+yoffset,k+zoffset) + gaa(i,j,k))
+ gxyh = 0.5d0 * (gab(i+xoffset,j+yoffset,k+zoffset) + gab(i,j,k))
+ gxzh = 0.5d0 * (gac(i+xoffset,j+yoffset,k+zoffset) + gac(i,j,k))
+ gyyh = 0.5d0 * (gbb(i+xoffset,j+yoffset,k+zoffset) + gbb(i,j,k))
+ gyzh = 0.5d0 * (gbc(i+xoffset,j+yoffset,k+zoffset) + gbc(i,j,k))
+ gzzh = 0.5d0 * (gcc(i+xoffset,j+yoffset,k+zoffset) + gcc(i,j,k))
avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\
gyyh,gyzh,gzzh)
diff --git a/src/GRHydro_Jacobian_state.c b/src/GRHydro_Jacobian_state.c
new file mode 100644
index 0000000..4dd6f47
--- /dev/null
+++ b/src/GRHydro_Jacobian_state.c
@@ -0,0 +1,32 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+
+
+
+/**
+ Checks states for Jacobians
+*/
+void
+GRHydro_check_Jacobian_state (CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ int ierr;
+
+ if (*jacobian_state == 0)
+ {
+ CCTK_WARN(0, "GRHydro requires storage for Jacobians. Please tell your Coordinates implementation to store the Jacobians!");
+ }
+
+ if (*inverse_jacobian_state == 0)
+ {
+ CCTK_WARN(0, "GRHydro requires storage for inverse Jacobians. Please tell your Coordinates implementation to store the inverse Jacobians!");
+ }
+
+}
+
+
+
diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90
index 963dbf6..a8058cc 100644
--- a/src/GRHydro_PPMMReconstruct_drv.F90
+++ b/src/GRHydro_PPMMReconstruct_drv.F90
@@ -14,15 +14,15 @@
#include "SpaceMask.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
-#define Bvecx(i,j,k) Bvec(i,j,k,1)
-#define Bvecy(i,j,k) Bvec(i,j,k,2)
-#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bvecx(i,j,k) lBvec(i,j,k,1)
+#define Bvecy(i,j,k) lBvec(i,j,k,2)
+#define Bvecz(i,j,k) lBvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)
@@ -153,7 +153,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS)
Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),&
1,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
- gzz(:,j,k), betax(:,j,k), alp(:,j,k),&
+ gzz(:,j,k), betaa(:,j,k), alp(:,j,k),&
w_lorentz(:,j,k), &
1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
GRHydro_mppm_eigenvalue_x_right, &
@@ -168,7 +168,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS)
Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),&
0,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
- gzz(:,j,k), betax(:,j,k), alp(:,j,k),&
+ gzz(:,j,k), betaa(:,j,k), alp(:,j,k),&
w_lorentz(:,j,k), &
1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
GRHydro_mppm_eigenvalue_x_right, &
@@ -224,7 +224,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS)
Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),&
1,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
- gxx(j,:,k), betay(j,:,k), alp(j,:,k),&
+ gxx(j,:,k), betab(j,:,k), alp(j,:,k),&
w_lorentz(j,:,k), &
2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
GRHydro_mppm_eigenvalue_y_right, &
@@ -239,7 +239,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS)
Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),&
0,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
- gxx(j,:,k), betay(j,:,k), alp(j,:,k),&
+ gxx(j,:,k), betab(j,:,k), alp(j,:,k),&
w_lorentz(j,:,k), &
2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
GRHydro_mppm_eigenvalue_y_right, &
@@ -294,7 +294,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS)
Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),&
1,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
- gyy(j,k,:), betaz(j,k,:), alp(j,k,:),&
+ gyy(j,k,:), betac(j,k,:), alp(j,k,:),&
w_lorentz(j,k,:), &
3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
GRHydro_mppm_eigenvalue_z_right, &
@@ -309,7 +309,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS)
Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),&
0,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
- gyy(j,k,:), betaz(j,k,:), alp(j,k,:),&
+ gyy(j,k,:), betac(j,k,:), alp(j,k,:),&
w_lorentz(j,k,:), &
3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
GRHydro_mppm_eigenvalue_z_right, &
diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90
index b514247..2ecbb71 100644
--- a/src/GRHydro_PPMReconstruct_drv.F90
+++ b/src/GRHydro_PPMReconstruct_drv.F90
@@ -14,9 +14,9 @@
#include "SpaceMask.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
@@ -125,8 +125,8 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),&
velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),&
trivial_rp(:,j,k), hydro_excision_mask(:,j,k),&
- gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), &
- gzz(:,j,k), betax(:,j,k), alp(:,j,k),&
+ gaa(:,j,k), gab(:,j,k), gac(:,j,k), gbb(:,j,k), gbc(:,j,k), &
+ gcc(:,j,k), betaa(:,j,k), alp(:,j,k),&
w_lorentz(:,j,k), &
1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, &
GRHydro_mppm_eigenvalue_x_right, &
@@ -176,8 +176,8 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),&
velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),&
trivial_rp(j,:,k), hydro_excision_mask(j,:,k),&
- gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), &
- gxx(j,:,k), betay(j,:,k), alp(j,:,k),&
+ gbb(j,:,k), gbc(j,:,k), gab(j,:,k), gcc(j,:,k), gac(j,:,k), &
+ gaa(j,:,k), betab(j,:,k), alp(j,:,k),&
w_lorentz(j,:,k), &
2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, &
GRHydro_mppm_eigenvalue_y_right, &
@@ -227,8 +227,8 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS)
velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),&
velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),&
trivial_rp(j,k,:), hydro_excision_mask(j,k,:),&
- gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), &
- gyy(j,k,:), betaz(j,k,:), alp(j,k,:),&
+ gcc(j,k,:), gac(j,k,:), gbc(j,k,:), gaa(j,k,:), gab(j,k,:), &
+ gbb(j,k,:), betac(j,k,:), alp(j,k,:),&
w_lorentz(j,k,:), &
3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, &
GRHydro_mppm_eigenvalue_z_right, &
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90
index fbe118b..2bd2d29 100644
--- a/src/GRHydro_Prim2Con.F90
+++ b/src/GRHydro_Prim2Con.F90
@@ -13,9 +13,9 @@
#include "cctk_Functions.h"
#include "GRHydro_Macros.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
@@ -44,44 +44,44 @@ subroutine primitive2conservative(CCTK_ARGUMENTS)
DECLARE_CCTK_PARAMETERS
integer :: i, j, k
- CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
- gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+ CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,&
+ gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr
CCTK_REAL :: xtemp(1)
if(evolve_temper.ne.1) then
!$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr,&
- !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
- !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
+ !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, &
+ !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr)
do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+ gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
+ avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl)
+ avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr)
- call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,&
- gyzl,gzzl, &
+ call prim2con(GRHydro_eos_handle, gaal,gabl,gacl,gbbl,&
+ gbcl,gccl, &
avg_detl,densminus(i,j,k),sxminus(i,j,k),&
syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),&
rhominus(i,j,k), &
velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
- call prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,&
- gyyr,gyzr,gzzr, &
+ call prim2con(GRHydro_eos_handle, gaar,gabr,gacr,&
+ gbbr,gbcr,gccr, &
avg_detr, densplus(i,j,k),sxplus(i,j,k),&
syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
@@ -94,27 +94,27 @@ subroutine primitive2conservative(CCTK_ARGUMENTS)
!$OMP END PARALLEL DO
else
!$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,&
- !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
- !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
+ !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, &
+ !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr)
do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+ gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
+ avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl)
+ avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr)
! we do not have plus/minus vars for temperature since
! eps is reconstructed. Hence, we do not update the
@@ -124,8 +124,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS)
xtemp(1) = 0.5d0*(temperature(i,j,k) + &
temperature(i-xoffset,j-yoffset,k-zoffset))
call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
- i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxl,gxyl,gxzl,gyyl,&
- gyzl,gzzl, &
+ i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaal,gabl,gacl,gbbl,&
+ gbcl,gccl, &
avg_detl,densminus(i,j,k),sxminus(i,j,k),&
syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),&
rhominus(i,j,k), &
@@ -141,8 +141,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS)
xtemp(1) = 0.5d0*(temperature(i,j,k) + &
temperature(i+xoffset,j+yoffset,k+zoffset))
call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, &
- i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxr,gxyr,gxzr,&
- gyyr,gyzr,gzzr, &
+ i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaar,gabr,gacr,&
+ gbbr,gbcr,gccr, &
avg_detr, densplus(i,j,k),sxplus(i,j,k),&
syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
@@ -353,12 +353,12 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS)
do j = 1,cctk_lsh(2)
do i = 1,cctk_lsh(1)
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
- call prim2con(GRHydro_eos_handle,gxx(i,j,k),&
- gxy(i,j,k),gxz(i,j,k),&
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ call prim2con(GRHydro_eos_handle,gaa(i,j,k),&
+ gab(i,j,k),gac(i,j,k),&
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
@@ -372,13 +372,13 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS)
do j = 1,cctk_lsh(2)
do i = 1,cctk_lsh(1)
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,i,j,k,&
- x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),&
- gxy(i,j,k),gxz(i,j,k),&
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),&
+ gab(i,j,k),gac(i,j,k),&
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), &
@@ -429,40 +429,40 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS)
DECLARE_CCTK_PARAMETERS
integer :: i, j, k
- CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
- gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+ CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,&
+ gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr
- !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
- !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr)
+ !$OMP PARALLEL DO PRIVATE(i, j, k, gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,&
+ !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr)
do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
-
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
-
- call prim2conpolytype(GRHydro_eos_handle, gxxl,gxyl,gxzl,&
- gyyl,gyzl,gzzl, &
+ gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
+
+ avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl)
+ avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr)
+
+ call prim2conpolytype(GRHydro_eos_handle, gaal,gabl,gacl,&
+ gbbl,gbcl,gccl, &
avg_detl,densminus(i,j,k),sxminus(i,j,k),&
syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), &
velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k))
- call prim2conpolytype(GRHydro_eos_handle, gxxr,gxyr,gxzr,&
- gyyr,gyzr,gzzr, &
+ call prim2conpolytype(GRHydro_eos_handle, gaar,gabr,gacr,&
+ gbbr,gbcr,gccr, &
avg_detr, densplus(i,j,k),sxplus(i,j,k),&
syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
@@ -586,11 +586,11 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS)
do j = 1,cctk_lsh(2)
do i = 1,cctk_lsh(1)
- det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\
- gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\
+ gbb(i,j,k),gbc(i,j,k),gcc(i,j,k))
- call prim2conpolytype(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
- gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+ call prim2conpolytype(GRHydro_eos_handle,gaa(i,j,k),gab(i,j,k),&
+ gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
@@ -637,46 +637,46 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS)
DECLARE_CCTK_PARAMETERS
integer :: i, j, k
- CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
- gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+ CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,&
+ gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr
do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
-
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
+ gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset))
+ gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset))
+ gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset))
+ gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset))
+ gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset))
+ gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset))
+ gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset))
+ gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset))
+ gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset))
+ gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset))
+ gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset))
+ gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset))
+
+ avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl)
+ avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr)
cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * &
sqrt(avg_detr) * rhoplus(i,j,k) / &
sqrt(1.d0 - &
- (gxxr * velxplus(i,j,k)**2 + &
- gyyr * velyplus(i,j,k)**2 + &
- gzzr * velzplus(i,j,k)**2 + &
- 2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + &
- gxzr * velxplus(i,j,k) * velzplus(i,j,k) + &
- gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) )
+ (gaar * velxplus(i,j,k)**2 + &
+ gbbr * velyplus(i,j,k)**2 + &
+ gccr * velzplus(i,j,k)**2 + &
+ 2.d0 * (gabr * velxplus(i,j,k) * velyplus(i,j,k) + &
+ gacr * velxplus(i,j,k) * velzplus(i,j,k) + &
+ gbcr * velyplus(i,j,k) * velzplus(i,j,k) ) ) )
cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * &
sqrt(avg_detl) * rhominus(i,j,k) / &
sqrt(1.d0 - &
- (gxxl * velxminus(i,j,k)**2 + &
- gyyl * velyminus(i,j,k)**2 + &
- gzzl * velzminus(i,j,k)**2 + &
- 2.d0 * (gxyl * velxminus(i,j,k) * velyminus(i,j,k) + &
- gxzl * velxminus(i,j,k) * velzminus(i,j,k) + &
- gyzl * velyminus(i,j,k) * velzminus(i,j,k) ) ) )
+ (gaal * velxminus(i,j,k)**2 + &
+ gbbl * velyminus(i,j,k)**2 + &
+ gccl * velzminus(i,j,k)**2 + &
+ 2.d0 * (gabl * velxminus(i,j,k) * velyminus(i,j,k) + &
+ gacl * velxminus(i,j,k) * velzminus(i,j,k) + &
+ gbcl * velyminus(i,j,k) * velzminus(i,j,k) ) ) )
end do
end do
diff --git a/src/GRHydro_RegisterVars.cc b/src/GRHydro_RegisterVars.cc
index 0d99955..5e91721 100644
--- a/src/GRHydro_RegisterVars.cc
+++ b/src/GRHydro_RegisterVars.cc
@@ -53,8 +53,13 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS)
register_constrained("HydroBase::rho");
register_constrained("HydroBase::press");
register_constrained("HydroBase::eps");
- register_constrained("HydroBase::vel");
register_constrained("HydroBase::w_lorentz");
+ register_constrained("HydroBase::vel");
+ register_constrained("GRHydro::lvel");
+
+ register_constrained("grhydro::local_shift");
+ register_constrained("grhydro::local_metric");
+ register_constrained("grhydro::local_extrinsic_curvature");
if (CCTK_EQUALS(evolution_method, "GRHydro"))
{
@@ -64,6 +69,7 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS)
if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) {
register_constrained("HydroBase::Bvec");
+ register_constrained("GRHydro::lBvec");
register_evolved("GRHydro::Bcons", "GRHydro::Bconsrhs");
if(clean_divergence) {
register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs");
diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90
index f8fd9eb..969e88f 100644
--- a/src/GRHydro_Source.F90
+++ b/src/GRHydro_Source.F90
@@ -9,27 +9,27 @@
! Second order f.d.
-#define DIFF_X_2(q) 0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx
-#define DIFF_Y_2(q) 0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy
-#define DIFF_Z_2(q) 0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz
+#define DIFF_X_2(q) 0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * ida
+#define DIFF_Y_2(q) 0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idb
+#define DIFF_Z_2(q) 0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idc
! Fourth order f.d.
#define DIFF_X_4(q) (-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \
- q(i-2,j,k)) / 12.d0 * idx
+ q(i-2,j,k)) / 12.d0 * ida
#define DIFF_Y_4(q) (-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \
- q(i,j-2,k)) / 12.d0 * idy
+ q(i,j-2,k)) / 12.d0 * idb
#define DIFF_Z_4(q) (-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \
- q(i,j,k-2)) / 12.d0 * idz
+ q(i,j,k-2)) / 12.d0 * idc
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "GRHydro_Macros.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define vela(i,j,k) lvel(i,j,k,1)
+#define velb(i,j,k) lvel(i,j,k,2)
+#define velc(i,j,k) lvel(i,j,k,3)
/*@@
@routine SourceTerms
@date Sat Jan 26 02:04:21 2002
@@ -52,26 +52,26 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
- CCTK_INT :: i, j, k, nx, ny, nz
+ CCTK_INT :: i, j, k, na, nb, nc
CCTK_REAL :: one, two, half
- CCTK_REAL :: t00, t0x, t0y, t0z, txx, txy, txz, tyy, tyz, tzz
- CCTK_REAL :: sqrtdet, det, uxx, uxy, uxz, uyy, uyz, uzz, rhoenthalpyW2
- CCTK_REAL :: shiftx, shifty, shiftz, velxshift, velyshift, velzshift
- CCTK_REAL :: vlowx, vlowy, vlowz
- CCTK_REAL :: dx_betax, dx_betay, dx_betaz, dy_betax, dy_betay,&
- dy_betaz, dz_betax, dz_betay, dz_betaz
- CCTK_REAL :: dx_alp, dy_alp, dz_alp
- CCTK_REAL :: tau_source, sx_source, sy_source, sz_source
- CCTK_REAL :: localgxx,localgxy,localgxz,localgyy,localgyz,localgzz
- CCTK_REAL :: dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz
- CCTK_REAL :: dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz
- CCTK_REAL :: dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz
- CCTK_REAL :: dx, dy, dz, idx, idy, idz
- CCTK_REAL :: psi4pt, dx_psi4, dy_psi4, dz_psi4
- CCTK_REAL :: shiftshiftk, shiftkx, shiftky, shiftkz
+ CCTK_REAL :: t00, t0a, t0b, t0c, taa, tab, tac, tbb, tbc, tcc
+ CCTK_REAL :: sqrtdet, det, uaa, uab, uac, ubb, ubc, ucc, rhoenthalpyW2
+ CCTK_REAL :: shifta, shiftb, shiftc, velashift, velbshift, velcshift
+ CCTK_REAL :: vlowa, vlowb, vlowc
+ CCTK_REAL :: da_betaa, da_betab, da_betac, db_betaa, db_betab,&
+ db_betac, dc_betaa, dc_betab, dc_betac
+ CCTK_REAL :: da_alp, db_alp, dc_alp
+ CCTK_REAL :: tau_source, sa_source, sb_source, sc_source
+ CCTK_REAL :: localgaa,localgab,localgac,localgbb,localgbc,localgcc
+ CCTK_REAL :: da_gaa, da_gab, da_gac, da_gbb, da_gbc, da_gcc
+ CCTK_REAL :: db_gaa, db_gab, db_gac, db_gbb, db_gbc, db_gcc
+ CCTK_REAL :: dc_gaa, dc_gab, dc_gac, dc_gbb, dc_gbc, dc_gcc
+ CCTK_REAL :: da, db, dc, ida, idb, idc
+ CCTK_REAL :: psi4pt, da_psi4, db_psi4, dc_psi4
+ CCTK_REAL :: shiftshiftk, shiftka, shiftkb, shiftkc
CCTK_REAL :: sumTK
- CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz
- CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz
+ CCTK_REAL :: halfshiftdga, halfshiftdgb, halfshiftdgc
+ CCTK_REAL :: halfTdga, halfTdgb, halfTdgc
CCTK_REAL :: invalp, invalp2
logical, allocatable, dimension (:,:,:) :: force_spatial_second_order
@@ -79,15 +79,15 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
one = 1.0d0
two = 2.0d0
half = 0.5d0
- nx = cctk_lsh(1)
- ny = cctk_lsh(2)
- nz = cctk_lsh(3)
- dx = CCTK_DELTA_SPACE(1)
- dy = CCTK_DELTA_SPACE(2)
- dz = CCTK_DELTA_SPACE(3)
- idx = 1.d0/dx
- idy = 1.d0/dy
- idz = 1.d0/dz
+ na = cctk_lsh(1)
+ nb = cctk_lsh(2)
+ nc = cctk_lsh(3)
+ da = CCTK_DELTA_SPACE(1)
+ db = CCTK_DELTA_SPACE(2)
+ dc = CCTK_DELTA_SPACE(3)
+ ida = 1.d0/da
+ idb = 1.d0/db
+ idc = 1.d0/dc
!!$ Initialize the update terms to be zero.
!!$ This will guarantee that no garbage in the boundaries is updated.
@@ -108,14 +108,14 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
!!$ differencing at boundaries and near excision regions.
!!$ Copied straight from BSSN.
- allocate (force_spatial_second_order(nx,ny,nz))
+ allocate (force_spatial_second_order(na,nb,nc))
force_spatial_second_order = .FALSE.
if (spatial_order > 2) then
!$OMP PARALLEL DO PRIVATE(i, j, k)
- do k = 1 + GRHydro_stencil, nz - GRHydro_stencil
- do j = 1 + GRHydro_stencil, ny - GRHydro_stencil
- do i = 1 + GRHydro_stencil, nx - GRHydro_stencil
+ do k = 1 + GRHydro_stencil, nc - GRHydro_stencil
+ do j = 1 + GRHydro_stencil, nb - GRHydro_stencil
+ do i = 1 + GRHydro_stencil, na - GRHydro_stencil
if ((i < 3).or.(i > cctk_lsh(1) - 2).or. &
(j < 3).or.(j > cctk_lsh(2) - 2).or. &
(k < 3).or.(k > cctk_lsh(3) - 2) ) then
@@ -132,22 +132,22 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
end if
!$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,&
- !$OMP localgxx,localgxy,localgxz,localgyy,localgyz,localgzz,&
- !$OMP psi4pt,det,sqrtdet,rhoenthalpyW2,shiftx,shifty,shiftz,&
- !$OMP dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz,&
- !$OMP dz_betax,dz_betay,dz_betaz,velxshift,velyshift,velzshift,&
- !$OMP vlowx,vlowy,vlowz,t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz,&
- !$OMP dx_alp,dy_alp,dz_alp,tau_source,sx_source,sy_source,sz_source,&
- !$OMP uxx, uxy, uxz, uyy, uyz, uzz,&
- !$OMP dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz,&
- !$OMP dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz,&
- !$OMP dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz,&
- !$OMP dx_psi4,dy_psi4,dz_psi4,shiftshiftk,shiftkx,shiftky,shiftkz,&
- !$OMP sumTK,halfshiftdgx,halfshiftdgy,halfshiftdgz,&
- !$OMP halfTdgx,halfTdgy,halfTdgz,invalp,invalp2)
- do k=1 + GRHydro_stencil,nz - GRHydro_stencil
- do j=1 + GRHydro_stencil,ny - GRHydro_stencil
- do i=1 + GRHydro_stencil,nx - GRHydro_stencil
+ !$OMP localgaa,localgab,localgac,localgbb,localgbc,localgcc,&
+ !$OMP psi4pt,det,sqrtdet,rhoenthalpyW2,shifta,shiftb,shiftc,&
+ !$OMP da_betaa,da_betab,da_betac,db_betaa,db_betab,db_betac,&
+ !$OMP dc_betaa,dc_betab,dc_betac,velashift,velbshift,velcshift,&
+ !$OMP vlowa,vlowb,vlowc,t00,t0a,t0b,t0c,taa,tab,tac,tbb,tbc,tcc,&
+ !$OMP da_alp,db_alp,dc_alp,tau_source,sa_source,sb_source,sc_source,&
+ !$OMP uaa, uab, uac, ubb, ubc, ucc,&
+ !$OMP da_gaa, da_gab, da_gac, da_gbb, da_gbc, da_gcc,&
+ !$OMP db_gaa, db_gab, db_gac, db_gbb, db_gbc, db_gcc,&
+ !$OMP dc_gaa, dc_gab, dc_gac, dc_gbb, dc_gbc, dc_gcc,&
+ !$OMP da_psi4,db_psi4,dc_psi4,shiftshiftk,shiftka,shiftkb,shiftkc,&
+ !$OMP sumTK,halfshiftdga,halfshiftdgb,halfshiftdgc,&
+ !$OMP halfTdga,halfTdgb,halfTdgc,invalp,invalp2)
+ do k=1 + GRHydro_stencil,nc - GRHydro_stencil
+ do j=1 + GRHydro_stencil,nb - GRHydro_stencil
+ do i=1 + GRHydro_stencil,na - GRHydro_stencil
local_spatial_order = spatial_order
if (force_spatial_second_order(i,j,k)) then
@@ -156,18 +156,18 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
!!$ Set the metric terms.
- localgxx = gxx(i,j,k)
- localgxy = gxy(i,j,k)
- localgxz = gxz(i,j,k)
- localgyy = gyy(i,j,k)
- localgyz = gyz(i,j,k)
- localgzz = gzz(i,j,k)
+ localgaa = gaa(i,j,k)
+ localgab = gab(i,j,k)
+ localgac = gac(i,j,k)
+ localgbb = gbb(i,j,k)
+ localgbc = gbc(i,j,k)
+ localgcc = gcc(i,j,k)
- det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\
- localgyy, localgyz, localgzz)
+ det = SPATIAL_DETERMINANT(localgaa, localgab, localgac,\
+ localgbb, localgbc, localgcc)
sqrtdet = sqrt(det)
- call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,&
- localgxy, localgxz, localgyy, localgyz, localgzz)
+ call UpperMetric(uaa, uab, uac, ubb, ubc, ucc, det, localgaa,&
+ localgab, localgac, localgbb, localgbc, localgcc)
!!$ All the matter variables (except velocity) always appear
!!$ together in this form
@@ -175,218 +175,218 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
rhoenthalpyW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k))*&
w_lorentz(i,j,k)**2
- shiftx = betax(i,j,k)
- shifty = betay(i,j,k)
- shiftz = betaz(i,j,k)
+ shifta = betaa(i,j,k)
+ shiftb = betab(i,j,k)
+ shiftc = betac(i,j,k)
if (local_spatial_order .eq. 2) then
- dx_betax = DIFF_X_2(betax)
- dx_betay = DIFF_X_2(betay)
- dx_betaz = DIFF_X_2(betaz)
+ da_betaa = DIFF_X_2(betaa)
+ da_betab = DIFF_X_2(betab)
+ da_betac = DIFF_X_2(betac)
- dy_betax = DIFF_Y_2(betax)
- dy_betay = DIFF_Y_2(betay)
- dy_betaz = DIFF_Y_2(betaz)
+ db_betaa = DIFF_Y_2(betaa)
+ db_betab = DIFF_Y_2(betab)
+ db_betac = DIFF_Y_2(betac)
- dz_betax = DIFF_Z_2(betax)
- dz_betay = DIFF_Z_2(betay)
- dz_betaz = DIFF_Z_2(betaz)
+ dc_betaa = DIFF_Z_2(betaa)
+ dc_betab = DIFF_Z_2(betab)
+ dc_betac = DIFF_Z_2(betac)
else
- dx_betax = DIFF_X_4(betax)
- dx_betay = DIFF_X_4(betay)
- dx_betaz = DIFF_X_4(betaz)
+ da_betaa = DIFF_X_4(betaa)
+ da_betab = DIFF_X_4(betab)
+ da_betac = DIFF_X_4(betac)
- dy_betax = DIFF_Y_4(betax)
- dy_betay = DIFF_Y_4(betay)
- dy_betaz = DIFF_Y_4(betaz)
+ db_betaa = DIFF_Y_4(betaa)
+ db_betab = DIFF_Y_4(betab)
+ db_betac = DIFF_Y_4(betac)
- dz_betax = DIFF_Z_4(betax)
- dz_betay = DIFF_Z_4(betay)
- dz_betaz = DIFF_Z_4(betaz)
+ dc_betaa = DIFF_Z_4(betaa)
+ dc_betab = DIFF_Z_4(betab)
+ dc_betac = DIFF_Z_4(betac)
end if
invalp = 1.0d0 / alp(i,j,k)
invalp2 = invalp**2
- velxshift = velx(i,j,k) - shiftx*invalp
- velyshift = vely(i,j,k) - shifty*invalp
- velzshift = velz(i,j,k) - shiftz*invalp
- vlowx = velx(i,j,k)*localgxx + vely(i,j,k)*localgxy +&
- velz(i,j,k)*localgxz
- vlowy = velx(i,j,k)*localgxy + vely(i,j,k)*localgyy +&
- velz(i,j,k)*localgyz
- vlowz = velx(i,j,k)*localgxz + vely(i,j,k)*localgyz +&
- velz(i,j,k)*localgzz
+ velashift = vela(i,j,k) - shifta*invalp
+ velbshift = velb(i,j,k) - shiftb*invalp
+ velcshift = velc(i,j,k) - shiftc*invalp
+ vlowa = vela(i,j,k)*localgaa + velb(i,j,k)*localgab +&
+ velc(i,j,k)*localgac
+ vlowb = vela(i,j,k)*localgab + velb(i,j,k)*localgbb +&
+ velc(i,j,k)*localgbc
+ vlowc = vela(i,j,k)*localgac + velb(i,j,k)*localgbc +&
+ velc(i,j,k)*localgcc
!!$ For a change, these are T^{ij}
t00 = (rhoenthalpyW2 - press(i,j,k))*invalp2
- t0x = rhoenthalpyW2*velxshift/alp(i,j,k) +&
- press(i,j,k)*shiftx*invalp2
- t0y = rhoenthalpyW2*velyshift/alp(i,j,k) +&
- press(i,j,k)*shifty*invalp2
- t0z = rhoenthalpyW2*velzshift/alp(i,j,k) +&
- press(i,j,k)*shiftz*invalp2
- txx = rhoenthalpyW2*velxshift*velxshift +&
- press(i,j,k)*(uxx - shiftx*shiftx*invalp2)
- txy = rhoenthalpyW2*velxshift*velyshift +&
- press(i,j,k)*(uxy - shiftx*shifty*invalp2)
- txz = rhoenthalpyW2*velxshift*velzshift +&
- press(i,j,k)*(uxz - shiftx*shiftz*invalp2)
- tyy = rhoenthalpyW2*velyshift*velyshift +&
- press(i,j,k)*(uyy - shifty*shifty*invalp2)
- tyz = rhoenthalpyW2*velyshift*velzshift +&
- press(i,j,k)*(uyz - shifty*shiftz*invalp2)
- tzz = rhoenthalpyW2*velzshift*velzshift +&
- press(i,j,k)*(uzz - shiftz*shiftz*invalp2)
+ t0a = rhoenthalpyW2*velashift/alp(i,j,k) +&
+ press(i,j,k)*shifta*invalp2
+ t0b = rhoenthalpyW2*velbshift/alp(i,j,k) +&
+ press(i,j,k)*shiftb*invalp2
+ t0c = rhoenthalpyW2*velcshift/alp(i,j,k) +&
+ press(i,j,k)*shiftc*invalp2
+ taa = rhoenthalpyW2*velashift*velashift +&
+ press(i,j,k)*(uaa - shifta*shifta*invalp2)
+ tab = rhoenthalpyW2*velashift*velbshift +&
+ press(i,j,k)*(uab - shifta*shiftb*invalp2)
+ tac = rhoenthalpyW2*velashift*velcshift +&
+ press(i,j,k)*(uac - shifta*shiftc*invalp2)
+ tbb = rhoenthalpyW2*velbshift*velbshift +&
+ press(i,j,k)*(ubb - shiftb*shiftb*invalp2)
+ tbc = rhoenthalpyW2*velbshift*velcshift +&
+ press(i,j,k)*(ubc - shiftb*shiftc*invalp2)
+ tcc = rhoenthalpyW2*velcshift*velcshift +&
+ press(i,j,k)*(ucc - shiftc*shiftc*invalp2)
!!$ Derivatives of the lapse, metric and shift
if (local_spatial_order .eq. 2) then
- dx_alp = DIFF_X_2(alp)
- dy_alp = DIFF_Y_2(alp)
- dz_alp = DIFF_Z_2(alp)
+ da_alp = DIFF_X_2(alp)
+ db_alp = DIFF_Y_2(alp)
+ dc_alp = DIFF_Z_2(alp)
else
- dx_alp = DIFF_X_4(alp)
- dy_alp = DIFF_Y_4(alp)
- dz_alp = DIFF_Z_4(alp)
+ da_alp = DIFF_X_4(alp)
+ db_alp = DIFF_Y_4(alp)
+ dc_alp = DIFF_Z_4(alp)
end if
if (local_spatial_order .eq. 2) then
- dx_gxx = DIFF_X_2(gxx)
- dx_gxy = DIFF_X_2(gxy)
- dx_gxz = DIFF_X_2(gxz)
- dx_gyy = DIFF_X_2(gyy)
- dx_gyz = DIFF_X_2(gyz)
- dx_gzz = DIFF_X_2(gzz)
- dy_gxx = DIFF_Y_2(gxx)
- dy_gxy = DIFF_Y_2(gxy)
- dy_gxz = DIFF_Y_2(gxz)
- dy_gyy = DIFF_Y_2(gyy)
- dy_gyz = DIFF_Y_2(gyz)
- dy_gzz = DIFF_Y_2(gzz)
- dz_gxx = DIFF_Z_2(gxx)
- dz_gxy = DIFF_Z_2(gxy)
- dz_gxz = DIFF_Z_2(gxz)
- dz_gyy = DIFF_Z_2(gyy)
- dz_gyz = DIFF_Z_2(gyz)
- dz_gzz = DIFF_Z_2(gzz)
+ da_gaa = DIFF_X_2(gaa)
+ da_gab = DIFF_X_2(gab)
+ da_gac = DIFF_X_2(gac)
+ da_gbb = DIFF_X_2(gbb)
+ da_gbc = DIFF_X_2(gbc)
+ da_gcc = DIFF_X_2(gcc)
+ db_gaa = DIFF_Y_2(gaa)
+ db_gab = DIFF_Y_2(gab)
+ db_gac = DIFF_Y_2(gac)
+ db_gbb = DIFF_Y_2(gbb)
+ db_gbc = DIFF_Y_2(gbc)
+ db_gcc = DIFF_Y_2(gcc)
+ dc_gaa = DIFF_Z_2(gaa)
+ dc_gab = DIFF_Z_2(gab)
+ dc_gac = DIFF_Z_2(gac)
+ dc_gbb = DIFF_Z_2(gbb)
+ dc_gbc = DIFF_Z_2(gbc)
+ dc_gcc = DIFF_Z_2(gcc)
else
- dx_gxx = DIFF_X_4(gxx)
- dx_gxy = DIFF_X_4(gxy)
- dx_gxz = DIFF_X_4(gxz)
- dx_gyy = DIFF_X_4(gyy)
- dx_gyz = DIFF_X_4(gyz)
- dx_gzz = DIFF_X_4(gzz)
- dy_gxx = DIFF_Y_4(gxx)
- dy_gxy = DIFF_Y_4(gxy)
- dy_gxz = DIFF_Y_4(gxz)
- dy_gyy = DIFF_Y_4(gyy)
- dy_gyz = DIFF_Y_4(gyz)
- dy_gzz = DIFF_Y_4(gzz)
- dz_gxx = DIFF_Z_4(gxx)
- dz_gxy = DIFF_Z_4(gxy)
- dz_gxz = DIFF_Z_4(gxz)
- dz_gyy = DIFF_Z_4(gyy)
- dz_gyz = DIFF_Z_4(gyz)
- dz_gzz = DIFF_Z_4(gzz)
+ da_gaa = DIFF_X_4(gaa)
+ da_gab = DIFF_X_4(gab)
+ da_gac = DIFF_X_4(gac)
+ da_gbb = DIFF_X_4(gbb)
+ da_gbc = DIFF_X_4(gbc)
+ da_gcc = DIFF_X_4(gcc)
+ db_gaa = DIFF_Y_4(gaa)
+ db_gab = DIFF_Y_4(gab)
+ db_gac = DIFF_Y_4(gac)
+ db_gbb = DIFF_Y_4(gbb)
+ db_gbc = DIFF_Y_4(gbc)
+ db_gcc = DIFF_Y_4(gcc)
+ dc_gaa = DIFF_Z_4(gaa)
+ dc_gab = DIFF_Z_4(gab)
+ dc_gac = DIFF_Z_4(gac)
+ dc_gbb = DIFF_Z_4(gbb)
+ dc_gbc = DIFF_Z_4(gbc)
+ dc_gcc = DIFF_Z_4(gcc)
end if
-!!$ Contract the shift with the extrinsic curvature
+!!$ Contract the shift with the eatrinsic curvature
- shiftshiftk = shiftx*shiftx*kxx(i,j,k) + &
- shifty*shifty*kyy(i,j,k) + &
- shiftz*shiftz*kzz(i,j,k) + &
- two*(shiftx*shifty*kxy(i,j,k) + &
- shiftx*shiftz*kxz(i,j,k) + &
- shifty*shiftz*kyz(i,j,k))
+ shiftshiftk = shifta*shifta*kaa(i,j,k) + &
+ shiftb*shiftb*kbb(i,j,k) + &
+ shiftc*shiftc*kcc(i,j,k) + &
+ two*(shifta*shiftb*kab(i,j,k) + &
+ shifta*shiftc*kac(i,j,k) + &
+ shiftb*shiftc*kbc(i,j,k))
- shiftkx = shiftx*kxx(i,j,k) + shifty*kxy(i,j,k) + shiftz*kxz(i,j,k)
- shiftky = shiftx*kxy(i,j,k) + shifty*kyy(i,j,k) + shiftz*kyz(i,j,k)
- shiftkz = shiftx*kxz(i,j,k) + shifty*kyz(i,j,k) + shiftz*kzz(i,j,k)
+ shiftka = shifta*kaa(i,j,k) + shiftb*kab(i,j,k) + shiftc*kac(i,j,k)
+ shiftkb = shifta*kab(i,j,k) + shiftb*kbb(i,j,k) + shiftc*kbc(i,j,k)
+ shiftkc = shifta*kac(i,j,k) + shiftb*kbc(i,j,k) + shiftc*kcc(i,j,k)
!!$ Contract the matter terms with the extrinsic curvature
- sumTK = txx*kxx(i,j,k) + tyy*kyy(i,j,k) + tzz*kzz(i,j,k) &
- + two*(txy*kxy(i,j,k) + txz*kxz(i,j,k) + tyz*kyz(i,j,k))
+ sumTK = taa*kaa(i,j,k) + tbb*kbb(i,j,k) + tcc*kcc(i,j,k) &
+ + two*(tab*kab(i,j,k) + tac*kac(i,j,k) + tbc*kbc(i,j,k))
!!$ Update term for tau
tau_source = t00* &
- (shiftshiftk - (shiftx*dx_alp + shifty*dy_alp + shiftz*dz_alp) )&
- + t0x*(-dx_alp + two*shiftkx) &
- + t0y*(-dy_alp + two*shiftky) &
- + t0z*(-dz_alp + two*shiftkz) &
+ (shiftshiftk - (shifta*da_alp + shiftb*db_alp + shiftc*dc_alp) )&
+ + t0a*(-da_alp + two*shiftka) &
+ + t0b*(-db_alp + two*shiftkb) &
+ + t0c*(-dc_alp + two*shiftkc) &
+ sumTK
-!!$ The following looks very little like the terms in the
+!!$ The following looks verb little like the terms in the
!!$ standard papers. Take a look in the ThornGuide to see why
!!$ it is really the same thing.
!!$ Contract the shift with derivatives of the metric
- halfshiftdgx = half*(shiftx*shiftx*dx_gxx + &
- shifty*shifty*dx_gyy + shiftz*shiftz*dx_gzz) + &
- shiftx*shifty*dx_gxy + shiftx*shiftz*dx_gxz + &
- shifty*shiftz*dx_gyz
- halfshiftdgy = half*(shiftx*shiftx*dy_gxx + &
- shifty*shifty*dy_gyy + shiftz*shiftz*dy_gzz) + &
- shiftx*shifty*dy_gxy + shiftx*shiftz*dy_gxz + &
- shifty*shiftz*dy_gyz
- halfshiftdgz = half*(shiftx*shiftx*dz_gxx + &
- shifty*shifty*dz_gyy + shiftz*shiftz*dz_gzz) + &
- shiftx*shifty*dz_gxy + shiftx*shiftz*dz_gxz + &
- shifty*shiftz*dz_gyz
+ halfshiftdga = half*(shifta*shifta*da_gaa + &
+ shiftb*shiftb*da_gbb + shiftc*shiftc*da_gcc) + &
+ shifta*shiftb*da_gab + shifta*shiftc*da_gac + &
+ shiftb*shiftc*da_gbc
+ halfshiftdgb = half*(shifta*shifta*db_gaa + &
+ shiftb*shiftb*db_gbb + shiftc*shiftc*db_gcc) + &
+ shifta*shiftb*db_gab + shifta*shiftc*db_gac + &
+ shiftb*shiftc*db_gbc
+ halfshiftdgc = half*(shifta*shifta*dc_gaa + &
+ shiftb*shiftb*dc_gbb + shiftc*shiftc*dc_gcc) + &
+ shifta*shiftb*dc_gab + shifta*shiftc*dc_gac + &
+ shiftb*shiftc*dc_gbc
!!$ Contract the matter with derivatives of the metric
- halfTdgx = half*(txx*dx_gxx + tyy*dx_gyy + tzz*dx_gzz) +&
- txy*dx_gxy + txz*dx_gxz + tyz*dx_gyz
- halfTdgy = half*(txx*dy_gxx + tyy*dy_gyy + tzz*dy_gzz) +&
- txy*dy_gxy + txz*dy_gxz + tyz*dy_gyz
- halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) +&
- txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz
-
- sx_source = t00*&
- (halfshiftdgx - alp(i,j,k)*dx_alp) +&
- t0x*(shiftx*dx_gxx + shifty*dx_gxy + shiftz*dx_gxz) +&
- t0y*(shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) +&
- t0z*(shiftx*dx_gxz + shifty*dx_gyz + shiftz*dx_gzz) +&
- halfTdgx + rhoenthalpyW2*&
- (vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz)*&
+ halfTdga = half*(taa*da_gaa + tbb*da_gbb + tcc*da_gcc) +&
+ tab*da_gab + tac*da_gac + tbc*da_gbc
+ halfTdgb = half*(taa*db_gaa + tbb*db_gbb + tcc*db_gcc) +&
+ tab*db_gab + tac*db_gac + tbc*db_gbc
+ halfTdgc = half*(taa*dc_gaa + tbb*dc_gbb + tcc*dc_gcc) +&
+ tab*dc_gab + tac*dc_gac + tbc*dc_gbc
+
+ sa_source = t00*&
+ (halfshiftdga - alp(i,j,k)*da_alp) +&
+ t0a*(shifta*da_gaa + shiftb*da_gab + shiftc*da_gac) +&
+ t0b*(shifta*da_gab + shiftb*da_gbb + shiftc*da_gbc) +&
+ t0c*(shifta*da_gac + shiftb*da_gbc + shiftc*da_gcc) +&
+ halfTdga + rhoenthalpyW2*&
+ (vlowa*da_betaa + vlowb*da_betab + vlowc*da_betac)*&
invalp
- sy_source = t00*&
- (halfshiftdgy - alp(i,j,k)*dy_alp) +&
- t0x*(shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) +&
- t0y*(shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) +&
- t0z*(shiftx*dy_gxz + shifty*dy_gyz + shiftz*dy_gzz) +&
- halfTdgy + rhoenthalpyW2*&
- (vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz)*&
+ sb_source = t00*&
+ (halfshiftdgb - alp(i,j,k)*db_alp) +&
+ t0a*(shifta*db_gaa + shiftb*db_gab + shiftc*db_gac) +&
+ t0b*(shifta*db_gab + shiftb*db_gbb + shiftc*db_gbc) +&
+ t0c*(shifta*db_gac + shiftb*db_gbc + shiftc*db_gcc) +&
+ halfTdgb + rhoenthalpyW2*&
+ (vlowa*db_betaa + vlowb*db_betab + vlowc*db_betac)*&
invalp
- sz_source = t00*&
- (halfshiftdgz - alp(i,j,k)*dz_alp) +&
- t0x*(shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) +&
- t0y*(shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) +&
- t0z*(shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) +&
- halfTdgz + rhoenthalpyW2*&
- (vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz)*&
+ sc_source = t00*&
+ (halfshiftdgc - alp(i,j,k)*dc_alp) +&
+ t0a*(shifta*dc_gaa + shiftb*dc_gab + shiftc*dc_gac) +&
+ t0b*(shifta*dc_gab + shiftb*dc_gbb + shiftc*dc_gbc) +&
+ t0c*(shifta*dc_gac + shiftb*dc_gbc + shiftc*dc_gcc) +&
+ halfTdgc + rhoenthalpyW2*&
+ (vlowa*dc_betaa + vlowb*dc_betab + vlowc*dc_betac)*&
invalp
densrhs(i,j,k) = 0.d0
- srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source
- srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source
- srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source
+ srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sa_source
+ srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sb_source
+ srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sc_source
taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source
enddo
diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90
index b314a62..96a71db 100644
--- a/src/GRHydro_TVDReconstruct_drv.F90
+++ b/src/GRHydro_TVDReconstruct_drv.F90
@@ -14,15 +14,15 @@
#include "SpaceMask.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
-#define Bvecx(i,j,k) Bvec(i,j,k,1)
-#define Bvecy(i,j,k) Bvec(i,j,k,2)
-#define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bvecx(i,j,k) lBvec(i,j,k,1)
+#define Bvecy(i,j,k) lBvec(i,j,k,2)
+#define Bvecz(i,j,k) lBvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)
@@ -155,20 +155,20 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS)
rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask)
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
- vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
+ lvel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
- vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
+ lvel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask)
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
- vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
+ lvel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
eps, epsplus, epsminus, trivial_rp, hydro_excision_mask)
if(evolve_mhd.ne.0) then
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
- Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
+ lBvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
- Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
+ lBvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
- Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
+ lBvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
endif
else if (CCTK_EQUALS(recon_vars,"conservative")) then
diff --git a/src/GRHydro_TransformTensorBasis.F90 b/src/GRHydro_TransformTensorBasis.F90
new file mode 100644
index 0000000..1aeda46
--- /dev/null
+++ b/src/GRHydro_TransformTensorBasis.F90
@@ -0,0 +1,211 @@
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+
+subroutine GRHydroTransformADMToLocalBasis(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: i,j,k
+
+ if (jacobian_state .eq. 0) then
+ call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!")
+ end if
+
+ if (inverse_jacobian_state .eq. 0) then
+ call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!")
+ end if
+
+ !$OMP PARALLEL DO PRIVATE(i,j,k)
+ do k=1,cctk_lsh(3)
+ do j=1,cctk_lsh(2)
+ do i=1,cctk_lsh(1)
+ gaa(i,j,k) = iJ11(i,j,k)**2 * gxx(i,j,k) + &
+ 2.0d0 * iJ11(i,j,k) * iJ21(i,j,k) * gxy(i,j,k) + &
+ 2.0d0 * iJ11(i,j,k) * iJ31(i,j,k) * gxz(i,j,k) + iJ21(i,j,k)**2 * gyy(i,j,k) + &
+ 2.0d0 * iJ21(i,j,k) * iJ31(i,j,k) * gyz(i,j,k) + iJ31(i,j,k)**2 * gzz(i,j,k)
+
+ gab(i,j,k) = iJ11(i,j,k) * iJ12(i,j,k) * gxx(i,j,k) + &
+ (iJ11(i,j,k) * iJ22(i,j,k) + iJ21(i,j,k) * iJ12(i,j,k)) * gxy(i,j,k) + &
+ (iJ11(i,j,k) * iJ32(i,j,k) + iJ31(i,j,k) * iJ12(i,j,k)) * gxz(i,j,k) + &
+ iJ21(i,j,k) * iJ22(i,j,k) * gyy(i,j,k) + (iJ21(i,j,k) * iJ32(i,j,k) + &
+ iJ31(i,j,k) * iJ22(i,j,k)) * gyz(i,j,k) + iJ31(i,j,k) * iJ32(i,j,k) * gzz(i,j,k)
+
+ gac(i,j,k) = iJ11(i,j,k) * iJ13(i,j,k) * gxx(i,j,k) + (iJ11(i,j,k) * iJ23(i,j,k) + &
+ iJ21(i,j,k) * iJ13(i,j,k)) * gxy(i,j,k) + &
+ (iJ11(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ13(i,j,k)) * gxz(i,j,k) + &
+ iJ21(i,j,k) * iJ23(i,j,k) * gyy(i,j,k) + &
+ (iJ21(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ23(i,j,k)) * gyz(i,j,k) + &
+ iJ31(i,j,k) * iJ33(i,j,k) * gzz(i,j,k)
+
+ gbb(i,j,k) = iJ12(i,j,k)**2 * gxx(i,j,k) + &
+ 2.0d0 * iJ12(i,j,k) * iJ22(i,j,k) * gxy(i,j,k) + &
+ 2.0d0 * iJ12(i,j,k) * iJ32(i,j,k) * gxz(i,j,k) + iJ22(i,j,k)**2 * gyy(i,j,k) + &
+ 2.0d0 * iJ22(i,j,k) * iJ32(i,j,k) * gyz(i,j,k) + iJ32(i,j,k)**2 * gzz(i,j,k)
+
+ gbc(i,j,k) = iJ12(i,j,k) * iJ13(i,j,k) * gxx(i,j,k) + &
+ (iJ12(i,j,k) * iJ23(i,j,k) + iJ22(i,j,k) * iJ13(i,j,k)) * gxy(i,j,k) + &
+ (iJ12(i,j,k) * iJ33(i,j,k) + iJ32(i,j,k) * iJ13(i,j,k)) * gxz(i,j,k) + &
+ iJ22(i,j,k) * iJ23(i,j,k) * gyy(i,j,k) + &
+ (iJ22(i,j,k)* iJ33(i,j,k) + iJ32(i,j,k) * iJ23(i,j,k)) * gyz(i,j,k) + &
+ iJ32(i,j,k) * iJ33(i,j,k) * gzz(i,j,k)
+
+ gcc(i,j,k) = iJ13(i,j,k)**2 * gxx(i,j,k) + &
+ 2.0d0 * iJ13(i,j,k) * iJ23(i,j,k) * gxy(i,j,k) + &
+ 2.0d0 * iJ13(i,j,k) * iJ33(i,j,k) * gxz(i,j,k) + &
+ iJ23(i,j,k)**2 * gyy(i,j,k) + &
+ 2.0d0 * iJ23(i,j,k) * iJ33(i,j,k) * gyz(i,j,k) + &
+ iJ33(i,j,k)**2 * gzz(i,j,k)
+
+ ! Transform extrinsic curvature from global to local basis.
+ ! Since extrinsic has covariant indices, use inverse Jacobian
+ kaa(i,j,k) = iJ11(i,j,k)**2 * kxx(i,j,k) + &
+ 2.0d0 * iJ11(i,j,k) * iJ21(i,j,k) * kxy(i,j,k) + &
+ 2.0d0 * iJ11(i,j,k) * iJ31(i,j,k) * kxz(i,j,k) + &
+ iJ21(i,j,k)**2 * kyy(i,j,k) + &
+ 2.0d0 * iJ21(i,j,k) * iJ31(i,j,k) * kyz(i,j,k) + &
+ iJ31(i,j,k)**2 * kzz(i,j,k)
+
+ kab(i,j,k) = iJ11(i,j,k) * iJ12(i,j,k) * kxx(i,j,k) + &
+ (iJ11(i,j,k) * iJ22(i,j,k) + iJ21(i,j,k) * iJ12(i,j,k)) * kxy(i,j,k) + &
+ (iJ11(i,j,k) * iJ32(i,j,k) + iJ31(i,j,k) * iJ12(i,j,k)) * kxz(i,j,k) + &
+ iJ21(i,j,k) * iJ22(i,j,k) * kyy(i,j,k) + &
+ (iJ21(i,j,k) * iJ32(i,j,k) + iJ31(i,j,k) * iJ22(i,j,k)) * kyz(i,j,k) + &
+ iJ31(i,j,k) * iJ32(i,j,k) * kzz(i,j,k)
+
+
+ kac(i,j,k) = iJ11(i,j,k) * iJ13(i,j,k) * kxx(i,j,k) + &
+ (iJ11(i,j,k) * iJ23(i,j,k) + iJ21(i,j,k) * iJ13(i,j,k)) * kxy(i,j,k) + &
+ (iJ11(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ13(i,j,k)) * kxz(i,j,k) + &
+ iJ21(i,j,k) * iJ23(i,j,k) * kyy(i,j,k) + &
+ (iJ21(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ23(i,j,k)) * kyz(i,j,k) + &
+ iJ31(i,j,k) * iJ33(i,j,k) * kzz(i,j,k)
+
+ kbb(i,j,k) = iJ12(i,j,k)**2 * kxx(i,j,k) + &
+ 2.0d0 * iJ12(i,j,k) * iJ22(i,j,k) * kxy(i,j,k) + &
+ 2.0d0 * iJ12(i,j,k) * iJ32(i,j,k) * kxz(i,j,k) + &
+ iJ22(i,j,k)**2 * kyy(i,j,k) + &
+ 2.0d0 * iJ22(i,j,k) * iJ32(i,j,k) * kyz(i,j,k) + &
+ iJ32(i,j,k)**2 * kzz(i,j,k)
+
+ kbc(i,j,k) = iJ12(i,j,k) * iJ13(i,j,k) * kxx(i,j,k) + &
+ (iJ12(i,j,k) * iJ23(i,j,k) + iJ22(i,j,k) * iJ13(i,j,k)) * kxy(i,j,k) + &
+ (iJ12(i,j,k) * iJ33(i,j,k) + iJ32(i,j,k) * iJ13(i,j,k)) * kxz(i,j,k) + &
+ iJ22(i,j,k) * iJ23(i,j,k) * kyy(i,j,k) + &
+ (iJ22(i,j,k) * iJ33(i,j,k) + iJ32(i,j,k) * iJ23(i,j,k)) * kyz(i,j,k) + &
+ iJ32(i,j,k) * iJ33(i,j,k) * kzz(i,j,k)
+
+ kcc(i,j,k) = iJ13(i,j,k)**2 * kxx(i,j,k) + &
+ 2.0d0 * iJ13(i,j,k) * iJ23(i,j,k) * kxy(i,j,k) + &
+ 2.0d0 * iJ13(i,j,k) * iJ33(i,j,k) * kxz(i,j,k) + &
+ iJ23(i,j,k)**2 * kyy(i,j,k) + &
+ 2.0d0 * iJ23(i,j,k) * iJ33(i,j,k) * kyz(i,j,k) + &
+ iJ33(i,j,k)**2 * kzz(i,j,k)
+
+ ! Transform shift from global to local basis.
+ ! Since shift has contravariant index, use Jacobian
+ betaa(i,j,k) = betax(i,j,k)*J11(i,j,k) + &
+ betay(i,j,k)*J12(i,j,k) + &
+ betaz(i,j,k)*J13(i,j,k)
+ betab(i,j,k) = betax(i,j,k)*J21(i,j,k) + &
+ betay(i,j,k)*J22(i,j,k) + &
+ betaz(i,j,k)*J23(i,j,k)
+ betac(i,j,k) = betax(i,j,k)*J31(i,j,k) + &
+ betay(i,j,k)*J32(i,j,k) + &
+ betaz(i,j,k)*J33(i,j,k)
+ enddo
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+end subroutine GRHydroTransformADMToLocalBasis
+
+
+
+subroutine GRHydroTransformPrimToLocalBasis(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_INT :: i,j,k
+
+ if (jacobian_state .eq. 0) then
+ call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!")
+ end if
+
+ if (inverse_jacobian_state .eq. 0) then
+ call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!")
+ end if
+
+ !call CCTK_INFO("Transforming primitives to local basis.")
+
+ !$OMP PARALLEL DO PRIVATE(i, j, k)
+ do k = 1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ ! Transform primitive velocity from global to local basis.
+ ! Since velocity has contravariant index, use Jacobian
+ lvel(i,j,k,1) = vel(i,j,k,1)*J11(i,j,k) + vel(i,j,k,2)*J12(i,j,k) + vel(i,j,k,3)*J13(i,j,k)
+ lvel(i,j,k,2) = vel(i,j,k,1)*J21(i,j,k) + vel(i,j,k,2)*J22(i,j,k) + vel(i,j,k,3)*J23(i,j,k)
+ lvel(i,j,k,3) = vel(i,j,k,1)*J31(i,j,k) + vel(i,j,k,2)*J32(i,j,k) + vel(i,j,k,3)*J33(i,j,k)
+
+ ! Transform primitive B-field from global to local basis.
+ ! Since B-field has contravariant index, use Jacobian
+! lBvec(i,j,k,1) = Bvec(i,j,k,1)*J11(i,j,k) + Bvec(i,j,k,2)*J12(i,j,k) + Bvec(i,j,k,3)*J13(i,j,k)
+! lBvec(i,j,k,2) = Bvec(i,j,k,1)*J21(i,j,k) + Bvec(i,j,k,2)*J22(i,j,k) + Bvec(i,j,k,3)*J23(i,j,k)
+! lBvec(i,j,k,3) = Bvec(i,j,k,1)*J31(i,j,k) + Bvec(i,j,k,2)*J32(i,j,k) + Bvec(i,j,k,3)*J33(i,j,k)
+ end do
+ end do
+ end do
+ !$OMP END PARALLEL DO
+
+end subroutine GRHydroTransformPrimToLocalBasis
+
+
+
+subroutine GRHydroTransformPrimToGlobalBasis(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_INT :: i,j,k
+
+ if (jacobian_state .eq. 0) then
+ call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!")
+ end if
+
+ if (inverse_jacobian_state .eq. 0) then
+ call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!")
+ end if
+
+ !call CCTK_INFO("Transforming primitives to global basis.")
+
+ !$OMP PARALLEL DO PRIVATE(i, j, k)
+ do k = 1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ ! Transform primitive velocity from local to global basis.
+ ! Since velocity has contravariant index, use inverse Jacobian
+ vel(i,j,k,1) = lvel(i,j,k,1)*iJ11(i,j,k) + lvel(i,j,k,2)*iJ12(i,j,k) + lvel(i,j,k,3)*iJ13(i,j,k)
+ vel(i,j,k,2) = lvel(i,j,k,1)*iJ21(i,j,k) + lvel(i,j,k,2)*iJ22(i,j,k) + lvel(i,j,k,3)*iJ23(i,j,k)
+ vel(i,j,k,3) = lvel(i,j,k,1)*iJ31(i,j,k) + lvel(i,j,k,2)*iJ32(i,j,k) + lvel(i,j,k,3)*iJ33(i,j,k)
+
+ end do
+ end do
+ end do
+ !$OMP END PARALLEL DO
+
+end subroutine GRHydroTransformPrimToGlobalBasis
+
+
+
+
+
diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90
index d6bf540..0bb092c 100644
--- a/src/GRHydro_UpdateMask.F90
+++ b/src/GRHydro_UpdateMask.F90
@@ -14,15 +14,15 @@
#include "GRHydro_Macros.h"
#include "SpaceMask.h"
-#define velx(i,j,k) vel(i,j,k,1)
-#define vely(i,j,k) vel(i,j,k,2)
-#define velz(i,j,k) vel(i,j,k,3)
-#define velx_p(i,j,k) vel_p(i,j,k,1)
-#define vely_p(i,j,k) vel_p(i,j,k,2)
-#define velz_p(i,j,k) vel_p(i,j,k,3)
-#define velx_p_p(i,j,k) vel_p_p(i,j,k,1)
-#define vely_p_p(i,j,k) vel_p_p(i,j,k,2)
-#define velz_p_p(i,j,k) vel_p_p(i,j,k,3)
+#define velx(i,j,k) lvel(i,j,k,1)
+#define vely(i,j,k) lvel(i,j,k,2)
+#define velz(i,j,k) lvel(i,j,k,3)
+#define velx_p(i,j,k) lvel_p(i,j,k,1)
+#define vely_p(i,j,k) lvel_p(i,j,k,2)
+#define velz_p(i,j,k) lvel_p(i,j,k,3)
+#define velx_p_p(i,j,k) lvel_p_p(i,j,k,1)
+#define vely_p_p(i,j,k) lvel_p_p(i,j,k,2)
+#define velz_p_p(i,j,k) lvel_p_p(i,j,k,3)
subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS)
@@ -246,9 +246,9 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
velx(i,j,k) = 0.0d0
vely(i,j,k) = 0.0d0
velz(i,j,k) = 0.0d0
- det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \
- gyy(i,j,k), gyz(i,j,k), gzz(i,j,k))
-
+ det = SPATIAL_DETERMINANT(gaa(i,j,k), gab(i,j,k), gac(i,j,k), \
+ gbb(i,j,k), gbc(i,j,k), gcc(i,j,k))
+
if(evolve_temper.ne.0) then
! ! set the temperature to be relatively low
temperature(i,j,k) = grhydro_hot_atmo_temp
@@ -259,8 +259,8 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
press(i,j,k),keyerr,anyerr)
call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
- i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),gxy(i,j,k),&
- gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+ i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),gab(i,j,k),&
+ gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k), &
det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k),&
@@ -268,19 +268,19 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS)
y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k)
else
call prim2conpolytype(GRHydro_polytrope_handle, &
- gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), &
- gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, &
+ gaa(i,j,k), gab(i,j,k), gac(i,j,k), &
+ gbb(i,j,k), gbc(i,j,k), gcc(i,j,k), det, &
dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k))
if (wk_atmosphere .eq. 0) then
- atmosphere_mask(i, j, k) = 0
- SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
- not_atmosphere)
+ atmosphere_mask(i, j, k) = 0
+ SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
+ not_atmosphere)
end if
endif
- end if
+ end if
end do
end do
@@ -330,8 +330,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
if (rho(i,j,k) .le. rho_min) then
rho(i,j,k) = rho_min
- det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \
- gyy(i,j,k), gyz(i,j,k), gzz(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa(i,j,k), gab(i,j,k), gac(i,j,k), \
+ gbb(i,j,k), gbc(i,j,k), gcc(i,j,k))
velx(i,j,k) = 0.0d0
vely(i,j,k) = 0.0d0
@@ -347,8 +347,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
press(i,j,k),keyerr,anyerr)
call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
- i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),gxy(i,j,k),&
- gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+ i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),gab(i,j,k),&
+ gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k), &
det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k),&
@@ -360,8 +360,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
call prim2conpolytype(eos_handle, &
- gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), &
- gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, &
+ gaa(i,j,k), gab(i,j,k), gac(i,j,k), &
+ gbb(i,j,k), gbc(i,j,k), gcc(i,j,k), det, &
dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), &
velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k))
@@ -374,8 +374,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
vely_p(i,j,k) = 0.0d0
velz_p(i,j,k) = 0.0d0
- det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \
- gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa_p(i,j,k), gab_p(i,j,k), gac_p(i,j,k), \
+ gbb_p(i,j,k), gbc_p(i,j,k), gcc_p(i,j,k))
if(evolve_temper.ne.0) then
! set the temperature to be relatively low
@@ -387,8 +387,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
press_p(i,j,k),keyerr,anyerr)
call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
- i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p(i,j,k),gxy_p(i,j,k),&
- gxz_p(i,j,k),gyy_p(i,j,k),gyz_p(i,j,k),gzz_p(i,j,k), &
+ i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa_p(i,j,k),gab_p(i,j,k),&
+ gac_p(i,j,k),gbb_p(i,j,k),gbc_p(i,j,k),gcc_p(i,j,k), &
det,dens_p(i,j,k),scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), &
tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), &
velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k),&
@@ -400,8 +400,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr)
call prim2conpolytype(eos_handle, &
- gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), &
- gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, &
+ gaa_p(i,j,k), gab_p(i,j,k), gac_p(i,j,k), &
+ gbb_p(i,j,k), gbc_p(i,j,k), gcc_p(i,j,k), det, &
dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), &
tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), &
velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k))
@@ -415,8 +415,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
velx_p_p(i,j,k) = 0.0d0
vely_p_p(i,j,k) = 0.0d0
velz_p_p(i,j,k) = 0.0d0
- det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \
- gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k))
+ det = SPATIAL_DETERMINANT(gaa_p_p(i,j,k), gab_p_p(i,j,k), gac_p_p(i,j,k), \
+ gbb_p_p(i,j,k), gbc_p_p(i,j,k), gcc_p_p(i,j,k))
if(evolve_temper.ne.0) then
! set the temperature to be relatively low
@@ -427,8 +427,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),&
press_p_p(i,j,k),keyerr,anyerr)
call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,&
- i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p_p(i,j,k),gxy_p_p(i,j,k),&
- gxz_p_p(i,j,k),gyy_p_p(i,j,k),gyz_p_p(i,j,k),gzz_p_p(i,j,k), &
+ i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa_p_p(i,j,k),gab_p_p(i,j,k),&
+ gac_p_p(i,j,k),gbb_p_p(i,j,k),gbc_p_p(i,j,k),gcc_p_p(i,j,k), &
det,dens_p_p(i,j,k),scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), &
tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &
velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k),&
@@ -440,8 +440,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr)
call prim2conpolytype(eos_handle, &
- gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), &
- gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, &
+ gaa_p_p(i,j,k), gab_p_p(i,j,k), gac_p_p(i,j,k), &
+ gbb_p_p(i,j,k), gbc_p_p(i,j,k), gcc_p_p(i,j,k), det, &
dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), &
tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &
velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k))
diff --git a/src/make.code.defn b/src/make.code.defn
index a4e0873..6ff530b 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -10,7 +10,6 @@ SRCS = Utils.F90 \
GRHydro_DivergenceClean.F90 \
GRHydro_Eigenproblem.F90 \
GRHydro_Eigenproblem_Marquina.F90 \
- GRHydro_ENOReconstruct_drv.F90 \
GRHydro_ENOReconstruct.F90 \
GRHydro_ENOScalars.F90 \
GRHydro_EOS.c \
@@ -20,7 +19,6 @@ SRCS = Utils.F90 \
GRHydro_Loop.F90 \
GRHydro_Minima.F90 \
GRHydro_Minima.cc \
- GRHydro_PPMReconstruct_drv.F90 \
GRHydro_PPM.F90 \
GRHydro_ParamCheck.F90 \
GRHydro_Particle.F90 \
@@ -38,7 +36,6 @@ SRCS = Utils.F90 \
GRHydro_SlopeLimiter.F90 \
GRHydro_Source.F90 \
GRHydro_Startup.F90 \
- GRHydro_TVDReconstruct_drv.F90 \
GRHydro_TVDReconstruct.F90 \
GRHydro_Marquina.F90 \
GRHydro_UpdateMask.F90 \
@@ -52,15 +49,19 @@ SRCS = Utils.F90 \
GRHydro_EigenproblemM.F90 \
GRHydro_FluxM.F90 \
GRHydro_HLLEM.F90 \
- GRHydro_PPMMReconstruct_drv.F90 \
GRHydro_PPMM.F90 \
GRHydro_Prim2ConM.F90 \
GRHydro_RiemannSolveM.F90 \
GRHydro_SourceM.F90 \
GRHydro_TmunuM.F90 \
GRHydro_UpdateMaskM.F90 \
- GRHydro_UtilsM.F90
-
+ GRHydro_UtilsM.F90 \
+ GRHydro_TransformTensorBasis.F90 \
+ GRHydro_Jacobian_state.c \
+ GRHydro_PPMMReconstruct_drv.F90\
+ GRHydro_ENOReconstruct_drv.F90\
+ GRHydro_PPMReconstruct_drv.F90\
+ GRHydro_TVDReconstruct_drv.F90
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-09 00:41:57 +0000