Implement a moving excision region.

darcs-hash:20040414201944-07bb3-3abc45a6f1034c083fa413d0b917d0f2791a787c.gz

5 files changed, 111 insertions, 25 deletions

diff --git a/Carpet/CarpetRegrid/param.ccl b/Carpet/CarpetRegrid/param.ccl
index eee4144b0..ca830a75f 100644
--- a/Carpet/CarpetRegrid/param.ccl
+++ b/Carpet/CarpetRegrid/param.ccl
@@ -1,5 +1,5 @@
 # Parameter definitions for thorn CarpetRegrid
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.16 2004/02/05 16:11:20 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.17 2004/04/14 22:19:44 schnetter Exp $
 
 
 
@@ -46,6 +46,7 @@ KEYWORD refined_regions "Regions where the grid is refined" STEERABLE=always
   "manual-coordinates"     :: "Refine the regions specified by coordinates l[123][xyz]{min,max}"
   "manual-gridpoint-list"  :: "Refine the regions specified by integer grid points in the parameter 'gridpoints'"
   "manual-coordinate-list" :: "Refine the regions specified by coordinates in the parameter 'coordinates'"
+  "moving"                 :: "Refine a moving region"
   "automatic"              :: "Refine automatically"
 } "centre"
 
@@ -65,6 +66,46 @@ BOOLEAN symmetry_z "Refine the lower half in z-direction"
 
 
 
+# Region specifications for moving boxes
+
+KEYWORD moving_trajectory "Type of trajectory" STEERABLE=always
+{
+  "point"  :: "Do not move"
+  "circle" :: "Move in a circle"
+} "point"
+
+CCTK_REAL moving_region_radius "Radius of the moving region" STEERABLE=always
+{
+  (0: :: ""
+} 1.0
+
+CCTK_REAL moving_centre_x "x-coordinate of the centre" STEERABLE=always
+{
+  : :: ""
+} 0.0
+
+CCTK_REAL moving_centre_y "y-coordinate of the centre" STEERABLE=always
+{
+  : :: ""
+} 0.0
+
+CCTK_REAL moving_centre_z "z-coordinate of the centre" STEERABLE=always
+{
+  : :: ""
+} 0.0
+
+CCTK_REAL moving_circle_radius "Radius of the circle" STEERABLE=always
+{
+  0: :: ""
+} 1.0
+
+CCTK_REAL moving_circle_frequency "Angular frequency on the circle" STEERABLE=always
+{
+  0: :: ""
+} 1.0
+
+
+
 # Region specifications for manual gridpoint refinement
 
 CCTK_INT l1ixmin "Lower boundary of level 1 box in x-direction" STEERABLE=always
diff --git a/Carpet/CarpetRegrid/src/make.code.defn b/Carpet/CarpetRegrid/src/make.code.defn
index 0b692a316..e48f9956e 100644
--- a/Carpet/CarpetRegrid/src/make.code.defn
+++ b/Carpet/CarpetRegrid/src/make.code.defn
@@ -1,5 +1,5 @@
 # Main make.code.defn file for thorn CarpetRegrid
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/make.code.defn,v 1.3 2004/01/25 14:57:30 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/make.code.defn,v 1.4 2004/04/14 22:19:44 schnetter Exp $
 
 # Source files in this directory
 SRCS =  automatic.cc            \
@@ -9,6 +9,7 @@ SRCS =  automatic.cc            \
         manualcoordinates.cc    \
         manualgridpointlist.cc  \
         manualgridpoints.cc     \
+        moving.cc               \
         regrid.cc               \
         paramcheck.cc
 
diff --git a/Carpet/CarpetRegrid/src/moving.cc b/Carpet/CarpetRegrid/src/moving.cc
index 5e5287b86..c4bf3abc6 100644
--- a/Carpet/CarpetRegrid/src/moving.cc
+++ b/Carpet/CarpetRegrid/src/moving.cc
@@ -1,5 +1,4 @@
-#include <cassert>
-#include <cmath>
+#include <assert.h>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -10,7 +9,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/moving.cc,v 1.5 2004/08/02 11:42:20 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/moving.cc,v 1.1 2004/04/14 22:19:44 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_moving_cc);
 }
 
