RIT MHD development:

 Add #include "cctk_Functions.h" and DECLARE_CCTK_FUNCTIONS 
to GRHydro_C2P2CM.F90.
 
 Solve several issues of rank mismatch in the EOS Omni 
interface.

 Initialize the velocities at GRHydro_C2P2CM.F90 since
GRHydro_Con2PrimM_pt.c uses past values as initial guess
for its Newton-Raphson procedure. These unitialized values
were causing C2P2CM to fail from time to time, since it
would access memory positions with very large numerical
values, i.e. v>>1. C2P2CM test always passes now.



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

4 files changed, 68 insertions, 35 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index 6383fd4..bee5e35 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -118,7 +118,7 @@ For testing purposes, this routine sets all the points to the values of the atmo
 \subsection{Simple wave}
 \label{sec:simple-wave}
 
-This routine stes initial data for a simple wave with sinusoidal initial function for the velocity,
+This routine testes initial data for a simple wave with sinusoidal initial function for the velocity,
 as described in Anile, Miller, Motta, {\it Formation and damping of relativistic strong
   shocks},Phys. Fluids {\bf 26}, 1450 (1983).
 
diff --git a/interface.ccl b/interface.ccl
index b8f7af7..4f3f7b5 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -79,12 +79,29 @@ SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \
 			CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed)
 
 
+void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \
+                           CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \
+                           CCTK_REAL INOUT tau, CCTK_REAL INOUT rho, \
+                           CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \
+                           CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \
+                           CCTK_REAL INOUT w_lorentz, \
+                           CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \
+                           CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \
+                           CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \
+                           CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \
+                           CCTK_REAL INOUT det, \
+                           CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \
+                           CCTK_REAL INOUT Bvecsq, \
+                           CCTK_INT OUT epsnegative, \
+                           CCTK_REAL OUT retval)
+
 USES FUNCTION SpatialDet
 USES FUNCTION Prim2ConPoly
 USES FUNCTION Prim2ConGen
 USES FUNCTION Prim2ConPolyM
 USES FUNCTION Prim2ConGenM
 USES FUNCTION Con2PrimPoly
+USES FUNCTION Con2PrimGenM
 
 CCTK_INT FUNCTION EOS_Omni_GetHandle(CCTK_STRING IN name)
 USES FUNCTION EOS_Omni_GetHandle
