diff --git a/changelog-entries/666.md b/changelog-entries/666.md new file mode 100644 index 000000000..069cca73c --- /dev/null +++ b/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) diff --git a/two-scale-heat-conduction/macro-dumux/appl/main.cc b/two-scale-heat-conduction/macro-dumux/appl/main.cc index bbcf54310..d36090d31 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/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("Precice.RunWithCoupling") == true) { + const auto runWithCoupling = getParam("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("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("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 temperatures; for (int solIdx = 0; solIdx < numberOfElements; ++solIdx) { temperatures.push_back(x[solIdx][problem->returnTemperatureIdx()]); }; - if (getParam("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(problem, gridGeometry); gridVariables->init(x); - // intialize the vtk output module + // initialize the vtk output module using IOFields = GetPropType; VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); @@ -220,7 +226,7 @@ int main(int argc, char **argv) // output every vtkOutputInterval time step const int vtkOutputInterval = getParam("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("Precice.RunWithCoupling") == true) { + if (runWithCoupling) { solverDt = getParam("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("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("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("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("TimeLoop.MaxDt")); + std::flush(std::cout); + couplingParticipant.advance(dt); - std::cout << "Using dt: " << dt << std::endl; - - if (getParam("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("Precice.RunWithCoupling") == true) { + if (runWithCoupling) { couplingParticipant.finalize(); } // print dumux end message diff --git a/two-scale-heat-conduction/macro-dumux/params.input b/two-scale-heat-conduction/macro-dumux/params.input index 660557e1c..37508f98f 100644 --- a/two-scale-heat-conduction/macro-dumux/params.input +++ b/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] diff --git a/two-scale-heat-conduction/precice-config.xml b/two-scale-heat-conduction/precice-config.xml index 4970b5168..a264649e8 100644 --- a/two-scale-heat-conduction/precice-config.xml +++ b/two-scale-heat-conduction/precice-config.xml @@ -63,7 +63,7 @@ - +