partial reassemble: only update the solution if the DOF gets reassembled

also make the MpNc model respect absolute convergence criteria. IMHO
an unweighted absolute criterion is pretty useless in most contexts,
though. that's because the criterion is about the amount of
conservation quantities that get lost per second but a losing or
gaining a joule of energy is much less of a problem than losing or
gaining a kilogram of mass...

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7442 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/boxmodels/MpNc/MpNcnewtoncontroller.hh

@@ -201,13 +201,13 @@ public:
                       const SolutionVector &uLastIter,
                       const SolutionVector &deltaU)
     {
-        this->writeConvergence_(uLastIter, deltaU);
-        this->newtonUpdateRelError(uLastIter, deltaU);
+        if (this->enableRelativeCriterion_ || this->enablePartialReassemble_)
+            this->newtonUpdateRelError(uLastIter, deltaU);
 
         // compute the vertex and element colors for partial
         // reassembly
         if (enablePartialReassemble) {
-            Scalar minReasmTol = 0.5*this->tolerance_;
+            Scalar minReasmTol = 0.01*this->tolerance_;
             Scalar reassembleTol = Dumux::geometricMean(this->error_, minReasmTol);
             reassembleTol = std::max(reassembleTol, minReasmTol);
 
@@ -215,28 +215,32 @@ public:
             this->model_().jacobianAssembler().computeColors(reassembleTol);
         }
 
+        this->writeConvergence_(uLastIter, deltaU);
+
         if (GET_PROP_VALUE(TypeTag, NewtonUseLineSearch)) {
             lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
         }
         else {
-            Scalar lambda = 1.0;
-            /*
-            if (this->error_ > 10) {
-                // Do a "poor man's line search"
-                lambda /= std::sqrt((this->numSteps_ + 1)*this->error_/10);
-                lambda = std::max(0.2, lambda);
+            for (int i = 0; i < uLastIter.size(); ++i) {
+                if (this->model_().jacobianAssembler().vertexColor(i) != JacobianAssembler::Red)
+                    continue;
+                
+                uCurrentIter[i] = uLastIter[i];
+                uCurrentIter[i] -= deltaU[i];
             }
-            */
-
-            uCurrentIter = deltaU;
-            uCurrentIter *= - lambda;
-            uCurrentIter += uLastIter;
-
-            if (this->numSteps_ < 4 && enableChop_) {
+            
+            if (this->numSteps_ < 2 && enableChop_) {
                 // put crash barriers along the update path at the
                 // beginning...
                 NewtonChop::chop(uCurrentIter, uLastIter);
             }
+            
+            if (this->enableAbsoluteCriterion_)
+            {
+                SolutionVector tmp(uLastIter);
+                this->absoluteError_ = this->method().model().globalResidual(tmp, uCurrentIter);
+                this->absoluteError_ /= this->initialAbsoluteError_;
+            }
         }
     }
 

dumux/nonlinear/newtoncontroller.hh

@@ -115,6 +115,9 @@ NEW_PROP_TAG(NewtonTargetSteps);
 //! Number of maximum iterations for the Newton method.
 NEW_PROP_TAG(NewtonMaxSteps);
 
+//! The assembler for the Jacobian matrix
+NEW_PROP_TAG(JacobianAssembler);
+
 // set default values
 SET_TYPE_PROP(NewtonMethod, NewtonController, Dumux::NewtonController<TypeTag>);
 SET_BOOL_PROP(NewtonMethod, NewtonWriteConvergence, false);
@@ -149,6 +152,7 @@ class NewtonController
     typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
     typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
     typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
+    typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
     typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
     typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
 
@@ -444,30 +448,33 @@ public:
                       const SolutionVector &uLastIter,
                       const SolutionVector &deltaU)
     {
-        writeConvergence_(uLastIter, deltaU);
-
         if (enableRelativeCriterion_ || enablePartialReassemble_)
             newtonUpdateRelError(uLastIter, deltaU);
 
         // compute the vertex and element colors for partial reassembly
         if (enablePartialReassemble_) {
-            Scalar minReasmTol, tmp, reassembleTol;
-            minReasmTol = 0.1*tolerance_;
-            tmp = Dumux::geometricMean(error_, minReasmTol);
-            reassembleTol = Dumux::geometricMean(error_, tmp);
+            Scalar minReasmTol = 0.01*tolerance_;
+            Scalar reassembleTol = Dumux::geometricMean(error_, minReasmTol);
             reassembleTol = std::max(reassembleTol, minReasmTol);
             this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
             this->model_().jacobianAssembler().computeColors(reassembleTol);
         }
 
+        writeConvergence_(uLastIter, deltaU);
+
         if (useLineSearch_)
         {
             lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
         }
         else {
-            uCurrentIter = uLastIter;
-            uCurrentIter -= deltaU;
-
+            for (int i = 0; i < uLastIter.size(); ++i) {
+                if (this->model_().jacobianAssembler().vertexColor(i) != JacobianAssembler::Red)
+                    continue;
+                
+                uCurrentIter[i] = uLastIter[i];
+                uCurrentIter[i] -= deltaU[i];
+            }
+            
             if (enableAbsoluteCriterion_)
             {
                 SolutionVector tmp(uLastIter);
@@ -712,7 +719,6 @@ protected:
     // the linear solver
     LinearSolver linearSolver_;
 
-private:
     bool enablePartialReassemble_;
     bool enableJacobianRecycling_;
     bool useLineSearch_;