src/sand-wedge-data/parset.cfg → src/one-body-problem.cfg

@@ -2,8 +2,12 @@
 gravity         = 9.81  # [m/s^2]
 
 [io]
+data.write      = true
 printProgress   = false
-writeVTK        = false
+restarts.first  = 0
+restarts.spacing= 20
+restarts.write  = true
+vtk.write       = false
 
 [problem]
 finalTime       = 1000  # [s]
@@ -35,11 +39,6 @@ b               = 0.017 # [ ]
 a               = 0.020 # [ ]
 b               = 0.005 # [ ]
 
-[restarts]
-template = restart_%06d
-first = 0
-spacing = 20
-
 [timeSteps]
 scheme = newmark
 

src/program_state.hh

@@ -39,7 +39,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
     using LocalVector = typename Vector::block_type;
     using LocalMatrix = typename Matrix::block_type;
-    auto const dims = LocalVector::dimension;
+    auto constexpr dims = LocalVector::dimension;
 
     // Solving a linear problem with a multigrid solver
     auto const solveLinearProblem = [&](

src/sand-wedge-data/boundaryconditions.py

-import math
-
-class neumannCondition:
-    def __call__(self, relativeTime):
-        return 0
-
-class velocityDirichletCondition:
-    def __call__(self, relativeTime):
-        # Assumed to vanish at time zero
-        finalVelocity = -5e-5
-        if relativeTime <= 0.1:
-            return finalVelocity * ( 1.0 - math.cos(relativeTime * math.pi / 0.1) ) / 2.0
-        return finalVelocity
-
-Functions = {
-    'neumannCondition' : neumannCondition(),
-    'velocityDirichletCondition' : velocityDirichletCondition()
-}

src/spatial-solving/fixedpointiterator.cc

@@ -83,7 +83,12 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
 
   updaters.rate_->postProcess(velocityIterate);
 
-  return { fixedPointIteration, multigridIterations };
+  // Cannot use return { fixedPointIteration, multigridIterations };
+  // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
+  FixedPointIterationCounter ret;
+  ret.iterations = fixedPointIteration;
+  ret.multigridIterations = multigridIterations;
+  return ret;
 }
 
 std::ostream &operator<<(std::ostream &stream,

src/spatial-solving/fixedpointiterator.hh

@@ -9,8 +9,8 @@
 #include <dune/solvers/solvers/solver.hh>
 
 struct FixedPointIterationCounter {
-  size_t iterations;
-  size_t multigridIterations;
+  size_t iterations = 0;
+  size_t multigridIterations = 0;
 
   void operator+=(FixedPointIterationCounter const &other);
 };

src/spatial-solving/solverfactory.hh

@@ -4,11 +4,7 @@
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/parametertree.hh>
 
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wsign-compare"
 #include <dune/solvers/iterationsteps/multigridstep.hh>
-#pragma clang diagnostic pop
-
 #include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>

src/time-stepping/adaptivetimestepper.cc

@@ -32,6 +32,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
       parset_(parset),
       globalFriction_(globalFriction),
       current_(current),
+      R1_(),
       externalForces_(externalForces),
       mustRefine_(mustRefine),
       errorNorm_(errorNorm) {}
@@ -42,73 +43,68 @@ bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::coarsen() {
-  bool didCoarsen = false;
-  while (relativeTime_ + relativeTau_ < 1.0) {
-    R2_ = { R1_.updaters.clone() };
-    step(R2_, relativeTime_ + relativeTau_, relativeTau_);
-
-    UpdatersWithCount C{ current_.clone() };
-    step(C, relativeTime_, 2.0 * relativeTau_);
+IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
+  /*
+    |     C     | We check here if making the step R1 of size tau is a
+    |  R1 | R2  | good idea. To check if we can coarsen, we compare
+    |F1|F2|     | the result of (R1+R2) with C, i.e. two steps of size
+                  tau with one of size 2*tau. To check if we need to
+    refine, we compare the result of (F1+F2) with R1, i.e. two steps
+    of size tau/2 with one of size tau. The method makes multiple
+    coarsening/refining attempts, with coarsening coming first. */
+  if (R1_.updaters == Updaters())
+    R1_ = step(current_, relativeTime_, relativeTau_);
 
-    if (!mustRefine_(C.updaters, R2_.updaters)) {
-      R2_ = {};
-      R1_ = C;
-      relativeTau_ *= 2.0;
-      didCoarsen = true;
-    } else {
+  bool didCoarsen = false;
+  iterationRegister_.reset();
+  UpdatersWithCount R2;
+  UpdatersWithCount C;
+  while (relativeTime_ + relativeTau_ <= 1.0) {
+    R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
+    C = step(current_, relativeTime_, 2 * relativeTau_);
+    if (mustRefine_(R2.updaters, C.updaters))
       break;
-    }
-  }
-  return didCoarsen;
-}
 
-template <class Factory, class Updaters, class ErrorNorm>
-void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::refine() {
-  while (true) {
-    UpdatersWithCount F1{ current_.clone() };
-    step(F1, relativeTime_, relativeTau_ / 2.0);
-    UpdatersWithCount F2{ F1.updaters.clone() };
-    step(F2, relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
+    didCoarsen = true;
+    relativeTau_ *= 2;
+    R1_ = C;
+  }
+  UpdatersWithCount F1;
+  UpdatersWithCount F2;
+  if (!didCoarsen) {
+    while (true) {
+      F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
+      F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
+                relativeTau_ / 2.0);
+      if (!mustRefine_(F2.updaters, R1_.updaters))
+        break;
 
-    if (!mustRefine_(R1_.updaters, F2.updaters)) {
-      break;
-    } else {
-      R1_ = F1;
-      R2_ = F2;
       relativeTau_ /= 2.0;
+      R1_ = F1;
+      R2 = F2;
     }
   }
-}
 
-template <class Factory, class Updaters, class ErrorNorm>
-IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
-  if (R2_.updaters == Updaters()) {
-    R1_ = { current_.clone() };
-    step(R1_, relativeTime_, relativeTau_); // counting again upon restart
-  } else {
-    R1_ = R2_;
-  }
-
-  iterationRegister_.reset();
-  if (!coarsen())
-    refine();
-
-  current_ = R1_.updaters;
   iterationRegister_.registerFinalCount(R1_.count);
   relativeTime_ += relativeTau_;
+  current_ = R1_.updaters;
+  R1_ = R2;
 
   return iterationRegister_;
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
-    UpdatersWithCount &updatersAndCount, double rTime, double rTau) {
-  updatersAndCount.count =
+typename AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
+    Updaters const &oldUpdaters, double rTime, double rTau) {
+  UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
+  newUpdatersAndCount.count =
       MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
-                           updatersAndCount.updaters, errorNorm_,
-                           externalForces_).step(rTime, rTau);
-  iterationRegister_.registerCount(updatersAndCount.count);
+                           newUpdatersAndCount.updaters, errorNorm_,
+                           externalForces_)
+          .step(rTime, rTau);
+  iterationRegister_.registerCount(newUpdatersAndCount.count);
+  return newUpdatersAndCount;
 }
 
 #include "adaptivetimestepper_tmpl.cc"

