multi threading steps in adaptivetimestepper

dune/tectonic/time-stepping/CMakeLists.txt

@@ -10,6 +10,8 @@ add_custom_target(tectonic_dune_time-stepping SOURCES
   rate.cc
   state.hh
   state.cc
+  step.hh
+  stepbase.hh
   updaters.hh
 )
 
@@ -19,5 +21,7 @@ install(FILES
   coupledtimestepper.hh
   rate.hh
   state.hh
+  step.hh
+  stepbase.hh
   updaters.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/time-stepping/adaptivetimestepper.hh

@@ -2,9 +2,14 @@
 #define SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
 
 #include <fstream>
+#include <future>
+#include <thread>
 
 //#include "../multi-threading/task.hh"
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
+
 #include "coupledtimestepper.hh"
+#include "stepbase.hh"
 
 struct IterationRegister {
   void registerCount(FixedPointIterationCounter count);
@@ -16,20 +21,23 @@ struct IterationRegister {
   FixedPointIterationCounter finalCount;
 };
 
-template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
-class AdaptiveTimeStepper {
+template <class Updaters>
 struct UpdatersWithCount {
   Updaters updaters;
   FixedPointIterationCounter count;
 };
 
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class AdaptiveTimeStepper {
+  using UpdatersWithCount = UpdatersWithCount<Updaters>;
+
+  using StepBase = StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>;
+
+  using NBodyAssembler = typename ContactNetwork::NBodyAssembler;
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
-  using IgnoreVector = typename Factory::BitVector;
-  //using ConvexProblem = typename Factory::ConvexProblem;
-  //using Nonlinearity = typename Factory::Nonlinearity;
 
-  using NBodyAssembler = typename ContactNetwork::NBodyAssembler;
+  using IgnoreVector = typename Factory::BitVector;
 
   using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>;
 
@@ -38,23 +46,15 @@ class AdaptiveTimeStepper {
   using ExternalForces = typename MyCoupledTimeStepper::ExternalForces;
 
 public:
-  AdaptiveTimeStepper(
-                      Dune::ParameterTree const &parset,
+  AdaptiveTimeStepper(const StepBase& stepBase,
                       ContactNetwork& contactNetwork,
-                      const IgnoreVector& ignoreNodes,
-                      GlobalFriction& globalFriction,
-                      const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
                       Updaters &current,
                       double relativeTime,
                       double relativeTau,
-                      ExternalForces& externalForces,
-                      const ErrorNorms& errorNorms,
                       std::function<bool(Updaters &, Updaters &)> mustRefine);
 
   bool reachedEnd();
 
-  auto makeLinearSolver() const;
-
   IterationRegister advance();
 
   double relativeTime_;
@@ -65,25 +65,19 @@ class AdaptiveTimeStepper {
 
   NBodyAssembler step(const NBodyAssembler& oldNBodyAssembler) const;
 
+  int refine(UpdatersWithCount& R2);
 
-  template <class LinearSolver>
-  UpdatersWithCount step(const Updaters& oldUpdaters, const NBodyAssembler& nBodyAssembler,
-                         std::shared_ptr<LinearSolver>& linearSolver,
-                         double rTime, double rTau);
+  int coarsen();
 
-  double finalTime_;
-  Dune::ParameterTree const &parset_;
-  ContactNetwork& contactNetwork_;
-  const IgnoreVector& ignoreNodes_;
+  int determineStrategy(UpdatersWithCount& R2);
 
-  GlobalFriction& globalFriction_;
-  const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+  const StepBase& stepBase_;
+  ContactNetwork& contactNetwork_;
 
   Updaters &current_;
   UpdatersWithCount R1_;
-  ExternalForces& externalForces_;
+
   std::function<bool(Updaters &, Updaters &)> mustRefine_;
-  const ErrorNorms& errorNorms_;
 
   IterationRegister iterationRegister_;
 };

dune/tectonic/time-stepping/adaptivetimestepper_tmpl.cc

@@ -2,6 +2,9 @@
 #error MY_DIM unset
 #endif
 
+#include <future>
+#include <thread>
+
 #include "../explicitgrid.hh"
 #include "../explicitvectors.hh"
 
@@ -37,6 +40,7 @@ using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
 
 template class AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
 
-template typename AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::UpdatersWithCount
+/*
+template std::packaged_task<typename AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::UpdatersWithCount()>
 AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::step<LinearSolver>(
-        const MyUpdaters&, const MyNBodyAssembler&, std::shared_ptr<LinearSolver>&, double, double);
+        const MyUpdaters&, const MyNBodyAssembler&, std::shared_ptr<LinearSolver>&, double, double); */

dune/tectonic/time-stepping/step.hh

+#ifndef DUNE_TECTONIC_TIME_STEPPING_STEP_HH
+#define DUNE_TECTONIC_TIME_STEPPING_STEP_HH
+
+
+#include <future>
+#include <thread>
+#include <chrono>
+#include <functional>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "stepbase.hh"
+#include "adaptivetimestepper.hh"
+
+template<class IterationStepType, class NormType, class ReductionFactorContainer>
+Dune::Solvers::Criterion reductionFactorCriterion(
+      IterationStepType& iterationStep,
+      const NormType& norm,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+        double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+
+        if (convRate>1.0)
+            std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
+
+template<class IterationStepType, class Functional, class ReductionFactorContainer>
+Dune::Solvers::Criterion energyCriterion(
+      const IterationStepType& iterationStep,
+      const Functional& f,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = std::abs(f(*lastIterate) - f(*iterationStep.getIterate())); //norm.diff(*lastIterate, *iterationStep.getIterate());
+
+        double convRate = (normOfOldCorrection != 0.0) ? 1.0 - (normOfCorrection / normOfOldCorrection) : 0.0;
+
+        if (convRate>1.0)
+            std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
+template <class ReductionFactorContainer>
+void updateReductionFactors(ReductionFactorContainer& reductionFactors) {
+    const size_t s = reductionFactors.size();
+
+    //print(reductionFactors, "reduction factors: ");
+
+    if (s>allReductionFactors.size()) {
+        allReductionFactors.resize(s);
+    }
+
+    for (size_t i=0; i<reductionFactors.size(); i++) {
+        allReductionFactors[i].push_back(reductionFactors[i]);
+    }
+
+    reductionFactors.clear();
+}
+
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class Step : protected StepBase<Factory, ContactNetwork, Updaters, ErrorNorms> {
+public:
+    enum Mode { sameThread, newThread };
+
+private:
+    using Base = StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>;
+    using NBodyAssembler = typename ContactNetwork::NBodyAssembler;
+    using UpdatersWithCount = UpdatersWithCount<Updaters>;
+
+    const Updaters& oldUpdaters_;
+    const NBodyAssembler& nBodyAssembler_;
+
+    double relativeTime_;
+    double relativeTau_;
+
+    IterationRegister& iterationRegister_;
+
+    std::packaged_task<UpdatersWithCount()> task_;
+    std::future<UpdatersWithCount> future_;
+    std::thread thread_;
+
+    Mode mode_;
+
+public:
+    Step(const Base& stepFactory, const Updaters& oldUpdaters, const NBodyAssembler& nBodyAssembler, double rTime, double rTau, IterationRegister& iterationRegister) :
+        Base(stepFactory.parset_, stepFactory.contactNetwork_, stepFactory.ignoreNodes_, stepFactory.globalFriction_,
+             stepFactory.bodywiseNonmortarBoundaries_, stepFactory.externalForces_, stepFactory.errorNorms_),
+        oldUpdaters_(oldUpdaters),
+        nBodyAssembler_(nBodyAssembler),
+        relativeTime_(rTime),
+        relativeTau_(rTau),
+        iterationRegister_(iterationRegister) {}
+
+    UpdatersWithCount doStep() {
+        // make linear solver for linear correction in TNNMGStep
+        using Vector = typename Factory::Vector;
+        using Matrix = typename Factory::Matrix;
+
+        /* old, pre multi threading, was unused
+        using Norm =  EnergyNorm<Matrix, Vector>;
+        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        const auto& preconditionerParset = parset_.sub("solver.tnnmg.linear.preconditioner");
+
+        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+        preconditioner.build();
+
+        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+        cgStep->setPreconditioner(preconditioner);
+
+        Norm norm(*cgStep);
+
+        return std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.linear.tolerance"), norm, Solver::QUIET);
+        */
+
+        // patch preconditioner only needs to be computed once per advance()
+        // make linear solver for linear correction in TNNMGStep
+        using Norm =  EnergyNorm<Matrix, Vector>;
+        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        /*const auto& preconditionerParset = parset_.sub("solver.tnnmg.preconditioner");
+
+        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+        preconditioner.build();
+
+        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+        cgStep->setPreconditioner(preconditioner);
+
+        Norm norm(*cgStep);
+
+        auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+        */
+
+        // set multigrid solver
+        auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+        // transfer operators need to be recomputed on change due to setDeformation()
+        using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+        using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+        TransferOperators transfer(this->contactNetwork_.nLevels()-1);
+        for (size_t i=0; i<transfer.size(); i++) {
+            transfer[i] = std::make_shared<TransferOperator>();
+            transfer[i]->setup(this->contactNetwork_, i, i+1);
+        }
+
+        // Remove any recompute filed so that initially the full transferoperator is assembled
+        for (size_t i=0; i<transfer.size(); i++)
+            std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+        auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+        linearMultigridStep->setMGType(1, 3, 3);
+        linearMultigridStep->setSmoother(smoother);
+        linearMultigridStep->setTransferOperators(transfer);
+
+        Norm norm(*linearMultigridStep);
+        auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, this->parset_.template get<int>("solver.tnnmg.main.multi"), this->parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+
+        Vector x;
+        x.resize(nBodyAssembler_.totalHasObstacle_.size());
+        x = 0;
+
+        linearSolver->getIterationStep().setProblem(x);
+      //dynamic_cast<Dune::Solvers::IterationStep<Vector>*>(linearMultigridStep.get())->setProblem(x);
+      //dynamic_cast<Dune::Solvers::IterationStep<Vector>*>(cgStep.get())->setProblem(x);
+
+      //std::vector<double> reductionFactors;
+      //linearSolver->addCriterion(reductionFactorCriterion(linearSolver->getIterationStep(), linearSolver->getErrorNorm(), reductionFactors));
+      //linearSolver->addCriterion(reductionFactorCriterion(*cgStep, norm, reductionFactors));
+
+      UpdatersWithCount newUpdatersAndCount = {oldUpdaters_.clone(), {}};
+
+      typename Base::MyCoupledTimeStepper coupledTimeStepper(
+        this->finalTime_, this->parset_, nBodyAssembler_,
+        this->ignoreNodes_, this->globalFriction_, this->bodywiseNonmortarBoundaries_,
+        newUpdatersAndCount.updaters, this->errorNorms_, this->externalForces_);
+
+      newUpdatersAndCount.count = coupledTimeStepper.step(linearSolver, relativeTime_, relativeTau_);
+      iterationRegister_.registerCount(newUpdatersAndCount.count);
+
+      //updateReductionFactors(reductionFactors);
+
+      return newUpdatersAndCount;
+    }
+
+   /* auto simple = [&]{
+        std::cout << "starting task... " << std::endl;
+        UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
+        std::this_thread::sleep_for(std::chrono::milliseconds(10000));
+        std::cout << "finishing task... " << std::endl;
+        return newUpdatersAndCount;
+    }*/
+
+    void run(Mode mode = sameThread) {
+        mode_ = mode;
+        task_ = std::packaged_task<UpdatersWithCount()>( [this]{ return doStep(); });
+        future_ = task_.get_future();
+
+        if (mode == Mode::sameThread) {
+            task_();
+        } else {
+            thread_ = std::thread(std::move(task_));
+        }
+    }
+
+    auto get() {
+        std::cout << "Step::get() called" << std::endl;
+        if (thread_.joinable()) {
+            thread_.join();
+            std::cout << "Thread joined " << std::endl;
+        }
+        return future_.get();
+    }
+
+    ~Step() {
+        if (thread_.joinable()) {
+            thread_.join();
+            std::cout << "Destructor:: Thread joined " << std::endl;
+        }
+    }
+};
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+Step<Factory, ContactNetwork, Updaters, ErrorNorms>
+buildStep(const StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>& base,
+          const Updaters& oldUpdaters,
+          const typename ContactNetwork::NBodyAssembler& nBodyAssembler,
+          double rTime, double rTau) {
+
+    return Step(base, oldUpdaters, nBodyAssembler, rTime, rTau);
+}
+
+#endif

dune/tectonic/time-stepping/stepbase.hh

+#ifndef DUNE_TECTONIC_TIME_STEPPING_STEPBASE_HH
+#define DUNE_TECTONIC_TIME_STEPPING_STEPBASE_HH
+
+#include "coupledtimestepper.hh"
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class StepBase {
+protected:
+  using NBodyAssembler = typename ContactNetwork::NBodyAssembler;
+  using IgnoreVector = typename Factory::BitVector;
+
+  using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+
+  using GlobalFriction = typename MyCoupledTimeStepper::GlobalFriction;
+  using BitVector = typename MyCoupledTimeStepper::BitVector;
+  using ExternalForces = typename MyCoupledTimeStepper::ExternalForces;
+
+public:
+  StepBase(
+    Dune::ParameterTree const &parset,
+    ContactNetwork& contactNetwork,
+    const IgnoreVector& ignoreNodes,
+    GlobalFriction& globalFriction,
+    const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+    ExternalForces& externalForces,
+    const ErrorNorms& errorNorms) :
+        parset_(parset),
+        finalTime_(parset_.get<double>("problem.finalTime")),
+        contactNetwork_(contactNetwork),
+        ignoreNodes_(ignoreNodes),
+        globalFriction_(globalFriction),
+        bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
+        externalForces_(externalForces),
+        errorNorms_(errorNorms) {}
+
+  Dune::ParameterTree const &parset_;
+  double finalTime_;
+
+  ContactNetwork& contactNetwork_;
+  const IgnoreVector& ignoreNodes_;
+
+  GlobalFriction& globalFriction_;
+  const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
+  ExternalForces& externalForces_;
+  const ErrorNorms& errorNorms_;
+};
+
+#endif

src/multi-body-problem/multi-body-problem.cc

@@ -74,6 +74,7 @@
 #include <dune/tectonic/time-stepping/adaptivetimestepper.hh>
 #include <dune/tectonic/time-stepping/rate.hh>
 #include <dune/tectonic/time-stepping/state.hh>
+#include <dune/tectonic/time-stepping/stepbase.hh>
 #include <dune/tectonic/time-stepping/updaters.hh>
 
 #include <dune/tectonic/utils/debugutils.hh>
@@ -341,10 +342,14 @@ int main(int argc, char *argv[]) {
     typename ContactNetwork::ExternalForces externalForces;
     contactNetwork.externalForces(externalForces);
 
+    StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
+                 externalForces, stateEnergyNorms);
+
     AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
-        adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
+        adaptiveTimeStepper(stepBase, contactNetwork, current,
                             programState.relativeTime, programState.relativeTau,
-                            externalForces, stateEnergyNorms, mustRefine);
+                            mustRefine);
 
     size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
 

src/strikeslip/strikeslip-2D.cfg

 # -*- mode:conf -*-
 [body0]
-smallestDiameter = 0.05 # 2e-3 [m]
+smallestDiameter = 0.02 # 2e-3 [m]
 
 [body1]
-smallestDiameter = 0.05  # 2e-3 [m]
+smallestDiameter = 0.02  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 5e-6 # 1e-5

src/strikeslip/strikeslip.cc

@@ -56,6 +56,8 @@
 
 #include <dune/tectonic/factories/strikeslipfactory.hh>
 
+
+#include <dune/tectonic/io/io-handler.hh>
 #include <dune/tectonic/io/hdf5-writer.hh>
 #include <dune/tectonic/io/hdf5/restart-io.hh>
 #include <dune/tectonic/io/vtk.hh>
@@ -167,29 +169,13 @@ int main(int argc, char *argv[]) {
 
     using MyProgramState = ProgramState<Vector, ScalarVector>;
     MyProgramState programState(nVertices);
-    auto const firstRestart = parset.get<size_t>("io.restarts.first");
-    auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
-    auto const writeRestarts = parset.get<bool>("io.restarts.write");
-    auto const writeData = parset.get<bool>("io.data.write");
-    bool const handleRestarts = writeRestarts or firstRestart > 0;
-
-    auto restartFile = handleRestarts
-                           ? std::make_unique<HDF5::File>(
-                                 "restarts.h5",
-                                 writeRestarts ? HDF5::Access::READWRITE
-                                               : HDF5::Access::READONLY)
-                           : nullptr;
-
-
-    auto restartIO = handleRestarts
-                         ? std::make_unique<RestartIO<MyProgramState>>(
-                               *restartFile, nVertices)
-                         : nullptr;
-
-    if (firstRestart > 0) // automatically adjusts the time and timestep
-      restartIO->read(firstRestart, programState);
-    else
+
+    IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
+
+    bool restartRead = ioHandler.read(programState);
+    if (!restartRead) {
       programState.setupInitialConditions(parset, contactNetwork);
+    }
 
     auto& nBodyAssembler = contactNetwork.nBodyAssembler();
     for (size_t i=0; i<bodyCount; i++) {
@@ -206,90 +192,9 @@ int main(int argc, char *argv[]) {
     auto& globalFriction = contactNetwork.globalFriction();
     globalFriction.updateAlpha(programState.alpha);
 
-    using MyVertexBasis = typename Assembler::VertexBasis;
-    using MyCellBasis = typename Assembler::CellBasis;
-    std::vector<Vector> vertexCoordinates(bodyCount);
-    std::vector<const MyVertexBasis* > vertexBases(bodyCount);
-    std::vector<const MyCellBasis* > cellBases(bodyCount);
-
-    auto& wPatches = blocksFactory.weakPatches();
-    std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
-
-
-    for (size_t i=0; i<bodyCount; i++) {
-      const auto& body = contactNetwork.body(i);
-      vertexBases[i] = &(body->assembler()->vertexBasis);
-      cellBases[i] = &(body->assembler()->cellBasis);
-
-      weakPatches[i].resize(1);
-      weakPatches[i][0] = wPatches[i].get();
-
-      auto& vertexCoords = vertexCoordinates[i];
-      vertexCoords.resize(nVertices[i]);
-
-      auto hostLeafView = body->grid()->hostGrid().leafGridView();
-      Dune::MultipleCodimMultipleGeomTypeMapper<
-          LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
-      for (auto &&v : vertices(hostLeafView))
-        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
-    }
-
-    typename ContactNetwork::BoundaryPatches frictionBoundaries;
-    contactNetwork.boundaryPatches("friction", frictionBoundaries);
-
-    std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
-    contactNetwork.frictionPatches(frictionPatches);
-
-
     IterationRegister iterationCount;
 
-    auto const report = [&](bool initial = false) {
-        auto dataFile =
-            writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
-
-        auto dataWriter =
-            writeData ? std::make_unique<
-                            HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
-                            *dataFile, vertexCoordinates, vertexBases,
-                            frictionPatches, weakPatches)
-                      : nullptr;
-
-      if (writeData) {
-        dataWriter->reportSolution(programState, contactNetwork, globalFriction);
-        if (!initial)
-          dataWriter->reportIterations(programState, iterationCount);
-        dataFile->flush();
-      }
-
-      if (writeRestarts and !initial and
-          programState.timeStep % restartSpacing == 0) {
-        restartIO->write(programState);
-        restartFile->flush();
-      }
-
-      if (parset.get<bool>("io.printProgress"))
-        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
-                  << ", time = " << std::setw(12) << programState.relativeTime
-                  << ", tau = " << std::setw(12) << programState.relativeTau
-                  << std::endl;
-
-      if (parset.get<bool>("io.vtk.write")) {
-        std::vector<ScalarVector> stress(bodyCount);
-
-        for (size_t i=0; i<bodyCount; i++) {
-          const auto& body = contactNetwork.body(i);
-          body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
-                                           body->data()->getPoissonRatio(),
-                                           programState.u[i], stress[i]);
-
-        }
-
-        const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "../debug_print/vtk/");
-        vtkWriter.write(programState.timeStep, programState.u, programState.v,
-                        programState.alpha, stress);
-      }
-    };
-    report(true);
+    ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, true);
 
     // -------------------
     // Set up TNNMG solver
@@ -428,10 +333,14 @@ int main(int argc, char *argv[]) {
     typename ContactNetwork::ExternalForces externalForces;
     contactNetwork.externalForces(externalForces);
 
+    StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
+                 externalForces, stateEnergyNorms);
+
     AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
-        adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
+        adaptiveTimeStepper(stepBase, contactNetwork, current,
                             programState.relativeTime, programState.relativeTau,
-                            externalForces, stateEnergyNorms, mustRefine);
+                            mustRefine);
 
     size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
 
@@ -456,7 +365,7 @@ int main(int argc, char *argv[]) {
 
       contactNetwork.setDeformation(programState.u);
 
-      report();
+      ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, false);
 
       if (programState.timeStep==timeSteps) {
         std::cout << "limit of timeSteps reached!" << std::endl;

src/strikeslip/strikeslip.cfg

@@ -55,11 +55,11 @@ sigmaN          = 200.0      # [Pa]
 finalVelocity   = 1e-4     # [m/s]
 
 [io]
-data.write      = true
+data.write      = false
 printProgress   = true
 restarts.first  = 0
-restarts.spacing= 20
-restarts.write  = false #true
+restarts.spacing= 1
+restarts.write  = true #true
 vtk.write       = true
 
 [problem]
@@ -73,7 +73,7 @@ relativeTau = 2e-4 # 1e-6
 
 [timeSteps]
 scheme = newmark
-timeSteps = 5000
+timeSteps = 5
 
 [u0.solver]
 maximumIterations = 100