From d467ef288fbe2f028193338146fe9ed4be978bcc Mon Sep 17 00:00:00 2001 From: Holger Class Date: Tue, 26 Jul 2016 12:19:32 +0200 Subject: [PATCH 1/2] the 3p2cni folder was moved into the stable part of dumux and renamed there to 3pwateroil. Necessary adaption has been made here. --- lecture/mm/heavyoil/sagd/problem.hh | 10 +++++----- lecture/mm/heavyoil/sagd/spatialparams.hh | 2 +- lecture/mm/heavyoil/sagdcyclic/problem.hh | 10 +++++----- lecture/mm/heavyoil/sagdcyclic/spatialparams.hh | 2 +- lecture/mm/heavyoil/sagdcyclichyst/problem.hh | 10 +++++----- lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh | 2 +- 6 files changed, 18 insertions(+), 18 deletions(-) diff --git a/lecture/mm/heavyoil/sagd/problem.hh b/lecture/mm/heavyoil/sagd/problem.hh index 2ea244d..b6154c1 100644 --- a/lecture/mm/heavyoil/sagd/problem.hh +++ b/lecture/mm/heavyoil/sagd/problem.hh @@ -27,8 +27,8 @@ #include -#include "../3p2cni/model.hh" -#include "../h2oheavyoilfluidsystem.hh" +#include +#include #include "spatialparams.hh" #define ISOTHERMAL 0 @@ -40,7 +40,7 @@ class SagdProblem; namespace Properties { -NEW_TYPE_TAG(SagdProblem, INHERITS_FROM(ThreePTwoCNI, SagdSpatialParams)); +NEW_TYPE_TAG(SagdProblem, INHERITS_FROM(ThreePWaterOil, SagdSpatialParams)); NEW_TYPE_TAG(SagdBoxProblem, INHERITS_FROM(BoxModel, SagdProblem)); NEW_TYPE_TAG(SagdCCProblem, INHERITS_FROM(CCModel, SagdProblem)); @@ -72,11 +72,11 @@ SET_BOOL_PROP(SagdProblem, UseMoles, false); /*! - * \ingroup ThreePTwoCNIBoxModel + * \ingroup ThreePWaterOilBoxModel * \ingroup ImplicitTestProblems * \brief Non-isothermal problem where ... * - * This problem uses the \ref ThreePTwoCNIModel. + * This problem uses the \ref ThreePWaterOilModel. * * */ template diff --git a/lecture/mm/heavyoil/sagd/spatialparams.hh b/lecture/mm/heavyoil/sagd/spatialparams.hh index 36c68c1..49e4585 100644 --- a/lecture/mm/heavyoil/sagd/spatialparams.hh +++ b/lecture/mm/heavyoil/sagd/spatialparams.hh @@ -27,7 +27,7 @@ #include -#include "../3p2cni/indices.hh" +#include #include "../3p/parkervangenuchtenzero.hh" #include "../3p/parkervangenuchtenzeroparams.hh" diff --git a/lecture/mm/heavyoil/sagdcyclic/problem.hh b/lecture/mm/heavyoil/sagdcyclic/problem.hh index e718b8e..4f9815d 100644 --- a/lecture/mm/heavyoil/sagdcyclic/problem.hh +++ b/lecture/mm/heavyoil/sagdcyclic/problem.hh @@ -28,8 +28,8 @@ #include -#include "../3p2cni/model.hh" -#include "../h2oheavyoilfluidsystem.hh" +#include +#include #include "spatialparams.hh" #define ISOTHERMAL 0 @@ -41,7 +41,7 @@ class SagdCyclicProblem; namespace Properties { -NEW_TYPE_TAG(SagdCyclicProblem, INHERITS_FROM(ThreePTwoCNI, SagdCyclicSpatialParams)); +NEW_TYPE_TAG(SagdCyclicProblem, INHERITS_FROM(ThreePWaterOil, SagdCyclicSpatialParams)); NEW_TYPE_TAG(SagdCyclicBoxProblem, INHERITS_FROM(BoxModel, SagdCyclicProblem)); @@ -74,11 +74,11 @@ SET_BOOL_PROP(SagdCyclicProblem, UseSimpleModel, true); /*! - * \ingroup ThreePTwoCNIBoxModel + * \ingroup ThreePWaterOilBoxModel * \ingroup ImplicitTestProblems * \brief Non-isothermal problem where ... * - * This problem uses the \ref ThreePTwoCNIModel. + * This problem uses the \ref ThreePWaterOilModel. * * */ template diff --git a/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh b/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh index fb6cdb0..75ca34a 100644 --- a/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh +++ b/lecture/mm/heavyoil/sagdcyclic/spatialparams.hh @@ -27,7 +27,7 @@ #include -#include "../3p2cni/indices.hh" +#include #include "../3p/parkervangenuchtenzero.hh" #include "../3p/parkervangenuchtenzeroparams.hh" diff --git a/lecture/mm/heavyoil/sagdcyclichyst/problem.hh b/lecture/mm/heavyoil/sagdcyclichyst/problem.hh index 8aad005..e038438 100644 --- a/lecture/mm/heavyoil/sagdcyclichyst/problem.hh +++ b/lecture/mm/heavyoil/sagdcyclichyst/problem.hh @@ -29,8 +29,8 @@ #include -#include "../3p2cni/model.hh" -#include "../h2oheavyoilfluidsystem.hh" +#include +#include #include "spatialparams.hh" #define ISOTHERMAL 0 @@ -42,7 +42,7 @@ class SagdCyclicHystProblem; namespace Properties { -NEW_TYPE_TAG(SagdCyclicHystProblem, INHERITS_FROM(ThreePTwoCNI, SagdCyclicHystSpatialParams)); +NEW_TYPE_TAG(SagdCyclicHystProblem, INHERITS_FROM(ThreePWaterOil, SagdCyclicHystSpatialParams)); NEW_TYPE_TAG(SagdCyclicHystBoxProblem, INHERITS_FROM(BoxModel, SagdCyclicHystProblem)); NEW_TYPE_TAG(SagdCyclicHystCCProblem, INHERITS_FROM(CCModel, SagdCyclicHystProblem)); @@ -73,11 +73,11 @@ SET_BOOL_PROP(SagdCyclicHystProblem, UseMoles, true); /*! - * \ingroup ThreePTwoCNIBoxModel + * \ingroup ThreePWaterOilBoxModel * \ingroup ImplicitTestProblems * \brief Non-isothermal problem where ... * - * This problem uses the \ref ThreePTwoCNIModel. + * This problem uses the \ref ThreePWaterOilModel. * * */ template diff --git a/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh b/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh index 62f6863..c4032db 100644 --- a/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh +++ b/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh @@ -26,7 +26,7 @@ #include -#include "../3p2cni/indices.hh" +#include #include "../3p/parkervangenuchtenzerohysteresis.hh" #include "../3p/parkervangenuchtenzerohysteresisparams.hh" -- GitLab From e57ddc0f86aa5c8f9aa6980a70282f827a484468 Mon Sep 17 00:00:00 2001 From: Holger Class Date: Tue, 26 Jul 2016 13:58:57 +0200 Subject: [PATCH 2/2] the old (local) files are deleted; no longer needed since it has been moved to dumux stable previously --- lecture/mm/heavyoil/3p2cni/fluxvariables.hh | 315 ----- lecture/mm/heavyoil/3p2cni/indices.hh | 80 -- lecture/mm/heavyoil/3p2cni/localresidual.hh | 367 ------ lecture/mm/heavyoil/3p2cni/model.hh | 1090 ----------------- .../mm/heavyoil/3p2cni/newtoncontroller.hh | 85 -- lecture/mm/heavyoil/3p2cni/properties.hh | 73 -- .../mm/heavyoil/3p2cni/propertydefaults.hh | 148 --- lecture/mm/heavyoil/3p2cni/volumevariables.hh | 876 ------------- lecture/mm/heavyoil/TODO.txt | 9 - lecture/mm/heavyoil/h2o_heavyoil.hh | 101 -- lecture/mm/heavyoil/h2oheavyoilfluidsystem.hh | 474 ------- lecture/mm/heavyoil/heavyoil.hh | 501 -------- 12 files changed, 4119 deletions(-) delete mode 100644 lecture/mm/heavyoil/3p2cni/fluxvariables.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/indices.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/localresidual.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/model.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/newtoncontroller.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/properties.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/propertydefaults.hh delete mode 100644 lecture/mm/heavyoil/3p2cni/volumevariables.hh delete mode 100644 lecture/mm/heavyoil/TODO.txt delete mode 100644 lecture/mm/heavyoil/h2o_heavyoil.hh delete mode 100644 lecture/mm/heavyoil/h2oheavyoilfluidsystem.hh delete mode 100644 lecture/mm/heavyoil/heavyoil.hh diff --git a/lecture/mm/heavyoil/3p2cni/fluxvariables.hh b/lecture/mm/heavyoil/3p2cni/fluxvariables.hh deleted file mode 100644 index 7be39b3..0000000 --- a/lecture/mm/heavyoil/3p2cni/fluxvariables.hh +++ /dev/null @@ -1,315 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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 -#include - -#include "properties.hh" - -namespace Dumux -{ - -/*! - * \ingroup ThreePTwoCNIModel - * \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 ThreePTwoCNIFluxVariables : 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 DimVector; - typedef Dune::FieldMatrix 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 diffusionCoefficientVector_i = volVarsI.diffusionCoefficient(); - Dune::FieldVector 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 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 porousDiffCoeff_; - - Scalar normalMatrixHeatFlux_; -}; - -} // end namespace - -#endif diff --git a/lecture/mm/heavyoil/3p2cni/indices.hh b/lecture/mm/heavyoil/3p2cni/indices.hh deleted file mode 100644 index 1d17e9f..0000000 --- a/lecture/mm/heavyoil/3p2cni/indices.hh +++ /dev/null @@ -1,80 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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 ThreePTwoCNIModel - * \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 ThreePTwoCNIIndices -{ - 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/lecture/mm/heavyoil/3p2cni/localresidual.hh b/lecture/mm/heavyoil/3p2cni/localresidual.hh deleted file mode 100644 index dda88d4..0000000 --- a/lecture/mm/heavyoil/3p2cni/localresidual.hh +++ /dev/null @@ -1,367 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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 ThreePTwoCNIModel - * \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 ThreePTwoCNILocalResidual: 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 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 (this); - } - - const Implementation *asImp_() const - { - return static_cast (this); - } -}; - -} // end namespace - -#endif diff --git a/lecture/mm/heavyoil/3p2cni/model.hh b/lecture/mm/heavyoil/3p2cni/model.hh deleted file mode 100644 index 1bce4a2..0000000 --- a/lecture/mm/heavyoil/3p2cni/model.hh +++ /dev/null @@ -1,1090 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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_3P2CNI_MODEL_HH -#define DUMUX_3P2CNI_MODEL_HH - -#include - -#include "properties.hh" - -namespace Dumux -{ -/*! - * \ingroup ThreePTwoCNIModel - * \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: - *
    - *
  • ... to be completed ...
  • - *