@@ -25,6 +24,13 @@ namespace CarpetRegrid {
   
   int Moving (cGH const * const cctkGH,
               gh<dim> const & hh,
+              int const reflevel,
+              int const map,
+              int const size,
+              jjvect const & nboundaryzones,
+              jjvect const & is_internal,
+              jjvect const & is_staggered,
+              jjvect const & shiftout,
               gh<dim>::rexts  & bbsss,
               gh<dim>::rbnds  & obss,
               gh<dim>::rprocs & pss)
@@ -32,37 +38,52 @@ namespace CarpetRegrid {
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     assert (refinement_levels >= 1);
     
+    // do nothing if the levels already exist
+    if (reflevel == refinement_levels) return 0;
+    
     assert (bbsss.size() >= 1);
     
     bbsss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
-    bvect const symmetric (symmetry_x, symmetry_y, symmetry_z);
-    bbvect const ob (false);
-    
-    assert (! smart_outer_boundaries);
+    ivect rstr = hh.baseextent.stride();
+    ivect rlb  = hh.baseextent.lower();
+    ivect rub  = hh.baseextent.upper();
     
     for (size_t rl=1; rl<bbsss.size(); ++rl) {
       
-      // calculate new extent
-      CCTK_REAL const argument = 2*M_PI * moving_circle_frequency * cctk_time;
-      rvect const pos
-        (moving_centre_x + moving_circle_radius * cos(argument),
-         moving_centre_y + moving_circle_radius * sin(argument),
-         moving_centre_z);
-      CCTK_REAL const radius = moving_region_radius / ipow(reffact, rl-1);
+      // save old values
+      ivect const oldrlb = rlb;
+      ivect const oldrub = rub;
+      
+      // refined boxes have smaller stride
+      assert (all(rstr%hh.reffact == 0));
+      rstr /= hh.reffact;
       
-      rvect const rlb (symmetric.ifthen (rvect(0),      pos - rvect(radius)));
-      rvect const rub (symmetric.ifthen (rvect(radius), pos + rvect(radius)));
+      // calculate new extent
+      rvect pos;
+      pos[0] = moving_centre_x + moving_circle_radius * cos(2*M_PI * moving_circle_frequency * cctk_time);
+      pos[1] = moving_centre_y + moving_circle_radius * sin(2*M_PI * moving_circle_frequency * cctk_time);
+      pos[2] = moving_centre_z;
+      ivect const centre = floor(rvect(rub - rlb) * pos / rstr + 0.5) * rstr;
+      ivect const radius = floor(rvect(rub - rlb) * moving_region_radius / rstr + 0.5) * rstr;
+      rlb = oldrlb + centre - radius;
+      rub = oldrlb + centre + radius;
+      assert (all(rlb >= oldrlb && rub <= oldrub));
       
-      vector<ibbox> bbs;
-      gh<dim>::cbnds obs;
+      ibbox const bb (rlb, rub, rstr);
+      vector<ibbox> bbs (1);
+      bbs.at(0) = bb;
       
-      ManualCoordinates_OneLevel
-        (cctkGH, hh, rl, refinement_levels, rlb, rub, ob, bbs, obs);
+      bbvect const ob (false);
+      gh<dim>::cbnds obs (1);
+      obs.at(0) = ob;
       
       // make multiprocessor aware
       gh<dim>::cprocs ps;
@@ -70,7 +91,10 @@ namespace CarpetRegrid {
       
       // make multigrid aware
       vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+      MakeMultigridBoxes
+        (cctkGH,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbs, obs, bbss);
       
       bbsss.at(rl) = bbss;
       obss.at(rl) = obs;
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 37d6b466c..279044927 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -13,7 +13,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.37 2004/02/09 16:48:43 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.38 2004/04/14 22:19:44 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_regrid_cc);
 }
 
@@ -175,6 +175,13 @@ namespace CarpetRegrid {
          size, nboundaryzones, is_internal, is_staggered, shiftout,
          bbsss, obss, pss);
                  
+    } else if (CCTK_EQUALS(refined_regions, "moving")) {
+      
+      do_recompose = Moving
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "automatic")) {
       
       do_recompose = Automatic
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index 27edc99c9..2de83f094 100644
--- a/Carpet/CarpetRegrid/src/regrid.hh
+++ b/Carpet/CarpetRegrid/src/regrid.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.11 2004/01/25 14:57:30 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.12 2004/04/14 22:19:44 schnetter Exp $
 
 #ifndef CARPETREGRID_HH
 #define CARPETREGRID_HH
@@ -151,6 +151,19 @@ namespace CarpetRegrid {
                             gh<dim>::rexts  & bbsss,
                             gh<dim>::rbnds  & obss,
                             gh<dim>::rprocs & pss);
+
+  int Moving (cGH const * const cctkGH,
+              gh<dim> const & hh,
+              int const reflevel,
+              int const map,
+              int const size,
+              jjvect const & nboundaryzones,
+              jjvect const & is_internal,
+              jjvect const & is_staggered,
+              jjvect const & shiftout,
+              gh<dim>::rexts  & bbsss,
+              gh<dim>::rbnds  & obss,
+              gh<dim>::rprocs & pss);
   
   int Automatic (cGH const * const cctkGH,
                  gh<dim> const & hh,