From cd39d3dea8df3abf2e60e89d79b9e94426cc1a8a Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 16 Dec 2011 13:23:27 +0000
Subject: [PATCH] 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               | 53 +++++++++++++++++----
 dumux/boxmodels/2p/2pvolumevariables.hh     | 33 ++-----------
 dumux/boxmodels/2pni/2pnimodel.hh           | 22 +++++++++
 dumux/boxmodels/2pni/2pnivolumevariables.hh | 31 +-----------
 4 files changed, 70 insertions(+), 69 deletions(-)

diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index 3035b1ca5e..62d4047d9a 100644
--- a/dumux/boxmodels/2p/2pmodel.hh
+++ b/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;
+    }
 };
 }
 
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index 50ec393383..4e15289f39 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/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,
diff --git a/dumux/boxmodels/2pni/2pnimodel.hh b/dumux/boxmodels/2pni/2pnimodel.hh
index f402b4b288..088c00b520 100644
--- a/dumux/boxmodels/2pni/2pnimodel.hh
+++ b/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);
+    }
 };
 
 }
diff --git a/dumux/boxmodels/2pni/2pnivolumevariables.hh b/dumux/boxmodels/2pni/2pnivolumevariables.hh
index 5399ac08ac..c3f0f19926 100644
--- a/dumux/boxmodels/2pni/2pnivolumevariables.hh
+++ b/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_);
-- 
GitLab