- */ -template -class ThreePTwoCNIModel: 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::Entity Vertex; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::template Codim<0>::Iterator ElementIterator; - typedef typename GridView::template Codim::Iterator VertexIterator; - - typedef Dune::FieldVector 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 (); - const VertexIterator &vEndIt = this->gridView_().template end (); - 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 (), - 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 - void addOutputVtkFields(const SolutionVector &sol, - MultiWriter &writer) - { - typedef Dune::BlockVector > ScalarField; - typedef Dune::BlockVector > 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(numDofs); - VectorField *velocityW = writer.template allocateManagedBuffer(numDofs); - VectorField *velocityG = writer.template allocateManagedBuffer(numDofs); - ImplicitVelocityOutput 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 - 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 - 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 staticDat_; - bool switchFlag_; -}; - -} - -#include "propertydefaults.hh" - -#endif diff --git a/lecture/mm/heavyoil/3p2cni/newtoncontroller.hh b/lecture/mm/heavyoil/3p2cni/newtoncontroller.hh deleted file mode 100644 index cfb7624..0000000 --- a/lecture/mm/heavyoil/3p2cni/newtoncontroller.hh +++ /dev/null @@ -1,85 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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 - -#include "properties.hh" - -namespace Dumux { -/*! - * \ingroup ThreePTwoCNIModel - * \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 ThreePTwoCNINewtonController : public NewtonController -{ - typedef NewtonController ParentType; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; - -public: - ThreePTwoCNINewtonController(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/lecture/mm/heavyoil/3p2cni/properties.hh b/lecture/mm/heavyoil/3p2cni/properties.hh deleted file mode 100644 index 318eb53..0000000 --- a/lecture/mm/heavyoil/3p2cni/properties.hh +++ /dev/null @@ -1,73 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \ingroup ThreePTwoCNIModel - */ -/*! - * \file - * - * \brief Defines the properties required for the 3p2cni fully implicit model. - */ -#ifndef DUMUX_3P2CNI_PROPERTIES_HH -#define DUMUX_3P2CNI_PROPERTIES_HH - -#include -#include - -namespace Dumux -{ - -namespace Properties -{ -////////////////////////////////////////////////////////////////// -// Type tags -////////////////////////////////////////////////////////////////// - -//! The type tags for the implicit three-phase three-component problems -NEW_TYPE_TAG(ThreePTwoCNI); -NEW_TYPE_TAG(BoxThreePTwoCNI, INHERITS_FROM(BoxModel, ThreePTwoCNI)); -NEW_TYPE_TAG(CCThreePTwoCNI, INHERITS_FROM(CCModel, ThreePTwoCNI)); - -////////////////////////////////////////////////////////////////// -// 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/lecture/mm/heavyoil/3p2cni/propertydefaults.hh b/lecture/mm/heavyoil/3p2cni/propertydefaults.hh deleted file mode 100644 index a21b6f8..0000000 --- a/lecture/mm/heavyoil/3p2cni/propertydefaults.hh +++ /dev/null @@ -1,148 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \ingroup ThreePTwoCNIModel - */ -/*! - * \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 -#include - -#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(ThreePTwoCNI, 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(ThreePTwoCNI, 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(ThreePTwoCNI, 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(ThreePTwoCNI, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params); - -//! The local residual function of the conservation equations -SET_TYPE_PROP(ThreePTwoCNI, LocalResidual, ThreePTwoCNILocalResidual); - -//! Use the 3p2cni specific newton controller for the 3p2cni model -SET_TYPE_PROP(ThreePTwoCNI, NewtonController, ThreePTwoCNINewtonController); - -//! the Model property -SET_TYPE_PROP(ThreePTwoCNI, Model, ThreePTwoCNIModel); - -//! the VolumeVariables property -SET_TYPE_PROP(ThreePTwoCNI, VolumeVariables, ThreePTwoCNIVolumeVariables); - -//! the FluxVariables property -SET_TYPE_PROP(ThreePTwoCNI, FluxVariables, ThreePTwoCNIFluxVariables); - -//! define the base flux variables to realize Darcy flow -SET_TYPE_PROP(ThreePTwoCNI, BaseFluxVariables, ImplicitDarcyFluxVariables); - -//! the upwind factor for the mobility. -SET_SCALAR_PROP(ThreePTwoCNI, ImplicitMassUpwindWeight, 1.0); - -//! set default mobility upwind weight to 1.0, i.e. fully upwind -SET_SCALAR_PROP(ThreePTwoCNI, ImplicitMobilityUpwindWeight, 1.0); - -//! Determines whether a constraint solver should be used explicitly -SET_BOOL_PROP(ThreePTwoCNI, UseSimpleModel, true); - -//! The indices required by the isothermal 3p2cni model -SET_TYPE_PROP(ThreePTwoCNI, Indices, ThreePTwoCNIIndices); - -//! The spatial parameters to be employed. -//! Use ImplicitSpatialParams by default. -SET_TYPE_PROP(ThreePTwoCNI, SpatialParams, ImplicitSpatialParams); - -// disable velocity output by default -SET_BOOL_PROP(ThreePTwoCNI, VtkAddVelocity, false); - -// enable gravity by default -SET_BOOL_PROP(ThreePTwoCNI, ProblemEnableGravity, true); - -SET_BOOL_PROP(ThreePTwoCNI, 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/lecture/mm/heavyoil/3p2cni/volumevariables.hh b/lecture/mm/heavyoil/3p2cni/volumevariables.hh deleted file mode 100644 index 0a2cd46..0000000 --- a/lecture/mm/heavyoil/3p2cni/volumevariables.hh +++ /dev/null @@ -1,876 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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 -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include "properties.hh" - -namespace Dumux -{ - -/*! - * \ingroup ThreePTwoCNIModel - * \brief Contains the quantities which are are constant within a - * finite volume in the two-phase, two-component model. - */ -template -class ThreePTwoCNIVolumeVariables : public ImplicitVolumeVariables -{ - typedef ImplicitVolumeVariables 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 MiscibleMultiPhaseComposition; - typedef Dumux::ComputeFromReferencePhase 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 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 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 diffusionCoefficient_; - FluidState fluidState_; - -private: - Implementation &asImp_() - { return *static_cast(this); } - - const Implementation &asImp_() const - { return *static_cast(this); } -}; - -template -const typename ThreePTwoCNIVolumeVariables::Scalar ThreePTwoCNIVolumeVariables::R = Constants::R; - -} // end namespace - -#endif diff --git a/lecture/mm/heavyoil/TODO.txt b/lecture/mm/heavyoil/TODO.txt deleted file mode 100644 index 2a5419d..0000000 --- a/lecture/mm/heavyoil/TODO.txt +++ /dev/null @@ -1,9 +0,0 @@ -This is an improvised version of the SAGD models. - -We plan to add the 3p2cni model to dumux-stable. It has as a special feature -that temperature is not always primary variable. If a gaseous phase exists, -temperature is no degree of freedom but linked to the vapor pressure(s) of the -boiling liquids. - -There is a 2p1cni model which follows the same philosophy. Thus, both the -3p2cni and the 2p1cni should be treated the same way. diff --git a/lecture/mm/heavyoil/h2o_heavyoil.hh b/lecture/mm/heavyoil/h2o_heavyoil.hh deleted file mode 100644 index 68ac6e8..0000000 --- a/lecture/mm/heavyoil/h2o_heavyoil.hh +++ /dev/null @@ -1,101 +0,0 @@ -// -*- 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: .@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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Binary coefficients for water and tce. - */ -#ifndef DUMUX_BINARY_COEFF_H2O_HEAVYOIL_HH -#define DUMUX_BINARY_COEFF_H2O_HEAVYOIL_HH - -#include -#include - -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 - 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 - 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 - 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 - 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/lecture/mm/heavyoil/h2oheavyoilfluidsystem.hh b/lecture/mm/heavyoil/h2oheavyoilfluidsystem.hh deleted file mode 100644 index 9e020d7..0000000 --- a/lecture/mm/heavyoil/h2oheavyoilfluidsystem.hh +++ /dev/null @@ -1,474 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \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 -#include -#include -#include -#include - -#include - -#include - -namespace Dumux -{ -namespace FluidSystems -{ - -/*! - * \brief A compositional fluid with water and heavy oil - * components in both, the liquid and the gas phase. - */ -template > > -class H2OHeavyOil - : public BaseFluidSystem > -{ - typedef H2OHeavyOil ThisType; - typedef BaseFluidSystem Base; - -public: - typedef Dumux::HeavyOil 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 - static Scalar density(const FluidState &fluidState, int phaseIdx) - { - if (phaseIdx == wPhaseIdx) { - // See: Ochs 2008 - // \todo: proper citation - 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 - 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 - 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 - 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 - 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 - 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 - 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 - static Scalar heatCapacity(const FluidState &fluidState, - int phaseIdx) - { - DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()"); - } - - using Base::thermalConductivity; - template - 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, 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 > H2O; - * - * // Apply e.g. untabulated water: - * typedef Dumux::H2O H2O; - * }; - * - * Also remember to initialize tabulated components (FluidSystem::init()), while this - * is not necessary for non-tabularized ones. - */ -template -class H2OHeavyOilFluidSystem -: public FluidSystems::H2OHeavyOil -{}; -#endif - -} // end namespace Dumux - -#endif diff --git a/lecture/mm/heavyoil/heavyoil.hh b/lecture/mm/heavyoil/heavyoil.hh deleted file mode 100644 index 12107f0..0000000 --- a/lecture/mm/heavyoil/heavyoil.hh +++ /dev/null @@ -1,501 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \file - * - * \ingroup Components - * - * \brief Properties of heavyoil. - * - */ -#ifndef DUMUX_HEAVYOIL_HH -#define DUMUX_HEAVYOIL_HH - -#include -#include -#include - -namespace Dumux -{ -/*! - * \ingroup Components - * \brief heavyoil - * - * \tparam Scalar The type used for scalar values - */ -template -class HeavyOil : public Component > -{ - typedef Dumux::Constants 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::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 -- GitLab