do not set conservatives, remove TOV_Atmosphere

remove TOV_Atmosphere[] from the code and no longer set the conservative
variables. Instead TOVSolver relies on GRHydro to set the conservatives in
InitialPrim2Con and to enforce atmosphere in initialatmospherereset. This also
fixes any asymmetry in atmosphere handling between the stars.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TOVSolver/trunk@143 1bdb13ef-5d69-4035-bb54-08abeb3aa7f1

6 files changed, 6 insertions, 112 deletions

diff --git a/interface.ccl b/interface.ccl
index b118756..3fd5ce0 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,18 +1,12 @@
 # Interface definition for thorn TOVSolver
 
 implements: TOVSolver
-inherits: GRHydro, Constants
+inherits: ADMBase, HydroBase, Constants, StaticConformal
 
 USES INCLUDE: constants.h
 
 public:
 
-void FUNCTION SpatialDet(CCTK_REAL IN gxx, CCTK_REAL IN gxy, \
-                         CCTK_REAL IN gxz, CCTK_REAL IN gyy, \
-                         CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
-                         CCTK_REAL OUT det)
-USES FUNCTION SpatialDet
-
 CCTK_INT FUNCTION Set_Rho_ADM(              \
     CCTK_POINTER_TO_CONST IN  cctkGH,       \
     CCTK_INT              IN  size,         \
diff --git a/schedule.ccl b/schedule.ccl
index d54c9c9..89b89dc 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -30,9 +30,6 @@ schedule GROUP TOV_Initial_Data IN HydroBase_Initial
   SYNC: press
   SYNC: eps
   SYNC: vel
-  SYNC: dens
-  SYNC: tau
-  SYNC: scon
   SYNC: w_lorentz
 } "Group for the TOV initial data"
 
diff --git a/src/external.inc b/src/external.inc
index 73010d9..73a8918 100644
--- a/src/external.inc
+++ b/src/external.inc
@@ -57,7 +57,6 @@ CCTK_INT TOV_Set_Rho_ADM(CCTK_POINTER_TO_CONST cctkGH,
 
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
-  assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
   assert(TOV_rbar_1d!=0);
@@ -130,7 +129,6 @@ CCTK_INT TOV_Set_Momentum_Source(
 
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
-  assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
   assert(TOV_rbar_1d!=0);
@@ -200,7 +198,6 @@ CCTK_INT TOV_Set_Initial_Guess_for_u(
 
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
-  assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
   assert(TOV_rbar_1d!=0);
@@ -587,7 +584,6 @@ CCTK_INT TOV_Rescale_Sources(
                 );
       }
 
-      /* FIXME: check for Atmosphere */
     }
   }
   if (debug)
diff --git a/src/tov.c b/src/tov.c
index a5d206d..ce7d1db 100644
--- a/src/tov.c
+++ b/src/tov.c
@@ -42,7 +42,6 @@
 CCTK_REAL * TOV_Surface=0;
 CCTK_REAL * TOV_R_Surface=0;
 CCTK_REAL * TOV_RProp_Surface=0;
-CCTK_REAL * TOV_Atmosphere=0;
 
 CCTK_REAL * TOV_r_1d=0;
 CCTK_REAL * TOV_rbar_1d=0;
@@ -200,7 +199,6 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
   assert(TOV_RProp_Surface!=0);
-  assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
   assert(TOV_rbar_1d!=0);
@@ -222,20 +220,6 @@ void TOV_C_Integrate_RHS(CCTK_ARGUMENTS)
     else
       rho_central=TOV_Rho_Central[star];
 
-    if (rho_abs_min > 0.0)
-      TOV_Atmosphere[star] = rho_abs_min;
-    else
-      TOV_Atmosphere[star] = rho_central * rho_rel_min;
-
-    if (initial_rho_abs_min > 0.0)
-      TOV_Atmosphere[star] = initial_rho_abs_min;
-    else
-      if (initial_rho_rel_min > 0.0)
-        TOV_Atmosphere[star] = rho_central * initial_rho_rel_min;
-
-    if (initial_atmosphere_factor > 0.0)
-      TOV_Atmosphere[star] *= initial_atmosphere_factor;
-
     /* Set conformal state like set in parameter file if we do not use
      * the old initial data. In this case we have to use what we get */
     if (!TOV_Use_Old_Initial_Data)
@@ -526,10 +510,9 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
 
   CCTK_INT i,j,k, i3D, star;
   CCTK_REAL *r_to_star;
-  CCTK_REAL g_diag, max_g_diag, max_rho, det, sqrt_det;
+  CCTK_REAL g_diag, max_g_diag, max_rho;
   CCTK_REAL my_velx, my_vely, my_velz, my_psi4;
-  CCTK_REAL vlowx, vlowy, vlowz;
-  CCTK_REAL D_h_w, PI, local_tiny;
+  CCTK_REAL PI, local_tiny;
 
   CCTK_INT tov_lapse;
 
