Fix subcycling and checkpointing for two-scale-heat-conduction participant macro-dumux (#666)

* [macro-dumux] Handle subcycling in macro-dumux. Fix checkpointing.

* [macro-dumux] Warn about different dt after solve

* [macro-dumux] Negotiate dt before solve

* [macro-dumux] Advance DuMux by default after solve

* [macro-dumux] Cleanup time loop

* Add changelog entry

* Formatting and making code comment more informative

---------

Co-authored-by: Ishaan Desai <[email protected]>

Author: mathiskelm

changelog-entries/666.md

@@ -0,0 +1 @@
+- Fix subcycling and checkpointing in macro-dumux of two-scale-heat-conduction [#666](https://github.com/precice/tutorials/pull/666)

two-scale-heat-conduction/macro-dumux/appl/main.cc

@@ -100,7 +100,9 @@ int main(int argc, char **argv)
 
   auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance();
 
-  if (getParam<bool>("Precice.RunWithCoupling") == true) {
+  const auto runWithCoupling = getParam<bool>("Precice.RunWithCoupling");
+
+  if (runWithCoupling) {
     couplingParticipant.announceSolver("macro-heat", preciceConfigFilename,
                                        mpiHelper.rank(), mpiHelper.size());
     // verify that dimensions match
@@ -132,13 +134,13 @@ int main(int argc, char **argv)
       std::cout << " ;" << std::endl;
     }
   }
+
   std::cout << "Number of Coupled Cells:" << coupledElementIdxs.size()
             << std::endl;
 
-  // initialize preCICE
   auto numberOfElements =
       coords.size() / couplingParticipant.getMeshDimensions(meshName);
-  if (getParam<bool>("Precice.RunWithCoupling") == true) {
+  if (runWithCoupling) {
     couplingParticipant.setMesh(meshName, coords);
 
     // couples between dumux element indices and preciceIndices;
@@ -154,7 +156,7 @@ int main(int argc, char **argv)
   const std::string writeDataConcentration = "concentration";
   // const std::string writeDataTemperature = "temperature";
 
-  if (getParam<bool>("Precice.RunWithCoupling") == true) {
+  if (runWithCoupling) {
     couplingParticipant.announceQuantity(meshName, readDatak00);
     couplingParticipant.announceQuantity(meshName, readDatak01);
     couplingParticipant.announceQuantity(meshName, readDatak10);
@@ -171,13 +173,17 @@ int main(int argc, char **argv)
   problem->applyInitialSolution(x);
   auto xOld = x;
 
+  auto   xCheckpoint        = x;
+  double timeCheckpoint     = 0.0;
+  int    timeStepCheckpoint = 0;
+
   // initialize the coupling data
   std::vector<double> temperatures;
   for (int solIdx = 0; solIdx < numberOfElements; ++solIdx) {
     temperatures.push_back(x[solIdx][problem->returnTemperatureIdx()]);
   };
 
-  if (getParam<bool>("Precice.RunWithCoupling") == true) {
+  if (runWithCoupling) {
     couplingParticipant.writeQuantityVector(meshName, writeDataConcentration,
                                             temperatures);
     if (couplingParticipant
@@ -203,7 +209,7 @@ int main(int argc, char **argv)
   auto gridVariables  = std::make_shared<GridVariables>(problem, gridGeometry);
   gridVariables->init(x);
 
-  // intialize the vtk output module
+  // initialize the vtk output module
   using IOFields = GetPropType<TypeTag, Properties::IOFields>;
   VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x,
                                                            problem->name());
@@ -220,7 +226,7 @@ int main(int argc, char **argv)
   // output every vtkOutputInterval time step
   const int vtkOutputInterval = getParam<int>("TimeLoop.OutputInterval");
 
-  // initialize
+  // initialize preCICE
   couplingParticipant.initialize();
 
   // time loop parameters
@@ -229,7 +235,7 @@ int main(int argc, char **argv)
   double     solverDt;
   double     dt;
 
-  if (getParam<bool>("Precice.RunWithCoupling") == true) {
+  if (runWithCoupling) {
     solverDt = getParam<Scalar>("TimeLoop.InitialDt");
     dt       = std::min(preciceDt, solverDt);
   } else {
@@ -261,16 +267,25 @@ int main(int argc, char **argv)
   std::cout << "Time Loop starts" << std::endl;
   timeLoop->start();
   do {
-    if (getParam<bool>("Precice.RunWithCoupling") == true) {
+    if (runWithCoupling) {
       if (couplingParticipant.isCouplingOngoing() == false)
         break;
 
       // write checkpoint
       if (couplingParticipant.requiresToWriteCheckpoint()) {
-        xOld = x;
+        xCheckpoint        = x;
+        timeCheckpoint     = timeLoop->time();
+        timeStepCheckpoint = timeLoop->timeStepIndex();
       }
 
+      preciceDt = couplingParticipant.getMaxTimeStepSize();
+      solverDt  = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()),
+                           timeLoop->maxTimeStepSize());
+      dt        = std::min(preciceDt, solverDt);
+
       // read porosity and conductivity data from other solver
+      // TODO: data needs to be updated if Newton solver adapts time-step size
+      // and coupling data is interpolated in time
       couplingParticipant.readQuantityFromOtherSolver(meshName, readDatak00,
                                                       dt);
       couplingParticipant.readQuantityFromOtherSolver(meshName, readDatak01,
@@ -281,89 +296,79 @@ int main(int argc, char **argv)
                                                       dt);
       couplingParticipant.readQuantityFromOtherSolver(meshName,
                                                       readDataPorosity, dt);
-      // store coupling data in problem
+      // store coupling data in spatial params, exchange with MPI
       problem->spatialParams().updateCouplingData();
+    } else {
+      dt = std::min(
+          nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()),
+          timeLoop->maxTimeStepSize());
     }
-    std::cout << "Solver starts" << std::endl;
+    // set new dt as suggested by the Newton solver or by preCICE
+    timeLoop->setTimeStepSize(dt);
+
+    std::cout << "Solver starts with target dt: " << dt << std::endl;
 
     // linearize & solve
     nonLinearSolver.solve(x, *timeLoop);
 
-    int solIdx = 0;
-    for (const auto &element : elements(leafGridView, Dune::Partitions::interior)) {
-      auto fvGeometry = localView(*gridGeometry);
-      fvGeometry.bindElement(element);
-      for (const auto &scv : scvs(fvGeometry)) {
-        temperatures[solIdx++] =
-            x[scv.elementIndex()][problem->returnTemperatureIdx()];
-      }
+    // save dt value that was actually used by the non-linear solver
+    dt = timeLoop->timeStepSize();
+
+    // DuMux advance + report
+    gridVariables->advanceTimeStep();
+    timeLoop->advanceTimeStep();
+    timeLoop->reportTimeStep();
+    xOld = x;
+
+    // Vtk output
+    // TODO: output interval does not work seamlessly when subcycling
+    n += 1;
+    if (n == vtkOutputInterval) {
+      problem->updateVtkOutput(x);
+      vtkWriter.write(timeLoop->time());
+      n = 0;
     }
 
-    if (getParam<bool>("Precice.RunWithCoupling") == true) {
+    if (runWithCoupling) {
+      int solIdx = 0;
+      for (const auto &element : elements(leafGridView, Dune::Partitions::interior)) {
+        auto fvGeometry = localView(*gridGeometry);
+        fvGeometry.bindElement(element);
+        for (const auto &scv : scvs(fvGeometry)) {
+          temperatures[solIdx++] =
+              x[scv.elementIndex()][problem->returnTemperatureIdx()];
+        }
+      }
+
       couplingParticipant.writeQuantityVector(meshName,
                                               writeDataConcentration, temperatures);
       couplingParticipant.writeQuantityToOtherSolver(meshName,
                                                      writeDataConcentration);
-    }
-
-    // advance precice
-    if (getParam<bool>("Precice.RunWithCoupling") == true) {
-      couplingParticipant.advance(timeLoop->timeStepSize());
-
-      preciceDt = couplingParticipant.getMaxTimeStepSize();
-      solverDt  = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), timeLoop->maxTimeStepSize());
-      dt        = std::min(preciceDt, solverDt);
 
+      // advance preCICE
       if ((!fabs(preciceDt - dt)) < 1e-14) {
-        std::cout << "dt from preCICE is different than dt from DuMuX. "
-                     "preCICE dt = "
-                  << preciceDt << " and DuMuX dt =" << solverDt << std::endl;
+        std::cout << "dt from preCICE is different than dt from DuMuX."
+                  << " preCICE dt = " << preciceDt
+                  << " and DuMuX dt = " << solverDt
+                  << " resulted in dt = " << dt
+                  << std::endl;
       }
-    } else
-      dt = std::min(
-          nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()),
-          getParam<Scalar>("TimeLoop.MaxDt"));
+      std::flush(std::cout);
+      couplingParticipant.advance(dt);
 
-    std::cout << "Using dt: " << dt << std::endl;
-
-    if (getParam<bool>("Precice.RunWithCoupling") == true) {
+      // reset to checkpoint if not converged
       if (couplingParticipant.requiresToReadCheckpoint()) {
-        // make the new solution the old solution
-        x = xOld;
+        x    = xCheckpoint;
+        xOld = x;
+        timeLoop->setTime(timeCheckpoint, timeStepCheckpoint);
+
+        // TODO: previousTimeStep might be more appropriate, last one could be small
+        timeLoop->setTimeStepSize(dt);
         gridVariables->update(x);
         gridVariables->advanceTimeStep();
-      } else // coupling successful
-      {
-        n += 1;
-        if (n == vtkOutputInterval) {
-          problem->updateVtkOutput(x);
-          vtkWriter.write(timeLoop->time());
-          n = 0;
-        }
-        gridVariables->advanceTimeStep();
-        // advance the time loop to the next step
-        timeLoop->advanceTimeStep();
-        // report statistics of this time step
-        timeLoop->reportTimeStep();
-      }
-    } else {
-      xOld = x;
-      gridVariables->advanceTimeStep();
-      // advance the time loop to the next step
-      timeLoop->advanceTimeStep();
-      // report statistics of this time step
-      timeLoop->reportTimeStep();
-
-      // output every outputinterval steps
-      n += 1;
-      if (n == vtkOutputInterval) {
-        problem->updateVtkOutput(x);
-        vtkWriter.write(timeLoop->time());
-        n = 0;
+        continue;
       }
     }
-    // set new dt as suggested by the newton solver or by preCICE
-    timeLoop->setTimeStepSize(dt);
 
     std::cout << "Time: " << timeLoop->time() << std::endl;
 
@@ -374,7 +379,7 @@ int main(int argc, char **argv)
   ////////////////////////////////////////////////////////////
   // finalize, print dumux message to say goodbye
   ////////////////////////////////////////////////////////////
-  if (getParam<bool>("Precice.RunWithCoupling") == true) {
+  if (runWithCoupling) {
     couplingParticipant.finalize();
   }
   // print dumux end message

two-scale-heat-conduction/macro-dumux/params.input

@@ -14,9 +14,9 @@ UpperRight = 1 0.5
 Cells = 8 4
 
 [TimeLoop]
-TEnd = 0.25 
+TEnd = 0.025
 OutputInterval = 1 # output after every n timesteps 
-MaxDt = 0.01 
+MaxDt = 0.004
 InitialDt = 0.01
 
 [Problem]

two-scale-heat-conduction/precice-config.xml

@@ -63,7 +63,7 @@
 
   <coupling-scheme:serial-implicit>
     <participants first="macro-heat" second="Micro-Manager" />
-    <max-time value="0.02" />
+    <max-time value="0.025" />
     <time-window-size value="1.0e-2" />
     <max-iterations value="30" />
     <exchange data="k_00" mesh="macro-mesh" from="Micro-Manager" to="macro-heat" />