Use boost::numeric::ublas::tensor as underlying data structure for DistributedFullGrid

examples/combi_example/TaskExample.hpp

@@ -87,7 +87,7 @@ class TaskExample : public Task {
     dfg_ = new DistributedFullGrid<CombiDataType>(dim, l, lcomm, this->getBoundary(), p);
 
     /* loop over local subgrid and set initial values */
-    std::vector<CombiDataType>& elements = dfg_->getElementVector();
+    auto& elements = dfg_->getElementVector();
 
     for (size_t i = 0; i < elements.size(); ++i) {
       IndexType globalLinearIndex = dfg_->getGlobalLinearIndex(i);
@@ -109,7 +109,7 @@ class TaskExample : public Task {
 
     int lrank = theMPISystem()->getLocalRank();
 
-    std::vector<CombiDataType>& elements = dfg_->getElementVector();
+    auto& elements = dfg_->getElementVector();
     // TODO if your Example uses another data structure, you need to copy
     // the data from elements to that data structure
 

examples/combi_example_faults/TaskExample.hpp

@@ -91,7 +91,7 @@ class TaskExample: public Task {
     dfg_ = new DistributedFullGrid<CombiDataType>(dim, l, lcomm, this->getBoundary(), p);
 
     /* loop over local subgrid and set initial values */
-    std::vector<CombiDataType>& elements = dfg_->getElementVector();
+    auto& elements = dfg_->getElementVector();
 
     for (size_t i = 0; i < elements.size(); ++i) {
       IndexType globalLinearIndex = dfg_->getGlobalLinearIndex(i);
@@ -115,7 +115,7 @@ class TaskExample: public Task {
 
     /* pseudo timestepping to demonstrate the behaviour of your typical
      * time-dependent simulation problem. */
-    std::vector<CombiDataType>& elements = dfg_->getElementVector();
+    auto& elements = dfg_->getElementVector();
 
     for (size_t step = stepsTotal_; step < stepsTotal_ + nsteps_; ++step) {
       real time = (step + 1)* dt_;

examples/distributed_advection/TaskAdvection.hpp

@@ -1,8 +1,9 @@
 #pragma once
 
+#include <boost/numeric/ublas/tensor.hpp>
+
 #include "fullgrid/DistributedFullGrid.hpp"
 #include "task/Task.hpp"
-
 namespace combigrid {
 
 // exact solution
@@ -21,17 +22,18 @@ class TestFn {
     // leave out normalization, such that maximum is always 1
     return std::exp(exponent);  // / std::sqrt( std::pow(2*pi*sigmaSquared, coords.size()));
   }
+
  private:
-  static constexpr double sigmaSquaredInv_ = 1. / ((1./3.)*(1./3.));
+  static constexpr double sigmaSquaredInv_ = 1. / ((1. / 3.) * (1. / 3.));
 };
 
 class TaskAdvection : public Task {
  public:
   /* if the constructor of the base task class is not sufficient we can provide an
    * own implementation. here, we add dt, nsteps, and p as a new parameters.
    */
-  TaskAdvection(DimType dim, const LevelVector& l, const std::vector<BoundaryType>& boundary,
-                real coeff, LoadModel* loadModel, real dt, size_t nsteps,
+  TaskAdvection(const LevelVector& l, const std::vector<BoundaryType>& boundary, real coeff,
+                LoadModel* loadModel, real dt, size_t nsteps,
                 std::vector<int> p = std::vector<int>(0),
                 FaultCriterion* faultCrit = (new StaticFaults({0, IndexVector(0), IndexVector(0)})))
       : Task(l, boundary, coeff, loadModel, faultCrit),
@@ -59,7 +61,7 @@ class TaskAdvection : public Task {
     dfg_ = new DistributedFullGrid<CombiDataType>(dim, l, lcomm, this->getBoundary(), p_, false,
                                                   decomposition);
     if (phi_ == nullptr) {
-      phi_ = new std::vector<CombiDataType>(dfg_->getNrLocalElements());
+      phi_ = new DistributedFullGrid<CombiDataType>::TensorType(dfg_->getLocalExtents());
     }
 
     std::vector<double> h = dfg_->getGridSpacing();
@@ -95,8 +97,7 @@ class TaskAdvection : public Task {
     std::vector<double> h = dfg_->getGridSpacing();
     const auto& fullOffsets = dfg_->getLocalOffsets();
 
-    phi_->resize(dfg_->getNrLocalElements());
-    std::fill(phi_->begin(), phi_->end(), 0.);
+    phi_->reshape(dfg_->getLocalExtents(), 0.);
 
     for (size_t i = 0; i < nsteps_; ++i) {
       // compute the gradient in the original dfg_, then update into phi_ and
@@ -122,7 +123,7 @@ class TaskAdvection : public Task {
           dfg_->getLocalVectorIndex(li, locAxisIndex);
           // TODO can be unrolled into ghost and other part, avoiding if-statement
           CombiDataType phi_neighbor = 0.;
-          if (locAxisIndex[d] == 0){
+          if (locAxisIndex[d] == 0) {
             assert(phi_ghost.size() > 0);
             // then use values from boundary exchange
             IndexType gli = 0;
@@ -132,10 +133,7 @@ class TaskAdvection : public Task {
             assert(gli > -1);
             assert(gli < phi_ghost.size());
             phi_neighbor = phi_ghost[gli];
-          } else{
-            // --locAxisIndex[d];
-            // IndexType lni = dfg_->getLocalLinearIndex(locAxisIndex);
-            // assert(lni == (li - fullOffsets[d]));
+          } else {
             phi_neighbor = dfg_->getElementVector()[li - fullOffsets[d]];
           }
           // calculate gradient of phi with backward differential quotient
@@ -144,10 +142,11 @@ class TaskAdvection : public Task {
           u_dot_dphi[li] += velocity[d] * dphi;
         }
       }
+      // update all values in phi_
       for (IndexType li = 0; li < dfg_->getNrLocalElements(); ++li) {
         (*phi_)[li] = dfg_->getElementVector()[li] - u_dot_dphi[li] * dt_;
       }
-      phi_->swap(dfg_->getElementVector());
+      std::swap(*phi_, dfg_->getElementVector());
     }
     stepsTotal_ += nsteps_;
 
@@ -161,7 +160,8 @@ class TaskAdvection : public Task {
    * solution on one process and then converting it to the full grid representation.
    * the DistributedFullGrid class offers a convenient function to do this.
    */
-  void getFullGrid(FullGrid<CombiDataType>& fg, RankType r, CommunicatorType lcomm, int n = 0) override {
+  void getFullGrid(FullGrid<CombiDataType>& fg, RankType r, CommunicatorType lcomm,
+                   int n = 0) override {
     assert(fg.getLevels() == dfg_->getLevels());
 
     dfg_->gatherFullGrid(fg, r);
@@ -181,15 +181,13 @@ class TaskAdvection : public Task {
     phi_ = nullptr;
   }
 
-  real getCurrentTime() const override {
-    return stepsTotal_ * dt_;
-  }
+  real getCurrentTime() const override { return stepsTotal_ * dt_; }
 
   CombiDataType analyticalSolution(const std::vector<real>& coords, int n = 0) const override {
     assert(n == 0);
     auto coordsCopy = coords;
     TestFn f;
-    return f(coordsCopy, stepsTotal_*dt_);
+    return f(coordsCopy, stepsTotal_ * dt_);
   }
 
  protected:
@@ -212,7 +210,7 @@ class TaskAdvection : public Task {
   bool initialized_;
   size_t stepsTotal_;
   DistributedFullGrid<CombiDataType>* dfg_;
-  static std::vector<CombiDataType>* phi_;
+  static typename DistributedFullGrid<CombiDataType>::TensorType* phi_;
 
   /**
    * The serialize function has to be extended by the new member variables.
@@ -234,6 +232,6 @@ class TaskAdvection : public Task {
   }
 };
 
-std::vector<CombiDataType>* TaskAdvection::phi_ = nullptr;
+typename DistributedFullGrid<CombiDataType>::TensorType* TaskAdvection::phi_ = nullptr;
 
 }  // namespace combigrid

examples/distributed_advection/combi_example.cpp

@@ -178,7 +178,7 @@ int main(int argc, char** argv) {
     TaskContainer tasks;
     std::vector<size_t> taskIDs;
     for (size_t i = 0; i < levels.size(); i++) {
-      Task* t = new TaskAdvection(dim, levels[i], boundary, coeffs[i], loadmodel.get(), dt,
+      Task* t = new TaskAdvection(levels[i], boundary, coeffs[i], loadmodel.get(), dt,
                                   nsteps, p);
 
       static_assert(!isGENE, "isGENE");

examples/distributed_third_level/TaskAdvection.hpp

@@ -1,8 +1,8 @@
 #pragma once
 
+
 #include "fullgrid/DistributedFullGrid.hpp"
 #include "task/Task.hpp"
-
 namespace combigrid {
 
 // exact solution
@@ -21,17 +21,18 @@ class TestFn {
     // leave out normalization, such that maximum is always 1
     return std::exp(exponent);  // / std::sqrt( std::pow(2*pi*sigmaSquared, coords.size()));
   }
+
  private:
-  static constexpr double sigmaSquaredInv_ = 1. / ((1./3.)*(1./3.));
+  static constexpr double sigmaSquaredInv_ = 1. / ((1. / 3.) * (1. / 3.));
 };
 
 class TaskAdvection : public Task {
  public:
   /* if the constructor of the base task class is not sufficient we can provide an
    * own implementation. here, we add dt, nsteps, and p as a new parameters.
    */
-  TaskAdvection(const LevelVector& l, const std::vector<BoundaryType>& boundary,
-                real coeff, LoadModel* loadModel, real dt, size_t nsteps,
+  TaskAdvection(const LevelVector& l, const std::vector<BoundaryType>& boundary, real coeff,
+                LoadModel* loadModel, real dt, size_t nsteps,
                 std::vector<int> p = std::vector<int>(0),
                 FaultCriterion* faultCrit = (new StaticFaults({0, IndexVector(0), IndexVector(0)})))
       : Task(l, boundary, coeff, loadModel, faultCrit),
@@ -59,7 +60,7 @@ class TaskAdvection : public Task {
     dfg_ = new DistributedFullGrid<CombiDataType>(dim, l, lcomm, this->getBoundary(), p_, false,
                                                   decomposition);
     if (phi_ == nullptr) {
-      phi_ = new std::vector<CombiDataType>(dfg_->getNrLocalElements());
+      phi_ = new DistributedFullGrid<CombiDataType>::TensorType(dfg_->getLocalExtents());
     }
 
     std::vector<double> h = dfg_->getGridSpacing();
@@ -95,8 +96,7 @@ class TaskAdvection : public Task {
     std::vector<double> h = dfg_->getGridSpacing();
     const auto& fullOffsets = dfg_->getLocalOffsets();
 
-    phi_->resize(dfg_->getNrLocalElements());
-    std::fill(phi_->begin(), phi_->end(), 0.);
+    phi_->reshape(dfg_->getLocalExtents(), 0.);
 
     for (size_t i = 0; i < nsteps_; ++i) {
       // compute the gradient in the original dfg_, then update into phi_ and
@@ -122,7 +122,7 @@ class TaskAdvection : public Task {
           dfg_->getLocalVectorIndex(li, locAxisIndex);
           // TODO can be unrolled into ghost and other part, avoiding if-statement
           CombiDataType phi_neighbor = 0.;
-          if (locAxisIndex[d] == 0){
+          if (locAxisIndex[d] == 0) {
             assert(phi_ghost.size() > 0);
             // then use values from boundary exchange
             IndexType gli = 0;
@@ -132,10 +132,7 @@ class TaskAdvection : public Task {
             assert(gli > -1);
             assert(gli < phi_ghost.size());
             phi_neighbor = phi_ghost[gli];
-          } else{
-            // --locAxisIndex[d];
-            // IndexType lni = dfg_->getLocalLinearIndex(locAxisIndex);
-            // assert(lni == (li - fullOffsets[d]));
+          } else {
             phi_neighbor = dfg_->getElementVector()[li - fullOffsets[d]];
           }
           // calculate gradient of phi with backward differential quotient
@@ -144,10 +141,11 @@ class TaskAdvection : public Task {
           u_dot_dphi[li] += velocity[d] * dphi;
         }
       }
+      // update all values in phi_
       for (IndexType li = 0; li < dfg_->getNrLocalElements(); ++li) {
         (*phi_)[li] = dfg_->getElementVector()[li] - u_dot_dphi[li] * dt_;
       }
-      phi_->swap(dfg_->getElementVector());
+      std::swap(*phi_, dfg_->getElementVector());
     }
     stepsTotal_ += nsteps_;
 
@@ -161,7 +159,8 @@ class TaskAdvection : public Task {
    * solution on one process and then converting it to the full grid representation.
    * the DistributedFullGrid class offers a convenient function to do this.
    */
-  void getFullGrid(FullGrid<CombiDataType>& fg, RankType r, CommunicatorType lcomm, int n = 0) override {
+  void getFullGrid(FullGrid<CombiDataType>& fg, RankType r, CommunicatorType lcomm,
+                   int n = 0) override {
     assert(fg.getLevels() == dfg_->getLevels());
 
     dfg_->gatherFullGrid(fg, r);
@@ -181,15 +180,13 @@ class TaskAdvection : public Task {
     phi_ = nullptr;
   }
 
-  real getCurrentTime() const override {
-    return stepsTotal_ * dt_;
-  }
+  real getCurrentTime() const override { return stepsTotal_ * dt_; }
 
   CombiDataType analyticalSolution(const std::vector<real>& coords, int n = 0) const override {
     assert(n == 0);
     auto coordsCopy = coords;
     TestFn f;
-    return f(coordsCopy, stepsTotal_*dt_);
+    return f(coordsCopy, stepsTotal_ * dt_);
   }
 
  protected:
@@ -212,7 +209,7 @@ class TaskAdvection : public Task {
   bool initialized_;
   size_t stepsTotal_;
   DistributedFullGrid<CombiDataType>* dfg_;
-  static std::vector<CombiDataType>* phi_;
+  static typename DistributedFullGrid<CombiDataType>::TensorType* phi_;
 
   /**
    * The serialize function has to be extended by the new member variables.
@@ -234,6 +231,6 @@ class TaskAdvection : public Task {
   }
 };
 
-std::vector<CombiDataType>* TaskAdvection::phi_ = nullptr;
+typename DistributedFullGrid<CombiDataType>::TensorType* TaskAdvection::phi_ = nullptr;
 
 }  // namespace combigrid

examples/distributed_third_level/TaskConst.hpp

@@ -38,7 +38,7 @@ class TaskConst : public combigrid::Task {
     dfg_ = new DistributedFullGrid<CombiDataType>(getDim(), getLevelVector(), lcomm, getBoundary(),
                                                   p, false, decomposition);
 
-    std::vector<CombiDataType>& elements = dfg_->getElementVector();
+    auto& elements = dfg_->getElementVector();
     for (auto& element : elements) {
       element = 10;
     }
@@ -49,7 +49,7 @@ class TaskConst : public combigrid::Task {
 
     // std::cout << "run " << getCommRank(lcomm) << std::endl;    
 
-    std::vector<CombiDataType>& elements = dfg_->getElementVector();
+    auto& elements = dfg_->getElementVector();
     for (auto& element : elements) {
       element = getLevelVector()[0] / (double)getLevelVector()[1];
       // std::cout << "e " << element << std::endl;

examples/selalib_distributed/src/SelalibTask.hpp

@@ -299,7 +299,7 @@ class SelalibTask : public combigrid::Task {
 
   void setDFGfromLocalDistribution() {
 #ifndef NDEBUG
-    auto& localDFGSize = dfg_->getLocalSizes();
+    auto& localDFGSize = dfg_->getLocalExtents();
     assert(dim_ == 6);
     // std::cout << localDFGSize << " " << std::vector<int>(localSize_.begin(), localSize_.end()) <<
     // std::endl;