@@ -543,7 +526,6 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
 
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
-  assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
   assert(TOV_rbar_1d!=0);
@@ -608,14 +590,9 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
   {
     CCTK_INFO("Not using old matter initial data");
     TOV_C_fill(rho,        LSH_MAX_I+1, 0.0);
-    TOV_C_fill(dens,       LSH_MAX_I+1, 0.0);
     TOV_C_fill(eps,        LSH_MAX_I+1, 0.0);
     TOV_C_fill(press,      LSH_MAX_I+1, 0.0);
-    TOV_C_fill(tau,        LSH_MAX_I+1, 0.0);
     TOV_C_fill(w_lorentz,  LSH_MAX_I+1, 0.0);
-    TOV_C_fill(sx,         LSH_MAX_I+1, 0.0);
-    TOV_C_fill(sy,         LSH_MAX_I+1, 0.0);
-    TOV_C_fill(sz,         LSH_MAX_I+1, 0.0);
     TOV_C_fill(velx,       LSH_MAX_I+1, 0.0);
     TOV_C_fill(vely,       LSH_MAX_I+1, 0.0);
     TOV_C_fill(velz,       LSH_MAX_I+1, 0.0);
@@ -789,12 +766,6 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
         rho[i3D] = max_rho;
         eps[i3D] = eps_point[star];
         press[i3D] = press_point[star];
-        /* Set atmosphere according to chosen star */
-        if(rho[i3D] <= TOV_Atmosphere[star]) {
-          rho[i3D] = TOV_Atmosphere[star];
-	  press[i3D] = TOV_K * pow(rho[i3D],TOV_Gamma);
-          eps[i3D] = press[i3D]/(TOV_Gamma-1.0)/rho[i3D];
-	}
 
       }
       else if (CCTK_EQUALS(TOV_Combine_Method, "average"))
@@ -832,33 +803,13 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
           eps[i3D] += eps_point[star_i];
           press[i3D] += press_point[star_i];
           /* we still have to know if we are inside one star - and which */
-          if ( (rho_point[star_i] > max_rho) &&
-               (rho_point[star_i] > TOV_Atmosphere[star_i]))
+          if (rho_point[star_i] > max_rho)
           {
             max_rho=rho_point[star_i];
             star=star_i;
           }
-
-
-          /* Reset atmosphere according to chosen star */
-	  /* It is absolutely idiotic to have different 
-             atmosphere thresholds for different stars that are placed
-             on the same goddamn grid. This also screws symmetry, so 
-             we get rid of it /*
-	     /* if(rho[i3D] <= TOV_Atmosphere[star_i]) {
-	     rho[i3D] = TOV_Atmosphere[star_i];
-	     press[i3D] = TOV_K[star_i] * pow(rho[i3D],TOV_Gamma[star_i]);
-             eps[i3D] = press[i3D]/(TOV_Gamma[star_i]-1.0)/rho[i3D];
-	    } */ 
-
         }
 
-	if(rho[i3D] <= TOV_Atmosphere[0]) {
-	  rho[i3D] = TOV_Atmosphere[0];
-	  press[i3D] = TOV_K * pow(rho[i3D],TOV_Gamma);
-	  eps[i3D] = press[i3D]/(TOV_Gamma-1.0)/rho[i3D];
-	} 
-
         if (TOV_Conformal_Flat_Three_Metric)
         {
           my_psi4 -= ((TOV_Num_TOVs+TOV_Use_Old_Initial_Data-1)/my_psi4);
@@ -888,34 +839,6 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
                                 2*gxz[i3D] * velx[i3D] * velz[i3D]+
                                 2*gyz[i3D] * vely[i3D] * velz[i3D]) * my_psi4);
       }
-      /* now all combine methods are doing the same */
-      if ((TOV_Num_TOVs==1) && (rho[i3D] <= TOV_Atmosphere[0]))
-      {
-        velx[i3D] = vely[i3D] = velz[i3D] = 0.0;
-        w_lorentz[i3D] = 1.0;
-      }
-
-      SpatialDet(gxx[i3D],gxy[i3D],gxz[i3D],gyy[i3D],gyz[i3D],gzz[i3D],&det);
-
-      sqrt_det=sqrt(det) * pow(my_psi4, 1.5);
-      dens[i3D] = sqrt_det * w_lorentz[i3D] * rho[i3D];
-
-      /* this variable is only for temporal storage */
-      D_h_w = dens[i3D] * (1 + eps[i3D] + press[i3D]/rho[i3D]) * w_lorentz[i3D];
-
-      vlowx = gxx[i3D]*velx[i3D] + gxy[i3D]*vely[i3D] + gxz[i3D]*velz[i3D];
-      vlowy = gxy[i3D]*velx[i3D] + gyy[i3D]*vely[i3D] + gyz[i3D]*velz[i3D];
-      vlowz = gxz[i3D]*velx[i3D] + gyz[i3D]*vely[i3D] + gzz[i3D]*velz[i3D];
-
-      sx[i3D] = D_h_w * vlowx;
-      sy[i3D] = D_h_w * vlowy;
-      sz[i3D] = D_h_w * vlowz;
-
-      /* One could use D_h_w here, but it would introduce more error */
-      tau[i3D] = sqrt_det * (rho[i3D]  *w_lorentz[i3D]*w_lorentz[i3D] *
-                                         (1.0 + eps[i3D]) +
-                              press[i3D]*(w_lorentz[i3D]*w_lorentz[i3D]-1.0) ) -
-                 dens[i3D];
     }
 
 
