updated completeFluidState for 2p(ni) models

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

dumux/boxmodels/2p/2pmodel.hh

@@ -76,6 +76,7 @@ class TwoPModel : public BoxModel<TypeTag>
     typedef TwoPModel<TypeTag> ThisType;
     typedef BoxModel<TypeTag> ParentType;
 
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Implementation;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
@@ -125,11 +126,20 @@ public:
         return 1;
     }
 
-    template <class PrimaryVariables, class MaterialParams, class FluidState>
+    template <class PrimaryVariables, class Problem, class Element, class ElementGeometry, class MaterialParams, class FluidState>
     static void completeFluidState(const PrimaryVariables& primaryVariables,
+                                   const Problem& problem,
+                                   const Element& element,
+                                   const ElementGeometry& elementGeometry,
+                                   int scvIdx,
                                    const MaterialParams& materialParams,
                                    FluidState& fluidState)
     {
+        Scalar t = Implementation::temperature_(primaryVariables, problem, element,
+                                                elementGeometry, scvIdx);
+        fluidState.setTemperature(t);
+
+        // material law parameters
         typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
 
         if (int(formulation) == pwSn) {
@@ -156,15 +166,20 @@ public:
         typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState);
-        fluidState.setViscosity(wPhaseIdx,
-                FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
-        fluidState.setViscosity(nPhaseIdx,
-                FluidSystem::viscosity(fluidState, paramCache, nPhaseIdx));
-
-        fluidState.setDensity(wPhaseIdx,
-                FluidSystem::density(fluidState, paramCache, wPhaseIdx));
-        fluidState.setDensity(nPhaseIdx,
-                FluidSystem::density(fluidState, paramCache, nPhaseIdx));
+
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            // compute and set the viscosity
+            Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
+            fluidState.setViscosity(phaseIdx, mu);
+
+            // compute and set the density
+            Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+            fluidState.setDensity(phaseIdx, rho);
+
+            // compute and set the enthalpy
+            Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+            fluidState.setEnthalpy(phaseIdx, h);
+        }
     }
 
     /*!
@@ -399,6 +414,24 @@ public:
         }
         writer.attachCellData(*rank, "process rank");
     }
+private:
+    template<class PrimaryVariables, class Problem, class Element, class FVElementGeometry>
+    static Scalar temperature_(const PrimaryVariables &priVars,
+                            const Problem& problem,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx)
+    {
+        return problem.boxTemperature(element, elemGeom, scvIdx);
+    }
+
+    template<class FluidState, class ParameterCache>
+    static Scalar enthalpy_(const FluidState& fluidState,
+                            const ParameterCache& paramCache,
+                            int phaseIdx)
+    {
+        return 0;
+    }
 };
 }
 

dumux/boxmodels/2p/2pvolumevariables.hh

@@ -102,17 +102,10 @@ public:
                            scvIdx,
                            isOldSol);
 
-        asImp_().updateTemperature_(priVars,
-                                    element,
-                                    elemGeom,
-                                    scvIdx,
-                                    problem);
-        
-        // material law parameters
         const MaterialLawParams &materialParams =
             problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
 
-        Model::completeFluidState(priVars, materialParams, fluidState_);
+        Model::completeFluidState(priVars, problem, element, elemGeom, scvIdx, materialParams, fluidState_);
 
         mobility_[wPhaseIdx] =
             MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx))
@@ -127,15 +120,8 @@ public:
                                                          elemGeom,
                                                          scvIdx);
 
-
-#warning "TODO: it would be nice to avoid creating the parameter cache a second time here," \
-    " but that requires to abandon/modify model.completeFluidState()."  \
-    " (alternatively a completeFluidState() method can be introduced in the 2pni model.)"
-        typename FluidSystem::ParameterCache paramCache;
-        paramCache.updateAll(fluidState_);
-
-        // energy related quantities
-        asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
+        // energy related quantities not belonging to the fluid state
+        asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol);
     }
 
     /*!
@@ -203,21 +189,10 @@ public:
     { return porosity_; }
 
 protected:
-    void updateTemperature_(const PrimaryVariables &priVars,
-                            const Element &element,
-                            const FVElementGeometry &elemGeom,
-                            int scvIdx,
-                            const Problem &problem)
-    {
-        fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
-    }
-
     /*!
      * \brief Called by update() to compute the energy related quantities
      */
-    template <class ParameterCache>
-    void updateEnergy_(ParameterCache &paramCache,
-                       const PrimaryVariables &sol,
+    void updateEnergy_(const PrimaryVariables &sol,
                        const Problem &problem,
                        const Element &element,
                        const FVElementGeometry &elemGeom,

dumux/boxmodels/2pni/2pnimodel.hh

@@ -69,6 +69,28 @@ namespace Dumux {
 template<class TypeTag>
 class TwoPNIModel: public TwoPModel<TypeTag>
 {
+    friend class TwoPModel<TypeTag>;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+    template<class PrimaryVariables, class Problem, class Element, class FVElementGeometry>
+    static Scalar temperature_(const PrimaryVariables &priVars,
+                            const Problem& problem,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx)
+    {
+        return priVars[Indices::temperatureIdx];
+    }
+
+    template<class FluidState, class ParameterCache>
+    static Scalar enthalpy_(const FluidState& fluidState,
+                            const ParameterCache& paramCache,
+                            int phaseIdx)
+    {
+        typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+        return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+    }
 };
 
 }

dumux/boxmodels/2pni/2pnivolumevariables.hh

@@ -95,45 +95,16 @@ protected:
     // is protected, we are friends with our parent..
     friend class TwoPVolumeVariables<TypeTag>;
 
-    /*!
-     * \brief Update the temperature for a given control volume.
-     *
-     * \param priVars The local primary variable vector
-     * \param element The current element
-     * \param elemGeom The finite-volume geometry in the box scheme
-     * \param scvIdx The local index of the SCV (sub-control volume)
-     * \param problem The problem object
-     *
-     */
-    void updateTemperature_(const PrimaryVariables &priVars,
-                            const Element &element,
-                            const FVElementGeometry &elemGeom,
-                            int scvIdx,
-                            const Problem &problem)
-    {
-        // retrieve temperature from primary variables
-        this->fluidState_.setTemperature(priVars[temperatureIdx]);
-    }
-    
     /*!
      * \brief Called by update() to compute the energy related quantities
      */
-    template <class ParameterCache>
-    void updateEnergy_(ParameterCache &paramCache,
-                       const PrimaryVariables &sol,
+    void updateEnergy_(const PrimaryVariables &sol,
                        const Problem &problem,
                        const Element &element,
                        const FVElementGeometry &elemGeom,
                        int scvIdx,
                        bool isOldSol)
     {
-        // copmute and set the internal energies of the fluid phases
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx);
-
-            this->fluidState_.setEnthalpy(phaseIdx, h);
-        }
-
         // copmute and set the heat capacity of the solid phase
         heatCapacity_ = problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
         Valgrind::CheckDefined(heatCapacity_);