Use SummationByParts to evaluate angular derivatives in order to

initialise the time derivatives of lapse and shift


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Meudon_Bin_NS/trunk@5 f2ea251b-07d6-4a45-8f2c-a162b3fa7596

3 files changed, 67 insertions, 13 deletions

diff --git a/interface.ccl b/interface.ccl
index c9b0207..e4ac277 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,4 +2,13 @@
 
 IMPLEMENTS: ID_Bin_NS
 
-INHERITS: grid ADMBase HydroBase
+INHERITS: grid SummationByParts ADMBase HydroBase
+
+
+
+SUBROUTINE Diff_gv (CCTK_POINTER_TO_CONST IN cctkGH, \
+                    CCTK_INT IN dir, \
+                    CCTK_REAL IN ARRAY var, \
+                    CCTK_REAL OUT ARRAY dvar, \
+                    CCTK_INT IN table_handle)
+USES FUNCTION Diff_gv
diff --git a/src/Bin_NS.cc b/src/Bin_NS.cc
index c04bf66..3714bc2 100644
--- a/src/Bin_NS.cc
+++ b/src/Bin_NS.cc
@@ -12,6 +12,30 @@ using namespace std;
 
 
 
+static void set_dt_from_domega (CCTK_ARGUMENTS,
+                                CCTK_REAL const* const var,
+                                CCTK_REAL      * const dtvar,
+                                CCTK_REAL const& omega)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  
+  int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+  vector<CCTK_REAL> dxvar(npoints), dyvar(npoints);
+  
+  Diff_gv (cctkGH, 0, var, &dxvar[0], -1);
+  Diff_gv (cctkGH, 1, var, &dyvar[0], -1);
+  
+#pragma omp parallel for
+  for (int i=0; i<npoints; ++i) {
+    CCTK_REAL const ephix = +y[i];
+    CCTK_REAL const ephiy = -x[i];
+    CCTK_REAL const dphi_var = ephix * dxvar[i] + ephiy * dyvar[i];
+    dtvar[i] = omega * dphi_var;
+  }
+}
+
+
+
 extern "C"
 void ID_Bin_NS_initialise (CCTK_ARGUMENTS)
 {
@@ -48,7 +72,6 @@ void ID_Bin_NS_initialise (CCTK_ARGUMENTS)
   CCTK_INFO ("Setting up coordinates");
   
   int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
-  
   vector<double> xx(npoints), yy(npoints), zz(npoints);
   
 #pragma omp parallel for
@@ -94,14 +117,6 @@ void ID_Bin_NS_initialise (CCTK_ARGUMENTS)
     betay[i] = bin_ns.beta_y[i];
     betaz[i] = bin_ns.beta_z[i];
     
-    // These initial data assume a helical Killing vector field
-    // TODO: calculate spatial derivatives to set this correctly
-    dtalp[i] = 0.0;
-    
-    dtbetax[i] = 0.0;
-    dtbetay[i] = 0.0;
-    dtbetaz[i] = 0.0;
-    
     CCTK_REAL g[3][3];
     g[0][0] = bin_ns.g_xx[i];
     g[0][1] = bin_ns.g_xy[i];
@@ -162,5 +177,33 @@ void ID_Bin_NS_initialise (CCTK_ARGUMENTS)
   
   
   
+  CCTK_INFO ("Calculating time derivatives of lapse and shift");
+  {
+    // Angular velocity
+    CCTK_REAL const omega = bin_ns.omega * cactusT;
+    
+    // These initial data assume a helical Killing vector field
+    
+    if (CCTK_EQUALS (initial_dtlapse, "ID_Bin_NS")) {
+      set_dt_from_domega (CCTK_PASS_CTOC, alp, dtalp, omega);
+    } else if (CCTK_EQUALS (initial_dtlapse, "none")) {
+      // do nothing
+    } else {
+      CCTK_WARN (CCTK_WARN_ABORT, "internal error");
+    }
+    
+    if (CCTK_EQUALS (initial_dtshift, "ID_Bin_NS")) {
+      set_dt_from_domega (CCTK_PASS_CTOC, betax, dtbetax, omega);
+      set_dt_from_domega (CCTK_PASS_CTOC, betay, dtbetay, omega);
+      set_dt_from_domega (CCTK_PASS_CTOC, betaz, dtbetaz, omega);
+    } else if (CCTK_EQUALS (initial_dtshift, "none")) {
+      // do nothing
+    } else {
+      CCTK_WARN (CCTK_WARN_ABORT, "internal error");
+    }
+  }
+  
+  
+  
   CCTK_INFO ("Done.");
 }
diff --git a/src/check_parameters.cc b/src/check_parameters.cc
index 40a760f..d96379d 100644
--- a/src/check_parameters.cc
+++ b/src/check_parameters.cc
@@ -13,10 +13,12 @@ void ID_Bin_NS_check_parameters (CCTK_ARGUMENTS)
   if (not CCTK_EQUALS (initial_data,    "ID_Bin_NS") or
       not CCTK_EQUALS (initial_lapse,   "ID_Bin_NS") or
       not CCTK_EQUALS (initial_shift,   "ID_Bin_NS") or
-      not CCTK_EQUALS (initial_dtlapse, "ID_Bin_NS") or
-      not CCTK_EQUALS (initial_dtshift, "ID_Bin_NS") or
+      not (CCTK_EQUALS (initial_dtlapse, "ID_Bin_NS") or
+           CCTK_EQUALS (initial_dtlapse, "none")) or
+      not (CCTK_EQUALS (initial_dtshift, "ID_Bin_NS") or
+           CCTK_EQUALS (initial_dtshift, "none")) or
       not CCTK_EQUALS (initial_hydro,   "ID_Bin_NS"))
   {
-    CCTK_PARAMWARN ("The parameters ADMBase::initial_data, ADMBase::initial_lapse, ADMBase::initial_shift, ADMBase::initial_dtlapse, ADMBase::initial_dtshift, and HydroBase::initial_hydro must all be set to the value \"ID_Bin_NS\"");
+    CCTK_PARAMWARN ("The parameters ADMBase::initial_data, ADMBase::initial_lapse, ADMBase::initial_shift, ADMBase::initial_dtlapse, ADMBase::initial_dtshift, and HydroBase::initial_hydro must all be set to the value \"ID_Bin_NS\" or \"none\"");
   }
 }