src/time-stepping/adaptivetimestepper.hh

@@ -38,15 +38,14 @@ class AdaptiveTimeStepper {
                       std::function<bool(Updaters &, Updaters &)> mustRefine);
 
   bool reachedEnd();
-  bool coarsen();
-  void refine();
   IterationRegister advance();
 
   double relativeTime_;
   double relativeTau_;
 
 private:
-  void step(UpdatersWithCount &updatersWithCount, double rTime, double rTau);
+  UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
+                         double rTau);
 
   double finalTime_;
   Factory &factory_;
@@ -54,7 +53,6 @@ class AdaptiveTimeStepper {
   std::shared_ptr<Nonlinearity> globalFriction_;
   Updaters &current_;
   UpdatersWithCount R1_;
-  UpdatersWithCount R2_;
   std::function<void(double, Vector &)> externalForces_;
   std::function<bool(Updaters &, Updaters &)> mustRefine_;
   ErrorNorm const &errorNorm_;

src/time-stepping/state/ageinglawstateupdater.cc

+#include <cmath>
+
 #include "ageinglawstateupdater.hh"
 #include "../../tobool.hh"
 

src/time-stepping/state/sliplawstateupdater.cc

+#include <cmath>
+
 #include "sliplawstateupdater.hh"
 #include "../../tobool.hh"
 

src/uniform-grid-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cmath>
+#include <exception>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+// #include <dune/common/parametertree.hh>
+// #include <dune/common/parametertreeparser.hh>
+
+#include "assemblers.hh"
+#include "diameter.hh"
+#include "gridselector.hh"
+#include "one-body-problem-data/mygrid.hh"
+#include "vtk.hh"
+
+size_t const dims = MY_DIM;
+size_t const refinements = 5; // FIXME?
+
+int main(int argc, char *argv[]) {
+  try {
+    // Dune::ParameterTree parset;
+    // Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+
+    using GridView = Grid::LeafGridView;
+    using MyAssembler = MyAssembler<GridView, dims>;
+
+    GridConstructor<Grid> gridConstructor;
+    auto grid = gridConstructor.getGrid();
+
+    // refine uniformly!
+    for (size_t refinement = 0; refinement < refinements; ++refinement)
+      grid->globalRefine(1);
+
+    double minDiameter = std::numeric_limits<double>::infinity();
+    double maxDiameter = 0.0;
+    for (auto &&e : elements(grid->leafGridView())) {
+      auto const geometry = e.geometry();
+      auto const diam = diameter(geometry);
+      minDiameter = std::min(minDiameter, diam);
+      maxDiameter = std::max(maxDiameter, diam);
+    }
+    std::cout << "min diameter: " << minDiameter << std::endl;
+    std::cout << "max diameter: " << maxDiameter << std::endl;
+
+    auto const leafView = grid->leafGridView();
+    auto const leafVertexCount = leafView.size(dims);
+
+    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
+
+    MyAssembler const myAssembler(leafView);
+    MyVTKWriter<typename MyAssembler::VertexBasis,
+                typename MyAssembler::CellBasis> const
+        vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+    vtkWriter.writeGrid();
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

src/vtk.cc

@@ -46,4 +46,13 @@ void MyVTKWriter<VertexBasis, CellBasis>::write(
   writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
+template <class VertexBasis, class CellBasis>
+void MyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
+  Dune::VTKWriter<typename VertexBasis::GridView> writer(
+      vertexBasis.getGridView());
+
+  std::string const filename = prefix + "_grid";
+  writer.write(filename.c_str(), Dune::VTK::appendedraw);
+}
+
 #include "vtk_tmpl.cc"

src/vtk.hh

@@ -15,5 +15,7 @@ template <class VertexBasis, class CellBasis> class MyVTKWriter {
   template <class Vector, class ScalarVector>
   void write(size_t record, Vector const &u, Vector const &v,
              ScalarVector const &alpha, ScalarVector const &stress) const;
+
+  void writeGrid() const;
 };
 #endif