box models: correct direction of flux vectors in base class as well as the individual models

also add a comment to the BoxLocalResidual to avoid confusion

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

dumux/boxmodels/1p2c/1p2clocalresidual.hh

@@ -164,10 +164,6 @@ public:
         fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
         flux[transEqIdx] +=
             c * (normalDisp * fluxVars.concentrationGrad(1));
-                
-        // we need to calculate the flux from i to j, not the other
-        // way round...
-        flux *= -1;
     }
 
     /*!

dumux/boxmodels/2p/2plocalresidual.hh

@@ -131,7 +131,6 @@ public:
         flux = 0;
         asImp_()->computeAdvectiveFlux(flux, vars);
         asImp_()->computeDiffusiveFlux(flux, vars);
-        flux *= -1;
     }
 
     /*!

dumux/boxmodels/2p2c/2p2clocalresidual.hh

@@ -203,10 +203,6 @@ public:
         flux = 0;
         asImp_()->computeAdvectiveFlux(flux, vars);
         asImp_()->computeDiffusiveFlux(flux, vars);
-
-        // the face normal points into the outward direction, so we
-        // have to multiply all fluxes with -1
-        flux *= -1;
     }
 
     /*!

dumux/boxmodels/common/boxlocalresidual.hh

@@ -414,12 +414,21 @@ protected:
             this->asImp_().computeFlux(flux, k);
             Valgrind::CheckDefined(flux);
 
-            // subtract fluxes from the local mass rates of the
-            // respective sub control volume adjacent to the face.
-            for (int eq = 0; eq < numEq; ++ eq) {
-                residual_[i][eq] -= flux[eq];
-                residual_[j][eq] += flux[eq];
-            }
+            // The balance equation for a finite volume is:
+            // 
+            // dStorage/dt = Flux + Source
+            //
+            // Re-arranging this, we get 
+            //
+            // dStorage/dt - Source - Flux = 0
+            //
+            // The residual already contains the "dStorage/dt -
+            // Source" term, so we have to subtract the flux
+            // term. Since the calculated flux goes from sub-control
+            // volume i to sub-control volume j, we need to add the
+            // flux to j and subtract it from i
+            residual_[i] += flux;
+            residual_[j] -= flux;
         }
     }
 
@@ -448,7 +457,7 @@ protected:
         // evaluate the volume terms (storage + source terms)
         for (int i=0; i < fvElemGeom_().numVertices; i++)
         {
-            PrimaryVariables massContrib(0), tmp(0);
+            PrimaryVariables dStorage_dt(0), tmp(0);
 
             // mass balance within the element. this is the
             // $\frac{m}{\partial t}$ term if using implicit
@@ -456,30 +465,25 @@ protected:
             //
             // TODO (?): we might need a more explicit way for
             // doing the time discretization...
-            this->asImp_().computeStorage(massContrib, i, false);
+            this->asImp_().computeStorage(dStorage_dt, i, false);
             this->asImp_().computeStorage(tmp, i, true);
 
-            massContrib -= tmp;
-            massContrib *=
+            dStorage_dt -= tmp;
+            dStorage_dt *=
                 fvElemGeom_().subContVol[i].volume
                 /
                 problem_().timeManager().timeStepSize();
-            
-            for (int j = 0; j < numEq; ++j)
-                residual_[i][j] += massContrib[j];
+            residual_[i] += dStorage_dt;
 
             // subtract the source term from the local rate
             PrimaryVariables source;
             this->asImp_().computeSource(source, i);
             source *= fvElemGeom_().subContVol[i].volume;
-
-            for (int j = 0; j < numEq; ++j) {
-                residual_[i][j] -= source[j];
-
-                // make sure that only defined quantities where used
-                // to calculate the residual.
-                Valgrind::CheckDefined(residual_[i][j]);
-            }
+            residual_[i] -= source;
+            
+            // make sure that only defined quantities where used
+            // to calculate the residual.
+            Valgrind::CheckDefined(residual_[i]);
         }
     }
 

dumux/boxmodels/richards/richardslocalresidual.hh

@@ -131,9 +131,6 @@ public:
             ((    mobilityUpwindAlpha)*up.density(wPhaseIdx)*up.mobility(wPhaseIdx) 
              +
              (1 - mobilityUpwindAlpha)*dn.density(wPhaseIdx)*dn.mobility(wPhaseIdx));
-
-        // we need the flux from i to j instead of the other way
-        flux *= -1;
     }
 
     /*!