From f6d1b7758610538dcc572480af74b08e732534c9 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Sun, 31 Aug 2025 21:14:12 +0200 Subject: [PATCH 1/7] [macro-dumux] Handle subcycling in macro-dumux. Fix checkpointing. --- .../macro-dumux/appl/main.cc | 19 ++++++++++++++++--- .../macro-dumux/params.input | 4 ++-- two-scale-heat-conduction/precice-config.xml | 2 +- 3 files changed, 19 insertions(+), 6 deletions(-) diff --git a/two-scale-heat-conduction/macro-dumux/appl/main.cc b/two-scale-heat-conduction/macro-dumux/appl/main.cc index 3491731a5..19675dbf2 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/two-scale-heat-conduction/macro-dumux/appl/main.cc @@ -171,6 +171,10 @@ 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) { @@ -267,7 +271,9 @@ int main(int argc, char **argv) // write checkpoint if (couplingParticipant.requiresToWriteCheckpoint()) { - xOld = x; + xCheckpoint = x; + timeCheckpoint = timeLoop->time(); + timeStepCheckpoint = timeLoop->timeStepIndex(); } // read porosity and conductivity data from other solver @@ -322,11 +328,18 @@ int main(int argc, char **argv) if (getParam("Precice.RunWithCoupling") == true) { if (couplingParticipant.requiresToReadCheckpoint()) { // make the new solution the old solution - x = xOld; + gridVariables->advanceTimeStep(); + timeLoop->advanceTimeStep(); + timeLoop->reportTimeStep(); + x = xCheckpoint; + xOld = x; + timeLoop->setTime(timeCheckpoint, timeStepCheckpoint); + timeLoop->setTimeStepSize(dt); gridVariables->update(x); gridVariables->advanceTimeStep(); - } else // coupling successful + } else // coupling successful OR not at time window!! { + xOld = x; n += 1; if (n == vtkOutputInterval) { problem->updateVtkOutput(x); 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 @@ - + From 694e8aee3bb9c2f18e5904ae0bdddbc2269947bc Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Sun, 31 Aug 2025 21:16:09 +0200 Subject: [PATCH 2/7] [macro-dumux] Warn about different dt after solve --- .../macro-dumux/appl/main.cc | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/two-scale-heat-conduction/macro-dumux/appl/main.cc b/two-scale-heat-conduction/macro-dumux/appl/main.cc index 19675dbf2..0ebfafb22 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/two-scale-heat-conduction/macro-dumux/appl/main.cc @@ -292,8 +292,10 @@ int main(int argc, char **argv) } std::cout << "Solver starts" << std::endl; + std::cout << "Using dt: " << dt << std::endl; // linearize & solve nonLinearSolver.solve(x, *timeLoop); + dt = timeLoop->timeStepSize(); for (int solIdx = 0; solIdx < numberOfElements; ++solIdx) temperatures[solIdx] = x[solIdx][problem->returnTemperatureIdx()]; @@ -309,22 +311,23 @@ int main(int argc, char **argv) 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); - 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; } + } + + if (getParam("Precice.RunWithCoupling") == true) { + preciceDt = couplingParticipant.getMaxTimeStepSize(); + solverDt = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), timeLoop->maxTimeStepSize()); + dt = std::min(preciceDt, solverDt); + } else dt = std::min( nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), getParam("TimeLoop.MaxDt")); - std::cout << "Using dt: " << dt << std::endl; - if (getParam("Precice.RunWithCoupling") == true) { if (couplingParticipant.requiresToReadCheckpoint()) { // make the new solution the old solution From d241f90ffc0d516752820ab06a386037b6b24c17 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Sun, 31 Aug 2025 21:20:23 +0200 Subject: [PATCH 3/7] [macro-dumux] Negotiate dt before solve --- .../macro-dumux/appl/main.cc | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/two-scale-heat-conduction/macro-dumux/appl/main.cc b/two-scale-heat-conduction/macro-dumux/appl/main.cc index 0ebfafb22..369a08df6 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/two-scale-heat-conduction/macro-dumux/appl/main.cc @@ -276,6 +276,11 @@ int main(int argc, char **argv) timeStepCheckpoint = timeLoop->timeStepIndex(); } + preciceDt = couplingParticipant.getMaxTimeStepSize(); + solverDt = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), timeLoop->maxTimeStepSize()); + dt = std::min(preciceDt, solverDt); + timeLoop->setTimeStepSize(dt); + // read porosity and conductivity data from other solver couplingParticipant.readQuantityFromOtherSolver(meshName, readDatak00, dt); @@ -289,7 +294,13 @@ int main(int argc, char **argv) readDataPorosity, dt); // store coupling data in problem problem->spatialParams().updateCouplingData(); + } else { + dt = std::min( + nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), + getParam("TimeLoop.MaxDt")); + timeLoop->setTimeStepSize(dt); } + std::cout << "Solver starts" << std::endl; std::cout << "Using dt: " << dt << std::endl; @@ -318,16 +329,6 @@ int main(int argc, char **argv) } } - if (getParam("Precice.RunWithCoupling") == true) { - preciceDt = couplingParticipant.getMaxTimeStepSize(); - solverDt = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), timeLoop->maxTimeStepSize()); - dt = std::min(preciceDt, solverDt); - - } else - dt = std::min( - nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), - getParam("TimeLoop.MaxDt")); - if (getParam("Precice.RunWithCoupling") == true) { if (couplingParticipant.requiresToReadCheckpoint()) { // make the new solution the old solution From 506e6703c62a1b7055444fdfa9ed33256569e2d7 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Sun, 31 Aug 2025 21:31:29 +0200 Subject: [PATCH 4/7] [macro-dumux] Advance DuMux by default after solve --- .../macro-dumux/appl/main.cc | 30 +++++++------------ 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/two-scale-heat-conduction/macro-dumux/appl/main.cc b/two-scale-heat-conduction/macro-dumux/appl/main.cc index 369a08df6..eb6d2de26 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/two-scale-heat-conduction/macro-dumux/appl/main.cc @@ -306,7 +306,12 @@ int main(int argc, char **argv) std::cout << "Using dt: " << dt << std::endl; // linearize & solve nonLinearSolver.solve(x, *timeLoop); + // save actual time-step size dt = timeLoop->timeStepSize(); + // DuMux advance + report + gridVariables->advanceTimeStep(); + timeLoop->advanceTimeStep(); + timeLoop->reportTimeStep(); for (int solIdx = 0; solIdx < numberOfElements; ++solIdx) temperatures[solIdx] = x[solIdx][problem->returnTemperatureIdx()]; @@ -320,21 +325,19 @@ int main(int argc, char **argv) // advance precice if (getParam("Precice.RunWithCoupling") == true) { - couplingParticipant.advance(timeLoop->timeStepSize()); 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; } + couplingParticipant.advance(dt); } if (getParam("Precice.RunWithCoupling") == true) { if (couplingParticipant.requiresToReadCheckpoint()) { - // make the new solution the old solution - gridVariables->advanceTimeStep(); - timeLoop->advanceTimeStep(); - timeLoop->reportTimeStep(); x = xCheckpoint; xOld = x; timeLoop->setTime(timeCheckpoint, timeStepCheckpoint); @@ -350,20 +353,9 @@ int main(int argc, char **argv) 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) { From 8357262d2104f10cb2e2948abee5be87c4ef3ed3 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Sun, 31 Aug 2025 22:37:38 +0200 Subject: [PATCH 5/7] [macro-dumux] Cleanup time loop --- .../macro-dumux/appl/main.cc | 82 +++++++++---------- 1 file changed, 38 insertions(+), 44 deletions(-) diff --git a/two-scale-heat-conduction/macro-dumux/appl/main.cc b/two-scale-heat-conduction/macro-dumux/appl/main.cc index eb6d2de26..fe8afbb5a 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/two-scale-heat-conduction/macro-dumux/appl/main.cc @@ -100,7 +100,8 @@ 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 @@ -138,7 +139,7 @@ int main(int argc, char **argv) // 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 +155,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); @@ -181,7 +182,7 @@ int main(int argc, char **argv) temperatures.push_back(x[solIdx][problem->returnTemperatureIdx()]); }; - if (getParam("Precice.RunWithCoupling") == true) { + if (runWithCoupling) { couplingParticipant.writeQuantityVector(meshName, writeDataConcentration, temperatures); if (couplingParticipant @@ -233,7 +234,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 { @@ -265,7 +266,7 @@ 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; @@ -277,11 +278,13 @@ int main(int argc, char **argv) } preciceDt = couplingParticipant.getMaxTimeStepSize(); - solverDt = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), timeLoop->maxTimeStepSize()); + solverDt = std::min(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), + timeLoop->maxTimeStepSize()); dt = std::min(preciceDt, solverDt); - timeLoop->setTimeStepSize(dt); // 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, @@ -292,80 +295,71 @@ 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()), - getParam("TimeLoop.MaxDt")); - timeLoop->setTimeStepSize(dt); + timeLoop->maxTimeStepSize()); } + // set new dt as suggested by the newton solver or by preCICE + timeLoop->setTimeStepSize(dt); - std::cout << "Solver starts" << std::endl; + std::cout << "Solver starts with Target dt: " << dt << std::endl; - std::cout << "Using dt: " << dt << std::endl; // linearize & solve nonLinearSolver.solve(x, *timeLoop); // save actual time-step size dt = timeLoop->timeStepSize(); + // DuMux advance + report gridVariables->advanceTimeStep(); timeLoop->advanceTimeStep(); timeLoop->reportTimeStep(); + xOld = x; + + // Vtk output + // TODO: output interval doesn't work easily with subcycling + n += 1; + if (n == vtkOutputInterval) { + problem->updateVtkOutput(x); + vtkWriter.write(timeLoop->time()); + n = 0; + } - for (int solIdx = 0; solIdx < numberOfElements; ++solIdx) - temperatures[solIdx] = x[solIdx][problem->returnTemperatureIdx()]; + if (runWithCoupling) { + // write coupling data to preCICE + for (int solIdx = 0; solIdx < numberOfElements; ++solIdx) + temperatures[solIdx] = x[solIdx][problem->returnTemperatureIdx()]; - if (getParam("Precice.RunWithCoupling") == true) { couplingParticipant.writeQuantityVector(meshName, writeDataConcentration, temperatures); couplingParticipant.writeQuantityToOtherSolver(meshName, writeDataConcentration); - } - - // advance precice - if (getParam("Precice.RunWithCoupling") == true) { + // 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 + << " and DuMuX dt = " << solverDt << " resulted in dt = " << dt << std::endl; } + std::flush(std::cout); couplingParticipant.advance(dt); - } - if (getParam("Precice.RunWithCoupling") == true) { + // reset to checkpoint if not converged if (couplingParticipant.requiresToReadCheckpoint()) { 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 OR not at time window!! - { - xOld = x; - n += 1; - if (n == vtkOutputInterval) { - problem->updateVtkOutput(x); - vtkWriter.write(timeLoop->time()); - n = 0; - } - } - } else { - xOld = x; - // 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; @@ -376,7 +370,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 From 80297a125d9008f0680ac681722fced9f1e07cbb Mon Sep 17 00:00:00 2001 From: Ishaan Desai Date: Sun, 7 Sep 2025 17:03:09 +0200 Subject: [PATCH 6/7] Add changelog entry --- changelog-entries/666.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 changelog-entries/666.md diff --git a/changelog-entries/666.md b/changelog-entries/666.md new file mode 100644 index 000000000..6aa0caa91 --- /dev/null +++ b/changelog-entries/666.md @@ -0,0 +1 @@ +- Fix subcycling and checkpointing in macro-dumux of tw0-scale-heat-conduction [#666](https://github.com/precice/tutorials/pull/666) From 676e941647fa6121a320e73ca200674e103befb5 Mon Sep 17 00:00:00 2001 From: Ishaan Desai Date: Sun, 7 Sep 2025 17:57:39 +0200 Subject: [PATCH 7/7] Formatting and making code comment more informative --- changelog-entries/666.md | 2 +- .../macro-dumux/appl/main.cc | 19 +++++++++++-------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/changelog-entries/666.md b/changelog-entries/666.md index 6aa0caa91..069cca73c 100644 --- a/changelog-entries/666.md +++ b/changelog-entries/666.md @@ -1 +1 @@ -- Fix subcycling and checkpointing in macro-dumux of tw0-scale-heat-conduction [#666](https://github.com/precice/tutorials/pull/666) +- 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 10ab987de..d36090d31 100644 --- a/two-scale-heat-conduction/macro-dumux/appl/main.cc +++ b/two-scale-heat-conduction/macro-dumux/appl/main.cc @@ -101,6 +101,7 @@ int main(int argc, char **argv) auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); const auto runWithCoupling = getParam("Precice.RunWithCoupling"); + if (runWithCoupling) { couplingParticipant.announceSolver("macro-heat", preciceConfigFilename, mpiHelper.rank(), mpiHelper.size()); @@ -133,10 +134,10 @@ 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 (runWithCoupling) { @@ -208,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()); @@ -225,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 @@ -302,14 +303,15 @@ int main(int argc, char **argv) nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()), timeLoop->maxTimeStepSize()); } - // set new dt as suggested by the newton solver or by preCICE + // 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; + std::cout << "Solver starts with target dt: " << dt << std::endl; // linearize & solve nonLinearSolver.solve(x, *timeLoop); - // save actual time-step size + + // save dt value that was actually used by the non-linear solver dt = timeLoop->timeStepSize(); // DuMux advance + report @@ -319,7 +321,7 @@ int main(int argc, char **argv) xOld = x; // Vtk output - // TODO: output interval doesn't work easily with subcycling + // TODO: output interval does not work seamlessly when subcycling n += 1; if (n == vtkOutputInterval) { problem->updateVtkOutput(x); @@ -343,7 +345,7 @@ int main(int argc, char **argv) couplingParticipant.writeQuantityToOtherSolver(meshName, writeDataConcentration); - // advance precice + // advance preCICE if ((!fabs(preciceDt - dt)) < 1e-14) { std::cout << "dt from preCICE is different than dt from DuMuX." << " preCICE dt = " << preciceDt @@ -359,6 +361,7 @@ int main(int argc, char **argv) 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);