From c67346644557106a0ceeca5a568f88b7d7784328 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Tue, 14 Sep 2021 00:27:28 +0200
Subject: [PATCH] add adaptivetimestepper
---
src/foam/foam.cc | 47 ++++++++++++++++++++++++++++++++++++++++++++---
1 file changed, 44 insertions(+), 3 deletions(-)
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index b606d65b..e9053748 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -445,6 +445,42 @@ int main(int argc, char *argv[]) {
const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+ auto const mustRefine = [&](Updaters &coarseUpdater,
+ Updaters &fineUpdater) {
+
+ //return false;
+ //std::cout << "Step 1" << std::endl;
+
+ std::vector<ScalarVector> coarseAlpha;
+ coarseAlpha.resize(bodyCount);
+ coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+ //print(coarseAlpha, "coarseAlpha:");
+
+ std::vector<ScalarVector> fineAlpha;
+ fineAlpha.resize(bodyCount);
+ fineUpdater.state_->extractAlpha(fineAlpha);
+
+ //print(fineAlpha, "fineAlpha:");
+
+ //std::cout << "Step 3" << std::endl;
+
+ ScalarVector::field_type energyNorm = 0;
+ for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+ //std::cout << "for " << i << std::endl;
+
+ //std::cout << not stateEnergyNorms[i] << std::endl;
+
+ if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
+ continue;
+
+ energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
+ }
+ //std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
+ //std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
+ return energyNorm > refinementTolerance;
+ };
+
std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
@@ -477,12 +513,17 @@ int main(int argc, char *argv[]) {
stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
externalForces, stateEnergyNorms);
- UniformTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+ const auto minTau = parset.get<double>("timeSteps.minTau");
+ AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
timeStepper(stepBase, contactNetwork, current,
- programState.relativeTime, programState.relativeTau);
+ programState.relativeTime, programState.relativeTau, minTau,
+ mustRefine);
- size_t timeSteps = std::round(parset.get<double>("timeSteps.timeSteps"));
+ /*UniformTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+ timeStepper(stepBase, contactNetwork, current,
+ programState.relativeTime, programState.relativeTau);*/
+ size_t timeSteps = std::round(parset.get<double>("timeSteps.timeSteps"));
while (!timeStepper.reachedEnd()) {
programState.timeStep++;
--
GitLab