diff --git a/dumux/boxmodels/2p/2pproperties.hh b/dumux/boxmodels/2p/2pproperties.hh index f32de6501913289c7db78182e278bcb1bd3abb3a..86a938c7e9b9df175b4efa96cdef021fb2f74a11 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 37ddd2b724ab593c948e200695931216e9b9a3ed..d42c6ed7f3adbde63b8e6f053ab78d71658fbb90 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 71ceef6371ec3c44da814af4a7e7951355d797fc..5517f9cc5ceaa20f3ec6bdb0ea1ad473a630c7c9 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 550965144ad05238090d5cf60cea3f29a2570434..7c915f049b4047ddbe889dceb68087764066e780 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