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

@@ -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);
 
 // \}
 }

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);
 
 // \}
 }

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);
 
 // \}
 

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