@@ -976,11 +899,6 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
         TOV_Copy(i3D, velx_p_p, velx);
         TOV_Copy(i3D, vely_p_p, vely);
         TOV_Copy(i3D, velz_p_p, velz);
-        TOV_Copy(i3D, dens_p_p, dens);
-        TOV_Copy(i3D, tau_p_p,  tau);
-        TOV_Copy(i3D, sx_p_p,   sx);
-        TOV_Copy(i3D, sy_p_p,   sy);
-        TOV_Copy(i3D, sz_p_p,   sz);
         TOV_Copy(i3D, w_lorentz_p_p, w_lorentz);
     case 2:
         TOV_Copy(i3D, gxx_p,  gxx);
@@ -991,11 +909,6 @@ void TOV_C_Exact(CCTK_ARGUMENTS)
         TOV_Copy(i3D, velx_p, velx);
         TOV_Copy(i3D, vely_p, vely);
         TOV_Copy(i3D, velz_p, velz);
-        TOV_Copy(i3D, dens_p, dens);
-        TOV_Copy(i3D, tau_p,  tau);
-        TOV_Copy(i3D, sx_p,   sx);
-        TOV_Copy(i3D, sy_p,   sy);
-        TOV_Copy(i3D, sz_p,   sz);
         TOV_Copy(i3D, w_lorentz_p, w_lorentz);
   }
   CCTK_INFO("Done interpolation.");
diff --git a/src/utils.inc b/src/utils.inc
index 20c379b..a5fa350 100644
--- a/src/utils.inc
+++ b/src/utils.inc
@@ -7,7 +7,6 @@ void TOV_C_AllocateMemory(CCTK_ARGUMENTS)
   assert(TOV_Surface==0);
   assert(TOV_R_Surface==0);
   assert(TOV_RProp_Surface==0);
-  assert(TOV_Atmosphere==0);
 
   assert(TOV_r_1d==0);
   assert(TOV_rbar_1d==0);
@@ -20,7 +19,6 @@ void TOV_C_AllocateMemory(CCTK_ARGUMENTS)
   TOV_Surface = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
   TOV_R_Surface = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
   TOV_RProp_Surface = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
-  TOV_Atmosphere = malloc(TOV_Num_TOVs * sizeof(CCTK_REAL));
 
   TOV_r_1d = malloc(TOV_Num_Radial * TOV_Num_TOVs * sizeof(CCTK_REAL));
   TOV_rbar_1d = malloc(TOV_Num_Radial * TOV_Num_TOVs * sizeof(CCTK_REAL));
@@ -39,7 +37,6 @@ void TOV_C_FreeMemory (CCTK_ARGUMENTS)
   assert(TOV_Surface!=0);
   assert(TOV_R_Surface!=0);
   assert(TOV_RProp_Surface!=0);
-  assert(TOV_Atmosphere!=0);
 
   assert(TOV_r_1d!=0);
   assert(TOV_rbar_1d!=0);
@@ -50,7 +47,6 @@ void TOV_C_FreeMemory (CCTK_ARGUMENTS)
   free(TOV_Surface);
   free(TOV_R_Surface);
   free(TOV_RProp_Surface);
-  free(TOV_Atmosphere);
 
   free(TOV_r_1d);
   free(TOV_rbar_1d);
@@ -63,7 +59,6 @@ void TOV_C_FreeMemory (CCTK_ARGUMENTS)
   TOV_Surface=0;
   TOV_R_Surface=0;
   TOV_RProp_Surface=0;
-  TOV_Atmosphere=0;
 
   TOV_r_1d=0;
   TOV_rbar_1d=0;
diff --git a/test/test_two_max.par b/test/test_two_max.par
index 8aed875..b41bf3e 100644
--- a/test/test_two_max.par
+++ b/test/test_two_max.par
@@ -43,9 +43,8 @@ ADMBase::initial_shift = "tov"
 tovsolver::TOV_Num_TOVs       = 2
 tovsolver::TOV_Num_Radial     = 200000
 tovsolver::TOV_Combine_Method = "maximum"
-
-tovsolver::TOV_Gamma       = 2.0
-tovsolver::TOV_K           = 100.0
+tovsolver::TOV_Gamma          = 2.0
+tovsolver::TOV_K              = 100.0
 
 tovsolver::TOV_Rho_Central[0] = 0.16e-3
 tovsolver::TOV_Position_x[0]  = -15.0