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_BH/trunk@5 0994f19a-7e1f-4f5a-9787-1e7a3bb7f33f

3 files changed, 68 insertions, 13 deletions

diff --git a/interface.ccl b/interface.ccl
index cb32a3c..d74bcc3 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,4 +2,13 @@
 
 IMPLEMENTS: ID_Bin_BH
 
-INHERITS: grid ADMBase
+INHERITS: grid SummationByParts ADMBase
+
+
+
+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_BH.cc b/src/Bin_BH.cc
index 87b0a7d..e501a38 100644
--- a/src/Bin_BH.cc
+++ b/src/Bin_BH.cc
@@ -1,4 +1,5 @@
 #include <cassert>
+#include <cmath>
 #include <cstdio>
 #include <vector>
 
@@ -12,6 +13,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_BH_initialise (CCTK_ARGUMENTS)
 {
@@ -25,7 +50,6 @@ void ID_Bin_BH_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
@@ -60,14 +84,6 @@ void ID_Bin_BH_initialise (CCTK_ARGUMENTS)
     betay[i] = bin_bh.beta_y[i];
     betaz[i] = bin_bh.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_bh.g_xx[i];
     g[0][1] = bin_bh.g_xy[i];
@@ -120,5 +136,33 @@ void ID_Bin_BH_initialise (CCTK_ARGUMENTS)
   
   
   
+  CCTK_INFO ("Calculating time derivatives of lapse and shift");
+  {
+    // Angular velocity
+    CCTK_REAL const omega = bin_bh.omega * bin_bh.dist;
+    
+    // These initial data assume a helical Killing vector field
+    
+    if (CCTK_EQUALS (initial_dtlapse, "ID_Bin_BH")) {
+      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_BH")) {
+      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 f3b626f..2775824 100644
--- a/src/check_parameters.cc
+++ b/src/check_parameters.cc
@@ -13,9 +13,11 @@ void ID_Bin_BH_check_parameters (CCTK_ARGUMENTS)
   if (not CCTK_EQUALS (initial_data,    "ID_Bin_BH") or
       not CCTK_EQUALS (initial_lapse,   "ID_Bin_BH") or
       not CCTK_EQUALS (initial_shift,   "ID_Bin_BH") or
-      not CCTK_EQUALS (initial_dtlapse, "ID_Bin_BH") or
-      not CCTK_EQUALS (initial_dtshift, "ID_Bin_BH"))
+      not (CCTK_EQUALS (initial_dtlapse, "ID_Bin_BH") or
+           CCTK_EQUALS (initial_dtlapse, "none")) or
+      not (CCTK_EQUALS (initial_dtshift, "ID_Bin_BH") or
+           CCTK_EQUALS (initial_dtshift, "none")))
   {
-    CCTK_PARAMWARN ("The parameters ADMBase::initial_data, ADMBase::initial_lapse, ADMBase::initial_shift , ADMBase::initial_dtlapse, and ADMBase::initial_dtshift must all be set to the value \"ID_Bin_BH\"");
+    CCTK_PARAMWARN ("The parameters ADMBase::initial_data, ADMBase::initial_lapse, ADMBase::initial_shift, ADMBase::initial_dtlapse, and ADMBase::initial_dtshift must all be set to the value \"ID_Bin_BH\" or \"none\"");
   }
 }