Implement vector groups correctly

darcs-hash:20040323114027-07bb3-1ac44b7b144e8936307942bc999d54ecb579719f.gz

7 files changed, 102 insertions, 40 deletions

diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index 5c631f0b5..8eec35113 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -10,7 +10,7 @@
 #include "carpet.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.31 2004/02/09 13:06:05 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.32 2004/03/23 12:40:27 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_Carpet_Storage_cc);
 }
 
@@ -79,21 +79,29 @@ namespace Carpet {
       }
     }
     
-    // VGF: allocate
+    cGroup gp;
+    const int ierr = CCTK_GroupData (group, &gp);
+    assert (!ierr);
+    const int vectorlength = gp.vectorgroup ? gp.vectorlength : 1;
     
     assert (arrdata.at(group).at(0).data.size()==0
 	    || arrdata.at(group).at(0).data.at(0) == 0);
     for (int m=0; m<maps; ++m) {
       for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
+        const int vectorindex = gp.vectorgroup ? var % vectorlength : 0;
+        ggf<dim>* vectorleader
+          = (gp.vectorgroup && vectorindex>0
+             ? arrdata.at(group).at(m).data.at(var - vectorindex)
+             : NULL);
         const int n = n0 + var;
         switch (CCTK_VarTypeI(n)) {
 #define TYPECASE(N,T)                                                   \
         case N:                                                         \
-          /* VGF */                                                     \
           arrdata.at(group).at(m).data.at(var) = new gf<T,dim>          \
             (n, groupdata.at(group).transport_operator,                 \
              *arrdata.at(group).at(m).tt, *arrdata.at(group).at(m).dd,  \
-             tmin, tmax, prolongation_order_time);                      \
+             tmin, tmax, prolongation_order_time,                       \
+             vectorlength, vectorindex, (gf<T,dim>*)vectorleader);      \
           break;
 #include "typecase"
 #undef TYPECASE
@@ -157,26 +165,24 @@ namespace Carpet {
         case N:                                                         \
           assert (arrdata.at(group).at(m).data.at(var));                \
           delete (gf<T,dim>*)arrdata.at(group).at(m).data.at(var);      \
-          arrdata.at(group).at(m).data.at(var) = 0;                     \
+          arrdata.at(group).at(m).data.at(var) = NULL;                  \
           break;
 #include "typecase"
 #undef TYPECASE
         default:
           UnsupportedVarType(n);
         } // switch
-        arrdata.at(group).at(m).data.at(var) = 0;
+        arrdata.at(group).at(m).data.at(var) = NULL;
         
         if (CCTK_GroupTypeI(group) != CCTK_GF) {
           for (int tl=0; tl<num_tl; ++tl) {
-            cgh->data[n][tl] = 0;
+            cgh->data[n][tl] = NULL;
           }
         }
         
       } // for
     } // for
     
-    // VGF: free
-    
     // storage was enabled previously
     return 1;
   }
@@ -198,7 +204,7 @@ namespace Carpet {
     
     assert (group<(int)arrdata.size());
     assert (var<(int)arrdata.at(group).at(0).data.size());
-    return arrdata.at(group).at(0).data.at(var) != 0;
+    return arrdata.at(group).at(0).data.at(var) != NULL;
   }
   
   
diff --git a/Carpet/Carpet/src/variables.hh b/Carpet/Carpet/src/variables.hh
index ba0b9aeed..82d4a41e5 100644
--- a/Carpet/Carpet/src/variables.hh
+++ b/Carpet/Carpet/src/variables.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.hh,v 1.1 2004/01/25 14:57:28 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.hh,v 1.2 2004/03/23 12:40:27 schnetter Exp $
 
 // It is assumed that each group has at least one map.  All arrays
 // have exactly one map.  All maps have the same number of refinement
@@ -118,7 +118,6 @@ namespace Carpet {
     dh<dim>* dd;
     th<dim>* tt;
     vector<ggf<dim>*> data;     // [var]
-    // VGF
   };
   extern vector<vector<arrdesc> > arrdata; // [group][map]
   
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 4c30b44cf..ec6c87858 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.46 2004/03/23 12:14:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.47 2004/03/23 12:40:27 schnetter Exp $
 
 #include <assert.h>
 #include <limits.h>
@@ -28,17 +28,34 @@ using namespace std;
 
 // Constructors
 template<class T, int D>
