Add a parameter "initial_hydro" to HydroBase with a default setting of

"zero".  This parameter, which can be extended by hydro initial data
thorns, determines which hydro thorn is setting up initial data.

The default setting "zero" is implemented by HydroBase itself, setting
all hydro variables to zero.  This corresponds to a vacuum, but
without atmosphere.

Reformat and/or clarify some comments in the ccl files.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinBase/HydroBase/trunk@17 57fe0bb3-ccba-405f-9b23-de0201f165b7

4 files changed, 188 insertions, 78 deletions

diff --git a/interface.ccl b/interface.ccl
index 84253b3..164a8ff 100755
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,14 +2,17 @@
 
 implements: HydroBase
 
+inherits: InitBase
+
 public:
 
 # These variables correspond to the Valencia formulation
 
-CCTK_REAL rho   type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "rest mass density"
-CCTK_REAL press type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "gas pressure"
-CCTK_REAL eps   type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "internal energy"
+CCTK_REAL rho     type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "rest mass density"
+CCTK_REAL press   type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "gas pressure"
+CCTK_REAL eps     type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "internal energy density"
+
+CCTK_REAL vel[3]  type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" interpolator="matter"' "velocity v^i"
 
-CCTK_REAL vel[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" interpolator="matter"' "velocity"
+CCTK_REAL Y_e     type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "electron fraction"
 
-CCTK_REAL Y_e    type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "electron fraction"
diff --git a/param.ccl b/param.ccl
index 3ccff24..5aca620 100755
--- a/param.ccl
+++ b/param.ccl
@@ -1,6 +1,15 @@
 # Parameter definitions for thorn HydroBase
 
-restricted:
+SHARES: InitBase
+
+USES KEYWORD initial_data_setup_method
+
+RESTRICTED:
+
+KEYWORD initial_hydro "The hydro initial data"
+{
+  "zero" :: "hydro variables are set to vacuum (without atmosphere)"
+} "zero"
 
 KEYWORD evolution_method "The hydro evolution method"
 {
@@ -9,7 +18,7 @@ KEYWORD evolution_method "The hydro evolution method"
 
 INT timelevels "Number of time levels in evolution scheme"
 {
-  1:3 :: "FMR/AMR runs probably need 3. Otherwise use 2"
+  1:3 :: ""
 } 1
 
 STRING prolongation_type "The prolongation operator used by Carpet for HydroBase variables"
@@ -19,6 +28,8 @@ STRING prolongation_type "The prolongation operator used by Carpet for HydroBase
   ".*"   :: "Anything else"
 } "ENO"
 
+# Settings for electron fraction Y_e
+
 KEYWORD initial_Y_e "Initial value for Y_e"
 {
   "none" :: "inactive"
diff --git a/schedule.ccl b/schedule.ccl
index 539efe5..57de555 100755
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,105 +1,119 @@
 # Schedule definitions for thorn HydroBase
 
-#######################################################################
-### We leave the storage on all the time for the variables required ###
-### by MoL. There is probably a better way of doing this.           ###
-#######################################################################
-
-if (timelevels == 3)
+if (timelevels == 1)
 {
-  STORAGE:rho[3]
-  STORAGE:press[3]
-  STORAGE:eps[3]
-  STORAGE:vel[3]
+  STORAGE: rho[1]
+  STORAGE: press[1]
+  STORAGE: eps[1]
+  STORAGE: vel[1]
   if (!CCTK_EQUALS(initial_Y_e, "none"))
   {
-    STORAGE:Y_e[3]
+    STORAGE: Y_e[1]
   }
 }
 else if (timelevels == 2)
 {
-  STORAGE:rho[2]
-  STORAGE:press[2]
-  STORAGE:eps[2]
-  STORAGE:vel[2]
+  STORAGE: rho[2]
+  STORAGE: press[2]
+  STORAGE: eps[2]
+  STORAGE: vel[2]
   if (!CCTK_EQUALS(initial_Y_e, "none"))
   {
-    STORAGE:Y_e[2]
+    STORAGE: Y_e[2]
   }
 }
-else if (timelevels == 1)
+else if (timelevels == 3)
 {
-  STORAGE:rho[1]
-  STORAGE:press[1]
-  STORAGE:eps[1]
-  STORAGE:vel[1]
+  STORAGE: rho[3]
+  STORAGE: press[3]
+  STORAGE: eps[3]
+  STORAGE: vel[3]
   if (!CCTK_EQUALS(initial_Y_e, "none"))
   {
-    STORAGE:Y_e[1]
+    STORAGE: Y_e[3]
   }
 }
 
-schedule group HydroBase_Initial AT Initial after ADMBase_InitialData after ADMBase_InitialGauge after IOUtil_RecoverIDFromDatafiles before ADMBase_PostInitial before SetTmunu
+schedule group HydroBase_Initial at CCTK_INITIAL    \
+    after (ADMBase_InitialData ADMBase_InitialGauge \
+           IOUtil_RecoverIDFromDatafiles)           \
+    before (ADMBase_PostInitial SetTmunu)
 {
 } "HydroBase initial data group"
 
-if(CCTK_Equals(initial_Y_e, "one"))
-{
-  SCHEDULE HydroBase_Y_e_one in HydroBase_Initial
-  {
-    LANG: C
-  } "Set electron fraction to 1 at all points"
-}
-
-schedule HydroBase_StartUp AT WRAGH
+schedule HydroBase_StartUp at CCTK_STARTUP
 {
   LANG: C
 } "Startup banner"
 
-schedule group HydroBase_RHS IN MoL_CalcRHS
+schedule group HydroBase_RHS in MoL_CalcRHS
 {
 } "Groups for scheduling tasks for calculating RHS of hydro variables"
 
-schedule group HydroBase_PostStep IN MoL_PostStep BEFORE SetTmunu
+schedule group HydroBase_PostStep in MoL_PostStep before SetTmunu
 {
 } "Post step tasks for hydro thorns"
 
-schedule group HydroBase_PostStep AT POSTRESTRICTINITIAL
+schedule group HydroBase_PostStep at CCTK_POSTRESTRICTINITIAL
 {
 } "Post step tasks for hydro thorns"
 
-# Hydro thorns should not schedule something in group HydroBase_Boundaries
-# (schedule them in HydroBase_Select_Boundaries), but they can refer to it
-# in order to schedule own routines before or after the boundary
-# treatment or they can schedule the whole group
-schedule group HydroBase_Boundaries IN HydroBase_PostStep BEFORE HydroBase_Con2Prim
+# Hydro thorns should not schedule anything in group
+# HydroBase_Boundaries; they should schedule in
+# HydroBase_Select_Boundaries instead.  However, they may refer to
+# HydroBase_Boundaries in order to schedule own routines before or
+# after the boundary treatment, or in order to schedule the whole
+# group.
+schedule group HydroBase_Boundaries in HydroBase_PostStep \
+    before HydroBase_Con2Prim
 {
 } "HydroBase-internal Boundary conditions group"
 
-# Here codes would have to schedule functions to specify which boundary conditions
-# should be applied to which variables
-schedule group HydroBase_Select_Boundaries IN HydroBase_Boundaries
+# This is where hydro thorns schedule functions that select which
+# boundary conditions should be applied to which variables.
+schedule group HydroBase_Select_Boundaries in HydroBase_Boundaries
 {
 } "Group to schedule the boundary condition functions"
 
-# This will then later actually apply all the selected boundary conditions. Thorns
-# using HydroBase should not need to change something here.
-schedule group ApplyBCs AS HydroBase_ApplyBCs IN HydroBase_Boundaries AFTER HydroBase_Select_Boundaries
+# This group actually applies all selected boundary conditions.  Hydro
+# thorns should not need to change anything here.
+schedule group ApplyBCs as HydroBase_ApplyBCs in HydroBase_Boundaries \
+    after HydroBase_Select_Boundaries
 {
 } "Apply the boundary conditions of HydroBase"
 
-# This is a group to schedule Con2Prim routines, scheduled at different times during
-# an evolution
-schedule group HydroBase_Con2Prim IN HydroBase_PostStep
+# Hydro thorns should schedule their Con2Prim routines in this group.
+# This group is then scheduled at different times during an evolution.
+schedule group HydroBase_Con2Prim in HydroBase_PostStep
 {
 } "Convert from conservative to primitive variables"
 
-schedule group HydroBase_Con2Prim AT CCTK_PostPostInitial AS Con2Prim BEFORE ADMConstraintsGroup
+schedule group HydroBase_Con2Prim at CCTK_PostPostInitial as Con2Prim \
+    before ADMConstraintsGroup
 {
 } "Convert from conservative to primitive variables (might be redundant)"
 
-# Schedule prim2con after initial data
-schedule group HydroBase_Prim2ConInitial AT Initial AFTER HydroBase_Initial
+# Hydro thorns should schedule their Prim2Con routines in this group.
+# This group is then scheduled at different times during an evolution.
+schedule group HydroBase_Prim2ConInitial at CCTK_INITIAL \
+    after HydroBase_Initial
 {
 } "Recover the conservative variables from the primitive variables"
 
+
+
+if (CCTK_EQUALS(initial_hydro, "zero"))
+{
+  schedule HydroBase_Zero in HydroBase_Initial
+  {
+    LANG: C
+  } "Set up zero (vacuum without atmosphere) hydro initial data"
+}
+
+if (CCTK_Equals(initial_Y_e, "one"))
+{
+  SCHEDULE HydroBase_Y_e_one in HydroBase_Initial
+  {
+    LANG: C
+  } "Set electron fraction to 1"
+}
diff --git a/src/Initialisation.c b/src/Initialisation.c
index de2411c..a739091 100644
--- a/src/Initialisation.c
+++ b/src/Initialisation.c
@@ -1,31 +1,113 @@
-#include "cctk.h"
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
 
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
 
-void HydroBase_Y_e_one(CCTK_ARGUMENTS)
+
+void HydroBase_Zero (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
-  int i;
-
-  for(i = 0; i < cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; i++)
-    Y_e[i] = 1.0;
   
-  if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::Y_e") > 1)
-    for(i = 0; i < cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; i++)
-      Y_e_p[i] = 1.0;
+  int const np = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+  
+#pragma omp parallel for
+  for (int i=0; i<np; ++i) {
+    rho[i     ] = 0.0;
+    vel[i     ] = 0.0;
+    vel[i+  np] = 0.0;
+    vel[i+2*np] = 0.0;
+    eps[i     ] = 0.0;
+  }
   
-  if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::Y_e") > 2)
-    for(i = 0; i < cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; i++)
-      Y_e_p_p[i] = 1.0;
+  if (CCTK_EQUALS (initial_data_setup_method, "init_some_levels") ||
+      CCTK_EQUALS (initial_data_setup_method, "init_single_levels"))
+  {
+    /* do nothing */
+  }
+  else if (CCTK_EQUALS (initial_data_setup_method, "init_all_levels"))
+  {
+    
+    if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::rho") >= 2) {
+#pragma omp parallel for
+      for (int i=0; i<np; ++i) {
+        rho_p[i     ] = 0.0;
+        vel_p[i     ] = 0.0;
+        vel_p[i+  np] = 0.0;
+        vel_p[i+2*np] = 0.0;
+        eps_p[i     ] = 0.0;
+      }
+    }
+    
+    if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::rho") >= 3) {
+#pragma omp parallel for
+      for (int i=0; i<np; ++i) {
+        rho_p_p[i     ] = 0.0;
+        vel_p_p[i     ] = 0.0;
+        vel_p_p[i+  np] = 0.0;
+        vel_p_p[i+2*np] = 0.0;
+        eps_p_p[i     ] = 0.0;
+      }
+    }
+    
+    if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::rho") >= 4) {
+      CCTK_WARN (CCTK_WARN_ABORT,
+                 "Too many active time levels for HydroBase variables");
+    }
+    
+  }
+  else
+  {
+    CCTK_WARN (CCTK_WARN_ABORT,
+               "Unsupported parameter value for InitBase::initial_data_setup_method");
+  }
 }
 
