From 93c17d1d78e060c429ce294d9d826811497ffdfe Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 17 May 2013 15:10:14 +0000
Subject: [PATCH] bugfix: correct calculation of effective thermal conductivity
 for 2p non-isothermal cell-centered models

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10703 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/implicit/2p2cni/2p2cnifluxvariables.hh | 37 +++++++++++++++++---
 dumux/implicit/2pni/2pnifluxvariables.hh     | 37 +++++++++++++++++---
 2 files changed, 66 insertions(+), 8 deletions(-)

diff --git a/dumux/implicit/2p2cni/2p2cnifluxvariables.hh b/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
index 4bae8dcdcf..6249ebf176 100644
--- a/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
+++ b/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
@@ -145,19 +145,48 @@ protected:
     {
         const unsigned i = this->face().i;
         const unsigned j = this->face().j;
+        Scalar lambdaI, lambdaJ;
 
-        const Scalar lambdaI =
-            ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            lambdaI =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, i),
                                                                    problem.spatialParams().porosity(element, this->fvGeometry_, i));
-        const Scalar lambdaJ =
-            ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
+            lambdaJ =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, j),
                                                                    problem.spatialParams().porosity(element, this->fvGeometry_, j));
+        }
+        else
+        {
+            const Element& elementI = *this->fvGeometry_.neighbors[i];
+            FVElementGeometry fvGeometryI;
+            fvGeometryI.subContVol[0].global = elementI.geometry().center();
+
+            lambdaI =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
+                                                                   elemVolVars[i].thermalConductivity(wPhaseIdx),
+                                                                   elemVolVars[i].thermalConductivity(nPhaseIdx),
+                                                                   problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
+                                                                   problem.spatialParams().porosity(elementI, fvGeometryI, 0));
+
+            const Element& elementJ = *this->fvGeometry_.neighbors[j];
+            FVElementGeometry fvGeometryJ;
+            fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+
+            lambdaJ =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
+                                                                   elemVolVars[j].thermalConductivity(wPhaseIdx),
+                                                                   elemVolVars[j].thermalConductivity(nPhaseIdx),
+                                                                   problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
+                                                                   problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
+        }
+
         // -> harmonic mean
         lambdaEff_ = harmonicMean(lambdaI, lambdaJ);
     }
diff --git a/dumux/implicit/2pni/2pnifluxvariables.hh b/dumux/implicit/2pni/2pnifluxvariables.hh
index 71f5763bfc..8e6025e484 100644
--- a/dumux/implicit/2pni/2pnifluxvariables.hh
+++ b/dumux/implicit/2pni/2pnifluxvariables.hh
@@ -151,19 +151,48 @@ protected:
     {
         const unsigned i = this->face().i;
         const unsigned j = this->face().j;
+        Scalar lambdaI, lambdaJ;
 
-        const Scalar lambdaI =
-            ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            lambdaI =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, i),
                                                                    problem.spatialParams().porosity(element, this->fvGeometry_, i));
-        const Scalar lambdaJ =
-            ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
+            lambdaJ =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, j),
                                                                    problem.spatialParams().porosity(element, this->fvGeometry_, j));
+        }
+        else
+        {
+            const Element& elementI = *this->fvGeometry_.neighbors[i];
+            FVElementGeometry fvGeometryI;
+            fvGeometryI.subContVol[0].global = elementI.geometry().center();
+
+            lambdaI =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
+                                                                   elemVolVars[i].thermalConductivity(wPhaseIdx),
+                                                                   elemVolVars[i].thermalConductivity(nPhaseIdx),
+                                                                   problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
+                                                                   problem.spatialParams().porosity(elementI, fvGeometryI, 0));
+
+            const Element& elementJ = *this->fvGeometry_.neighbors[j];
+            FVElementGeometry fvGeometryJ;
+            fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
+
+            lambdaJ =
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
+                                                                   elemVolVars[j].thermalConductivity(wPhaseIdx),
+                                                                   elemVolVars[j].thermalConductivity(nPhaseIdx),
+                                                                   problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
+                                                                   problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
+        }
+
         // -> harmonic mean
         lambdaEff_ = harmonicMean(lambdaI, lambdaJ);
     }
-- 
GitLab