diff --git a/src/GRHydro_C2P2CM.F90 b/src/GRHydro_C2P2CM.F90
index 3bd61b6..c382037 100644
--- a/src/GRHydro_C2P2CM.F90
+++ b/src/GRHydro_C2P2CM.F90
@@ -10,6 +10,7 @@
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
+#include "cctk_Functions.h"
 
  /*@@
    @routine    c2p2cM
@@ -29,10 +30,12 @@
 
 subroutine c2p2cM(CCTK_ARGUMENTS)
 
+
   implicit none
 
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
 
   integer didit,i,j,k,nx,ny,nz
   CCTK_REAL det
@@ -42,19 +45,18 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   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 bvcx_send, bvcy_send, bvcz_send, b2_send
-  CCTK_REAL pmin, epsmin, local_gam
-  CCTK_INT C2P_failed
-  logical epsnegative
+  CCTK_REAL :: pmin(1), epsmin(1), local_gam(1)
+  CCTK_REAL C2P_failed
+  CCTK_INT epsnegative
+
+  CCTK_REAL :: rhoval(1), epsval(1)
+  CCTK_REAL :: gamval
 
 ! begin EOS Omni vars
-  integer :: n = 1
-  integer :: keytemp = 0
-  integer :: anyerr = 0
-  integer :: keyerr(1) = 0
-  real*8  :: xpress = 0.0d0
-  real*8  :: xeps = 0.0d0
-  real*8  :: xtemp = 0.0d0
-  real*8  :: xye = 0.0d0
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1)
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0;xeps(1)=0.0d0
 ! end EOS Omni vars
 
   call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS")
@@ -83,6 +85,12 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   uxy = 0.0d0
   uxz = 0.0d0
   uyz = 0.0d0
+
+! 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
@@ -98,22 +106,27 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   press_send = 6.666666666666667d-7
   w_lorentz_send = 1.0d0
 
-  epsnegative = .false.
+!  epsnegative = .false.
+  epsnegative = 0
   
-  GRHydro_rho_min = 1.e-10
+  GRHydro_rho_min = 1.0d-10
+  rhoval(1) = GRHydro_rho_min
+  epsval(1) = 1.0d0
 
   call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-       GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr)
+       rhoval,epsval,xtemp,xye,pmin,keyerr,anyerr)
 
   call EOS_Omni_EpsFromPress(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-       GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
+       rhoval,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
 
-  local_gam=0.0
+  local_gam(1)=0.0d0
+  rhoval(1) = 1.0d0
   call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
-       1.d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr)
-  local_gam = local_gam+1.0
+       rhoval,epsval,xtemp,xye,local_gam,keyerr,anyerr)
+  local_gam(1) = local_gam(1)+1.0
+  gamval = local_gam(1)
 
-  C2P_failed = 0
+  C2P_failed = 0.0d0
 
   write(*,*) 'C2P2CM test: initial values.'
   write(*,*) '   conservative variables: '
@@ -129,7 +142,8 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   write(*,*) '   Bvecz  : ',bvcz_send
   
   write(*,*) 'C2P2CM test: getting the associated primitive variables.'
-  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,dens_send,sx_send,sy_send,sz_send, &
+!  call GRHydro_Con2PrimM_pt(GRHydro_eos_handle,local_gam,dens_send,sx_send,sy_send,sz_send, &
+  call Con2PrimGenM(GRHydro_eos_handle,gamval,dens_send,sx_send,sy_send,sz_send, &
        tau_send,rho_send,velx_send,vely_send,velz_send, &
        eps_send,press_send,w_lorentz_send, &
        gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,&
@@ -149,9 +163,11 @@ subroutine c2p2cM(CCTK_ARGUMENTS)
   write(*,*) '   Bvecx  : ',bvcx_send
   write(*,*) '   Bvecy  : ',bvcy_send
   write(*,*) '   Bvecz  : ',bvcz_send
+  write(*,*) '   C2P_failed : ',C2P_failed
   
   write(*,*) 'C2P2CM test: converting back to conserved variables.'
-  call prim2conM(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+!  call prim2conM(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
+  call Prim2ConGenM(GRHydro_eos_handle,gxx_send, gxy_send, gxz_send, gyy_send, gyz_send, gzz_send, det, &
        dens_send, sx_send, sy_send, sz_send, tau_send, bvcx_send, bvcy_send, bvcz_send, rho_send, &
        velx_send, vely_send, velz_send, eps_send, press_send, w_lorentz_send) 
   
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90
index 96b0e7c..186cf69 100644
--- a/src/GRHydro_ShockTubeM.F90
+++ b/src/GRHydro_ShockTubeM.F90
@@ -118,22 +118,22 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           epsl = 1500.d0
           epsr = 1.5d-2
        else if (CCTK_EQUALS(shock_case,"Balsara1")) then
-          rhol = 1.d0
+          rhol = 1.0d0
           rhor = 0.125d0
-          velxl = 0.d0
-          velxr = 0.d0
-          velyl = 0.d0
-          velyr = 0.d0
-          velzl = 0.d0
-          velzr = 0.d0
+          velxl = 0.0d0
+          velxr = 0.0d0
+          velyl = 0.0d0
+          velyr = 0.0d0
+          velzl = 0.0d0
+          velzr = 0.0d0
           bvcxl=0.5d0
           bvcxr=0.5d0
-          bvcyl=1.d0
-          bvcyr=-1.d0
-          bvczl=0.d0
-          bvczr=0.d0
-          epsl = 1.0/rhol
-          epsr = 0.1/rhor
+          bvcyl=1.0d0
+          bvcyr=-1.0d0
+          bvczl=0.0d0
+          bvczr=0.0d0
+          epsl = 1.0d0/rhol
+          epsr = 0.1d0/rhor
        else if (CCTK_EQUALS(shock_case,"Balsara2")) then
           rhol = 1.d0
           rhor = 1.d0