-data<T,D>::data (const int varindex_, const operator_type transport_operator_)
+data<T,D>::data (const int varindex_, const operator_type transport_operator_,
+                 const int vectorlength, const int vectorindex,
+                 data* const vectorleader)
   : gdata<D>(varindex_, transport_operator_),
-    _storage(0)
-{ }
+    _storage(0),
+    vectorlength(vectorlength), vectorindex(vectorindex),
+    vectorleader(vectorleader)
+{
+  assert (vectorlength>=1);
+  assert (vectorindex>=0 && vectorindex<vectorlength);
+  assert ((vectorindex==0 && !vectorleader)
+          || (vectorindex!=0 && vectorleader));
+}
 
 template<class T, int D>
 data<T,D>::data (const int varindex_, const operator_type transport_operator_,
+                 const int vectorlength, const int vectorindex,
+                 data* const vectorleader,
                  const ibbox& extent_, const int proc_)
   : gdata<D>(varindex_, transport_operator_),
-    _storage(0)
+    _storage(0),
+    vectorlength(vectorlength), vectorindex(vectorindex),
+    vectorleader(vectorleader)
 {
+  assert (vectorlength>=1);
+  assert (vectorindex>=0 && vectorindex<vectorlength);
+  assert ((vectorindex==0 && !vectorleader)
+          || (vectorindex!=0 && vectorleader));
   allocate(extent_, proc_);
 }
 
@@ -80,7 +97,13 @@ void data<T,D>::allocate (const ibbox& extent_, const int proc_,
   if (rank==this->_proc) {
     this->_owns_storage = !mem;
     if (this->_owns_storage) {
-      this->_storage = new T[this->_size];
+      if (this->vectorindex == 0) {
+        assert (! this->vectorleader);
+        this->_storage = new T[this->vectorlength * this->_size];
+      } else {
+        assert (this->vectorleader);
+        this->_storage = this->vectorleader->vectordata (this->vectorindex);
+      }
     } else {
       this->_storage = (T*)mem;
     }
@@ -91,19 +114,33 @@ void data<T,D>::allocate (const ibbox& extent_, const int proc_,
 
 template<class T, int D>
 void data<T,D>::free () {
-  if (this->_storage && this->_owns_storage) delete [] _storage;
+  if (this->_storage && this->_owns_storage && this->vectorindex==0) {
+    delete [] _storage;
+  }
   _storage = 0;
   this->_has_storage = false;
 }
 
 template<class T, int D>
 void data<T,D>::transfer_from (gdata<D>* gsrc) {
+  assert (this->vectorlength==1);
   data* src = (data*)gsrc;
+  assert (src->vectorlength==1);
   assert (!_storage);
   *this = *src;
   *src = data(this->varindex, this->transport_operator);
 }
 
+template<class T, int D>
+T* data<T,D>::vectordata (const int vectorindex) const
+{
+  assert (this->vectorindex==0);
+  assert (! this->vectorleader);
+  assert (vectorindex>=0 && vectorindex<this->vectorlength);
+  assert (this->_storage && this->_owns_storage);
+  return this->_storage + vectorindex * this->_size;
+}
+
 
 
 // Processor management
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index a6658e23e..85e8a28fe 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.17 2004/03/23 12:14:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.18 2004/03/23 12:40:27 schnetter Exp $
 
 #ifndef DATA_HH
 #define DATA_HH
@@ -30,13 +30,21 @@ class data: public gdata<D> {
 
   // Fields
   T* _storage;			// the data (if located on this processor)
+   
+  int vectorlength;
+  int vectorindex;
+  data* vectorleader;
   
 public:
   
   // Constructors
   data (const int varindex = -1,
-        const operator_type transport_operator = op_error);
+        const operator_type transport_operator = op_error,
+        const int vectorlength = 1, const int vectorindex = 0,
+        data* const vectorleader = NULL);
   data (const int varindex, const operator_type transport_operator,
+        const int vectorlength, const int vectorindex,
+        data* const vectorleader,
         const ibbox& extent, const int proc);
 
   // Destructors
@@ -52,6 +60,10 @@ public:
   virtual void free ();
   virtual void transfer_from (gdata<D>* gsrc);
 
+private:
+  T* vectordata (const int vectorindex) const;
+public:
+
   // Processor management
   virtual void change_processor (comm_state<D>& state,
                                  const int newproc, void* const mem=0);
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 2fcec0d2f..b2b5198d4 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.15 2004/01/25 14:57:30 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.16 2004/03/23 12:40:27 schnetter Exp $
 
 #include <assert.h>
 
@@ -13,15 +13,16 @@ using namespace std;
 
 
 // Constructors
-// VGF
 template<class T,int D>
 gf<T,D>::gf (const int varindex, const operator_type transport_operator,
              th<D>& t, dh<D>& d,
-	     const int tmin, const int tmax, const int prolongation_order_time)
+	     const int tmin, const int tmax, const int prolongation_order_time,
+             const int vectorlength, const int vectorindex,
+             gf* const vectorleader)
   : ggf<D>(varindex, transport_operator,
-           t, d, tmin, tmax, prolongation_order_time)
+           t, d, tmin, tmax, prolongation_order_time,
+           vectorlength, vectorindex, vectorleader)
 {
-  // VGF
   this->recompose (0, true);
 }
 
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index a1da9d3d6..877308389 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.34 2004/03/23 09:26:49 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.35 2004/03/23 12:40:27 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -23,12 +23,16 @@ template<int D>
 ggf<D>::ggf (const int varindex, const operator_type transport_operator,
              th<D>& t, dh<D>& d,
              const int tmin, const int tmax,
-             const int prolongation_order_time)
+             const int prolongation_order_time,
+             const int vectorlength, const int vectorindex,
+             ggf* const vectorleader)
   : varindex(varindex), transport_operator(transport_operator), t(t),
     tmin(tmin), tmax(tmax),
     prolongation_order_time(prolongation_order_time),
     h(d.h), d(d),
