Calculate the weight function correctly.

darcs-hash:20040802094335-07bb3-8cf89dfef6eab7c25a50842c68b62d8a7c230d57.gz

5 files changed, 66 insertions, 28 deletions

diff --git a/Carpet/CarpetReduce/schedule.ccl b/Carpet/CarpetReduce/schedule.ccl
index 048fa5a62..db1530ae7 100644
--- a/Carpet/CarpetReduce/schedule.ccl
+++ b/Carpet/CarpetReduce/schedule.ccl
@@ -1,5 +1,5 @@
 # Schedule definitions for thorn CarpetReduce
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/schedule.ccl,v 1.3 2004/06/14 07:01:21 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/schedule.ccl,v 1.4 2004/08/02 11:43:35 schnetter Exp $
 
 schedule CarpetReduceStartup at STARTUP
 {
@@ -22,6 +22,7 @@ SCHEDULE GROUP MaskBase_SetupMask AT postregrid
 SCHEDULE MaskBase_InitMask IN MaskBase_SetupMask
 {
   LANG: C
+  OPTIONS: global loop-local
 } "Initialise the weight function"
 
 SCHEDULE GROUP SetupMask IN MaskBase_SetupMask AFTER MaskBase_InitMask
@@ -32,11 +33,12 @@ SCHEDULE GROUP SetupMask IN MaskBase_SetupMask AFTER MaskBase_InitMask
 SCHEDULE CoordBase_SetupMask IN SetupMask
 {
   LANG: C
+  OPTIONS: global loop-local
 } "Set up the outer boundaries of the weight function"
 
 # This might move to CarpetMask
 SCHEDULE CarpetMaskSetup IN SetupMask
 {
   LANG: C
-  OPTIONS: singlemap
+  OPTIONS: global loop-singlemap
 } "Set up the weight function for the restriction regions"
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc
index 383d3cb6e..2e6f558ae 100644
--- a/Carpet/CarpetReduce/src/mask_carpet.cc
+++ b/Carpet/CarpetReduce/src/mask_carpet.cc
@@ -10,7 +10,7 @@
 #include "mask_carpet.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_carpet.cc,v 1.2 2004/06/15 16:45:29 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_carpet.cc,v 1.3 2004/08/02 11:43:35 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetMask_Mask_cc);
 }
 
