[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 205)

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2013/03/27 08:49 PM

Modified:
 /trunk/
  param.ccl, schedule.ccl
 /trunk/src/
  GRHydro_BondiM.c

Log:
 GRHydro_Init_Data: add routine to force boundary data to Bondi
 
 From: Roland Haas <rhaas at tapir.caltech.edu>

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_BondiM.c
Delta lines: +46 -8
===================================================================
--- trunk/src/GRHydro_BondiM.c	2013-03-28 01:49:03 UTC (rev 204)
+++ trunk/src/GRHydro_BondiM.c	2013-03-28 01:49:05 UTC (rev 205)
@@ -890,22 +890,41 @@
    -- driver routine for the Bondi solution;
 
 ***********************************************************************************/
-static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max);
+static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, 
+  CCTK_REAL range_min, CCTK_REAL range_max, 
+  const int range_imin[], const int range_imax[]);
 
 void GRHydro_BondiM(CCTK_ARGUMENTS)
 {
-  GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100);
+  DECLARE_CCTK_ARGUMENTS;
+  const int null[3] = {0,0,0};
+  GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100, cctk_lsh, null);
 }
 
 void GRHydro_BondiM_Range(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_PARAMETERS;
+  DECLARE_CCTK_ARGUMENTS;
+  const int null[3] = {0,0,0};
   GRHydro_BondiM_Internal(CCTK_PASS_CTOC, bondi_freeze_inner_radius,
-                          bondi_freeze_outer_radius);
+                          bondi_freeze_outer_radius,
+                          cctk_lsh, null);
 }
 
+void GRHydro_BondiM_Boundary(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_PARAMETERS;
+  DECLARE_CCTK_ARGUMENTS;
+  const int imax[3] = {cctk_lsh[0]-cctk_nghostzones[0],
+                       cctk_lsh[1]-cctk_nghostzones[1],
+                       cctk_lsh[2]-cctk_nghostzones[2]};
+  GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100, 
+                          cctk_nghostzones, imax);
+}
 
-static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max)
+static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, 
+  CCTK_REAL range_min, CCTK_REAL range_max,
+  const int range_imin[], const int range_imax[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
@@ -920,7 +939,7 @@
 
   int nbondi; 
   
-  const int size = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+  const int size = cctk_ash[0] * cctk_ash[1] * cctk_ash[2];
   int i, imin, j;
 
 #    define velx(i_)  vel[i_ + 0*size]
@@ -952,7 +971,7 @@
   } else if( CCTK_EQUALS(bondi_Bvec_method, "transform")) {
     bvec_method = BVEC_TRANSFORM;
   } else {
-    CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+    CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING,
                "Unknown bvec setup method '%s'", bondi_Bvec_method);
   }
 
@@ -1109,8 +1128,12 @@
   *********************************************************************************/
   xpos[TT] = 0.;
   
-  for(i=0; i < size; i++) {
+  for(int kk = 0 ; kk < cctk_lsh[2] ; kk++) {
+  for(int jj = 0 ; jj < cctk_lsh[1] ; jj++) {
+  for(int ii = 0 ; ii < cctk_lsh[0] ; ii++) {
 
+    i = CCTK_GFINDEX3D(cctkGH, ii,jj,kk);
+
     xpos[XX] = x[i] ;
     xpos[YY] = y[i] ;
     xpos[ZZ] = z[i] ;
@@ -1152,8 +1175,21 @@
 
     rspher = x_spher[RR]; 
 
+    // the conditions here are a bit strange. The logic goes like this:
+    // range_min,range_max (and range_imin, range_imax) are radii (and indices)
+    // bounding the volume where we want to evolve the fields. So during
+    // evolution we only reset to ID data if we are NOT inside this volume. For
+    // ID we pass special values for min/max such that these conditions never
+    // trigger. 
+    // TODO: do this more nicely. The fact that a comment is needed should tell
+    // me that this is getting too complex.
     if(rspher > range_min && rspher < range_max)
       continue;
+    if(ii >= range_imin[0] && ii < range_imax[0] &&
+       jj >= range_imin[1] && jj < range_imax[1] &&
+       kk >= range_imin[2] && kk < range_imax[2]) {
+      continue;
+    }
 
     /* Find nearest point in the Bondi solution : */
     j = (int)  ( 0.5 +  (log10(rspher) - logrmin) / dlogr ) ; 
@@ -1266,7 +1302,9 @@
     if(psidc)
       psidc[i] = 0.;
                  
-  }
+  } // ii
+  } // jj
+  } // kk
 
   free(r_bondi);      free(logr_bondi);      free(rho_bondi);     free(u_bondi);   free(v_bondi);
 

Directory: /trunk/
==================

File [modified]: param.ccl
Delta lines: +4 -0
===================================================================
--- trunk/param.ccl	2013-03-28 01:49:03 UTC (rev 204)
+++ trunk/param.ccl	2013-03-28 01:49:05 UTC (rev 205)
@@ -420,6 +420,10 @@
   *:* :: "any value"
 } 1e300
 
+CCTK_BOOLEAN bondi_overwrite_boundary "reset data to initial data in outer boundary in boundary condition"
+{
+} "no"
+
 # For Poloidal Magnetic field test:
 
 CCTK_REAL poloidal_A_b  "Vector potential strength"

File [modified]: schedule.ccl
Delta lines: +10 -0
===================================================================
--- trunk/schedule.ccl	2013-03-28 01:49:03 UTC (rev 204)
+++ trunk/schedule.ccl	2013-03-28 01:49:05 UTC (rev 205)
@@ -257,6 +257,16 @@
 } "force analytic solution outside anulus"
 }
 
+if(bondi_overwrite_boundary ) {
+# this is a terrible HACK. We really should at least schedule ourselves in
+# Boundaries and look at cctk_bbox and the symmetry flags
+# by now Cactus is convenient CCTK_LOOP macros for this...
+SCHEDULE GRHydro_BondiM_Boundary IN HydroBase_Boundaries BEFORE HydroBase_Select_Boundaries
+{
+  LANG: C
+} "force analytic solution in boundaries"
+}
+
 if (CCTK_EQUALS(initial_Bvec, "poloidalmagfield") || CCTK_EQUALS(initial_Avec, "poloidalmagfield"))
 {
 #  SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_INITIAL AFTER IN HydroBase_Initial AFTER rnsid_init AFTER TOV_Initial_Data after CCCC_StarMapper_InitialData