-    storage(tmax-tmin+1)
+    storage(tmax-tmin+1),
+    vectorlength(vectorlength), vectorindex(vectorindex),
+    vectorleader(vectorleader)
 {
   assert (&t.h == &d.h);
 
@@ -50,7 +54,6 @@ bool ggf<D>::operator== (const ggf<D>& f) const {
 
 
 // Modifiers
-// VGF
 template<int D>
 void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
   
@@ -68,7 +71,7 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
       for (int c=0; c<h.components(rl); ++c) {
       	storage[tl-tmin][rl][c].resize(h.mglevels(rl,c));
 	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-          storage[tl-tmin][rl][c][ml] = 0;
+          storage[tl-tmin][rl][c][ml] = NULL;
 	} // for ml
       } // for c
     } // for rl
@@ -80,10 +83,9 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
       for (int c=0; c<h.components(rl); ++c) {
       	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
 	  
-	  storage[tl-tmin][rl][c][ml] = typed_data();
+	  storage[tl-tmin][rl][c][ml] = typed_data(tl,rl,c,ml);
 	  
       	  // Allocate storage
-          // VGF
       	  storage[tl-tmin][rl][c][ml]->allocate
       	    (d.boxes[rl][c][ml].exterior, h.proc(rl,c));
 	  
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 7b7c5c159..64998c8a2 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.19 2004/03/23 09:26:49 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.20 2004/03/23 12:40:27 schnetter Exp $
 
 #ifndef GGF_HH
 #define GGF_HH
@@ -47,9 +47,9 @@ class ggf {
   typedef vector<mdata> cdata;  // ... for each component
   typedef vector<cdata> rdata;  // ... for each refinement level
   typedef vector<rdata> fdata;  // ... for each time level
-
+  
 public:				// should be readonly
-
+  
   // Fields
   int varindex;                 // Cactus variable index
   operator_type transport_operator;
@@ -63,14 +63,20 @@ public:				// should be readonly
 
 protected:
   fdata storage;		// storage
-
+  
+  int vectorlength;             // vector length
+  int vectorindex;              // index of *this
+  ggf* vectorleader;            // first vector element
+  
 public:
 
   // Constructors
   ggf (const int varindex, const operator_type transport_operator,
        th<D>& t, dh<D>& d,
        const int tmin, const int tmax,
-       const int prolongation_order_time);
+       const int prolongation_order_time,
+       const int vectorlength, const int vectorindex,
+       ggf* const vectorleader);
 
   // Destructors
   virtual ~ggf ();
@@ -81,7 +87,6 @@ public:
 
 
   // Modifiers
-  // VGF
   void recompose (const int initialise_from, const bool do_prolongate);
 
   // Cycle the time levels by rotating the data sets
@@ -96,7 +101,7 @@ public:
   
 protected:
   
-  virtual gdata<D>* typed_data() = 0;
+  virtual gdata<D>* typed_data (int tl, int rl, int c, int ml) = 0;