From 4ba12f4fe2cc4c7d1b82249f5f60e80bccf4b035 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 13 Aug 2010 08:48:54 +0000
Subject: [PATCH] box models: implement time step ramp up
the biggest open issue is how to combine this with partial jacobian
reassembly. (it works as it is implemented now but the percentage of
elements which needs to be reassembled is quite high between the
newton iterations.) one idea is to just reassemble a fixed fraction of
vertices, e.g. 20%. Here the problem is how to determine this set of
vertices efficiently...
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4066 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/2p/2pproperties.hh | 2 +
dumux/boxmodels/2p2c/2p2cproperties.hh | 2 +
dumux/boxmodels/common/boxproperties.hh | 6 +++
dumux/nonlinear/newtoncontroller.hh | 57 ++++++++++++++++++++++---
4 files changed, 61 insertions(+), 6 deletions(-)
diff --git a/dumux/boxmodels/2p/2pproperties.hh b/dumux/boxmodels/2p/2pproperties.hh
index f32de65019..86a938c7e9 100644
--- a/dumux/boxmodels/2p/2pproperties.hh
+++ b/dumux/boxmodels/2p/2pproperties.hh
@@ -225,6 +225,8 @@ SET_TYPE_PROP(BoxTwoP, FluidState, TwoPFluidState<TypeTag>);
SET_BOOL_PROP(BoxTwoP, EnableJacobianRecycling, true);
// enable partial reassembling by default
SET_BOOL_PROP(BoxTwoP, EnablePartialReassemble, true);
+// enable time-step ramp up by default
+SET_BOOL_PROP(BoxTwoP, EnableTimeStepRampUp, true);
// \}
}
diff --git a/dumux/boxmodels/2p2c/2p2cproperties.hh b/dumux/boxmodels/2p2c/2p2cproperties.hh
index 37ddd2b724..d42c6ed7f3 100644
--- a/dumux/boxmodels/2p2c/2p2cproperties.hh
+++ b/dumux/boxmodels/2p2c/2p2cproperties.hh
@@ -274,6 +274,8 @@ public:
SET_BOOL_PROP(BoxTwoPTwoC, EnableJacobianRecycling, true);
// enable partial reassembling by default
SET_BOOL_PROP(BoxTwoPTwoC, EnablePartialReassemble, true);
+// enable time-step ramp up by default
+SET_BOOL_PROP(BoxTwoPTwoC, EnableTimeStepRampUp, true);
// \}
}
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index 71ceef6371..5517f9cc5c 100644
--- a/dumux/boxmodels/common/boxproperties.hh
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -97,6 +97,10 @@ NEW_PROP_TAG(EnableJacobianRecycling);
//! tolerance
NEW_PROP_TAG(EnablePartialReassemble);
+//! Specify whether the time step should be increased in between
+//! newton iterations to achive larger time step sizes
+NEW_PROP_TAG(EnableTimeStepRampUp);
+
// mappers from local to global indices
NEW_PROP_TAG(VertexMapper);
NEW_PROP_TAG(ElementMapper);
@@ -372,6 +376,8 @@ SET_PROP(BoxModel, LocalOperator)
SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);
// disable partial reassembling by default
SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);
+// disable time-step ramp up by default
+SET_BOOL_PROP(BoxModel, EnableTimeStepRampUp, false);
// \}
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 550965144a..7c915f049b 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -183,6 +183,7 @@ class NewtonController
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
@@ -190,18 +191,21 @@ class NewtonController
typedef NewtonConvergenceWriter<TypeTag, GET_PROP_VALUE(TypeTag, PTAG(NewtonWriteConvergence))> ConvergenceWriter;
+ enum { enableTimeStepRampUp = GET_PROP_VALUE(TypeTag, PTAG(EnableTimeStepRampUp)) };
+
public:
NewtonController()
: endIterMsgStream_(std::ostringstream::out),
convergenceWriter_(asImp_())
{
numSteps_ = 0;
+ rampUpSteps_ = 4;
// maximum acceptable difference of any primary variable
// between two iterations for convergence
setRelTolerance(1e-7);
- setTargetSteps(8);
- setMaxSteps(12);
+ setTargetSteps(9);
+ setMaxSteps(15);
};
/*!
@@ -230,7 +234,7 @@ public:
*/
bool newtonProceed(const SolutionVector &u)
{
- if (numSteps_ < 2)
+ if (numSteps_ < rampUpSteps_)
return true; // we always do at least two iterations
else if (numSteps_ >= maxSteps_) {
// we have exceeded the allowed number of steps. if the
@@ -263,6 +267,18 @@ public:
method_ = &method;
numSteps_ = 0;
+ if (enableTimeStepRampUp) {
+ destTimeStepSize_ = timeManager_().timeStepSize();
+
+ // reduce initial time step size for ramp-up
+ timeManager_().setTimeStepSize(destTimeStepSize_ / rampUpSteps_);
+
+ // reduce the tolerance for partial reassembly during the
+ // ramp up
+ finalTolerance_ = tolerance_;
+ tolerance_ = tolerance_*1e8;
+ }
+
convergenceWriter_.beginTimestep();
}
@@ -397,6 +413,17 @@ public:
void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
{
++numSteps_;
+
+ if (enableTimeStepRampUp) {
+ if (numSteps_ < rampUpSteps_) {
+ // increase time step size
+ Scalar dt = destTimeStepSize_ * (numSteps_ + 1) / rampUpSteps_;
+ timeManager_().setTimeStepSize(dt);
+ }
+ else // if we're not in the ramp-up reduce the target
+ // tolerance
+ tolerance_ = finalTolerance_;
+ }
if (verbose())
std::cout << "\rNewton iteration " << numSteps_ << " done: "
@@ -409,6 +436,8 @@ public:
*/
void newtonEnd()
{
+ if (enableTimeStepRampUp)
+ timeManager_().setTimeStepSize(destTimeStepSize_);
convergenceWriter_.endTimestep();
}
@@ -439,20 +468,24 @@ public:
*/
Scalar suggestTimeStepSize(Scalar oldTimeStep) const
{
+ Scalar n = numSteps_;
+ if (enableTimeStepRampUp) {
+ n -= rampUpSteps_/2;
+ }
// be agressive reducing the timestep size but
// conservative when increasing it. the rationale is
// that we want to avoid failing in the next newton
// iteration which would require another linearization
// of the problem.
- if (numSteps_ > targetSteps_) {
- Scalar percent = ((Scalar) numSteps_ - targetSteps_)/targetSteps_;
+ if (n > targetSteps_) {
+ Scalar percent = (n - targetSteps_)/targetSteps_;
return oldTimeStep/(1.0 + percent);
}
else {
/*Scalar percent = (Scalar(1))/targetSteps_;
return oldTimeStep*(1 + percent);
*/
- Scalar percent = ((Scalar) targetSteps_ - numSteps_)/targetSteps_;
+ Scalar percent = (targetSteps_ - n)/targetSteps_;
return oldTimeStep*(1.0 + percent/1.2);
}
}
@@ -499,6 +532,12 @@ protected:
const Problem &problem_() const
{ return method_->problem(); }
+ /*!
+ * \brief Returns a reference to the time manager.
+ */
+ TimeManager &timeManager_()
+ { return problem_().timeManager(); }
+
/*!
* \brief Returns a reference to the problem.
*/
@@ -579,10 +618,16 @@ protected:
ConvergenceWriter convergenceWriter_;
Scalar tolerance_;
+ Scalar finalTolerance_;
Scalar error_;
Scalar lastError_;
+ // number of iterations for the time-step ramp-up
+ Scalar rampUpSteps_;
+ // time step size after the ramp-up
+ Scalar destTimeStepSize_;
+
// optimal number of iterations we want to achive
int targetSteps_;
// maximum number of iterations we do before giving up
--
GitLab