GRHydro_InitData: initialize c2p2c values as parameters and update routine to test inversions in the gr context.

From: Bruno Coutinho Mundim <bcmsma@astro.rit.edu>

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@192 ac85fae7-cede-4708-beff-ae01c7fa1c26

2 files changed, 145 insertions, 38 deletions

diff --git a/param.ccl b/param.ccl
index 67e51e1..bb9dad6 100644
--- a/param.ccl
+++ b/param.ccl
@@ -126,6 +126,8 @@ BOOLEAN attenuate_atmosphere "Attenuate the velocity in the atmosphere"
 {
 } "no"
 
+# Initial magnetic field used by different tests
+
 REAL Bx_init "Initial B-field in the x-dir"
 {
   *:*     :: "Anything"
@@ -141,6 +143,95 @@ REAL Bz_init "Initial B-field in the z-dir"
   *:*     :: "Anything"
 } 0.0
 
+# Initial primitive values:
+REAL rho_init "Initial rest mass density"
+{
+ (0:*    :: "Anything positive."
+} 1.0d-6
+
+REAL velx_init "Initial x velocity" 
+{
+ *:*    :: "Anything."
+} 1.0d-1
+
+REAL vely_init "Initial y velocity" 
+{
+ *:*    :: "Anything."
+} 1.0d-1
+
+REAL velz_init "Initial z velocity" 
+{
+ *:*    :: "Anything."
+} 1.0d-1
+
+REAL eps_init "Initial specific internal energy"
+{
+ (0:*    :: "Anything positive."
+} 1.0d-6
+
+REAL press_init "Initial pressure"
+{
+ (0:*    :: "Anything positive."
+} 6.666666666666667d-7
+
+# Initial conservative values used by c2p tests:
+REAL dens_init "Initial conserved mass density"
+{
+ (0:*    :: "Anything positive."
+} 1.29047362
+
+REAL sx_init "Initial x component of conserved momentum density"
+{
+ *:*    :: "Anything."
+} 0.166666658
+
+REAL sy_init "Initial y component of conserved momentum density"
+{
+ *:*    :: "Anything."
+} 0.166666658
+
+REAL sz_init "Initial z component of conserved momentum density"
+{
+ *:*    :: "Anything."
+} 0.166666658
+
+REAL tau_init "Initial conserved total energy density"
+{
+ (0:*    :: "Anything positive."
+} 0.484123939
+
+# Initial values for 3-metric components. Default to euclidian 3-metric
+REAL gxx_init "Initial xx metric componenent"
+{
+ *:*    :: "Anything, but be carefull to set a positive definite 3-metric!"
+} 1.0
+
+REAL gxy_init "Initial xy metric componenent"
+{
+ *:*    :: "Anything, but be carefull to set a positive definite 3-metric!"
+} 0.0
+
+REAL gxz_init "Initial xz metric componenent"
+{
+ *:*    :: "Anything, but be carefull to set a positive definite 3-metric!"
+} 0.0
+
+REAL gyy_init "Initial yy metric componenent"
+{
+ *:*    :: "Anything, but be carefull to set a positive definite 3-metric!"
+} 1.0
+
+REAL gyz_init "Initial yz metric componenent"
+{
+ *:*    :: "Anything, but be carefull to set a positive definite 3-metric!"
+} 0.0
+
+REAL gzz_init "Initial zz metric componenent"
+{
+ *:*    :: "Anything, but be carefull to set a positive definite 3-metric!"
+} 1.0
+
+
 KEYWORD monopole_type "Which kind of monopole?"
 {
   "Point" :: "Single point with Bx /= 0"
diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90
index c167504..9394005 100644
--- a/src/GRHydro_C2P2CM.F90
+++ b/src/GRHydro_C2P2CM.F90
@@ -11,6 +11,8 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 #include "cctk_Functions.h"
+#include "GRHydro_Macros.h"
+
 
  /*@@
    @routine    c2p2cM
@@ -37,13 +39,13 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_REAL :: det
+  CCTK_REAL :: det, sdet, invdet
   CCTK_REAL :: uxx,uxy,uxz,uyy,uyz,uzz
   CCTK_REAL :: gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send
   CCTK_REAL :: dens_send,sx_send,sy_send,sz_send,tau_send
   CCTK_REAL :: bconsx_send, bconsy_send, bconsz_send
   CCTK_REAL :: rho_send,velx_send,vely_send,velz_send,eps_send
-  CCTK_REAL :: press_send,w_lorentz_send,x_send,y_send,z_send,r_send
+  CCTK_REAL :: press_send,w_lorentz_send
   CCTK_REAL :: bvcx_send, bvcy_send, bvcz_send, b2_send
   CCTK_REAL :: C2P_failed
   CCTK_INT :: epsnegative
@@ -58,50 +60,64 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
 
   call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
   
-  x_send = 0.0d0
-  y_send = 0.0d0
-  z_send = 0.0d0
-  r_send = 0.0d0
-  
-  gxx_send = 1.0d0
-  gyy_send = 1.0d0 
-  gzz_send = 1.0d0
-  gxy_send = 0.0d0
-  gxz_send = 0.0d0
-  gyz_send = 0.0d0
+  gxx_send = gxx_init
+  gxy_send = gxy_init
+  gxz_send = gxz_init
+  gyy_send = gyy_init
+  gyz_send = gyz_init
+  gzz_send = gzz_init
   
-  det = 1.0d0
-  
-  uxx = 1.0d0
-  uyy = 1.0d0 
-  uzz = 1.0d0
-  uxy = 0.0d0
-  uxz = 0.0d0
-  uyz = 0.0d0
+  det = SPATIAL_DETERMINANT(gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send)
+  sdet = sqrt(det)
+  invdet = 1.d0 / det
+  uxx = (-gyz_send**2 + gyy_send*gzz_send)*invdet
+  uxy = (gxz_send*gyz_send - gxy_send*gzz_send)*invdet
+  uyy = (-gxz_send**2 + gxx_send*gzz_send)*invdet
+  uxz = (-gxz_send*gyy_send + gxy_send*gyz_send)*invdet
+  uyz = (gxy_send*gxz_send - gxx_send*gyz_send)*invdet
+  uzz = (-gxy_send**2 + gxx_send*gyy_send)*invdet
 
 ! Initialize the velocity as GRHydro_Con2PrimM_pt may use these
 ! values as initial guess for its Newton-Raphson procedure.
-  velx_send = 0.1d0
-  vely_send = 0.1d0 
-  velz_send = 0.1d0
-  
-  dens_send = 1.29047362d0
-  sx_send   = 0.166666658d0
-  sy_send   = 0.166666658d0
-  sz_send   = 0.166666658d0
-  tau_send  = 0.484123939d0
+!  velx_send = 0.1d0
+!  vely_send = 0.1d0 
+!  velz_send = 0.1d0
+!  
+!  dens_send = 1.29047362d0
+!  sx_send   = 0.166666658d0
+!  sy_send   = 0.166666658d0
+!  sz_send   = 0.166666658d0
+!  tau_send  = 0.484123939d0
+
+  velx_send = velx_init
+  vely_send = vely_init
+  velz_send = velz_init
+
+  dens_send = dens_init
+  sx_send   = sx_init
+  sy_send   = sy_init
+  sz_send   = sz_init
+  tau_send  = tau_init
 
   bvcx_send   = Bx_init
   bvcy_send   = By_init
   bvcz_send   = Bz_init
-  bconsx_send   = Bx_init
-  bconsy_send   = By_init
-  bconsz_send   = Bz_init
-
-  rho_send=1.0d-6
-  eps_send = 1.0d-6
-  press_send = 6.666666666666667d-7
-  w_lorentz_send = 1.0d0
+  bconsx_send   = sdet*Bx_init
+  bconsy_send   = sdet*By_init
+  bconsz_send   = sdet*Bz_init
+
+  rho_send = rho_init
+  eps_send = eps_init
+  press_send = press_init
+  w_lorentz_send = sqrt(1.0d0-(gxx_send*velx_send**2+gyy_send*vely_send**2+ &
+                               gzz_send*velz_send**2 &
+                               +2.0d0*(gxy_send*velx_send*vely_send+ &
+                                       gxz_send*velx_send*velz_send+ &
+                                       gyz_send*vely_send*velz_send &
+                                      ) &
+                               ) &
+                        ) 
+  w_lorentz_send = 1.0d0/w_lorentz_send
 
   epsnegative = 0