@@ -45,6 +45,7 @@ namespace CarpetMask {
     if (reflevel > 0) {
       
       ivect const izero = ivect(0);
+      ivect const ione  = ivect(1);
       
       gh<dim> const & hh = *vhh.at(Carpet::map);
       dh<dim> const & dd = *vdd.at(Carpet::map);
@@ -90,6 +91,12 @@ namespace CarpetMask {
         boundaries[d].normalize();
       }
       
+      // Subtract the boundaries from the refined region
+      for (int d=0; d<dim; ++d) {
+        refined -= boundaries[d];
+      }
+      refined.normalize();
+      
       
       
       // Set prolongation boundaries of this level
@@ -116,8 +123,8 @@ namespace CarpetMask {
               
               if (verbose) {
                 ostringstream buf;
-                buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax;
-                CCTK_VInfo (CCTK_THORNSTRING, buf.str().c_str());
+                buf << "Setting prolongation boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione;
+                CCTK_INFO (buf.str().c_str());
               }
               
               // Set weight on the boundary to 1/2
@@ -169,8 +176,8 @@ namespace CarpetMask {
             
             if (verbose) {
               ostringstream buf;
-              buf << "Setting restricted region on level " << reflevel << " to weight 0: " << imin << ":" << imax;
-              CCTK_VInfo (CCTK_THORNSTRING, buf.str().c_str());
+              buf << "Setting restricted region on level " << reflevel << " to weight 0: " << imin << ":" << imax-ione;
+              CCTK_INFO (buf.str().c_str());
             }
             
             // Set weight in the restricted region to 0
@@ -201,8 +208,8 @@ namespace CarpetMask {
               
               if (verbose) {
                 ostringstream buf;
-                buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax;
-                CCTK_VInfo (CCTK_THORNSTRING, buf.str().c_str());
+                buf << "Setting restriction boundary on level " << reflevel << " direction " << d << " to weight 1/2: " << imin << ":" << imax-ione;
+                CCTK_INFO (buf.str().c_str());
               }
               
               // Set weight on the boundary to 1/2
diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c
index 0f54bc6a3..90084d7f5 100644
--- a/Carpet/CarpetReduce/src/mask_coords.c
+++ b/Carpet/CarpetReduce/src/mask_coords.c
@@ -1,7 +1,8 @@
-/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_coords.c,v 1.1 2004/06/14 07:01:21 schnetter Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_coords.c,v 1.2 2004/08/02 11:43:35 schnetter Exp $ */
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
 
 
 
@@ -9,6 +10,7 @@ void
 CoordBase_SetupMask (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
   
   CCTK_INT nboundaryzones[6];
   CCTK_INT is_internal[6];
@@ -69,6 +71,10 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
           }
           
           /* Loop over the boundary */
+          if (verbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Setting symmetry boundary in direction %d face %d to weight 0", d, f);
+          }
           for (k=imin[2]; k<imax[2]; ++k) {
             for (j=imin[1]; j<imax[1]; ++j) {
               for (i=imin[0]; i<imax[0]; ++i) {
@@ -95,17 +101,21 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
               imin[d] = imax[d] - 1;
             }
           
-          /* Loop over the points next to boundary */
-          for (k=imin[2]; k<imax[2]; ++k) {
-            for (j=imin[1]; j<imax[1]; ++j) {
-              for (i=imin[0]; i<imax[0]; ++i) {
-                
-                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
-                weight[ind] *= 0.5;
-                
+            /* Loop over the points next to boundary */
+            if (verbose) {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "Setting symmetry axis in direction %d face %d to weight 1/2", d, f);
+            }
+            for (k=imin[2]; k<imax[2]; ++k) {
+              for (j=imin[1]; j<imax[1]; ++j) {
+                for (i=imin[0]; i<imax[0]; ++i) {
+                  
+                  int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                  weight[ind] *= 0.5;
+                  
+                }
               }
             }
-          }
             
           } /* if the boundary is not staggered */
           
diff --git a/Carpet/CarpetReduce/src/mask_init.c b/Carpet/CarpetReduce/src/mask_init.c
index cd3d0a965..70768b973 100644
--- a/Carpet/CarpetReduce/src/mask_init.c
+++ b/Carpet/CarpetReduce/src/mask_init.c
@@ -1,7 +1,9 @@
-/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_init.c,v 1.1 2004/06/14 07:01:21 schnetter Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_init.c,v 1.2 2004/08/02 11:43:35 schnetter Exp $ */
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
 #include "util_Table.h"
 
 
@@ -10,6 +12,7 @@ void
 MaskBase_InitMask (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
   
 #if 0
   CCTK_INT symtable;
@@ -30,6 +33,9 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
   
   
   /* Initialise the weight to 1 everywhere */
+  if (verbose) {
+    CCTK_INFO ("Initialising to weight 1");
+  }
   
   for (k=0; k<cctk_lsh[2]; ++k) {
     for (j=0; j<cctk_lsh[1]; ++j) {
@@ -45,6 +51,9 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
   
   
   /* Set the weight to 0 on inter-processor boundaries */
+  if (verbose) {
+    CCTK_INFO ("Setting inter-processor boundaries to weight 0");
+  }
   
   /* Loop over all dimensions and faces */
   for (d=0; d<3; ++d) {
@@ -86,6 +95,9 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
 #if 0
   
   /* Take the symmetry boundaries into account */
+  if (verbose) {
+    CCTK_INFO ("Setting symmetry boundaries to weight 0");
+  }
   
   /* Get symmetry information */
   symtable = SymmetryTableHandleForGrid (cctkGH);
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index 9d6a8453a..622d63291 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.42 2004/07/06 14:44:32 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.43 2004/08/02 11:43:35 schnetter Exp $
 
 #include <assert.h>
 #include <float.h>
@@ -23,7 +23,7 @@
 #include "reduce.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.42 2004/07/06 14:44:32 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.43 2004/08/02 11:43:35 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetReduce_reduce_cc);
 }
 
@@ -954,11 +954,18 @@ namespace CarpetReduce {
         if (need_time_interp) {
           
           if (CCTK_ActiveTimeLevelsVI(cgh, vi) < num_tl) {
-            char * const fullname = CCTK_FullName(vi);
-            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                        "Grid function \"%s\" has only %d active time levels on refinement level %d; this is not enough for time interpolation.  Using the current time level instead",
-                        fullname, CCTK_ActiveTimeLevelsVI(cgh, vi), reflevel);
-            free (fullname);
+            static vector<bool> have_warned;
+            if (have_warned.empty()) {
+              have_warned.resize (CCTK_NumVars(), false);
+            }
+            if (! have_warned.at(vi)) {
+              char * const fullname = CCTK_FullName(vi);
+              CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                          "Grid function \"%s\" has only %d active time levels on refinement level %d; this is not enough for time interpolation.  Using the current time level instead",
+                          fullname, CCTK_ActiveTimeLevelsVI(cgh, vi), reflevel);
+              free (fullname);
+              have_warned.at(vi) = true;
+            }
             
             // fall back
             need_time_interp = false;
@@ -1052,7 +1059,7 @@ namespace CarpetReduce {
               weight = (static_cast<CCTK_REAL const *>
                         (CCTK_VarDataPtr (cgh, 0, "CarpetReduce::weight")));
               assert (weight);
-              levfac = pow(1.0*reflevelfact, -grpdim);
+              levfac = pow(CCTK_REAL(reflevelfact), -grpdim);
             } else {
               weight = NULL;
               levfac = 1.0;