implicit 3p3c: add velocity calculation and output. Reviewed by Christoph.

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

dumux/implicit/3p3c/3p3cmodel.hh

@@ -27,6 +27,7 @@
 #ifndef DUMUX_3P3C_MODEL_HH
 #define DUMUX_3P3C_MODEL_HH
 
+#include <dumux/implicit/common/implicitvelocityoutput.hh>
 #include "3p3cproperties.hh"
 
 namespace Dumux
@@ -97,6 +98,7 @@ class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
 
@@ -284,6 +286,7 @@ public:
                             MultiWriter &writer)
     {
         typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+        typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
 
         // get the number of degrees of freedom
         unsigned numDofs = this->numDofs();
@@ -307,13 +310,24 @@ public:
         ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
         ScalarField *poro = writer.allocateManagedBuffer(numDofs);
         ScalarField *perm = writer.allocateManagedBuffer(numDofs);
+        VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
+        VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
+        VectorField *velocityG = writer.template allocateManagedBuffer<double, dim>(numDofs);
+        ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
 
-        unsigned numElements = this->gridView_().size(0);
-        ScalarField *rank =
-            writer.allocateManagedBuffer (numElements);
+        if (velocityOutput.enableOutput()) // check if velocity output is demanded
+        {
+            // initialize velocity fields
+            for (unsigned int i = 0; i < numDofs; ++i)
+            {
+                (*velocityN)[i] = Scalar(0);
+                (*velocityW)[i] = Scalar(0);
+                (*velocityG)[i] = Scalar(0);
+            }
+        }
 
-        FVElementGeometry fvGeometry;
-        VolumeVariables volVars;
+        unsigned numElements = this->gridView_().size(0);
+        ScalarField *rank = writer.allocateManagedBuffer (numElements);
 
         ElementIterator elemIt = this->gridView_().template begin<0>();
         ElementIterator elemEndIt = this->gridView_().template end<0>();
@@ -321,43 +335,54 @@ public:
         {
             int idx = this->problem_().elementMapper().map(*elemIt);
             (*rank)[idx] = this->gridView_().comm().rank();
+
+            FVElementGeometry fvGeometry;
             fvGeometry.update(this->gridView_(), *elemIt);
 
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
 
-                volVars.update(sol[globalIdx],
-                               this->problem_(),
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
                                *elemIt,
                                fvGeometry,
-                               scvIdx,
-                               false);
+                               false /* oldSol? */);
+
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
 
                 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                    (*saturation[phaseIdx])[globalIdx] = volVars.fluidState().saturation(phaseIdx);
-                    (*pressure[phaseIdx])[globalIdx] = volVars.fluidState().pressure(phaseIdx);
-                    (*density[phaseIdx])[globalIdx] = volVars.fluidState().density(phaseIdx);
+                    (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().saturation(phaseIdx);
+                    (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().pressure(phaseIdx);
+                    (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().density(phaseIdx);
                 }
 
                 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                     for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                         (*moleFraction[phaseIdx][compIdx])[globalIdx] =
-                            volVars.fluidState().moleFraction(phaseIdx,
+                            elemVolVars[scvIdx].fluidState().moleFraction(phaseIdx,
                                                               compIdx);
 
                         Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
                     }
                 }
 
-                (*poro)[globalIdx] = volVars.porosity();
-                (*perm)[globalIdx] = volVars.permeability();
-                (*temperature)[globalIdx] = volVars.temperature();
+                (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
+                (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
+                (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
                 (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
             }
 
+            // velocity output
+            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *elemIt, wPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *elemIt, nPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *elemIt, gPhaseIdx);
+
         }
 
+        velocityOutput.completeVelocityCalculation(*velocityW);
+        velocityOutput.completeVelocityCalculation(*velocityN);
+        velocityOutput.completeVelocityCalculation(*velocityG);
+
         writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
         writer.attachDofData(*saturation[nPhaseIdx], "Sn", isBox);
         writer.attachDofData(*saturation[gPhaseIdx], "Sg", isBox);
@@ -384,7 +409,14 @@ public:
         writer.attachDofData(*perm, "permeability", isBox);
         writer.attachDofData(*temperature, "temperature", isBox);
         writer.attachDofData(*phasePresence, "phase presence", isBox);
-        
+
+        if (velocityOutput.enableOutput()) // check if velocity output is demanded
+        {
+            writer.attachDofData(*velocityW,  "velocityW", isBox, dim);
+            writer.attachDofData(*velocityN,  "velocityN", isBox, dim);
+            writer.attachDofData(*velocityG,  "velocityG", isBox, dim);
+        }
+
         writer.attachCellData(*rank, "process rank");
     }
 

dumux/implicit/3p3c/3p3cproperties.hh

@@ -63,6 +63,7 @@ NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility
 NEW_PROP_TAG(UseConstraintSolver); //!< Determines whether a constraint solver should be used explicitly
 NEW_PROP_TAG(BaseFluxVariables); //! The base flux variables
 NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
 }
 }
 

dumux/implicit/3p3c/3p3cpropertydefaults.hh

@@ -126,6 +126,9 @@ SET_TYPE_PROP(ThreePThreeC, Indices, ThreePThreeCIndices<TypeTag, /*PVOffset=*/0
 //! Use ImplicitSpatialParams by default.
 SET_TYPE_PROP(ThreePThreeC, SpatialParams, ImplicitSpatialParams<TypeTag>);
 
+// disable velocity output by default
+SET_BOOL_PROP(ThreePThreeC, VtkAddVelocity, false);
+
 // enable gravity by default
 SET_BOOL_PROP(ThreePThreeC, ProblemEnableGravity, true);