diff --git a/dumux/material/binarycoefficients/h2o_heavyoil.hh b/dumux/material/binarycoefficients/h2o_heavyoil.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a6df04b32ae57c3282bbe3075495c107d9eee147
--- /dev/null
+++ b/dumux/material/binarycoefficients/h2o_heavyoil.hh
@@ -0,0 +1,101 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * Copyright (C) 2010 by Andreas Lauser *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Binary coefficients for water and tce.
+ */
+#ifndef DUMUX_BINARY_COEFF_H2O_HEAVYOIL_HH
+#define DUMUX_BINARY_COEFF_H2O_HEAVYOIL_HH
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/heavyoil.hh>
+
+namespace Dumux
+{
+namespace BinaryCoeff
+{
+
+/*!
+ * \brief Binary coefficients for water and heavy oil as in SAGD processes
+ */
+class H2O_HeavyOil
+{
+public:
+ /*!
+ * \brief Henry coefficent \f$[N/m^2]\f$ for heavy oil in liquid water.
+ *
+ * See:
+ *
+ */
+
+ template <class Scalar>
+ static Scalar henryOilInWater(Scalar temperature)
+ {
+ // values copied from TCE, TODO: improve this!!
+ Scalar dumuxH = 1.5e-1 / 101.325; // unit [(mol/m^3)/Pa]
+ dumuxH *= 18.02e-6; //multiplied by molar volume of reference phase = water
+ return 1.0/dumuxH; // [Pa]
+ }
+
+ /*!
+ * \brief Henry coefficent \f$[N/m^2]\f$ for water in liquid heavy oil.
+ *
+ * See:
+ *
+ */
+
+ template <class Scalar>
+ static Scalar henryWaterInOil(Scalar temperature)
+ {
+ // arbitrary, TODO: improve it!!
+ return 1.0e8; // [Pa]
+ }
+
+
+ /*!
+ * \brief Binary diffusion coefficent [m^2/s] for molecular water and heavy oil.
+ *
+ */
+ template <class Scalar>
+ static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ return 1e-6; // [m^2/s] This is just an order of magnitude. Please improve it!
+ }
+
+ /*!
+ * \brief Diffusion coefficent [m^2/s] for tce in liquid water.
+ *
+ * \todo
+ */
+ template <class Scalar>
+ static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ return 1.e-9; // This is just an order of magnitude. Please improve it!
+ }
+};
+
+}
+} // end namespace
+
+#endif
diff --git a/dumux/material/components/heavyoil.hh b/dumux/material/components/heavyoil.hh
new file mode 100644
index 0000000000000000000000000000000000000000..12107f094923d857a247fde66bbe7a6a04c0d4ce
--- /dev/null
+++ b/dumux/material/components/heavyoil.hh
@@ -0,0 +1,501 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \ingroup Components
+ *
+ * \brief Properties of heavyoil.
+ *
+ */
+#ifndef DUMUX_HEAVYOIL_HH
+#define DUMUX_HEAVYOIL_HH
+
+#include <dumux/material/idealgas.hh>
+#include <dumux/material/components/component.hh>
+#include <dumux/material/constants.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief heavyoil
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class HeavyOil : public Component<Scalar, HeavyOil<Scalar> >
+{
+ typedef Dumux::Constants<Scalar> Consts;
+
+public:
+ /*!
+ * \brief A human readable name for heavyoil
+ */
+ static const char *name()
+ { return "heavyoil"; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of heavyoil
+ */
+ constexpr static Scalar molarMass()
+ { return .350; }
+
+ /*!
+ * \brief The MolecularWeight in \f$\mathrm{[kg/mol]}\f$ of refComponent
+ */
+ constexpr static Scalar refComponentMolecularWeight()
+ { return .400; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of heavyoil
+ */
+
+ constexpr static Scalar molecularWeight()
+ { return .350; }
+
+ /*!
+ * \brief The Specific Gravity \f$\mathrm{[ ]}\f$ of heavyoil
+ */
+ constexpr static Scalar specificGravity()
+ { return 0.91; }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ at heavyoil's triple point.
+ */
+ static Scalar tripleTemperature()
+ {
+ DUNE_THROW(Dune::NotImplemented, "tripleTemperature for heavyoil");
+ }
+
+ /*!
+ * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at heavyoil's triple point.
+ */
+ static Scalar triplePressure()
+ {
+ DUNE_THROW(Dune::NotImplemented, "triplePressure for heavyoil");
+ }
+
+ static Scalar refComponentSpecificGravity()
+ {
+ const Scalar A = 0.83;
+ const Scalar B = 89.9513;
+ const Scalar C = 139.6612;
+ const Scalar D = 3.2033;
+ const Scalar E = 1.0564;
+
+ const Scalar mW = refComponentMolecularWeight() *1000. ; // in [g/mol];
+
+ return A+(B/mW)-(C/std::pow((mW+D),E));
+
+ }
+
+ static Scalar perbutationFactorBoilingTemperature()
+ {
+ const Scalar A = -7.4120e-2; //All factors for 1 atm / 101325 pascals [760 mmHg]
+ const Scalar B = -7.5041e-3;
+ const Scalar C = -2.6031;
+ const Scalar D = 9.0180e-2;
+ const Scalar E = -1.0482;
+
+ Scalar deltaSpecificGravity = std::log(refComponentSpecificGravity()/specificGravity());
+ Scalar deltaMolecularWeight = std::log(refComponentMolecularWeight()/molecularWeight());
+
+ return A*std::pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*std::pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+ + E*deltaSpecificGravity*deltaMolecularWeight;
+
+ }
+
+ static Scalar perbutationFactorCriticalTemperature()
+ {
+ const Scalar A = -6.1294e-2;
+ const Scalar B = -7.0862e-2;
+ const Scalar C = 6.1976e-1;
+ const Scalar D = -5.7090e-2;
+ const Scalar E = -8.4583e-2;
+
+ Scalar deltaSpecificGravity = std::log(refComponentSpecificGravity()/specificGravity());
+ Scalar deltaMolecularWeight = std::log(refComponentMolecularWeight()/molecularWeight());
+
+ return A*std::pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*std::pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+ + E*deltaSpecificGravity*deltaMolecularWeight;
+
+ }
+
+ static Scalar perbutationFactorCriticalPressure()
+ {
+ const Scalar A = 1.8270e-1;
+ const Scalar B = -2.4864e-1;
+ const Scalar C = 8.3611;
+ const Scalar D = -2.2389e-1;
+ const Scalar E = 2.6984;
+
+ Scalar deltaSpecificGravity = std::log(refComponentSpecificGravity()/specificGravity());
+ Scalar deltaMolecularWeight = std::log(refComponentMolecularWeight()/molecularWeight());
+
+ return A*std::pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*std::pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+ + E*deltaSpecificGravity*deltaMolecularWeight;
+
+ }
+
+ static Scalar refComponentBoilingTemperature()
+ {
+ const Scalar A = 477.63; //All factors for 1 atm / 101325 pascals [760 mmHg]
+ const Scalar B = 88.51;
+ const Scalar C = 1007;
+ const Scalar D = 1214.40;
+
+ return A*std::log((1000.*refComponentMolecularWeight() + B)/(1000.*refComponentMolecularWeight()+C)) + D;
+
+ }
+
+ static Scalar refComponentCriticalTemperature()
+ {
+ const Scalar A = 226.50;
+ const Scalar B = 6.78;
+ const Scalar C = 1.282e6;
+ const Scalar D = 2668;
+
+ return A*std::log((1000.*refComponentMolecularWeight() + B)/(1000.*refComponentMolecularWeight()+C)) + D ;
+
+ }
+
+ static Scalar refComponentCriticalPressure()
+ {
+ const Scalar A = 141.20;
+ const Scalar B = 45.66e-2;
+ const Scalar C = 16.59e-3;
+ const Scalar D = 2.19;
+
+ return (A*1000.*molecularWeight())/(std::pow(B + (C*1000.*molecularWeight()),D)) ;
+
+ }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ at heavyoil's boiling point (1 atm)
+ */
+ static Scalar boilingTemperature()
+ {
+
+ return refComponentBoilingTemperature() * std::pow((1 + 2*perbutationFactorBoilingTemperature())/(1 - 2*perbutationFactorBoilingTemperature()),2);
+
+ }
+
+ /*!
+ * \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of heavyoil
+ */
+ static Scalar criticalTemperature()
+ {
+
+ return refComponentCriticalTemperature() * std::pow((1 + 2*perbutationFactorCriticalTemperature())/(1 - 2*perbutationFactorCriticalTemperature()),2);
+
+ }
+
+ /*!
+ * \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of heavyoil
+ */
+ static Scalar criticalPressure()
+ {
+
+ return refComponentCriticalPressure() * std::pow((1 + 2*perbutationFactorCriticalPressure())/(1 - 2*perbutationFactorCriticalPressure()),2);
+
+ }
+
+ /*!
+ * \brief The saturation vapor pressure in \f$\mathrm{[Pa]}\f$ of
+ *
+ *
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ */
+ static Scalar vaporPressure(Scalar temperature)
+ {
+ const Scalar A = 8.25990;
+ const Scalar B = 2830.065;
+ const Scalar C = 42.95101;
+ /*const Scalar A = 7.00909;;
+ const Scalar B = 1462.266;;
+ const Scalar C = 215.110;;*/
+ //const Scalar A = 6.90027; //n-Heptane http://tinyurl.com/lpo2h3s
+ //const Scalar B = 1266.871;
+ //const Scalar C = 216.757;
+
+ Scalar T = temperature - 273.15;
+ return 100*1.334*std::pow(10.0, (A - (B/(T + C)))); // in [Pa]
+
+ // return value2;
+
+ }
+
+
+ /* P = 100 * 1.334 * (10.0)^(A - (B / (T + C)) =>
+ * P/(100*1.334)=(10.0)^(A - (B / (T + C)) =>
+ * P/133.4=(10.0)^(A - (B / (T + C)) =>
+ * log(10)(P/133.4)=A - B / (T + C) =>
+ * B / (T + C) =A-log(10) (P/133.4) =>
+ * B/(A-log(10) (P/133.4))=T+C =>
+ * T=B/(A-log(10) (P/133.4))-C
+ */
+
+ static Scalar vaporTemperature(Scalar pressure)
+ {
+ const Scalar A = 8.25990;
+ const Scalar B = 2830.065;
+ const Scalar C = 42.95101;
+
+ const Scalar P = pressure;
+
+ return Scalar ((B/(A-std::log10(P/100*1.334)))-C); //T=B/(A-log(10) (P/133.4))-C
+ //std::log(arg) / std::log(base);
+ }
+
+ /*!
+ * \brief Specific enthalpy of liquid heavyoil \f$\mathrm{[J/kg]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidEnthalpy(const Scalar temperature,
+ const Scalar pressure)
+ {
+ // Gauss quadrature rule:
+ // Interval: [0K; temperature (K)]
+ // Gauss-Legendre-Integration with variable transformation:
+ // \int_a^b f(T) dT \approx (b-a)/2 \sum_i=1^n \alpha_i f( (b-a)/2 x_i + (a+b)/2 )
+ // with: n=2, legendre -> x_i = +/- \sqrt(1/3), \apha_i=1
+ // here: a=273.15K, b=actual temperature in Kelvin
+ // \leadsto h(T) = \int_273.15^T c_p(T) dT
+ // \approx 0.5 (T-273.15) * (cp( 0.5(temperature-273.15)sqrt(1/3) ) + cp(0.5(temperature-273.15)(-1)sqrt(1/3))
+
+ // Enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input and is
+ // fit over a certain temperature range. This suggests choosing an interval of integration being in the actual fit range.
+ // I.e. choosing T=273.15K as reference point for liquid enthalpy.
+
+ const Scalar sqrt1over3 = std::sqrt(1./3.);
+ const Scalar TEval1 = 0.5*(temperature-273.15)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration
+ const Scalar TEval2 = 0.5*(temperature-273.15)* (-1)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration
+
+ const Scalar h_n = 0.5 * (temperature-273.15) * ( liquidHeatCapacity(TEval1, pressure) + liquidHeatCapacity(TEval2, pressure) ) ;
+
+ return h_n;
+ }
+
+ /*!
+ * \brief Latent heat of vaporization for heavyoil \f$\mathrm{[J/kg]}\f$.
+ *
+ * source : Reid et al. (fourth edition): Chen method (chap. 7-11, Delta H_v = Delta H_v (T) according to chap. 7-12)
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar heatVap(Scalar temperature,
+ const Scalar pressure)
+ {
+ temperature = std::min(temperature, criticalTemperature()); // regularization
+ temperature = std::max(temperature, 0.0); // regularization
+
+ Scalar T_crit = criticalTemperature();
+ Scalar Tr1 = boilingTemperature()/criticalTemperature();
+ Scalar p_crit = criticalPressure();
+
+ // Chen method, eq. 7-11.4 (at boiling)
+ const Scalar DH_v_boil = Consts::R * T_crit * Tr1
+ * (3.978 * Tr1 - 3.958 + 1.555*std::log(p_crit * 1e-5 /*Pa->bar*/ ) )
+ / (1.07 - Tr1); /* [J/mol] */
+
+ /* Variation with temp according to Watson relation eq 7-12.1*/
+ const Scalar Tr2 = temperature/criticalTemperature();
+ const Scalar n = 0.375;
+ const Scalar DH_vap = DH_v_boil * std::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
+
+ return (DH_vap/molarMass()); // we need [J/kg]
+ }
+
+
+ /*!
+ * \brief Specific enthalpy of heavyoil vapor \f$\mathrm{[J/kg]}\f$.
+ *
+ * This relation is true on the vapor pressure curve, i.e. as long
+ * as there is a liquid phase present.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
+ {
+ return liquidEnthalpy(temperature,pressure) + heatVap(temperature, pressure);
+ }
+
+ /*!
+ * \brief
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
+ {
+ return IdealGas<Scalar>::density(molarMass(),
+ temperature,
+ pressure);
+ }
+
+ /*!
+ * \brief The density of pure heavyoil at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ // return molarLiquidDensity_(temperature)*molarMass(); // [kg/m^3]
+
+ /* according to Lashanizadegan et al (2008) in Chemical Engineering Communications: */
+ /* Simultaneous Heat and Fluid Flow in Porous Media: Case Study: Steam Injection for Tertiary Oil Recovery */
+ Scalar rhoReference = 906.; // [kg/m^3] at reference pressure and temperature
+ Scalar compressCoeff = 1.e-8; // just a value without justification
+ Scalar expansCoeff = 1.e-7; // also just a value
+ Scalar rho = rhoReference * (1. + (pressure - 1.e5)*compressCoeff) * (1. - (temperature - 293.)*expansCoeff);
+
+ return rho; // [kg/m^3]
+ }
+
+ /*!
+ * \brief Returns true iff the gas phase is assumed to be compressible
+ */
+ static bool gasIsCompressible()
+ { return true; }
+
+ /*!
+ * \brief Returns true iff the gas phase is assumed to be ideal
+ */
+ static bool gasIsIdeal()
+ { return true; }
+
+ /*!
+ * \brief Returns true iff the liquid phase is assumed to be compressible
+ */
+ static bool liquidIsCompressible()
+ { return true; }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of heavyoil vapor
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ * \param regularize defines, if the functions is regularized or not, set to true by default
+ */
+ static Scalar gasViscosity(Scalar temperature, Scalar pressure, bool regularize=true)
+ {
+ temperature = std::min(temperature, 500.0); // regularization
+ temperature = std::max(temperature, 250.0);
+
+ // reduced temperature
+ Scalar Tr = temperature/criticalTemperature();
+
+ Scalar Fp0 = 1.0;
+ Scalar xi = 0.00474;
+ Scalar eta_xi =
+ Fp0*(0.807*std::pow(Tr,0.618)
+ - 0.357*std::exp(-0.449*Tr)
+ + 0.34*std::exp(-4.058*Tr)
+ + 0.018);
+
+ return eta_xi/xi/1e7; // [Pa s]
+ }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure heavyoil.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ // static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ // {
+
+ /* according to Lashanizadegan et al (2008) in Chemical Engineering Communications: */
+ /* Simultaneous Heat and Fluid Flow in Porous Media: Case Study: Steam Injection for Tertiary Oil Recovery */
+
+ //return 1027919.422*std::exp(-0.04862*temperature); // [Pa s]
+
+ //according to http://www.ecltechnology.com/subsur/reports/pvt_tgb.pdf[Page 10]
+
+ // Scalar temperatureFahrenheit = ((temperature-273.15)*1.8)+32;
+ // Scalar API = 15;
+ // return (std::pow(10,0.052*std::pow(API,2)-2.2704*API-5.7567)*std::pow(temperatureFahrenheit,-0.0222*std::pow(API,2)+0.9415*API-12.839))*0.001; // [Pa s]
+
+ // }
+
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ {
+
+ /* according to Lashanizadegan et al (2008) in Chemical Engineering Communications: */
+ /* Simultaneous Heat and Fluid Flow in Porous Media: Case Study: Steam Injection for Tertiary Oil Recovery */
+
+ //return 1027919.422*std::exp(-0.04862*temperature); // [Pa s]
+
+ //according to http://www.ecltechnology.com/subsur/reports/pvt_tgb.pdf [Page 10]
+ Scalar temperatureFahrenheit = (9/5)*(temperature-273.15)+32;
+ Scalar API = 9;
+ return ((std::pow(10,0.10231*std::pow(API,2)-3.9464*API+46.5037))*(std::pow(temperatureFahrenheit,-0.04542*std::pow(API,2)+1.70405*API-19.18)))*0.001;
+
+ }
+ /*!
+ * \brief Specific heat cap of liquid heavyoil \f$\mathrm{[J/kg]}\f$.
+ *
+ * source : Reid et al. (fourth edition): Missenard group contrib. method (chap 5-7, Table 5-11, s. example 5-8)
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar liquidHeatCapacity(const Scalar temperature,
+ const Scalar pressure)
+ {
+ /* according to Lashanizadegan et al (2008) in Chemical Engineering Communications: */
+ /* Simultaneous Heat and Fluid Flow in Porous Media: Case Study: Steam Injection for Tertiary Oil Recovery */
+
+ return 618.; // J/(kg K)
+ }
+
+protected:
+ /*!
+ * \brief The molar density of pure heavyoil at a given pressure and temperature
+ * \f$\mathrm{[mol/m^3]}\f$.
+ *
+ * source : Reid et al. (fourth edition): Modified Racket technique (chap. 3-11, eq. 3-11.9)
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ */
+ static Scalar molarLiquidDensity_(Scalar temperature)
+ {
+ temperature = std::min(temperature, 500.0); // regularization
+ temperature = std::max(temperature, 250.0);
+
+ const Scalar Z_RA = 0.2556; // from equation
+ const Scalar expo = 1.0 + std::pow(1.0 - temperature/criticalTemperature(), 2.0/7.0);
+ Scalar V = Consts::R*criticalTemperature()/criticalPressure()*std::pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
+
+ return 1.0/V; // molar density [mol/m^3]
+ }
+
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh b/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9a53594d4bf2ddcb1597211281f84c99d03bc303
--- /dev/null
+++ b/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh
@@ -0,0 +1,475 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A fluid system with water, gas and NAPL as phases and
+ * \f$H_2O\f$ and \f$NAPL (contaminant)\f$ as components.
+ * It uses heavyoil Properties, but allows for a scaling of density
+ * thus enabling an artifical DNAPL if desired
+ */
+#ifndef DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
+#define DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
+
+#include <dumux/material/idealgas.hh>
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/heavyoil.hh>
+
+#include <dumux/material/binarycoefficients/h2o_heavyoil.hh>
+
+#include <dumux/material/fluidsystems/base.hh>
+
+namespace Dumux
+{
+namespace FluidSystems
+{
+
+/*!
+ * \brief A compositional fluid with water and heavy oil
+ * components in both, the liquid and the gas phase.
+ */
+template <class TypeTag, class Scalar,
+ class H2OType = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > >
+class H2OHeavyOil
+ : public BaseFluidSystem<Scalar, H2OHeavyOil<TypeTag, Scalar, H2OType> >
+{
+ typedef H2OHeavyOil<TypeTag, Scalar, H2OType> ThisType;
+ typedef BaseFluidSystem<Scalar, ThisType> Base;
+
+public:
+ typedef Dumux::HeavyOil<Scalar> HeavyOil;
+ typedef H2OType H2O;
+
+
+ static const int numPhases = 3;
+ static const int numComponents = 2;
+
+ static const int wPhaseIdx = 0; // index of the water phase
+ static const int nPhaseIdx = 1; // index of the NAPL phase
+ static const int gPhaseIdx = 2; // index of the gas phase
+
+ static const int H2OIdx = 0;
+ static const int NAPLIdx = 1;
+
+ // export component indices to indicate the main component
+ // of the corresponding phase at atmospheric pressure 1 bar
+ // and room temperature 20°C:
+ static const int wCompIdx = H2OIdx;
+ static const int nCompIdx = NAPLIdx;
+
+ /*!
+ * \brief Initialize the fluid system's static parameters generically
+ *
+ * If a tabulated H2O component is used, we do our best to create
+ * tables that always work.
+ */
+ static void init()
+ {
+ init(/*tempMin=*/273.15,
+ /*tempMax=*/623.15,
+ /*numTemp=*/100,
+ /*pMin=*/0.0,
+ /*pMax=*/20e6,
+ /*numP=*/200);
+ }
+
+ /*!
+ * \brief Initialize the fluid system's static parameters using
+ * problem specific temperature and pressure ranges
+ *
+ * \param tempMin The minimum temperature used for tabulation of water [K]
+ * \param tempMax The maximum temperature used for tabulation of water [K]
+ * \param nTemp The number of ticks on the temperature axis of the table of water
+ * \param pressMin The minimum pressure used for tabulation of water [Pa]
+ * \param pressMax The maximum pressure used for tabulation of water [Pa]
+ * \param nPress The number of ticks on the pressure axis of the table of water
+ */
+ static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
+ Scalar pressMin, Scalar pressMax, unsigned nPress)
+ {
+ if (H2O::isTabulated) {
+ std::cout << "Initializing tables for the H2O fluid properties ("
+ << nTemp*nPress
+ << " entries).\n";
+
+ H2O::init(tempMin, tempMax, nTemp,
+ pressMin, pressMax, nPress);
+ }
+ }
+
+
+ /*!
+ * \brief Return whether a phase is liquid
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static bool isLiquid(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return phaseIdx != gPhaseIdx;
+ }
+
+ static bool isIdealGas(int phaseIdx)
+ { return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && HeavyOil::gasIsIdeal(); }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be an ideal mixture.
+ *
+ * We define an ideal mixture as a fluid phase where the fugacity
+ * coefficients of all components times the pressure of the phase
+ * are indepent on the fluid composition. This assumtion is true
+ * if Henry's law and Rault's law apply. If you are unsure what
+ * this function should return, it is safe to return false. The
+ * only damage done will be (slightly) increased computation times
+ * in some cases.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static bool isIdealMixture(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ // we assume Henry's and Rault's laws for the water phase and
+ // and no interaction between gas molecules of different
+ // components, so all phases are ideal mixtures!
+ return true;
+ }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be compressible.
+ *
+ * Compressible means that the partial derivative of the density
+ * to the fluid pressure is always larger than zero.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static bool isCompressible(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ // gases are always compressible
+ if (phaseIdx == gPhaseIdx)
+ return true;
+ else if (phaseIdx == wPhaseIdx)
+ // the water component decides for the water phase...
+ return H2O::liquidIsCompressible();
+
+ // the NAPL component decides for the napl phase...
+ return HeavyOil::liquidIsCompressible();
+ }
+
+ /*!
+ * \brief Return the human readable name of a phase (used in indices)
+ */
+ static const char *phaseName(int phaseIdx)
+ {
+ switch (phaseIdx) {
+ case wPhaseIdx: return "w";
+ case nPhaseIdx: return "n";
+ case gPhaseIdx: return "g";;
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Return the human readable name of a component (used in indices)
+ */
+ static const char *componentName(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::name();
+ case NAPLIdx: return HeavyOil::name();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Return the molar mass of a component in [kg/mol].
+ */
+ static Scalar molarMass(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::molarMass();
+ case NAPLIdx: return HeavyOil::molarMass();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the phase
+ * density [kg/m^3].
+ */
+ using Base::density;
+ template <class FluidState>
+ static Scalar density(const FluidState &fluidState, int phaseIdx)
+ {
+ if (phaseIdx == wPhaseIdx) {
+ // See: doctoral thesis of Steffen Ochs 2007
+ // Steam injection into saturated porous media : process analysis including experimental and numerical investigations
+ // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf
+ Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ Scalar clH2O = rholH2O/H2O::molarMass();
+
+ // this assumes each dissolved molecule displaces exactly one
+ // water molecule in the liquid
+ return
+ clH2O*(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
+ +
+ HeavyOil::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
+ }
+ else if (phaseIdx == nPhaseIdx) {
+ // assume pure NAPL for the NAPL phase
+ Scalar pressure = HeavyOil::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
+ return HeavyOil::liquidDensity(fluidState.temperature(phaseIdx), pressure);
+ }
+
+ assert (phaseIdx == gPhaseIdx);
+ Scalar pH2O =
+ fluidState.moleFraction(gPhaseIdx, H2OIdx) *
+ fluidState.pressure(gPhaseIdx);
+ Scalar pNAPL =
+ fluidState.moleFraction(gPhaseIdx, NAPLIdx) *
+ fluidState.pressure(gPhaseIdx);
+ return
+ H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) +
+ HeavyOil::gasDensity(fluidState.temperature(phaseIdx), pNAPL);
+ }
+
+ /*!
+ * \brief Return the viscosity of a phase.
+ */
+ using Base::viscosity;
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ if (phaseIdx == wPhaseIdx) {
+ // assume pure water viscosity
+ return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
+ fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == nPhaseIdx) {
+ // assume pure NAPL viscosity
+ return HeavyOil::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+
+ assert (phaseIdx == gPhaseIdx);
+
+ /* Wilke method. See:
+ *
+ * See: R. Reid, et al.: The Properties of Gases and Liquids,
+ * 4th edition, McGraw-Hill, 1987, 407-410
+ * 5th edition, McGraw-Hill, 20001, p. 9.21/22
+ *
+ * in this case, we use a simplified version in order to avoid
+ * computationally costly evaluation of sqrt and pow functions and
+ * divisions
+ * -- compare e.g. with Promo Class p. 32/33
+ */
+ const Scalar mu[numComponents] = {
+ H2O::gasViscosity(fluidState.temperature(phaseIdx), H2O::vaporPressure(fluidState.temperature(phaseIdx))),
+ HeavyOil::gasViscosity(fluidState.temperature(phaseIdx), HeavyOil::vaporPressure(fluidState.temperature(phaseIdx)))
+ };
+
+ return mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
+ + mu[NAPLIdx]*fluidState.moleFraction(gPhaseIdx, NAPLIdx);
+ }
+
+
+ using Base::binaryDiffusionCoefficient;
+ template <class FluidState>
+ static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx)
+ {
+ const Scalar T = fluidState.temperature(phaseIdx);
+ const Scalar p = fluidState.pressure(phaseIdx);
+
+ // liquid phase
+ if (phaseIdx == wPhaseIdx)
+ return BinaryCoeff::H2O_HeavyOil::liquidDiffCoeff(T, p);
+
+ // gas phase
+ else if (phaseIdx == gPhaseIdx)
+ return BinaryCoeff::H2O_HeavyOil::gasDiffCoeff(T, p);
+
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent diffusion coefficient for phase index " << phaseIdx);
+ }
+
+ /* Henry coefficients
+ */
+ template <class FluidState>
+ static Scalar henryCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ const Scalar p = fluidState.pressure(phaseIdx);
+
+ if (compIdx == NAPLIdx && phaseIdx == wPhaseIdx)
+ return Dumux::BinaryCoeff::H2O_HeavyOil::henryOilInWater(T)/p;
+
+ else if (phaseIdx == nPhaseIdx && compIdx == H2OIdx)
+ return Dumux::BinaryCoeff::H2O_HeavyOil::henryWaterInOil(T)/p;
+
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx
+ << " and component index " << compIdx);
+ }
+
+ /* partial pressures in the gas phase, taken from saturation vapor pressures
+ */
+ template <class FluidState>
+ static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ if (compIdx == NAPLIdx)
+ return HeavyOil::vaporPressure(T);
+ else if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+ }
+
+ /* inverse vapor pressures, taken from inverse saturation vapor pressures
+ */
+ template <class FluidState>
+ static Scalar inverseVaporPressureCurve(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar pressure = fluidState.pressure(phaseIdx);
+ if (compIdx == NAPLIdx)
+ return HeavyOil::vaporTemperature(pressure);
+ else if (compIdx == H2OIdx)
+ return H2O::vaporTemperature(pressure);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+ }
+
+
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the specific
+ * phase enthalpy [J/kg].
+ */
+ /*!
+ * \todo This system neglects the contribution of gas-molecules in the liquid phase.
+ * This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ...
+ */
+ using Base::enthalpy;
+ template <class FluidState>
+ static Scalar enthalpy(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ if (phaseIdx == wPhaseIdx) {
+ return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == nPhaseIdx) {
+ return HeavyOil::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+ else if (phaseIdx == gPhaseIdx) { // gas phase enthalpy depends strongly on composition
+ Scalar hgc = HeavyOil::gasEnthalpy(fluidState.temperature(phaseIdx),
+ fluidState.pressure(phaseIdx));
+ Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx),
+ fluidState.pressure(phaseIdx));
+
+ Scalar result = 0;
+ result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx);
+ result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx);
+
+ return result;
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ using Base::heatCapacity;
+ template <class FluidState>
+ static Scalar heatCapacity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()");
+ }
+
+ using Base::thermalConductivity;
+ template <class FluidState>
+ static Scalar thermalConductivity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::thermalConductivity()");
+ }
+
+private:
+
+};
+} // end namespace FluidSystems
+
+#ifdef DUMUX_PROPERTIES_HH
+// forward defintions of the property tags
+namespace Properties {
+ NEW_PROP_TAG(Scalar);
+ NEW_PROP_TAG(Components);
+}
+
+/*!
+ * \brief A threephase fluid system with water, heavyoil as components.
+ *
+ * This is an adapter to use Dumux::H2OheavyoilFluidSystem<TypeTag>, as is
+ * done with most other classes in Dumux.
+ * This fluidsystem is applied by default with the tabulated version of
+ * water of the IAPWS-formulation.
+ *
+ * To change the component formulation (ie to use nontabulated or
+ * incompressible water), or to switch on verbosity of tabulation,
+ * use the property system and the property "Components":
+ *
+ * // Select desired version of the component
+ * SET_PROP(myApplicationProperty, Components) : public GET_PROP(TypeTag, DefaultComponents)
+ * {
+ * typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ *
+ * // Do not use the defaults !
+ * // typedef Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > H2O;
+ *
+ * // Apply e.g. untabulated water:
+ * typedef Dumux::H2O<Scalar> H2O;
+ * };
+ *
+ * Also remember to initialize tabulated components (FluidSystem::init()), while this
+ * is not necessary for non-tabularized ones.
+ */
+template<class TypeTag>
+class H2OHeavyOilFluidSystem
+: public FluidSystems::H2OHeavyOil<TypeTag, typename GET_PROP_TYPE(TypeTag, Scalar),
+ typename GET_PROP(TypeTag, Components)::H2O>
+{};
+#endif
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/fluxvariables.hh b/dumux/porousmediumflow/3pwateroil/fluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f63df3bbafad2f16d663bf3ecc7147ca36a4fbd8
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/fluxvariables.hh
@@ -0,0 +1,315 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ * all fluxes of components over a face of a finite volume for
+ * the three-phase, two-component model.
+ */
+#ifndef DUMUX_3P2CNI_FLUX_VARIABLES_HH
+#define DUMUX_3P2CNI_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/spline.hh>
+
+#include "properties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \ingroup ImplicitFluxVariables
+ * \brief This template class contains the data which is required to
+ * calculate all fluxes of components over a face of a finite
+ * volume for the three-phase, two-component model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+template <class TypeTag>
+class ThreePWaterOilFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
+{
+ friend typename GET_PROP_TYPE(TypeTag, BaseFluxVariables); // be friends with base class
+ typedef typename GET_PROP_TYPE(TypeTag, BaseFluxVariables) BaseFluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+ typedef Dune::FieldVector<Scalar, dim> DimVector;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx,
+ };
+
+public:
+ /*!
+ * \brief Compute / update the flux variables
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param fIdx The local index of the SCV (sub-control-volume) face
+ * \param elemVolVars The volume variables of the current element
+ * \param onBoundary A boolean variable to specify whether the flux variables
+ * are calculated for interior SCV faces or boundary faces, default=false
+ */
+ void update(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int fIdx,
+ const ElementVolumeVariables &elemVolVars,
+ const bool onBoundary = false)
+ {
+ BaseFluxVariables::update(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary);
+ calculatePorousDiffCoeff_(problem, element, elemVolVars);
+
+ // The spatial parameters calculates the actual heat flux vector
+ DimVector temperatureGrad(0);
+ DimVector tmp(0.0);
+ problem.spatialParams().matrixHeatFlux(tmp, *this,
+ elemVolVars,
+ temperatureGrad,
+ element,
+ fvGeometry,
+ fIdx);
+ // project the heat flux vector on the face's normal vector
+ normalMatrixHeatFlux_ = tmp*this->face().normal;
+ }
+
+private:
+ void calculateGradients_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ // initialize to zero
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ density_[phaseIdx] = Scalar(0);
+ molarDensity_[phaseIdx] = Scalar(0);
+ massFractionCompWGrad_[phaseIdx] = Scalar(0);
+ massFractionCompNGrad_[phaseIdx] = Scalar(0);
+ moleFractionCompWGrad_[phaseIdx] = Scalar(0);
+ moleFractionCompNGrad_[phaseIdx] = Scalar(0);
+ }
+
+ BaseFluxVariables::calculateGradients_(problem, element, elemVolVars);
+
+ // calculate gradients
+ DimVector tmp(0.0);
+ for (int idx = 0; idx < this->face().numFap; idx++) // loop over adjacent vertices
+ {
+ // FE gradient at vertex idx
+ const DimVector &feGrad = this->face().grad[idx];
+
+ // index for the element volume variables
+ int volVarsIdx = this->face().fapIndices[idx];
+
+ // the concentration gradient of the components
+ // component in the phases
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, wCompIdx);
+ massFractionCompWGrad_[wPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, wCompIdx);
+ massFractionCompWGrad_[nPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, wCompIdx);
+ massFractionCompWGrad_[gPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, nCompIdx);
+ massFractionCompNGrad_[wPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, nCompIdx);
+ massFractionCompNGrad_[nPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, nCompIdx);
+ massFractionCompNGrad_[gPhaseIdx] += tmp;
+
+ // the molar concentration gradients of the components
+ // in the phases
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ moleFractionCompWGrad_[wPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
+ moleFractionCompWGrad_[nPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ moleFractionCompWGrad_[gPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ moleFractionCompNGrad_[wPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, nCompIdx);
+ moleFractionCompNGrad_[nPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, nCompIdx);
+ moleFractionCompNGrad_[gPhaseIdx] += tmp;
+
+ }
+ // calculate temperature gradient using finite element
+ // gradients
+ DimVector temperatureGrad(0);
+ for (int idx = 0; idx < this->face().numFap; idx++)
+ {
+ tmp = this->face().grad[idx];
+
+ // index for the element volume variables
+ int volVarsIdx = this->face().fapIndices[idx];
+
+ tmp *= elemVolVars[volVarsIdx].temperature();
+ temperatureGrad += tmp;
+ }
+ }
+
+ void calculatePorousDiffCoeff_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+
+ const VolumeVariables &volVarsI = elemVolVars[this->face().i];
+ const VolumeVariables &volVarsJ = elemVolVars[this->face().j];
+
+ Dune::FieldVector<Scalar, numPhases> diffusionCoefficientVector_i = volVarsI.diffusionCoefficient();
+ Dune::FieldVector<Scalar, numPhases> diffusionCoefficientVector_j = volVarsJ.diffusionCoefficient();
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // make sure to calculate only diffusion coefficents
+ // for phases which exist in both finite volumes
+ /* \todo take care: This should be discussed once again
+ * as long as a meaningful value can be found for the required mole fraction
+ * diffusion should work even without this one here */
+ if (volVarsI.saturation(phaseIdx) <= 0 ||
+ volVarsJ.saturation(phaseIdx) <= 0)
+ {
+ porousDiffCoeff_[phaseIdx] = 0.0;
+ continue;
+ }
+
+ // calculate tortuosity at the nodes i and j needed
+ // for porous media diffusion coefficient
+
+ Scalar tauI =
+ 1.0/(volVarsI.porosity() * volVarsI.porosity()) *
+ pow(volVarsI.porosity() * volVarsI.saturation(phaseIdx), 7.0/3);
+ Scalar tauJ =
+ 1.0/(volVarsJ.porosity() * volVarsJ.porosity()) *
+ pow(volVarsJ.porosity() * volVarsJ.saturation(phaseIdx), 7.0/3);
+ // Diffusion coefficient in the porous medium
+
+ // -> harmonic mean
+ porousDiffCoeff_[phaseIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientVector_i[phaseIdx],
+ volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientVector_j[phaseIdx]);
+
+ }
+ }
+
+public:
+ /*!
+ * \brief The diffusivity vector
+ */
+ Dune::FieldVector<Scalar, numPhases> porousDiffCoeff() const
+ { return porousDiffCoeff_; };
+
+ /*!
+ * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase.
+ */
+ Scalar density(int phaseIdx) const
+ { return density_[phaseIdx]; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase.
+ */
+ Scalar molarDensity(int phaseIdx) const
+ { return molarDensity_[phaseIdx]; }
+
+ const DimVector &massFractionCompWGrad(int phaseIdx) const
+ {return massFractionCompWGrad_[phaseIdx];}
+
+ const DimVector &massFractionCompNGrad(int phaseIdx) const
+ { return massFractionCompNGrad_[phaseIdx]; };
+
+ const DimVector &moleFractionCompWGrad(int phaseIdx) const
+ { return moleFractionCompWGrad_[phaseIdx]; };
+
+ const DimVector &moleFractionCompNGrad(int phaseIdx) const
+ { return moleFractionCompNGrad_[phaseIdx]; };
+
+ /*!
+ * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
+ * of the rock matrix over the sub-control volume's face in
+ * direction of the face normal.
+ */
+ Scalar normalMatrixHeatFlux() const
+ { return normalMatrixHeatFlux_; }
+
+
+protected:
+ // gradients
+ DimVector massFractionCompWGrad_[numPhases];
+ DimVector massFractionCompNGrad_[numPhases];
+ DimVector moleFractionCompWGrad_[numPhases];
+ DimVector moleFractionCompNGrad_[numPhases];
+
+ // density of each face at the integration point
+ Scalar density_[numPhases], molarDensity_[numPhases];
+
+ // the diffusivity matrix for the porous medium
+ Dune::FieldVector<Scalar, numPhases> porousDiffCoeff_;
+
+ Scalar normalMatrixHeatFlux_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/indices.hh b/dumux/porousmediumflow/3pwateroil/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cba433d9c67e9bb6d5795bb0f913d2ec58b25383
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/indices.hh
@@ -0,0 +1,80 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Defines the indices required for the 3p2cni model.
+ */
+#ifndef DUMUX_3P2CNI_INDICES_HH
+#define DUMUX_3P2CNI_INDICES_HH
+
+#include "properties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal 3p2cni model.
+ *
+ * \tparam formulation The formulation, only pgSwSn
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+class ThreePWaterOilIndices
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+ // Phase indices
+ static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< index of the wetting liquid phase
+ static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the nonwetting liquid phase
+ static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< index of the gas phase
+
+ // Component indices to indicate the main component
+ // of the corresponding phase at atmospheric pressure 1 bar
+ // and room temperature 20°C:
+ static const int wCompIdx = FluidSystem::wCompIdx;
+ static const int nCompIdx = FluidSystem::nCompIdx;
+
+ // present phases (-> 'pseudo' primary variable)
+ static const int threePhases = 1; //!< All three phases are present
+ static const int wPhaseOnly = 2; //!< Only the water phase is present
+ static const int gnPhaseOnly = 3; //!< Only gas and NAPL phases are present
+ static const int wnPhaseOnly = 4; //!< Only water and NAPL phases are present
+ static const int gPhaseOnly = 5; //!< Only gas phase is present
+ static const int wgPhaseOnly = 6; //!< Only water and gas phases are present
+
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for phase pressure in a solution vector
+ static const int switch1Idx = PVOffset + 1; //!< Index 1 of saturation, mole fraction or temperature
+ static const int switch2Idx = PVOffset + 2; //!< Index 2 of saturation, mole fraction or temperature
+
+ // equation indices
+ static const int conti0EqIdx = PVOffset + wCompIdx; //!< Index of the mass conservation equation for the water component
+ static const int conti1EqIdx = conti0EqIdx + nCompIdx; //!< Index of the mass conservation equation for the contaminant component
+ static const int energyEqIdx = PVOffset + 2; //! The index for energy in equation vectors.
+
+ static const int contiWEqIdx = conti0EqIdx + wCompIdx; //!< index of the mass conservation equation for the water component
+ static const int contiNEqIdx = conti0EqIdx + nCompIdx; //!< index of the mass conservation equation for the contaminant component
+};
+
+}
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/localresidual.hh b/dumux/porousmediumflow/3pwateroil/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..060b90ce42b72ffaa11e4a1b6b815e768e63cdf4
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/localresidual.hh
@@ -0,0 +1,367 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the three-phase three-component fully implicit model.
+ */
+#ifndef DUMUX_3P2CNI_LOCAL_RESIDUAL_HH
+#define DUMUX_3P2CNI_LOCAL_RESIDUAL_HH
+
+#include "properties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the three-phase three-component fully implicit model.
+ *
+ * This class is used to fill the gaps in BoxLocalResidual for the 3P2C flow.
+ */
+template<class TypeTag>
+class ThreePWaterOilLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
+{
+protected:
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+
+
+ enum {
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+
+ conti0EqIdx = Indices::conti0EqIdx,//!< Index of the mass conservation equation for the water component
+ conti1EqIdx = Indices::conti1EqIdx,//!< Index of the mass conservation equation for the contaminant component
+ energyEqIdx = Indices::energyEqIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ //! property that defines whether mole or mass fractions are used
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+public:
+
+ /*!
+ * \brief Evaluate the storage term of the current solution in a
+ * single phase.
+ *
+ * \param element The element
+ * \param phaseIdx The index of the fluid phase
+ */
+ void evalPhaseStorage(const Element &element, const int phaseIdx)
+ {
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(this->gridView_(), element);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeometry);
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(), element, fvGeometry, false);
+
+ this->storageTerm_.resize(fvGeometry.numScv);
+ this->storageTerm_ = 0;
+
+ this->elemPtr_ = &element;
+ this->fvElemGeomPtr_ = &fvGeometry;
+ this->bcTypesPtr_ = &bcTypes;
+ this->prevVolVarsPtr_ = 0;
+ this->curVolVarsPtr_ = &elemVolVars;
+ evalPhaseStorage_(phaseIdx);
+ }
+
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \param storage The mass of the component within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \param usePrevSol Evaluate function with solution of current or previous time step
+ */
+ void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const
+ {
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemVolVars =
+ usePrevSol
+ ? this->prevVolVars_()
+ : this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
+
+ // compute storage term of all components within all phases
+ storage = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ storage[conti0EqIdx + compIdx] +=
+ volVars.porosity()
+ * volVars.saturation(phaseIdx)
+ * volVars.molarDensity(phaseIdx)
+ * volVars.fluidState().moleFraction(phaseIdx, compIdx);
+ }
+ }
+
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+
+ // compute the energy storage
+ storage[energyEqIdx] = volVars.porosity()
+ *(
+ volVars.density(wPhaseIdx)
+ *volVars.internalEnergy(wPhaseIdx)
+ *volVars.saturation(wPhaseIdx)
+ +
+ volVars.density(nPhaseIdx)
+ *volVars.internalEnergy(nPhaseIdx)
+ *volVars.saturation(nPhaseIdx)
+ +
+ volVars.density(gPhaseIdx)
+ *volVars.internalEnergy(gPhaseIdx)
+ *volVars.saturation(gPhaseIdx)
+ )
+ + volVars.temperature()*volVars.heatCapacity();
+ }
+
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param faceIdx The index of the SCV face
+ * \param onBoundary A boolean variable to specify whether the flux variables
+ * are calculated for interior SCV faces or boundary faces, default=false
+ */
+ void computeFlux(PrimaryVariables &flux, const int fIdx, const bool onBoundary=false) const
+ {
+ FluxVariables fluxVars;
+ fluxVars.update(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ fIdx,
+ this->curVolVars_(),
+ onBoundary);
+
+ flux = 0;
+ asImp_()->computeAdvectiveFlux(flux, fluxVars);
+ asImp_()->computeDiffusiveFlux(flux, fluxVars);
+ }
+
+ /*!
+ * \brief Evaluates the advective mass flux of all components over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The advective flux over the sub-control-volume face for each component
+ * \param fluxVars The flux variables at the current SCV
+ */
+
+ void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ {
+ Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+
+ ////////
+ // advective fluxes of all components in all phases
+ ////////
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ // add advective flux of current component in current
+ // phase
+ // if alpha > 0 und alpha < 1 then both upstream and downstream
+ // nodes need their contribution
+ // if alpha == 1 (which is mostly the case) then, the downstream
+ // node is not evaluated
+ int eqIdx = conti0EqIdx + compIdx;
+ flux[eqIdx] += fluxVars.volumeFlux(phaseIdx)
+ * (massUpwindWeight
+ * up.fluidState().molarDensity(phaseIdx)
+ * up.fluidState().moleFraction(phaseIdx, compIdx)
+ +
+ (1.0 - massUpwindWeight)
+ * dn.fluidState().molarDensity(phaseIdx)
+ * dn.fluidState().moleFraction(phaseIdx, compIdx));
+ }
+ }
+
+ // advective heat flux in all phases
+ flux[energyEqIdx] = 0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // vertex data of the upstream and the downstream vertices
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
+
+ flux[energyEqIdx] +=
+ fluxVars.volumeFlux(phaseIdx) * (
+ massUpwindWeight * // upstream vertex
+ ( up.density(phaseIdx) *
+ up.enthalpy(phaseIdx))
+ +
+ (1-massUpwindWeight) * // downstream vertex
+ ( dn.density(phaseIdx) *
+ dn.enthalpy(phaseIdx)) );
+ }
+
+ }
+
+ /*!
+ * \brief Adds the diffusive mass flux of all components over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The diffusive flux over the sub-control-volume face for each component
+ * \param fluxVars The flux variables at the current SCV
+ */
+
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ {
+ // TODO: reference!? Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = fluxVars.porousDiffCoeff();
+ // add diffusive flux of gas component in liquid phase
+ Scalar tmp;
+ tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx] * fluxVars.molarDensity(wPhaseIdx);
+ tmp *= (fluxVars.moleFractionCompNGrad(wPhaseIdx) * fluxVars.face().normal);
+ Scalar jNW = tmp;
+
+ Scalar jWW = -jNW;
+
+ tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx] * fluxVars.molarDensity(gPhaseIdx);
+ tmp *= (fluxVars.moleFractionCompWGrad(gPhaseIdx) * fluxVars.face().normal);
+ Scalar jWG = tmp;
+
+ Scalar jNG = -jWG;
+
+ tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx] * fluxVars.molarDensity(nPhaseIdx);
+ tmp *= (fluxVars.moleFractionCompWGrad(nPhaseIdx) * fluxVars.face().normal);
+ Scalar jWN = tmp;
+
+ Scalar jNN = -jWN;
+
+ flux[conti0EqIdx] += jWW+jWG+jWN;
+ flux[conti1EqIdx] += jNW+jNG+jNN;
+ }
+
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param source The source/sink in the SCV for each component
+ * \param scvIdx The index of the SCV
+ */
+ void computeSource(PrimaryVariables &source, const int scvIdx)
+ {
+ this->problem_().solDependentSource(source,
+ this->element_(),
+ this->fvGeometry_(),
+ scvIdx,
+ this->curVolVars_());
+ }
+
+protected:
+ void evalPhaseStorage_(const int phaseIdx)
+ {
+ if(!useMoles) //mass-fraction formulation
+ {
+ // evaluate the storage terms of a single phase
+ for (int i=0; i < this->fvGeometry_().numScv; i++) {
+ PrimaryVariables &storage = this->storageTerm_[i];
+ const ElementVolumeVariables &elemVolVars = this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[i];
+
+ // compute storage term of all components within all phases
+ storage = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == wCompIdx) ? conti0EqIdx : conti1EqIdx;
+ storage[eqIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFraction(phaseIdx, compIdx);
+ }
+
+ storage *= volVars.porosity();
+ storage *= this->fvGeometry_().subContVol[i].volume;
+ }
+ }
+ else //mole-fraction formulation
+ {
+ // evaluate the storage terms of a single phase
+ for (int i=0; i < this->fvGeometry_().numScv; i++) {
+ PrimaryVariables &storage = this->storageTerm_[i];
+ const ElementVolumeVariables &elemVolVars = this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[i];
+
+ // compute storage term of all components within all phases
+ storage = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == wCompIdx) ? conti0EqIdx : conti1EqIdx;
+ storage[eqIdx] += volVars.molarDensity(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().moleFraction(phaseIdx, compIdx);
+ }
+
+ storage *= volVars.porosity();
+ storage *= this->fvGeometry_().subContVol[i].volume;
+ }
+ }
+ }
+ Implementation *asImp_()
+ {
+ return static_cast<Implementation *> (this);
+ }
+
+ const Implementation *asImp_() const
+ {
+ return static_cast<const Implementation *> (this);
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/model.hh b/dumux/porousmediumflow/3pwateroil/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..11e374b31b6f04e5ad917d67848cc1bc1fcff119
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/model.hh
@@ -0,0 +1,1090 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Adaption of the fully implicit scheme to the three-phase three-component flow model.
+ *
+ * The model is designed for simulating three fluid phases with water, gas, and
+ * a liquid contaminant (NAPL - non-aqueous phase liquid)
+ */
+#ifndef DUMUX_3PWATEROIL_MODEL_HH
+#define DUMUX_3PWATEROIL_MODEL_HH
+
+#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+
+#include "properties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \brief Adaption of the fully implicit scheme to the three-phase two-component flow model.
+ *
+ * This model implements three-phase two-component flow of three fluid phases
+ * \f$\alpha \in \{ water, gas, NAPL \}\f$ each composed of up to two components
+ * \f$\kappa \in \{ water, contaminant \}\f$. The standard multiphase Darcy
+ * approach is used as the equation for the conservation of momentum:
+ * \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ * \f]
+ *
+ * By inserting this into the equations for the conservation of the
+ * components, one transport equation for each component is obtained as
+ * \f{eqnarray*}
+ && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa
+ S_\alpha )}{\partial t}
+ - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha}
+ \varrho_\alpha x_\alpha^\kappa \mbox{\bf K}
+ (\textbf{grad}\, p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\}
+ \nonumber \\
+ \nonumber \\
+ && - \sum\limits_\alpha \text{div} \left\{ D_\text{pm}^\kappa \varrho_\alpha \frac{M^\kappa}{M_\alpha}
+ \textbf{grad} x^\kappa_{\alpha} \right\}
+ - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha
+ \f}
+ *
+ * Note that these balance equations are molar.
+ *
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme as spatial
+ * and the implicit Euler method as time discretization.
+ *
+ * The model uses commonly applied auxiliary conditions like
+ * \f$S_w + S_n + S_g = 1\f$ for the saturations and
+ * \f$x^w_\alpha + x^c_\alpha = 1\f$ for the mole fractions.
+ * Furthermore, the phase pressures are related to each other via
+ * capillary pressures between the fluid phases, which are functions of
+ * the saturation, e.g. according to the approach of Parker et al.
+ *
+ * The used primary variables are dependent on the locally present fluid phases
+ * An adaptive primary variable switch is included. The phase state is stored for all nodes
+ * of the system. Different cases can be distinguished:
+ * <ul>
+ * <li> ... to be completed ... </li>
+ * </ul>
+ */
+template<class TypeTag>
+class ThreePWaterOilModel: public GET_PROP_TYPE(TypeTag, BaseModel)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, BaseModel) ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ 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;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+
+ pressureIdx = Indices::pressureIdx,
+ switch1Idx = Indices::switch1Idx,
+ switch2Idx = Indices::switch2Idx,
+
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx,
+
+ threePhases = Indices::threePhases,
+ wPhaseOnly = Indices::wPhaseOnly,
+ gnPhaseOnly = Indices::gnPhaseOnly,
+ wnPhaseOnly = Indices::wnPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ wgPhaseOnly = Indices::wgPhaseOnly
+
+ };
+
+
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+public:
+ /*!
+ * \brief Initialize the static data with the initial solution.
+ *
+ * \param problem The problem to be solved
+ */
+ void init(Problem &problem)
+ {
+ ParentType::init(problem);
+
+ staticDat_.resize(this->numDofs());
+
+ setSwitched_(false);
+
+ if (isBox)
+ {
+ VertexIterator vIt = this->gridView_().template begin<dim> ();
+ const VertexIterator &vEndIt = this->gridView_().template end<dim> ();
+ for (; vIt != vEndIt; ++vIt)
+ {
+ int globalIdx = this->dofMapper().index(*vIt);
+ const GlobalPosition &globalPos = vIt->geometry().corner(0);
+
+ // initialize phase presence
+ staticDat_[globalIdx].phasePresence
+ = this->problem_().initialPhasePresence(*vIt, globalIdx,
+ globalPos);
+ staticDat_[globalIdx].wasSwitched = false;
+
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
+ }
+ }
+ else
+ {
+ ElementIterator eIt = this->gridView_().template begin<0> ();
+ const ElementIterator &eEndIt = this->gridView_().template end<0> ();
+ for (; eIt != eEndIt; ++eIt)
+ {
+ int globalIdx = this->dofMapper().index(*eIt);
+ const GlobalPosition &globalPos = eIt->geometry().center();
+
+ // initialize phase presence
+ staticDat_[globalIdx].phasePresence
+ = this->problem_().initialPhasePresence(*this->gridView_().template begin<dim> (),
+ globalIdx, globalPos);
+ staticDat_[globalIdx].wasSwitched = false;
+
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
+ }
+ }
+ }
+
+ /*!
+ * \brief Compute the total storage inside one phase of all
+ * conservation quantities.
+ *
+ * \param storage Contains the storage of each component for one phase
+ * \param phaseIdx The phase index
+ */
+ void globalPhaseStorage(PrimaryVariables &storage, const int phaseIdx)
+ {
+ storage = 0;
+
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ const ElementIterator eEndIt = this->gridView_().template end<0>();
+ for (; eIt != eEndIt; ++eIt) {
+ if(eIt->partitionType() == Dune::InteriorEntity)
+ {
+ this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
+
+ for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+ storage += this->localResidual().storageTerm()[i];
+ }
+ }
+ if (this->gridView_().comm().size() > 1)
+ storage = this->gridView_().comm().sum(storage);
+ }
+
+
+ /*!
+ * \brief Called by the update() method if applying the newton
+ * method was unsuccessful.
+ */
+ void updateFailed()
+ {
+ ParentType::updateFailed();
+
+ setSwitched_(false);
+ resetPhasePresence_();
+ };
+
+ /*!
+ * \brief Called by the problem if a time integration was
+ * successful, post processing of the solution is done and the
+ * result has been written to disk.
+ *
+ * This should prepare the model for the next time integration.
+ */
+ void advanceTimeLevel()
+ {
+ ParentType::advanceTimeLevel();
+
+ // update the phase state
+ updateOldPhasePresence_();
+ setSwitched_(false);
+ }
+
+ /*!
+ * \brief Return true if the primary variables were switched for
+ * at least one vertex after the last timestep.
+ */
+ bool switched() const
+ {
+ return switchFlag_;
+ }
+
+ /*!
+ * \brief Returns the phase presence of the current or the old solution of a degree of freedom.
+ *
+ * \param globalIdx The global index of the degree of freedom
+ * \param oldSol Evaluate function with solution of current or previous time step
+ */
+ int phasePresence(int globalIdx, bool oldSol) const
+ {
+ return
+ oldSol
+ ? staticDat_[globalIdx].oldPhasePresence
+ : staticDat_[globalIdx].phasePresence;
+ }
+
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ *
+ * \param sol The solution vector
+ * \param writer The writer for multi-file VTK datasets
+ */
+ template<class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol,
+ 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();
+
+ // create the required scalar fields
+ ScalarField *saturation[numPhases];
+ ScalarField *pressure[numPhases];
+ ScalarField *density[numPhases];
+
+ ScalarField *mobility[numPhases];
+ ScalarField *viscosity[numPhases];
+
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+
+ mobility[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ viscosity[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ }
+
+ ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
+ ScalarField *moleFraction[numPhases][numComponents];
+ for (int i = 0; i < numPhases; ++i)
+ for (int j = 0; j < numComponents; ++j)
+ moleFraction[i][j] = writer.allocateManagedBuffer (numDofs);
+ 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_());
+
+ 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);
+ }
+ }
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField *rank = writer.allocateManagedBuffer (numElements);
+
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ ElementIterator eEndIt = this->gridView_().template end<0>();
+ for (; eIt != eEndIt; ++eIt)
+ {
+ int idx = this->problem_().elementMapper().index(*eIt);
+ (*rank)[idx] = this->gridView_().comm().rank();
+
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(this->gridView_(), *eIt);
+
+
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(),
+ *eIt,
+ fvGeometry,
+ false /* oldSol? */);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ int globalIdx = this->dofMapper().subIndex(*eIt, scvIdx, dofCodim);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ 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);
+
+ (*mobility[phaseIdx])[globalIdx] = elemVolVars[scvIdx].mobility(phaseIdx);
+ (*viscosity[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().viscosity(phaseIdx);
+ }
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+ (*moleFraction[phaseIdx][compIdx])[globalIdx] =
+ elemVolVars[scvIdx].fluidState().moleFraction(phaseIdx,
+ compIdx);
+
+ Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
+ }
+ }
+
+ (*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, *eIt, wPhaseIdx);
+ velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
+ velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
+
+ }
+
+ writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
+ writer.attachDofData(*saturation[nPhaseIdx], "sn", isBox);
+ writer.attachDofData(*saturation[gPhaseIdx], "sg", isBox);
+ writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
+ writer.attachDofData(*pressure[nPhaseIdx], "pn", isBox);
+ writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
+ writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
+ writer.attachDofData(*density[nPhaseIdx], "rhon", isBox);
+ writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
+
+ writer.attachDofData(*mobility[wPhaseIdx], "MobW", isBox);
+ writer.attachDofData(*mobility[nPhaseIdx], "MobN", isBox);
+ writer.attachDofData(*mobility[gPhaseIdx], "MobG", isBox);
+
+ writer.attachDofData(*viscosity[wPhaseIdx], "ViscosW", isBox);
+ writer.attachDofData(*viscosity[nPhaseIdx], "ViscosN", isBox);
+ writer.attachDofData(*viscosity[gPhaseIdx], "ViscosG", isBox);
+ for (int i = 0; i < numPhases; ++i)
+ {
+ for (int j = 0; j < numComponents; ++j)
+ {
+ std::ostringstream oss;
+ oss << "x^"
+ << FluidSystem::phaseName(i)
+ << "_"
+ << FluidSystem::componentName(j);
+ writer.attachDofData(*moleFraction[i][j], oss.str().c_str(), isBox);
+ }
+ }
+ writer.attachDofData(*poro, "porosity", isBox);
+ 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");
+ }
+
+ /*!
+ * \brief Write the current solution to a restart file.
+ *
+ * \param outStream The output stream of one entity for the restart file
+ * \param entity The entity, either a vertex or an element
+ */
+ template<class Entity>
+ void serializeEntity(std::ostream &outStream, const Entity &entity)
+ {
+ // write primary variables
+ ParentType::serializeEntity(outStream, entity);
+
+ int globalIdx = this->dofMapper().index(entity);
+ if (!outStream.good())
+ DUNE_THROW(Dune::IOError, "Could not serialize entity " << globalIdx);
+
+ outStream << staticDat_[globalIdx].phasePresence << " ";
+ }
+
+ /*!
+ * \brief Reads the current solution from a restart file.
+ *
+ * \param inStream The input stream of one entity from the restart file
+ * \param entity The entity, either a vertex or an element
+ */
+ template<class Entity>
+ void deserializeEntity(std::istream &inStream, const Entity &entity)
+ {
+ // read primary variables
+ ParentType::deserializeEntity(inStream, entity);
+
+ // read phase presence
+ int globalIdx = this->dofMapper().index(entity);
+ if (!inStream.good())
+ DUNE_THROW(Dune::IOError,
+ "Could not deserialize entity " << globalIdx);
+
+ inStream >> staticDat_[globalIdx].phasePresence;
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
+
+ }
+
+ /*!
+ * \brief Update the static data of all vertices in the grid.
+ *
+ * \param curGlobalSol The current global solution
+ * \param oldGlobalSol The previous global solution
+ */
+ void updateStaticData(SolutionVector &curGlobalSol,
+ const SolutionVector &oldGlobalSol)
+ {
+ bool wasSwitched = false;
+
+ bool useSimpleModel = GET_PROP_VALUE(TypeTag, UseSimpleModel);
+
+ for (unsigned i = 0; i < staticDat_.size(); ++i)
+ staticDat_[i].visited = false;
+
+ FVElementGeometry fvGeometry;
+ static VolumeVariables volVars;
+ ElementIterator eIt = this->gridView_().template begin<0> ();
+ const ElementIterator &eEndIt = this->gridView_().template end<0> ();
+ for (; eIt != eEndIt; ++eIt)
+ {
+ fvGeometry.update(this->gridView_(), *eIt);
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ int globalIdx = this->dofMapper().subIndex(*eIt, scvIdx, dofCodim);
+
+ if (staticDat_[globalIdx].visited)
+ continue;
+
+ staticDat_[globalIdx].visited = true;
+ volVars.update(curGlobalSol[globalIdx],
+ this->problem_(),
+ *eIt,
+ fvGeometry,
+ scvIdx,
+ false);
+ const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+
+ if(useSimpleModel)
+ {
+ if (primaryVarSwitchSimple_(curGlobalSol,
+ volVars,
+ globalIdx,
+ globalPos))
+ {
+ this->jacobianAssembler().markDofRed(globalIdx);
+ wasSwitched = true;
+ }
+ }
+ else
+ {
+ if (primaryVarSwitch_(curGlobalSol,
+ volVars,
+ globalIdx,
+ globalPos))
+ {
+ this->jacobianAssembler().markDofRed(globalIdx);
+ wasSwitched = true;
+ }
+ }
+ }
+ }
+
+ // make sure that if there was a variable switch in an
+ // other partition we will also set the switch flag
+ // for our partition.
+ if (this->gridView_().comm().size() > 1)
+ wasSwitched = this->gridView_().comm().max(wasSwitched);
+
+ setSwitched_(wasSwitched);
+ }
+
+protected:
+ /*!
+ * \brief Data which is attached to each vertex and is not only
+ * stored locally.
+ */
+ struct StaticVars
+ {
+ int phasePresence;
+ bool wasSwitched;
+
+ int oldPhasePresence;
+ bool visited;
+ };
+
+ /*!
+ * \brief Reset the current phase presence of all vertices to the old one.
+ *
+ * This is done after an update failed.
+ */
+ void resetPhasePresence_()
+ {
+ for (unsigned int i = 0; i < this->numDofs(); ++i)
+ {
+ staticDat_[i].phasePresence
+ = staticDat_[i].oldPhasePresence;
+ staticDat_[i].wasSwitched = false;
+ }
+ }
+
+ /*!
+ * \brief Set the old phase of all verts state to the current one.
+ */
+ void updateOldPhasePresence_()
+ {
+ for (unsigned int i = 0; i < this->numDofs(); ++i)
+ {
+ staticDat_[i].oldPhasePresence
+ = staticDat_[i].phasePresence;
+ staticDat_[i].wasSwitched = false;
+ }
+ }
+
+ /*!
+ * \brief Set whether there was a primary variable switch after in
+ * the last timestep.
+ */
+ void setSwitched_(bool yesno)
+ {
+ switchFlag_ = yesno;
+ }
+
+ // perform variable switch at a vertex; Returns true if a
+ // variable switch was performed.
+ bool primaryVarSwitch_(SolutionVector &globalSol,
+ const VolumeVariables &volVars,
+ int globalIdx,
+ const GlobalPosition &globalPos)
+ {
+ // evaluate primary variable switch
+ bool wouldSwitch = false;
+ int phasePresence = staticDat_[globalIdx].phasePresence;
+ int newPhasePresence = phasePresence;
+
+ // check if a primary var switch is necessary
+ if (phasePresence == threePhases)
+ {
+ Scalar Smin = 0;
+ if (staticDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ newPhasePresence = wnPhaseOnly;
+
+ globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
+ }
+ else if (volVars.saturation(wPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // water phase disappears
+ std::cout << "Water phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sw: "
+ << volVars.saturation(wPhaseIdx) << std::endl;
+ newPhasePresence = gnPhaseOnly;
+
+ globalSol[globalIdx][switch1Idx] = volVars.fluidState().saturation(nPhaseIdx);
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
+ }
+ else if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // NAPL phase disappears
+ std::cout << "NAPL phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sn: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ newPhasePresence = wgPhaseOnly;
+
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == wPhaseOnly)
+ {
+ bool gasFlag = 0;
+ bool nonwettingFlag = 0;
+ // calculate fractions of the partial pressures in the
+ // hypothetical gas phase
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xwg + xng > xgMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xwg + xng > xgMax)
+ {
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xwg + xng: "
+ << xwg + xng << std::endl;
+ gasFlag = 1;
+ }
+
+ // calculate fractions in the hypothetical NAPL phase
+ Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+ /* take care:
+ for xnn in case wPhaseOnly we compute xnn=henry_mesitylene*x1w,
+ where a hypothetical gas pressure is assumed for the Henry
+ x0n is set to NULL (all NAPL phase is dirty)
+ x2n is set to NULL (all NAPL phase is dirty)
+ */
+
+ Scalar xnMax = 1.0;
+ if (xnn > xnMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xnMax *= 1.02;
+
+ // if the sum of the hypothetical mole fractions would be larger than
+ // 100%, NAPL phase appears
+ if (xnn > xnMax)
+ {
+ // NAPL phase appears
+ std::cout << "NAPL phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xnn: "
+ << xnn << std::endl;
+ nonwettingFlag = 1;
+ }
+
+ if ((gasFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = wgPhaseOnly;
+ globalSol[globalIdx][switch1Idx] = 0.9999;
+ globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
+ }
+ else if ((gasFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = threePhases;
+ globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = 0.999;
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wnPhaseOnly;
+ globalSol[globalIdx][switch2Idx] = 0.0001;
+ }
+ }
+ else if (phasePresence == gnPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool wettingFlag = 0;
+
+ Scalar Smin = 0.0;
+ if (staticDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // NAPL phase disappears
+ std::cout << "NAPL phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sn: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ nonwettingFlag = 1;
+ }
+
+
+ // calculate fractions of the hypothetical water phase
+ Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ /*
+ take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
+ If this is larger than 1, then water appears
+ */
+ Scalar xwMax = 1.0;
+ if (xww > xwMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xwMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, water phase appears
+ if (xww > xwMax)
+ {
+ // water phase appears
+ std::cout << "water phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
+ << xww << std::endl;
+ wettingFlag = 1;
+ }
+
+ if ((wettingFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = threePhases;
+ globalSol[globalIdx][switch1Idx] = 0.0001;
+ globalSol[globalIdx][switch2Idx] = volVars.saturation(nPhaseIdx);
+ }
+ else if ((wettingFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wgPhaseOnly;
+ globalSol[globalIdx][switch1Idx] = 0.0001;
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((wettingFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = gPhaseOnly;
+ globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool gasFlag = 0;
+
+ Scalar Smin = 0.0;
+ if (staticDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // NAPL phase disappears
+ std::cout << "NAPL phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sn: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ nonwettingFlag = 1;
+ }
+
+ // calculate fractions of the partial pressures in the
+ // hypothetical gas phase
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xwg + xng > xgMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xwg + xng > xgMax)
+ {
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xwg + xng: "
+ << xwg + xng << std::endl;
+ gasFlag = 1;
+ }
+
+ if ((gasFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = threePhases;
+ globalSol[globalIdx][pressureIdx] = volVars.pressure(gPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
+ }
+ else if ((gasFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wgPhaseOnly;
+ globalSol[globalIdx][pressureIdx] = volVars.pressure(gPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = volVars.temperature();
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wPhaseOnly;
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == gPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool wettingFlag = 0;
+
+ // calculate fractions in the hypothetical NAPL phase
+ Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+ /*
+ take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
+ if this is larger than 1, then NAPL appears
+ */
+
+ Scalar xnMax = 1.0;
+ if (xnn > xnMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xnMax *= 1.02;
+
+ // if the sum of the hypothetical mole fraction would be larger than
+ // 100%, NAPL phase appears
+ if (xnn > xnMax)
+ {
+ // NAPL phase appears
+ std::cout << "NAPL phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xnn: "
+ << xnn << std::endl;
+ nonwettingFlag = 1;
+ }
+ // calculate fractions of the hypothetical water phase
+ Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ /*
+ take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
+ If this is larger than 1, then water appears
+ */
+ Scalar xwMax = 1.0;
+ if (xww > xwMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xwMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xww > xwMax)
+ {
+ // water phase appears
+ std::cout << "water phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
+ << xww << std::endl;
+ wettingFlag = 1;
+ }
+ if ((wettingFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = wgPhaseOnly;
+ globalSol[globalIdx][switch1Idx] = 0.0001;
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((wettingFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = threePhases;
+ globalSol[globalIdx][switch1Idx] = 0.0001;
+ globalSol[globalIdx][switch2Idx] = 0.0001;
+ }
+ else if ((wettingFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = gnPhaseOnly;
+ globalSol[globalIdx][switch2Idx]
+ = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == wgPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool gasFlag = 0;
+ bool wettingFlag = 0;
+
+ // get the fractions in the hypothetical NAPL phase
+ Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+
+ // take care: if the NAPL phase is not present, take
+ // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
+ // appears
+ Scalar xnMax = 1.0;
+ if (xnn > xnMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xnMax *= 1.02;
+
+ // if the sum of the hypothetical mole fraction would be larger than
+ // 100%, NAPL phase appears
+ if (xnn > xnMax)
+ {
+ // NAPL phase appears
+ std::cout << "NAPL phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xnn: "
+ << xnn << std::endl;
+ nonwettingFlag = 1;
+ }
+
+ Scalar Smin = -1.e-6;
+ if (staticDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ gasFlag = 1;
+ }
+
+ Smin = 0.0;
+ if (staticDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(wPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Water phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sw: "
+ << volVars.saturation(wPhaseIdx) << std::endl;
+ wettingFlag = 1;
+ }
+
+ if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 1))
+ {
+ newPhasePresence = gnPhaseOnly;
+ globalSol[globalIdx][switch1Idx] = 0.0001;
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
+;
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 0))
+ {
+ newPhasePresence = threePhases;
+ globalSol[globalIdx][switch2Idx] = 0.0001;
+ }
+ else if ((gasFlag == 1) && (nonwettingFlag == 0) && (wettingFlag == 0))
+ {
+ newPhasePresence = wPhaseOnly;
+ globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
+ globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 0) && (wettingFlag == 1))
+ {
+ newPhasePresence = gPhaseOnly;
+ globalSol[globalIdx][switch1Idx]
+ = volVars.fluidState().temperature();
+ globalSol[globalIdx][switch2Idx]
+ = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+ }
+ }
+
+ staticDat_[globalIdx].phasePresence = newPhasePresence;
+ staticDat_[globalIdx].wasSwitched = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+
+ // perform variable switch at a vertex; Returns true if a
+ // variable switch was performed.
+ // Note that this one is the simplified version with only two phase states
+ bool primaryVarSwitchSimple_(SolutionVector &globalSol,
+ const VolumeVariables &volVars,
+ int globalIdx,
+ const GlobalPosition &globalPos)
+ {
+ // evaluate primary variable switch
+ bool wouldSwitch = false;
+ int phasePresence = staticDat_[globalIdx].phasePresence;
+ int newPhasePresence = phasePresence;
+
+ // check if a primary var switch is necessary
+ if (phasePresence == threePhases)
+ {
+ Scalar Smin = 0;
+ if (staticDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ newPhasePresence = wnPhaseOnly;
+
+ globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
+ }
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ bool gasFlag = 0;
+
+ // calculate fractions of the partial pressures in the
+ // hypothetical gas phase
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xwg + xng > xgMax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xwg + xng > xgMax)
+ {
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xwg + xng: "
+ << xwg + xng << std::endl;
+ gasFlag = 1;
+ }
+
+ if (gasFlag == 1)
+ {
+ newPhasePresence = threePhases;
+ globalSol[globalIdx][pressureIdx] = volVars.pressure(gPhaseIdx);
+ globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
+ }
+ }
+
+ staticDat_[globalIdx].phasePresence = newPhasePresence;
+ staticDat_[globalIdx].wasSwitched = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+
+ // parameters given in constructor
+ std::vector<StaticVars> staticDat_;
+ bool switchFlag_;
+};
+
+}
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/newtoncontroller.hh b/dumux/porousmediumflow/3pwateroil/newtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..033fc13a275c8eabc714a47c6422da1e0083e7ba
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/newtoncontroller.hh
@@ -0,0 +1,85 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A 3p2cni specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+#ifndef DUMUX_3P2CNI_NEWTON_CONTROLLER_HH
+#define DUMUX_3P2CNI_NEWTON_CONTROLLER_HH
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include "properties.hh"
+
+namespace Dumux {
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \brief A 3p2cni specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class ThreePWaterOilNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+public:
+ ThreePWaterOilNewtonController(const Problem &problem)
+ : ParentType(problem)
+ {};
+
+
+ /*!
+ * \brief Called after each Newton update
+ *
+ * \param uCurrentIter The current global solution vector
+ * \param uLastIter The previous global solution vector
+ */
+ void newtonEndStep(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter)
+ {
+ // call the method of the base class
+ this->method().model().updateStaticData(uCurrentIter, uLastIter);
+ ParentType::newtonEndStep(uCurrentIter, uLastIter);
+ }
+
+
+ /*!
+ * \brief Returns true if the current solution can be considered to
+ * be accurate enough
+ */
+ bool newtonConverged()
+ {
+ if (this->method().model().switched())
+ return false;
+
+ return ParentType::newtonConverged();
+ };
+};
+}
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/properties.hh b/dumux/porousmediumflow/3pwateroil/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7c42af3d694e188e2b03c7625a3bc77cc5274211
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/properties.hh
@@ -0,0 +1,73 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup ThreePWaterOilModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the 3p2cni fully implicit model.
+ */
+#ifndef DUMUX_3P2CNI_PROPERTIES_HH
+#define DUMUX_3P2CNI_PROPERTIES_HH
+
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/implicit/cellcentered/properties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit three-phase three-component problems
+NEW_TYPE_TAG(ThreePWaterOil);
+NEW_TYPE_TAG(BoxThreePWaterOil, INHERITS_FROM(BoxModel, ThreePWaterOil));
+NEW_TYPE_TAG(CCThreePWaterOil, INHERITS_FROM(CCModel, ThreePWaterOil));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
+NEW_PROP_TAG(Indices); //!< Enumerations for the model
+NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
+NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
+
+NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
+
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+
+NEW_PROP_TAG(UseMoles); //!Defines whether mole (true) or mass (false) fractions are used
+
+NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind parameter for the mobility
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
+NEW_PROP_TAG(UseSimpleModel); //!< Determines whether a simple model with only two phase states (wng) and (wn) should be used
+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
+}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/propertydefaults.hh b/dumux/porousmediumflow/3pwateroil/propertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..48f0a337e73bc08a4ddcb686687e12574d08bc58
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/propertydefaults.hh
@@ -0,0 +1,148 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup ThreePWaterOilModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ * 3p2cni fully implicit model.
+ */
+#ifndef DUMUX_3P2CNI_PROPERTY_DEFAULTS_HH
+#define DUMUX_3P2CNI_PROPERTY_DEFAULTS_HH
+
+#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
+#include <dumux/material/spatialparams/implicit.hh>
+
+#include "indices.hh"
+#include "model.hh"
+#include "fluxvariables.hh"
+#include "volumevariables.hh"
+#include "properties.hh"
+#include "newtoncontroller.hh"
+#include "localresidual.hh"
+
+namespace Dumux
+{
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \brief Set the property for the number of components.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(ThreePWaterOil, NumComponents)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ public:
+ static const int value = FluidSystem::numComponents;
+
+ static_assert(value == 2,
+ "Only fluid systems with 2 components are supported by the 3p2cni model!");
+};
+
+/*!
+ * \brief Set the property for the number of fluid phases.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(ThreePWaterOil, NumPhases)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 3,
+ "Only fluid systems with 3 phases are supported by the 3p2cni model!");
+};
+
+SET_INT_PROP(ThreePWaterOil, NumEq, 3); //!< set the number of equations to 2
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_TYPE_PROP(ThreePWaterOil, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+//! The local residual function of the conservation equations
+SET_TYPE_PROP(ThreePWaterOil, LocalResidual, ThreePWaterOilLocalResidual<TypeTag>);
+
+//! Use the 3p2cni specific newton controller for the 3p2cni model
+SET_TYPE_PROP(ThreePWaterOil, NewtonController, ThreePWaterOilNewtonController<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(ThreePWaterOil, Model, ThreePWaterOilModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(ThreePWaterOil, VolumeVariables, ThreePWaterOilVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(ThreePWaterOil, FluxVariables, ThreePWaterOilFluxVariables<TypeTag>);
+
+//! define the base flux variables to realize Darcy flow
+SET_TYPE_PROP(ThreePWaterOil, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
+
+//! the upwind factor for the mobility.
+SET_SCALAR_PROP(ThreePWaterOil, ImplicitMassUpwindWeight, 1.0);
+
+//! set default mobility upwind weight to 1.0, i.e. fully upwind
+SET_SCALAR_PROP(ThreePWaterOil, ImplicitMobilityUpwindWeight, 1.0);
+
+//! Determines whether a constraint solver should be used explicitly
+SET_BOOL_PROP(ThreePWaterOil, UseSimpleModel, true);
+
+//! The indices required by the isothermal 3p2cni model
+SET_TYPE_PROP(ThreePWaterOil, Indices, ThreePWaterOilIndices<TypeTag, /*PVOffset=*/0>);
+
+//! The spatial parameters to be employed.
+//! Use ImplicitSpatialParams by default.
+SET_TYPE_PROP(ThreePWaterOil, SpatialParams, ImplicitSpatialParams<TypeTag>);
+
+// disable velocity output by default
+SET_BOOL_PROP(ThreePWaterOil, VtkAddVelocity, false);
+
+// enable gravity by default
+SET_BOOL_PROP(ThreePWaterOil, ProblemEnableGravity, true);
+
+SET_BOOL_PROP(ThreePWaterOil, UseMoles, true); //!< Define that mole fractions are used in the balance equations per default
+
+
+
+//! default value for the forchheimer coefficient
+// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
+// Actually the Forchheimer coefficient is also a function of the dimensions of the
+// porous medium. Taking it as a constant is only a first approximation
+// (Nield, Bejan, Convection in porous media, 2006, p. 10)
+SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
+
+}
+
+}
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/volumevariables.hh b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..016924f699244f573cda1ab6d8197dbde327f498
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
@@ -0,0 +1,876 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the three-phase, two-component model.
+ */
+#ifndef DUMUX_3P2CNI_VOLUME_VARIABLES_HH
+#define DUMUX_3P2CNI_VOLUME_VARIABLES_HH
+
+#include <dumux/implicit/model.hh>
+#include <dumux/common/math.hh>
+
+#include <dune/common/parallel/collectivecommunication.hh>
+#include <vector>
+#include <iostream>
+
+#include <dumux/material/constants.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
+#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+
+#include "properties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the two-phase, two-component model.
+ */
+template <class TypeTag>
+class ThreePWaterOilVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+ typedef ImplicitVolumeVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+
+ // constraint solvers
+ typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
+ typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ dim = GridView::dimension,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ switch1Idx = Indices::switch1Idx,
+ switch2Idx = Indices::switch2Idx,
+ pressureIdx = Indices::pressureIdx
+ };
+
+ // present phases
+ enum {
+ threePhases = Indices::threePhases,
+ wPhaseOnly = Indices::wPhaseOnly,
+ gnPhaseOnly = Indices::gnPhaseOnly,
+ wnPhaseOnly = Indices::wnPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ wgPhaseOnly = Indices::wgPhaseOnly
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ static const Scalar R; // universial gas constant
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+public:
+ //! The type of the object returned by the fluidState() method
+ typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;
+
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx,
+ bool isOldSol)
+ {
+ ParentType::update(priVars,
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ isOldSol);
+
+ bool useSimpleModel = GET_PROP_VALUE(TypeTag, UseSimpleModel);
+
+ // capillary pressure parameters
+ const MaterialLawParams &materialParams =
+ problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+
+ int globalIdx = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
+
+ int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
+
+ if(!useSimpleModel)
+ {
+ /* first the saturations */
+ if (phasePresence == threePhases)
+ {
+ sw_ = priVars[switch1Idx];
+ sn_ = priVars[switch2Idx];
+ sg_ = 1. - sw_ - sn_;
+ }
+ else if (phasePresence == wPhaseOnly)
+ {
+ sw_ = 1.;
+ sn_ = 0.;
+ sg_ = 0.;
+ }
+ else if (phasePresence == gnPhaseOnly)
+ {
+ sw_ = 0.;
+ sn_ = priVars[switch1Idx];
+ sg_ = 1. - sn_;
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ sn_ = priVars[switch2Idx];
+ sw_ = 1. - sn_;
+ sg_ = 0.;
+ }
+ else if (phasePresence == gPhaseOnly)
+ {
+ sw_ = 0.;
+ sn_ = 0.;
+ sg_ = 1.;
+ }
+ else if (phasePresence == wgPhaseOnly)
+ {
+ sw_ = priVars[switch1Idx];
+ sn_ = 0.;
+ sg_ = 1. - sw_;
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(sg_);
+
+ fluidState_.setSaturation(wPhaseIdx, sw_);
+ fluidState_.setSaturation(gPhaseIdx, sg_);
+ fluidState_.setSaturation(nPhaseIdx, sn_);
+
+ /* now the pressures */
+ if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
+ {
+ pg_ = priVars[pressureIdx];
+
+ // calculate capillary pressures
+ Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
+ Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
+ Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
+
+ Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
+ Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
+
+ pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
+ pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
+ }
+ else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
+ {
+ pw_ = priVars[pressureIdx];
+
+ // calculate capillary pressures
+ Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
+ Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
+ Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
+
+ Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
+ Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
+
+ pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
+ pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(pw_);
+
+ fluidState_.setPressure(wPhaseIdx, pw_);
+ fluidState_.setPressure(gPhaseIdx, pg_);
+ fluidState_.setPressure(nPhaseIdx, pn_);
+
+ /* now the temperature */
+ if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
+ {
+ temp_ = priVars[switch1Idx];
+ }
+ else if (phasePresence == threePhases)
+ {
+ // temp from inverse pwsat and pnsat which have to sum up to pg
+ Scalar temp = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); // initial guess
+ fluidState_.setTemperature(temp);
+ Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+
+ while(std::abs(defect) > 0.01) // simply a small number chosen ...
+ {
+ Scalar deltaT = 1.e-8 * temp;
+ fluidState_.setTemperature(temp+deltaT);
+ Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+
+ fluidState_.setTemperature(temp-deltaT);
+ Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+
+ temp = temp - defect * 2. * deltaT / (fUp - fDown);
+
+ fluidState_.setTemperature(temp);
+ defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+ }
+ temp_ = temp;
+ }
+ else if (phasePresence == wgPhaseOnly)
+ {
+ // temp from inverse pwsat
+ temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
+ }
+ else if (phasePresence == gnPhaseOnly)
+ {
+ // temp from inverse pnsat
+ temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(temp_);
+
+ fluidState_.setTemperature(temp_);
+
+ // now comes the tricky part: calculate phase composition
+ if (phasePresence == threePhases) {
+
+ // all phases are present, phase compositions are a
+ // result of the the gas <-> liquid equilibrium.
+ Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
+ Scalar partPressNAPL = pg_ - partPressH2O;
+
+ Scalar xgn = partPressNAPL/pg_;
+ Scalar xgw = partPressH2O/pg_;
+
+ // Henry
+ Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
+ Scalar xww = 1.-xwn;
+
+ // Not yet filled with real numbers for the NAPL phase
+ Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
+ Scalar xnn = 1.-xnw;
+
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+ fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else if (phasePresence == wPhaseOnly) {
+ // only the water phase is present, water phase composition is
+ // stored explicitly.
+
+ // extract mole fractions in the water phase
+ Scalar xwn = priVars[switch2Idx];
+ Scalar xww = 1 - xwn;
+
+ // write water mole fractions in the fluid state
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
+
+ // note that the gas phase is actually not existing!
+ // thus, this is used as phase switch criterion
+ Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
+ Scalar xgw = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
+
+
+ // note that the NAPL phase is actually not existing!
+ // thus, this is used as phase switch criterion
+ // maybe solubility would be better than this approach via Henry
+ Scalar xnn = xwn * FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
+ Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
+
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+ fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else if (phasePresence == gnPhaseOnly) {
+
+ // only gas and NAPL phases are present
+
+ Scalar xnw = priVars[switch2Idx];
+ Scalar xnn = 1.-xnw;
+ Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
+ Scalar xgw = 1.-xgn;
+
+ // note that the water phase is actually not present
+ // the values are used as switching criteria
+ Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
+
+ // write mole fractions in the fluid state
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+ fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+
+ }
+ else if (phasePresence == wnPhaseOnly) {
+ // water and NAPL are present, phase compositions are a
+ // mole fractions of non-existing gas phase are used as switching criteria
+ Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
+ Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);;
+
+ Scalar xgn = partPressNAPL/pg_;
+ Scalar xgw = partPressH2O/pg_;
+
+ // Henry
+ Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
+ Scalar xww = 1.-xwn;
+
+ Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
+ Scalar xnn = 1.-xnw;
+
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+ fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else if (phasePresence == gPhaseOnly) {
+ // only the gas phase is present, gas phase composition is
+ // stored explicitly here below.
+
+ const Scalar xgn = priVars[switch2Idx];
+ Scalar xgw = 1 - xgn;
+
+ // write mole fractions in the fluid state
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+
+ // note that the water and NAPL phase is actually not present
+ // the values are used as switching criteria
+ Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
+ Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
+
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else if (phasePresence == wgPhaseOnly) {
+ // only water and gas phases are present
+ const Scalar xgn = priVars[switch2Idx];
+ Scalar xgw = 1 - xgn;
+
+ // write mole fractions in the fluid state
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+
+
+ Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
+ Scalar xww = 1.-xwn;
+
+ // write mole fractions in the fluid state
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
+
+ // note that the NAPL phase is actually not existing!
+ // thus, this is used as phase switch criterion
+ Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
+
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else
+ assert(false); // unhandled phase state
+ } // end of if(!UseSimpleModel), i.e. the more complex version with six phase states
+ else // use the simpler model with only two phase states
+ {
+ /* first the saturations */
+ if (phasePresence == threePhases)
+ {
+ sw_ = priVars[switch1Idx];
+ sn_ = priVars[switch2Idx];
+ sg_ = 1. - sw_ - sn_;
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ sn_ = priVars[switch2Idx];
+ sw_ = 1. - sn_;
+ sg_ = 0.;
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(sg_);
+
+ fluidState_.setSaturation(wPhaseIdx, sw_);
+ fluidState_.setSaturation(gPhaseIdx, sg_);
+ fluidState_.setSaturation(nPhaseIdx, sn_);
+
+ /* now the pressures */
+ if (phasePresence == threePhases)
+ {
+ pg_ = priVars[pressureIdx];
+
+ // calculate capillary pressures
+ Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
+ Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
+ Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
+
+ Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
+ Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
+
+ pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
+ pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ pw_ = priVars[pressureIdx];
+
+ // calculate capillary pressures
+ Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
+ Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
+ Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
+
+ Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
+ Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
+
+ pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
+ pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(pw_);
+
+ fluidState_.setPressure(wPhaseIdx, pw_);
+ fluidState_.setPressure(gPhaseIdx, pg_);
+ fluidState_.setPressure(nPhaseIdx, pn_);
+
+ /* now the temperature */
+ if (phasePresence == wnPhaseOnly)
+ {
+ temp_ = priVars[switch1Idx];
+ }
+ else if (phasePresence == threePhases)
+ {
+ if(sn_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
+ {
+ Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
+ temp_ = tempOnlyWater;
+ }
+ if(sw_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
+ {
+ Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
+ temp_ = tempOnlyNAPL;
+ }
+ else
+ {
+ // temp from inverse pwsat and pnsat which have to sum up to pg
+ Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
+ Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
+ fluidState_.setTemperature(tempOnlyWater);
+ Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+
+ Scalar temp = tempOnlyWater; // initial guess
+ int counter = 0;
+ while(std::abs(defect) > 0.01) // simply a small number chosen ...
+ {
+ Scalar deltaT = 1.e-6; // fixed number, but T should always be in the order of a few hundred Kelvin
+ fluidState_.setTemperature(temp+deltaT);
+ Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+
+ fluidState_.setTemperature(temp-deltaT);
+ Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+
+ temp = temp - defect * 2. * deltaT / (fUp - fDown);
+
+ fluidState_.setTemperature(temp);
+ defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
+ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
+ counter +=1;
+ if (counter>10) break;
+ }
+ if ((sw_>1.e-10)&&(sw_<0.01))
+ temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10);
+ if ((sn_>1.e-10)&&(sn_<0.01))
+ temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10);
+ temp_ = temp;
+ }
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(temp_);
+
+ fluidState_.setTemperature(temp_);
+
+ // now comes the tricky part: calculate phase composition
+ if (phasePresence == threePhases) {
+
+ // all phases are present, phase compositions are a
+ // result of the the gas <-> liquid equilibrium.
+ Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
+ Scalar partPressNAPL = pg_ - partPressH2O;
+ // regularized evaporation for small liquid phase saturations
+ // avoids negative saturations of liquid phases
+ if (sw_<0.02) partPressH2O *= sw_/0.02;
+ if (partPressH2O < 0.) partPressH2O = 0;
+ if (sn_<0.02) partPressNAPL *= sn_ / 0.02;
+ if (partPressNAPL < 0.) partPressNAPL = 0;
+
+ Scalar xgn = partPressNAPL/pg_;
+ Scalar xgw = partPressH2O/pg_;
+
+ // Immiscible liquid phases, mole fractions are just dummy values
+ Scalar xwn = 0;
+ Scalar xww = 1.-xwn;
+
+ // Not yet filled with real numbers for the NAPL phase
+ Scalar xnw = 0;
+ Scalar xnn = 1.-xnw;
+
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+ fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else if (phasePresence == wnPhaseOnly) {
+ // mole fractions of non-existing gas phase are used as switching criteria
+ Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
+ Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);;
+
+ Scalar xgn = partPressNAPL/pg_;
+ Scalar xgw = partPressH2O/pg_;
+
+ // immiscible liquid phases, mole fractions are just dummy values
+ Scalar xwn = 0;
+ Scalar xww = 1.-xwn;
+
+ Scalar xnw = 0;
+ Scalar xnn = 1.-xnw;
+
+ fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+ fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
+ fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+ fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
+ fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
+ fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+ }
+ else
+ assert(false); // unhandled phase state
+ }
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // Mobilities
+ const Scalar mu =
+ FluidSystem::viscosity(fluidState_,
+ phaseIdx);
+ fluidState_.setViscosity(phaseIdx,mu);
+
+ Scalar kr;
+ kr = MaterialLaw::kr(materialParams, phaseIdx,
+ fluidState_.saturation(wPhaseIdx),
+ fluidState_.saturation(nPhaseIdx),
+ fluidState_.saturation(gPhaseIdx));
+ mobility_[phaseIdx] = kr / mu;
+ Valgrind::CheckDefined(mobility_[phaseIdx]);
+ }
+
+ // material dependent parameters for NAPL adsorption
+ bulkDensTimesAdsorpCoeff_ =
+ MaterialLaw::bulkDensTimesAdsorpCoeff(materialParams);
+
+ /* ATTENTION: The conversion to effective diffusion parameters
+ * for the porous media happens at another place!
+ */
+
+ // diffusivity coefficents
+ diffusionCoefficient_[gPhaseIdx] = FluidSystem::diffusionCoefficient(fluidState_, gPhaseIdx);
+
+ diffusionCoefficient_[wPhaseIdx] = FluidSystem::diffusionCoefficient(fluidState_, wPhaseIdx);
+
+ /* no diffusion in NAPL phase considered at the moment, dummy values */
+ diffusionCoefficient_[nPhaseIdx] = 1.e-10;
+
+ Valgrind::CheckDefined(diffusionCoefficient_);
+
+ // porosity
+ porosity_ = problem.spatialParams().porosity(element,
+ fvGeometry,
+ scvIdx);
+ Valgrind::CheckDefined(porosity_);
+
+ // permeability
+ permeability_ = problem.spatialParams().intrinsicPermeability(element,
+ fvGeometry,
+ scvIdx);
+ Valgrind::CheckDefined(permeability_);
+
+ // the enthalpies (internal energies are directly calculated in the fluidstate
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
+ fluidState_.setEnthalpy(phaseIdx, h);
+ }
+
+
+ // energy related quantities not contained in the fluid state
+ asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
+ }
+
+ /*!
+ * \brief Returns the phase state for the control-volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the effective saturation of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturation(const int phaseIdx) const
+ { return fluidState_.saturation(phaseIdx); }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density(const int phaseIdx) const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Returns the molar density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(const int phaseIdx) const
+ { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar pressure(const int phaseIdx) const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperatures of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(/*phaseIdx=*/0); }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(const int phaseIdx) const
+ {
+ return mobility_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the effective capillary pressure within the control volume.
+ */
+ Scalar capillaryPressure() const
+ { return fluidState_.capillaryPressure(); }
+
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+ /*!
+ * \brief Returns the permeability within the control volume.
+ */
+ Scalar permeability() const
+ { return permeability_; }
+
+ /*!
+ * \brief Returns the diffusivity coefficient matrix
+ */
+ Dune::FieldVector<Scalar, numPhases> diffusionCoefficient() const
+ { return diffusionCoefficient_; }
+
+ /*!
+ * \brief Returns the adsorption information
+ */
+ Scalar bulkDensTimesAdsorpCoeff() const
+ { return bulkDensTimesAdsorpCoeff_; }
+
+/*!
+ * \brief Returns the total internal energy of a phase in the
+ * sub-control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar internalEnergy(int phaseIdx) const
+ { return fluidState_.internalEnergy(phaseIdx); };
+
+ /*!
+ * \brief Returns the total enthalpy of a phase in the sub-control
+ * volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar enthalpy(int phaseIdx) const
+ { return fluidState_.enthalpy(phaseIdx); };
+
+ /*!
+ * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
+ * the sub-control volume.
+ */
+ Scalar heatCapacity() const
+ { return heatCapacity_; };
+
+
+
+protected:
+
+ /*!
+ * \brief Called by update() to compute the energy related quantities
+ */
+ void updateEnergy_(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx,
+ bool isOldSol)
+ {
+ // compute and set the heat capacity of the solid phase
+ heatCapacity_ = problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
+ Valgrind::CheckDefined(heatCapacity_);
+ }
+
+ Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_;
+
+ Scalar moleFrac_[numPhases][numComponents];
+ Scalar massFrac_[numPhases][numComponents];
+
+ Scalar heatCapacity_;
+
+ Scalar porosity_; //!< Effective porosity within the control volume
+ Scalar permeability_; //!< Effective porosity within the control volume
+ Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
+ Scalar bulkDensTimesAdsorpCoeff_; //!< the basis for calculating adsorbed NAPL
+ /* We need a tensor here !! */
+ //!< Binary diffusion coefficients of the 3 components in the phases
+ Dune::FieldVector<Scalar, numPhases> diffusionCoefficient_;
+ FluidState fluidState_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+template <class TypeTag>
+const typename ThreePWaterOilVolumeVariables<TypeTag>::Scalar ThreePWaterOilVolumeVariables<TypeTag>::R = Constants<typename GET_PROP_TYPE(TypeTag, Scalar)>::R;
+
+} // end namespace
+
+#endif