Merge branch '403-add-stdist-warnings' into 'feature'

Resolve "add stdist warnings"

Closes #403

See merge request mpp/mpp!394

src/lib2_mesh/distribution/Distribution.cpp

@@ -69,7 +69,7 @@ void Distribution::DistributeMesh(Mesh &mesh) {
 
   if (distribute) {
 
-    for(int i = startDist; i<cells.size(); ++i){
+    for (int i = startDist; i < cells.size(); ++i) {
       distributedCells.emplace_back(cells[i]());
     }
 
@@ -85,8 +85,7 @@ void Distribution::DistributeMesh(Mesh &mesh) {
            << endl;
 
   if (PPM->Min(cellCount, commSplit) == 0)
-    VerboseWarning(
-        "At least one processor has no cells!", 1);
+    VerboseWarning("At least one processor has no cells!", 1);
 
   emptyData(clearDistributed);
   mout.EndBlock(true);
@@ -97,7 +96,7 @@ void Distribution::Communicate(Mesh &mesh) {
   for (int q = 0; q < PPM->Size(commSplit); ++q) {
     for (auto c = markedCells[q].begin(); c != markedCells[q].end(); ++c) {
       for (int i = 0; i < c->Corners(); ++i) {
-          mesh.procSets.Add(c->Corner(i), q);
+        mesh.procSets.Add(c->Corner(i), q);
       }
       for (int i = 0; i < c->SpaceCell().Edges(); ++i)
         mesh.procSets.Add(c->SpaceCell().Edge(i), q);
@@ -327,6 +326,21 @@ bool TLess_t(const cell &c0, const cell &c1) {
   return false;
 }
 
+std::function<bool(const cell &, const cell &)> GetLess(char var) {
+  switch (var) {
+    case 'x':
+      return Less_x;
+    case 'y':
+      return Less_y;
+    case 'z':
+      return Less_z;
+    case 't':
+      return Less_t;
+    default:
+      THROW("Dimension " + std::to_string(var) + " not known");
+  }
+}
+
 void Distribution::RCB_x(int rekDepth, int begin, int end) {
   if (rekDepth == 0) return;
   sort(cells.begin() + begin, cells.begin() + end, Less_x);
@@ -794,19 +808,7 @@ void Distribution::timeOpt(int tProc, int level, int begin, int end) {
 
 void Distribution::timercb(int level, int begin, int end, char d) {
   int N = end - begin;
-  switch (d) {
-    case 'x':
-      sort(cells.begin() + begin, cells.begin() + end, TLess_x);
-      break;
-    case 'y':
-      sort(cells.begin() + begin, cells.begin() + end, TLess_y);
-      break;
-    case 't':
-      sort(cells.begin() + begin, cells.begin() + end, TLess_t);
-      break;
-    default:
-      break;
-  }
+  sort(cells.begin() + begin, cells.begin() + end, GetLess(d));
   if (level == 0) {
     return;
   } else {
@@ -865,23 +867,8 @@ void Distribution::timercb_weighted(int tlevel, int xlevel, int ylevel,
                                     int default_weight) {
   int N = end - begin;
   int level = tlevel + xlevel + ylevel;
-
-  switch (d) {
-    case 'x':
-      sort(cells.begin() + begin, cells.begin() + end, TLess_x);
-      break;
-    case 'y':
-      sort(cells.begin() + begin, cells.begin() + end, TLess_y);
-      break;
-    case 't':
-      sort(cells.begin() + begin, cells.begin() + end, TLess_t);
-      break;
-    default:
-      break;
-  }
-
+  sort(cells.begin() + begin, cells.begin() + end, GetLess(d));
   if (PPM->proc() == 0) {
-
     int W_sum = 0;
     for (int i = begin; i < end; ++i) {
       if (W.find(cells[i]()) != W.end())
@@ -919,7 +906,7 @@ void Distribution::timercb_weighted(int tlevel, int xlevel, int ylevel,
 
     while (cells[mid]().t() == cells[mid - 1]().t() || a == b) {
       --mid;
-      if (mid < begin) break;
+      if (mid < begin || mid == 0) break;
       a = cells[mid]().CopyWithT(0.0);
       b = cells[mid + 1]().CopyWithT(0.0);
     }
@@ -1046,8 +1033,13 @@ void Distribution::DistributeSTMesh(Mesh &mesh) {
   } else if (_distName == "deformed_optimized") {
     int S = mesh.steps() - 1;
     int Slice_Maximum = int(log(double(S)) / log(2.0) + 1e-10);
-    int X_Maximum = int(log(double(sqrt(N / (S)))) / log(2.0) + 1e-10);
+    int X_Maximum = int(log(double(sqrt(N / S))) / log(2.0) + 1e-10);
     int Y_Maximum = X_Maximum;
+
+    if (Slice_Maximum == 0 && X_Maximum == 0 && Y_Maximum == 0) {
+      THROW("Not enough cells to distribute. Level: " + mesh.Level().str())
+    }
+
     int tL = 0;
     int yL = 0;
     int xL = 0;

src/lib2_mesh/mesh/CoarseGeometry.cpp

@@ -77,7 +77,7 @@ std::vector<double> CoarseGeometry::refineTimeSteps(const std::vector<double> &t
       double dt = init_timesteps[S - 1][k + 2] - init_timesteps[S - 1][k + 1];
       for (unsigned int l = 1; l < SF_inv; ++l)
         init_timesteps[S][k * SF_inv + l + 1] =
-          init_timesteps[S - 1][k + 1] + dt * 1 / double(SF_inv) * l;
+            init_timesteps[S - 1][k + 1] + dt * 1 / double(SF_inv) * l;
     }
   }
   return init_timesteps[refineTimeStepCount];
@@ -136,7 +136,7 @@ void CoarseGeometry::CommunicateGeometry() {
     ExchangeBuffer exBuffer;
     if (PPM->Master()) {
       for (int q = 1; q < PPM->Size(); ++q) {
-        exBuffer.Send(q) << timeSteps.size();
+        exBuffer.Send(q) << int(timeSteps.size());
         for (double timeStep : timeSteps) {
           exBuffer.Send(q) << double(timeStep);
         }
@@ -147,8 +147,9 @@ void CoarseGeometry::CommunicateGeometry() {
       int m;
       exBuffer.Receive(0) >> m;
       timeSteps.resize(m);
-      for (double & timeStep : timeSteps) {
-        exBuffer.Receive(0) >> timeStep;
+      for (int i = 0; i < m; i++) {
+        exBuffer.Receive(0) >> timeSteps[i];
+
       }
     }
     exBuffer.ClearBuffers();
@@ -157,7 +158,7 @@ void CoarseGeometry::CommunicateGeometry() {
 
 CoarseGeometry::CoarseGeometry(const MeshPoints &points, const MeshCells &cells,
                                const MeshFaces &faces, const TimeSteps &timeSteps) :
-  timeSteps(timeSteps) {
+    timeSteps(timeSteps) {
   for (const auto &point : points) {
     double z[]{0.0, 0.0, 0.0, 0.0};
     for (int j = 0; j < min(3, int(point.size())); ++j)
@@ -175,8 +176,8 @@ void CoarseGeometry::Scale(double scalingFactor) {
 
 Logging &operator<<(Logging &s, const CoarseGeometry &M) {
   return s << "POINTS: " << M.coordinates.size() << endl << M.coordinates
-  << "CELLS: " << M.cellIds.size() << endl << M.cellIds
-  << "FACES: " << M.faceIds.size() << endl << M.faceIds
+           << "CELLS: " << M.cellIds.size() << endl << M.cellIds
+           << "FACES: " << M.faceIds.size() << endl << M.faceIds
            #ifdef USE_DATAMESH
            << "VDATA: " << M.vdataList.size() << endl << M.vdataList
            << "CDATA: " << M.cdataList.size() << endl << M.cdataList

src/lib2_mesh/mesh/Mesh.hpp

@@ -24,6 +24,13 @@ struct LevelPair {
   LevelPair CoarserInTime() const { return LevelPair{space, time - 1}; }
 
   LevelPair Next() const { return LevelPair{space + 1, time + 1}; }
+
+  std::string str() const {
+    std::stringstream ss;
+    ss << "[" << space << ", " << time << "]";
+    return ss.str();
+  }
+
 };
 
 inline bool operator==(const LevelPair &l1, const LevelPair &l2) {

src/lib5_solvers/solver/linearsolver/BiCGStab.hpp

@@ -164,8 +164,7 @@ public:
       if (d < epsilon) break;
       vout(1) << "d(" << iter << ")= " << d << endl;
       Scalar rho_new = r0 * r;
-      if (((abs(rho) < 1e-8 * abs(rho_new)) && (iter)) ||
-      (abs(omega) < 1e-8 * abs(alpha))) {
+      if ((abs(rho) < 1e-8 * abs(rho_new) && iter) || abs(omega) < 1e-8 * abs(alpha)) {
         vout(2) << "breakdown occured (rho = " << rho
                 << " , omega = " << omega << ")" << endl;
         Restart(cnt, rho, alpha, omega, p, v);
@@ -201,6 +200,7 @@ public:
       rho = rho_new;
       d = norm(r);
     }
+    mout.EndBlock();
     //    vout(0) << "d(" << iter << ")= " << d << " rate " << rate(iter, d_0, d) << endl;
     //        vout(0) << IterationEnd();
   }