-/*void HydroBase_SetY_eStateOn(CCTK_ARGUMENTS)
+
+
+void HydroBase_Y_e_one (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
-
-  *Y_e_state = 1;
-}*/
-
+  
+  int const np = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+  
+#pragma omp parallel for
+  for (int i=0; i<np; ++i) {
+    Y_e[i] = 1.0;
+  }
+  
+  if (CCTK_EQUALS (initial_data_setup_method, "init_some_levels") ||
+      CCTK_EQUALS (initial_data_setup_method, "init_single_levels"))
+  {
+    /* do nothing */
+  }
+  else if (CCTK_EQUALS (initial_data_setup_method, "init_all_levels"))
+  {
+    
+    if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::Y_e") >= 2) {
+#pragma omp parallel for
+      for (int i=0; i<np; ++i) {
+        Y_e_p[i] = 1.0;
+      }
+    }
+    
+    if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::Y_e") >= 3) {
+#pragma omp parallel for
+      for (int i=0; i<np; ++i) {
+        Y_e_p_p[i] = 1.0;
+      }
+    }
+    
+    if (CCTK_ActiveTimeLevels(cctkGH, "HydroBase::Y_e") >= 4) {
+      CCTK_WARN (CCTK_WARN_ABORT,
+                 "Too many active time levels for HydroBase::Y_e");
+    }
+    
+  }
+  else
+  {
+    CCTK_WARN (CCTK_WARN_ABORT,
+               "Unsupported parameter value for InitBase::initial_data_setup_method");
+  }
+}