From 5c03f62f94d9a550cb37c0382e3bc178977fe292 Mon Sep 17 00:00:00 2001
From: vishal jambhekar <vishal.jambhekar@iws.uni-stuttgart.de>
Date: Mon, 7 May 2012 13:27:50 +0000
Subject: [PATCH] New naminig convention FIRST SESSION 3p3c model
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8235 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/3p3c/3p3cfluxvariables.hh | 279 ++++++++++----------
dumux/boxmodels/3p3c/3p3cindices.hh | 13 +-
dumux/boxmodels/3p3c/3p3clocalresidual.hh | 73 +++--
dumux/boxmodels/3p3c/3p3cmodel.hh | 79 +++---
dumux/boxmodels/3p3c/3p3cvolumevariables.hh | 106 ++++----
5 files changed, 280 insertions(+), 270 deletions(-)
diff --git a/dumux/boxmodels/3p3c/3p3cfluxvariables.hh b/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
index 4b948d810e..15f728595f 100644
--- a/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
+++ b/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
@@ -74,8 +74,8 @@ class ThreePThreeCFluxVariables
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
- typedef Dune::FieldVector<Scalar, dimWorld> Vector;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
@@ -83,9 +83,9 @@ class ThreePThreeCFluxVariables
nPhaseIdx = Indices::nPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx,
- wCompIdx = Indices::wCompIdx,
- cCompIdx = Indices::cCompIdx,
- aCompIdx = Indices::aCompIdx
+ comp0Idx = Indices::comp0Idx,
+ comp1Idx = Indices::comp1Idx,
+ comp2Idx = Indices::comp2Idx
};
public:
@@ -94,52 +94,51 @@ public:
*
* \param problem The problem
* \param element The finite element
- * \param elemGeom The finite-volume geometry in the box scheme
+ * \param fvGeometry The finite-volume geometry in the box scheme
* \param faceIdx The local index of the SCV (sub-control-volume) face
- * \param elemDat The volume variables of the current element
+ * \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
*/
ThreePThreeCFluxVariables(const Problem &problem,
const Element &element,
- const FVElementGeometry &elemGeom,
- int faceIdx,
- const ElementVolumeVariables &elemDat,
+ const FVElementGeometry &fvGeometry,
+ const int faceIdx,
+ const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
- : fvElemGeom_(elemGeom), onBoundary_(onBoundary)
+ : fvGeometry_(fvGeometry),scvfIdx_(faceIdx), onBoundary_(onBoundary)
{
- scvfIdx_ = faceIdx;
-
+ //scvfIdx_ = faceIdx;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- densityAtIP_[phaseIdx] = Scalar(0);
- molarDensityAtIP_[phaseIdx] = Scalar(0);
+ density_[phaseIdx] = Scalar(0);
+ molarDensity_[phaseIdx] = Scalar(0);
potentialGrad_[phaseIdx] = Scalar(0);
- wConcentrationGrad_[phaseIdx] = Scalar(0);
- cConcentrationGrad_[phaseIdx] = Scalar(0);
- aConcentrationGrad_[phaseIdx] = Scalar(0);
- molarWConcGrad_[phaseIdx] = Scalar(0);
- molarCConcGrad_[phaseIdx] = Scalar(0);
- molarAConcGrad_[phaseIdx] = Scalar(0);
+ massFractionGradComp0_[phaseIdx] = Scalar(0);
+ massFractionGradComp1_[phaseIdx] = Scalar(0);
+ massFractionGradComp2_[phaseIdx] = Scalar(0);
+ moleFractionGradComp0_[phaseIdx] = Scalar(0);
+ moleFractionGradComp1_[phaseIdx] = Scalar(0);
+ moleFractionGradComp2_[phaseIdx] = Scalar(0);
}
- calculateGradients_(problem, element, elemDat);
- calculateVelocities_(problem, element, elemDat);
- calculateDiffCoeffPM_(problem, element, elemDat);
+ calculateGradients_(problem, element, elemVolVars);
+ calculateVelocities_(problem, element, elemVolVars);
+ calculateporousDiffCoeff_(problem, element, elemVolVars);
};
private:
void calculateGradients_(const Problem &problem,
const Element &element,
- const ElementVolumeVariables &elemDat)
+ const ElementVolumeVariables &elemVolVars)
{
// calculate gradients
- Vector tmp(0.0);
+ GlobalPosition tmp(0.0);
for (int idx = 0;
- idx < fvElemGeom_.numFAP;
+ idx < fvGeometry_.numVertices;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
- const Vector &feGrad = face().grad[idx];
+ const GlobalPosition &feGrad = face().grad[idx];
// index for the element volume variables
int volVarsIdx = face().fapIndices[idx];
@@ -149,85 +148,85 @@ private:
{
// the pressure gradient
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].pressure(phaseIdx);
+ tmp *= elemVolVars[idx].pressure(phaseIdx);
potentialGrad_[phaseIdx] += tmp;
}
// the concentration gradient of the components
// component in the phases
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
- wConcentrationGrad_[wPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(wPhaseIdx, comp0Idx);
+ massFractionGradComp0_[wPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
- wConcentrationGrad_[nPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(nPhaseIdx, comp0Idx);
+ massFractionGradComp0_[nPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
- wConcentrationGrad_[gPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(gPhaseIdx, comp0Idx);
+ massFractionGradComp0_[gPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, cCompIdx);
- cConcentrationGrad_[wPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(wPhaseIdx, comp1Idx);
+ massFractionGradComp1_[wPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, cCompIdx);
- cConcentrationGrad_[nPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(nPhaseIdx, comp1Idx);
+ massFractionGradComp1_[nPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, cCompIdx);
- cConcentrationGrad_[gPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(gPhaseIdx, comp1Idx);
+ massFractionGradComp1_[gPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, aCompIdx);
- aConcentrationGrad_[wPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(wPhaseIdx, comp2Idx);
+ massFractionGradComp2_[wPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, aCompIdx);
- aConcentrationGrad_[nPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(nPhaseIdx, comp2Idx);
+ massFractionGradComp2_[nPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, aCompIdx);
- aConcentrationGrad_[gPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().massFraction(gPhaseIdx, comp2Idx);
+ massFractionGradComp2_[gPhaseIdx] += tmp;
// the molar concentration gradients of the components
// in the phases
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
- molarWConcGrad_[wPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(wPhaseIdx, comp0Idx);
+ moleFractionGradComp0_[wPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
- molarWConcGrad_[nPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(nPhaseIdx, comp0Idx);
+ moleFractionGradComp0_[nPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
- molarWConcGrad_[gPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(gPhaseIdx, comp0Idx);
+ moleFractionGradComp0_[gPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, cCompIdx);
- molarCConcGrad_[wPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(wPhaseIdx, comp1Idx);
+ moleFractionGradComp1_[wPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, cCompIdx);
- molarCConcGrad_[nPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(nPhaseIdx, comp1Idx);
+ moleFractionGradComp1_[nPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, cCompIdx);
- molarCConcGrad_[gPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(gPhaseIdx, comp1Idx);
+ moleFractionGradComp1_[gPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, aCompIdx);
- molarAConcGrad_[wPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(wPhaseIdx, comp2Idx);
+ moleFractionGradComp2_[wPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, aCompIdx);
- molarAConcGrad_[nPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(nPhaseIdx, comp2Idx);
+ moleFractionGradComp2_[nPhaseIdx] += tmp;
tmp = feGrad;
- tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, aCompIdx);
- molarAConcGrad_[gPhaseIdx] += tmp;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(gPhaseIdx, comp2Idx);
+ moleFractionGradComp2_[gPhaseIdx] += tmp;
}
// correct the pressure gradients by the hydrostatic
@@ -236,26 +235,26 @@ private:
{
int i = face().i;
int j = face().j;
- Scalar fI = rhoFactor_(phaseIdx, i, elemDat);
- Scalar fJ = rhoFactor_(phaseIdx, j, elemDat);
+ Scalar fI = rhoFactor_(phaseIdx, i, elemVolVars);
+ Scalar fJ = rhoFactor_(phaseIdx, j, elemVolVars);
if (fI + fJ <= 0)
fI = fJ = 0.5; // doesn't matter because no phase is
// present in both cells!
- densityAtIP_[phaseIdx] =
- (fI*elemDat[i].density(phaseIdx) +
- fJ*elemDat[j].density(phaseIdx))
+ density_[phaseIdx] =
+ (fI*elemVolVars[i].density(phaseIdx) +
+ fJ*elemVolVars[j].density(phaseIdx))
/
(fI + fJ);
// phase density
- molarDensityAtIP_[phaseIdx]
+ molarDensity_[phaseIdx]
=
- (fI*elemDat[i].molarDensity(phaseIdx) +
- fJ*elemDat[j].molarDensity(phaseIdx))
+ (fI*elemVolVars[i].molarDensity(phaseIdx) +
+ fJ*elemVolVars[j].molarDensity(phaseIdx))
/
(fI + fJ);
tmp = problem.gravity();
- tmp *= densityAtIP_[phaseIdx];
+ tmp *= density_[phaseIdx];
potentialGrad_[phaseIdx] -= tmp;
}
@@ -278,20 +277,20 @@ private:
void calculateVelocities_(const Problem &problem,
const Element &element,
- const ElementVolumeVariables &elemDat)
+ const ElementVolumeVariables &elemVolVars)
{
- const SpatialParameters &spatialParams = problem.spatialParameters();
+ const SpatialParameters &spatialParams = problem.spatialParams();
// multiply the pressure potential with the intrinsic
// permeability
- Tensor K;
+ DimMatrix K;
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
spatialParams.meanK(K,
spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
+ fvGeometry_,
face().i),
spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
+ fvGeometry_,
face().j));
K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
KmvpNormal_[phaseIdx] = - (Kmvp_[phaseIdx] * face().normal);
@@ -310,16 +309,16 @@ private:
}
}
- void calculateDiffCoeffPM_(const Problem &problem,
+ void calculateporousDiffCoeff_(const Problem &problem,
const Element &element,
- const ElementVolumeVariables &elemDat)
+ const ElementVolumeVariables &elemVolVars)
{
- const VolumeVariables &vDat_i = elemDat[face().i];
- const VolumeVariables &vDat_j = elemDat[face().j];
+ const VolumeVariables &volVarsI = elemVolVars[face().i];
+ const VolumeVariables &volVarsJ = elemVolVars[face().j];
- Dune::FieldMatrix<Scalar, numPhases, numComponents> diffusionCoefficientMatrix_i = vDat_i.diffusionCoefficient();
- Dune::FieldMatrix<Scalar, numPhases, numComponents> diffusionCoefficientMatrix_j = vDat_j.diffusionCoefficient();
+ Dune::FieldMatrix<Scalar, numPhases, numComponents> diffusionCoefficientMatrix_i = volVarsI.diffusionCoefficient();
+ Dune::FieldMatrix<Scalar, numPhases, numComponents> diffusionCoefficientMatrix_j = volVarsJ.diffusionCoefficient();
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
@@ -328,33 +327,33 @@ private:
/* \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 (vDat_i.saturation(phaseIdx) <= 0 ||
- vDat_j.saturation(phaseIdx) <= 0)
+ if (volVarsI.saturation(phaseIdx) <= 0 ||
+ volVarsJ.saturation(phaseIdx) <= 0)
{
- porousDiffCoeff_[phaseIdx][wCompIdx] = 0.0;
- porousDiffCoeff_[phaseIdx][cCompIdx] = 0.0;
- porousDiffCoeff_[phaseIdx][aCompIdx] = 0.0;
+ porousDiffCoeff_[phaseIdx][comp0Idx] = 0.0;
+ porousDiffCoeff_[phaseIdx][comp1Idx] = 0.0;
+ porousDiffCoeff_[phaseIdx][comp2Idx] = 0.0;
continue;
}
// calculate tortuosity at the nodes i and j needed
// for porous media diffusion coefficient
- Scalar tau_i =
- 1.0/(vDat_i.porosity() * vDat_i.porosity()) *
- pow(vDat_i.porosity() * vDat_i.saturation(phaseIdx), 7.0/3);
- Scalar tau_j =
- 1.0/(vDat_j.porosity() * vDat_j.porosity()) *
- pow(vDat_j.porosity() * vDat_j.saturation(phaseIdx), 7.0/3);
+ 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][wCompIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * diffusionCoefficientMatrix_i[phaseIdx][wCompIdx],
- vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * diffusionCoefficientMatrix_j[phaseIdx][wCompIdx]);
- porousDiffCoeff_[phaseIdx][cCompIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * diffusionCoefficientMatrix_i[phaseIdx][cCompIdx],
- vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * diffusionCoefficientMatrix_j[phaseIdx][cCompIdx]);
- porousDiffCoeff_[phaseIdx][aCompIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * diffusionCoefficientMatrix_i[phaseIdx][aCompIdx],
- vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * diffusionCoefficientMatrix_j[phaseIdx][aCompIdx]);
+ porousDiffCoeff_[phaseIdx][comp0Idx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][comp0Idx],
+ volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][comp0Idx]);
+ porousDiffCoeff_[phaseIdx][comp1Idx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][comp1Idx],
+ volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][comp1Idx]);
+ porousDiffCoeff_[phaseIdx][comp2Idx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][comp2Idx],
+ volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][comp2Idx]);
}
}
@@ -377,9 +376,9 @@ public:
/*!
* \brief Return the pressure potential multiplied with the
- * intrinsic permeability as vector (for velocity output)
+ * intrinsic permeability as GlobalPosition (for velocity output)
*/
- Vector Kmvp(int phaseIdx) const
+ GlobalPosition Kmvp(int phaseIdx) const
{ return Kmvp_[phaseIdx]; }
/*!
@@ -406,67 +405,79 @@ public:
* \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
* point.
*/
+ DUMUX_DEPRECATED_MSG("use density instead")
Scalar densityAtIP(int phaseIdx) const
- { return densityAtIP_[phaseIdx]; }
+ {
+ density(phaseIdx);
+ }
+
+ /*!
+ * \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 at the integration
* point.
*/
+ DUMUX_DEPRECATED_MSG("use molarDensity instead")
Scalar molarDensityAtIP(int phaseIdx) const
- { return molarDensityAtIP_[phaseIdx]; }
+ {
+ molarDensity(phaseIdx);
+ }
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase.
+ */
+ Scalar molarDensity(int phaseIdx) const
+ { return molarDensity_[phaseIdx]; }
/*!
* \brief The concentration gradients of the components in a phase.
*/
- const Vector &wConcentrationGrad(int phaseIdx) const
- { return wConcentrationGrad_[phaseIdx]; };
+ const GlobalPosition &wmassFractionGrad(int phaseIdx) const
+ { return massFractionGradComp0_[phaseIdx]; };
- const Vector &cConcentrationGrad(int phaseIdx) const
- { return cConcentrationGrad_[phaseIdx]; };
+ const GlobalPosition &nmassFractionGrad(int phaseIdx) const
+ { return massFractionGradComp1_[phaseIdx]; };
- const Vector &aConcentrationGrad(int phaseIdx) const
- { return aConcentrationGrad_[phaseIdx]; };
+ const GlobalPosition &amassFractionGrad(int phaseIdx) const
+ { return massFractionGradComp2_[phaseIdx]; };
/*!
* \brief The molar concentration gradients of the components in a phase.
*/
- const Vector &molarWConcGrad(int phaseIdx) const
- { return molarWConcGrad_[phaseIdx]; };
+ const GlobalPosition &wmoleFractionGrad(int phaseIdx) const
+ { return moleFractionGradComp0_[phaseIdx]; };
- const Vector &molarCConcGrad(int phaseIdx) const
- { return molarCConcGrad_[phaseIdx]; };
+ const GlobalPosition &nmoleFractionGrad(int phaseIdx) const
+ { return moleFractionGradComp1_[phaseIdx]; };
- const Vector &molarAConcGrad(int phaseIdx) const
- { return molarAConcGrad_[phaseIdx]; };
+ const GlobalPosition &amoleFractionGrad(int phaseIdx) const
+ { return moleFractionGradComp2_[phaseIdx]; };
const SCVFace &face() const
- {
- if (onBoundary_)
- return fvElemGeom_.boundaryFace[scvfIdx_];
- else
- return fvElemGeom_.subContVolFace[scvfIdx_];
- }
+ { return fvGeometry_.subContVolFace[scvfIdx_]; }
protected:
- const FVElementGeometry &fvElemGeom_;
+ const FVElementGeometry &fvGeometry_;
int scvfIdx_;
const bool onBoundary_;
// gradients
- Vector potentialGrad_[numPhases];
- Vector wConcentrationGrad_[numPhases];
- Vector cConcentrationGrad_[numPhases];
- Vector aConcentrationGrad_[numPhases];
- Vector molarWConcGrad_[numPhases];
- Vector molarCConcGrad_[numPhases];
- Vector molarAConcGrad_[numPhases];
+ GlobalPosition potentialGrad_[numPhases];
+ GlobalPosition massFractionGradComp0_[numPhases];
+ GlobalPosition massFractionGradComp1_[numPhases];
+ GlobalPosition massFractionGradComp2_[numPhases];
+ GlobalPosition moleFractionGradComp0_[numPhases];
+ GlobalPosition moleFractionGradComp1_[numPhases];
+ GlobalPosition moleFractionGradComp2_[numPhases];
// density of each face at the integration point
- Scalar densityAtIP_[numPhases], molarDensityAtIP_[numPhases];
+ Scalar density_[numPhases], molarDensity_[numPhases];
// intrinsic permeability times pressure potential gradient
- Vector Kmvp_[numPhases];
+ GlobalPosition Kmvp_[numPhases];
// projected on the face normal
Scalar KmvpNormal_[numPhases];
diff --git a/dumux/boxmodels/3p3c/3p3cindices.hh b/dumux/boxmodels/3p3c/3p3cindices.hh
index b2cac440d2..73083c0491 100644
--- a/dumux/boxmodels/3p3c/3p3cindices.hh
+++ b/dumux/boxmodels/3p3c/3p3cindices.hh
@@ -53,9 +53,9 @@ public:
static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the gas phase
// Component indices
- static const int wCompIdx = 0; //!< Index of the water component
- static const int cCompIdx = 1; //!< Index of the NAPL component
- static const int aCompIdx = 2; //!< Index of the air component
+ static const int comp0Idx = 0; //!< Index of the water component
+ static const int comp1Idx = 1; //!< Index of the NAPL component
+ static const int comp2Idx = 2; //!< Index of the gas component
// present phases (-> 'pseudo' primary variable)
static const int threePhases = 1; //!< All three phases are present
@@ -75,10 +75,9 @@ public:
static const int SOrX2Idx = switch2Idx; //!< Index of the either the saturation of the gas phase or the mass fraction secondary component if a phase is not present
// equation indices
- static const int conti0EqIdx = PVOffset; //!< Index of the mass conservation equation for the first component
- static const int contiWEqIdx = conti0EqIdx + wCompIdx; //!< Index of the mass conservation equation for the water component
- static const int contiCEqIdx = conti0EqIdx + cCompIdx; //!< Index of the mass conservation equation for the contaminant component
- static const int contiAEqIdx = conti0EqIdx + aCompIdx; //!< Index of the mass conservation equation for the air component
+ static const int conti0EqIdx = PVOffset + comp0Idx; //!< Index of the mass conservation equation for the water component
+ static const int conti1EqIdx = conti0EqIdx + comp1Idx; //!< Index of the mass conservation equation for the contaminant component
+ static const int conti2EqIdx = conti0EqIdx + comp2Idx; //!< Index of the mass conservation equation for the gas component
};
}
diff --git a/dumux/boxmodels/3p3c/3p3clocalresidual.hh b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
index 4cf1092bea..594c43a909 100644
--- a/dumux/boxmodels/3p3c/3p3clocalresidual.hh
+++ b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
@@ -68,18 +68,17 @@ protected:
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
- conti0EqIdx = Indices::conti0EqIdx,
- contiWEqIdx = Indices::contiWEqIdx,
- contiCEqIdx = Indices::contiCEqIdx,
- contiAEqIdx = Indices::contiAEqIdx,
+ 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
+ conti2EqIdx = Indices::conti2EqIdx,//!< Index of the mass conservation equation for the gas component
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx,
- wCompIdx = Indices::wCompIdx,
- cCompIdx = Indices::cCompIdx,
- aCompIdx = Indices::aCompIdx
+ comp0Idx = Indices::comp0Idx,
+ comp1Idx = Indices::comp1Idx,
+ comp2Idx = Indices::comp2Idx
};
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
@@ -135,9 +134,9 @@ public:
* \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, int faceIdx, const bool onBoundary=false) const
+ void computeFlux(PrimaryVariables &flux, const int faceIdx,bool onBoundary=false) const
{
- FluxVariables vars(this->problem_(),
+ FluxVariables fluxVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceIdx,
@@ -145,8 +144,8 @@ public:
onBoundary);
flux = 0;
- asImp_()->computeAdvectiveFlux(flux, vars);
- asImp_()->computeDiffusiveFlux(flux, vars);
+ asImp_()->computeAdvectiveFlux(flux, fluxVars);
+ asImp_()->computeDiffusiveFlux(flux, fluxVars);
}
/*!
@@ -154,10 +153,10 @@ public:
* a face of a subcontrol volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
- * \param vars The flux variables at the current SCV
+ * \param fluxVars The flux variables at the current SCV
*/
- void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
+ void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
Scalar massUpwindWeight = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
@@ -168,8 +167,8 @@ public:
{
// data attached to upstream and the downstream vertices
// of the current phase
- const VolumeVariables &up = this->curVolVars_(vars.upstreamIdx(phaseIdx));
- const VolumeVariables &dn = this->curVolVars_(vars.downstreamIdx(phaseIdx));
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
@@ -180,7 +179,7 @@ public:
// if alpha == 1 (which is mostly the case) then, the downstream
// node is not evaluated
int eqIdx = conti0EqIdx + compIdx;
- flux[eqIdx] += vars.KmvpNormal(phaseIdx)
+ flux[eqIdx] += fluxVars.KmvpNormal(phaseIdx)
* (massUpwindWeight
* up.mobility(phaseIdx)
* up.fluidState().molarDensity(phaseIdx)
@@ -199,39 +198,39 @@ public:
* a face of a subcontrol volume.
*
* \param flux The diffusive flux over the sub-control-volume face for each component
- * \param vars The flux variables at the current SCV
+ * \param fluxVars The flux variables at the current SCV
*/
- void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
- // TODO: reference!? Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = vars.porousDiffCoeff();
+ // TODO: reference!? Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = fluxVars.porousDiffCoeff();
// add diffusive flux of gas component in liquid phase
- Scalar tmp = - vars.porousDiffCoeff()[wPhaseIdx][aCompIdx] * vars.molarDensityAtIP(wPhaseIdx);
- tmp *= (vars.molarAConcGrad(wPhaseIdx) * vars.face().normal);
+ Scalar tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx][comp2Idx] * fluxVars.molarDensity(wPhaseIdx);
+ tmp *= (fluxVars.amoleFractionGrad(wPhaseIdx) * fluxVars.face().normal);
Scalar jAW = tmp;
- tmp = - vars.porousDiffCoeff()[wPhaseIdx][cCompIdx] * vars.molarDensityAtIP(wPhaseIdx);
- tmp *= (vars.molarCConcGrad(wPhaseIdx) * vars.face().normal);
+ tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx][comp1Idx] * fluxVars.molarDensity(wPhaseIdx);
+ tmp *= (fluxVars.nmoleFractionGrad(wPhaseIdx) * fluxVars.face().normal);
Scalar jCW = tmp;
Scalar jWW = -(jAW+jCW);
- tmp = - vars.porousDiffCoeff()[gPhaseIdx][wCompIdx] * vars.molarDensityAtIP(gPhaseIdx);
- tmp *= (vars.molarWConcGrad(gPhaseIdx) * vars.face().normal);
+ tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx][comp0Idx] * fluxVars.molarDensity(gPhaseIdx);
+ tmp *= (fluxVars.wmoleFractionGrad(gPhaseIdx) * fluxVars.face().normal);
Scalar jWG = tmp;
- tmp = - vars.porousDiffCoeff()[gPhaseIdx][cCompIdx] * vars.molarDensityAtIP(gPhaseIdx);
- tmp *= (vars.molarCConcGrad(gPhaseIdx) * vars.face().normal);
+ tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx][comp1Idx] * fluxVars.molarDensity(gPhaseIdx);
+ tmp *= (fluxVars.nmoleFractionGrad(gPhaseIdx) * fluxVars.face().normal);
Scalar jCG = tmp;
Scalar jAG = -(jWG+jCG);
- tmp = - vars.porousDiffCoeff()[nPhaseIdx][wCompIdx] * vars.molarDensityAtIP(nPhaseIdx);
- tmp *= (vars.molarWConcGrad(nPhaseIdx) * vars.face().normal);
+ tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx][comp0Idx] * fluxVars.molarDensity(nPhaseIdx);
+ tmp *= (fluxVars.wmoleFractionGrad(nPhaseIdx) * fluxVars.face().normal);
Scalar jWN = tmp;
- tmp = - vars.porousDiffCoeff()[nPhaseIdx][aCompIdx] * vars.molarDensityAtIP(nPhaseIdx);
- tmp *= (vars.molarAConcGrad(nPhaseIdx) * vars.face().normal);
+ tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx][comp2Idx] * fluxVars.molarDensity(nPhaseIdx);
+ tmp *= (fluxVars.amoleFractionGrad(nPhaseIdx) * fluxVars.face().normal);
Scalar jAN = tmp;
Scalar jCN = -(jAN+jWN);
@@ -248,23 +247,23 @@ public:
jAN = 0;
*/
- flux[contiWEqIdx] += jWW+jWG+jWN;
- flux[contiCEqIdx] += jCW+jCG+jCN;
- flux[contiAEqIdx] += jAW+jAG+jAN;
+ flux[conti0EqIdx] += jWW+jWG+jWN;
+ flux[conti1EqIdx] += jCW+jCG+jCN;
+ flux[conti2EqIdx] += jAW+jAG+jAN;
}
/*!
* \brief Calculate the source term of the equation
*
* \param q The source/sink in the SCV for each component
- * \param localVertexIdx The index of the SCV
+ * \param scvIdx The index of the SCV
*/
- void computeSource(PrimaryVariables &q, int localVertexIdx)
+ void computeSource(PrimaryVariables &q, int scvIdx)
{
this->problem_().boxSDSource(q,
this->elem_(),
this->fvElemGeom_(),
- localVertexIdx,
+ scvIdx,
this->curVolVars_());
}
diff --git a/dumux/boxmodels/3p3c/3p3cmodel.hh b/dumux/boxmodels/3p3c/3p3cmodel.hh
index 3b454ca9c6..f40ef057cd 100644
--- a/dumux/boxmodels/3p3c/3p3cmodel.hh
+++ b/dumux/boxmodels/3p3c/3p3cmodel.hh
@@ -126,9 +126,9 @@ class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
nPhaseIdx = Indices::nPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx,
- wCompIdx = Indices::wCompIdx,
- cCompIdx = Indices::cCompIdx,
- aCompIdx = Indices::aCompIdx,
+ comp0Idx = Indices::comp0Idx,
+ comp1Idx = Indices::comp1Idx,
+ comp2Idx = Indices::comp2Idx,
threePhases = Indices::threePhases,
wPhaseOnly = Indices::wPhaseOnly,
@@ -307,7 +307,7 @@ public:
ScalarField *rank =
writer.allocateManagedBuffer (numElements);
- FVElementGeometry fvElemGeom;
+ FVElementGeometry fvGeometry;
VolumeVariables volVars;
ElementIterator elemIt = this->gridView_().template begin<0>();
@@ -316,7 +316,7 @@ public:
{
int idx = this->problem_().elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
- fvElemGeom.update(this->gridView_(), *elemIt);
+ fvGeometry.update(this->gridView_(), *elemIt);
int numVerts = elemIt->template count<dim> ();
for (int i = 0; i < numVerts; ++i)
@@ -325,7 +325,7 @@ public:
volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
- fvElemGeom,
+ fvGeometry,
i,
false);
@@ -438,14 +438,14 @@ public:
for (unsigned i = 0; i < staticVertexDat_.size(); ++i)
staticVertexDat_[i].visited = false;
- FVElementGeometry fvElemGeom;
+ FVElementGeometry fvGeometry;
static VolumeVariables volVars;
ElementIterator it = this->gridView_().template begin<0> ();
const ElementIterator &endit = this->gridView_().template end<0> ();
for (; it != endit; ++it)
{
- fvElemGeom.update(this->gridView_(), *it);
- for (int i = 0; i < fvElemGeom.numSCV; ++i)
+ fvGeometry.update(this->gridView_(), *it);
+ for (int i = 0; i < fvGeometry.numVertices; ++i)
{
int globalIdx = this->vertexMapper().map(*it, i, dim);
@@ -456,7 +456,7 @@ public:
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
*it,
- fvElemGeom,
+ fvGeometry,
i,
false);
const GlobalPosition &global = it->geometry().corner(i);
@@ -533,7 +533,8 @@ protected:
// perform variable switch at a vertex; Returns true if a
// variable switch was performed.
bool primaryVarSwitch_(SolutionVector &globalSol,
- const VolumeVariables &volVars, int globalIdx,
+ const VolumeVariables &volVars,
+ int globalIdx,
const GlobalPosition &globalPos)
{
// evaluate primary variable switch
@@ -558,7 +559,7 @@ protected:
newPhasePresence = wnPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, aCompIdx);
+ = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
}
else if (volVars.saturation(wPhaseIdx) <= Smin)
{
@@ -570,7 +571,7 @@ protected:
newPhasePresence = gnPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
}
else if (volVars.saturation(nPhaseIdx) <= Smin)
{
@@ -582,7 +583,7 @@ protected:
newPhasePresence = wgPhaseOnly;
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
}
}
else if (phasePresence == wPhaseOnly)
@@ -591,9 +592,9 @@ protected:
bool NAPLFlag = 0;
// calculate fractions of the partial pressures in the
// hypothetical gas phase
- Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
- Scalar xag = volVars.fluidState().moleFraction(gPhaseIdx, aCompIdx);
- Scalar xcg = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+ Scalar xag = volVars.fluidState().moleFraction(gPhaseIdx, comp2Idx);
+ Scalar xcg = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
/* take care:
for xag in case wPhaseOnly we compute xag=henry_air*xaw
for xwg in case wPhaseOnly we compute xwg=pwsat
@@ -618,7 +619,7 @@ protected:
}
// calculate fractions in the hypothetical NAPL phase
- Scalar xnc = volVars.fluidState().moleFraction(nPhaseIdx, cCompIdx);
+ Scalar xnc = volVars.fluidState().moleFraction(nPhaseIdx, comp1Idx);
/* take care:
for xnc in case wPhaseOnly we compute xnc=henry_mesitylene*xcw,
where a hypothetical gas pressure is assumed for the Henry
@@ -659,7 +660,7 @@ protected:
{
newPhasePresence = wnPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, aCompIdx);
+ = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
globalSol[globalIdx][switch2Idx] = 0.0001;
}
}
@@ -684,7 +685,7 @@ protected:
// calculate fractions of the hypothetical water phase
- Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, comp0Idx);
/*
take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
If this is larger than 1, then water appears
@@ -717,15 +718,15 @@ protected:
newPhasePresence = wgPhaseOnly;
globalSol[globalIdx][switch1Idx] = 0.0001;
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
}
else if ((waterFlag == 0) && (NAPLFlag == 1))
{
newPhasePresence = gPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
}
}
else if (phasePresence == wnPhaseOnly)
@@ -749,9 +750,9 @@ protected:
// calculate fractions of the partial pressures in the
// hypothetical gas phase
- Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
- Scalar xag = volVars.fluidState().moleFraction(gPhaseIdx, aCompIdx);
- Scalar xcg = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+ Scalar xag = volVars.fluidState().moleFraction(gPhaseIdx, comp2Idx);
+ Scalar xcg = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
/* take care:
for xag in case wPhaseOnly we compute xag=henry_air*xaw
for xwg in case wPhaseOnly we compute xwg=pwsat
@@ -785,15 +786,15 @@ protected:
newPhasePresence = wgPhaseOnly;
globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
}
else if ((gasFlag == 0) && (NAPLFlag == 1))
{
newPhasePresence = wPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, aCompIdx);
+ = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(wPhaseIdx, comp1Idx);
}
}
else if (phasePresence == gPhaseOnly)
@@ -802,7 +803,7 @@ protected:
bool waterFlag = 0;
// calculate fractions in the hypothetical NAPL phase
- Scalar xnc = volVars.fluidState().moleFraction(nPhaseIdx, cCompIdx);
+ Scalar xnc = volVars.fluidState().moleFraction(nPhaseIdx, comp1Idx);
/*
take care:, xnc, if no NAPL phase is there, take xnc=xcg*pg/pcsat
if this is larger than 1, then NAPL appears
@@ -825,7 +826,7 @@ protected:
NAPLFlag = 1;
}
// calculate fractions of the hypothetical water phase
- Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, comp0Idx);
/*
take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
If this is larger than 1, then water appears
@@ -851,7 +852,7 @@ protected:
newPhasePresence = wgPhaseOnly;
globalSol[globalIdx][switch1Idx] = 0.0001;
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
}
else if ((waterFlag == 1) && (NAPLFlag == 1))
{
@@ -863,7 +864,7 @@ protected:
{
newPhasePresence = gnPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
globalSol[globalIdx][switch2Idx] = 0.0001;
}
}
@@ -874,7 +875,7 @@ protected:
bool waterFlag = 0;
// get the fractions in the hypothetical NAPL phase
- Scalar xnc = volVars.fluidState().moleFraction(nPhaseIdx, cCompIdx);
+ Scalar xnc = volVars.fluidState().moleFraction(nPhaseIdx, comp1Idx);
// take care: if the NAPL phase is not present, take
// xnc=xcg*pg/pcsat if this is larger than 1, then NAPL
@@ -928,7 +929,7 @@ protected:
{
newPhasePresence = gnPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
globalSol[globalIdx][switch2Idx] = 0.0001;
}
else if ((gasFlag == 0) && (NAPLFlag == 1) && (waterFlag == 0))
@@ -941,17 +942,17 @@ protected:
{
newPhasePresence = wPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, aCompIdx);
+ = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(wPhaseIdx, comp1Idx);
}
else if ((gasFlag == 0) && (NAPLFlag == 0) && (waterFlag == 1))
{
newPhasePresence = gPhaseOnly;
globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, cCompIdx);
+ = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
}
}
diff --git a/dumux/boxmodels/3p3c/3p3cvolumevariables.hh b/dumux/boxmodels/3p3c/3p3cvolumevariables.hh
index 95d516ab91..497e2e7be4 100644
--- a/dumux/boxmodels/3p3c/3p3cvolumevariables.hh
+++ b/dumux/boxmodels/3p3c/3p3cvolumevariables.hh
@@ -79,9 +79,9 @@ class ThreePThreeCVolumeVariables : public BoxVolumeVariables<TypeTag>
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
- wCompIdx = Indices::wCompIdx,
- aCompIdx = Indices::aCompIdx,
- cCompIdx = Indices::cCompIdx,
+ comp0Idx = Indices::comp0Idx,
+ comp2Idx = Indices::comp2Idx,
+ comp1Idx = Indices::comp1Idx,
wPhaseIdx = Indices::wPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx,
@@ -137,7 +137,7 @@ public:
// capillary pressure parameters
const MaterialLawParams &materialParams =
- problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
+ problem.spatialParams().materialLawParams(element, elemGeom, scvIdx);
int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
@@ -237,14 +237,14 @@ public:
// stored explicitly.
// extract mole fractions in the water phase
- Scalar xwa = primaryVars[switch1Idx];
- Scalar xwc = primaryVars[switch2Idx];
- Scalar xww = 1 - xwa - xwc;
+ Scalar xw2 = primaryVars[switch1Idx];
+ Scalar xw1 = primaryVars[switch2Idx];
+ Scalar xw0 = 1 - xw2 - xw1;
// write water mole fractions in the fluid state
- fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
- fluidState_.setMoleFraction(wPhaseIdx, aCompIdx, xwa);
- fluidState_.setMoleFraction(wPhaseIdx, cCompIdx, xwc);
+ fluidState_.setMoleFraction(wPhaseIdx, comp0Idx, xw0);
+ fluidState_.setMoleFraction(wPhaseIdx, comp2Idx, xw2);
+ fluidState_.setMoleFraction(wPhaseIdx, comp1Idx, xw1);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
@@ -264,16 +264,16 @@ public:
// we have all (partly hypothetical) phase pressures
// and temperature and the mole fraction of water in
// the gas phase
- Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx, cCompIdx)*pn_;
+ Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx, comp1Idx)*pn_;
- Scalar xgw = primaryVars[switch1Idx];
- Scalar xgc = partPressNAPL/pg_;
- Scalar xga = 1.-xgw-xgc;
+ Scalar xg0 = primaryVars[switch1Idx];
+ Scalar xg1 = partPressNAPL/pg_;
+ Scalar xg2 = 1.-xg0-xg1;
// write mole fractions in the fluid state
- fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
- fluidState_.setMoleFraction(gPhaseIdx, aCompIdx, xga);
- fluidState_.setMoleFraction(gPhaseIdx, cCompIdx, xgc);
+ fluidState_.setMoleFraction(gPhaseIdx, comp0Idx, xg0);
+ fluidState_.setMoleFraction(gPhaseIdx, comp2Idx, xg2);
+ fluidState_.setMoleFraction(gPhaseIdx, comp1Idx, xg1);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
@@ -286,17 +286,17 @@ public:
}
else if (phasePresence == wnPhaseOnly) {
// only water and NAPL phases are present
- Scalar pPartialC = fluidState_.fugacityCoefficient(nPhaseIdx,cCompIdx)*pn_;
- Scalar henryC = fluidState_.fugacityCoefficient(wPhaseIdx,cCompIdx)*pw_;
+ Scalar pPartialC = fluidState_.fugacityCoefficient(nPhaseIdx,comp1Idx)*pn_;
+ Scalar henryC = fluidState_.fugacityCoefficient(wPhaseIdx,comp1Idx)*pw_;
- Scalar xwa = primaryVars[switch1Idx];
- Scalar xwc = pPartialC/henryC;
- Scalar xww = 1.-xwa-xwc;
+ Scalar xw2 = primaryVars[switch1Idx];
+ Scalar xw1 = pPartialC/henryC;
+ Scalar xw0 = 1.-xw2-xw1;
// write mole fractions in the fluid state
- fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
- fluidState_.setMoleFraction(wPhaseIdx, aCompIdx, xwa);
- fluidState_.setMoleFraction(wPhaseIdx, cCompIdx, xwc);
+ fluidState_.setMoleFraction(wPhaseIdx, comp0Idx, xw0);
+ fluidState_.setMoleFraction(wPhaseIdx, comp2Idx, xw2);
+ fluidState_.setMoleFraction(wPhaseIdx, comp1Idx, xw1);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
@@ -311,14 +311,14 @@ public:
// only the gas phase is present, gas phase composition is
// stored explicitly here below.
- const Scalar xgw = primaryVars[switch1Idx];
- const Scalar xgc = primaryVars[switch2Idx];
- Scalar xga = 1 - xgw - xgc;
+ const Scalar xg0 = primaryVars[switch1Idx];
+ const Scalar xg1 = primaryVars[switch2Idx];
+ Scalar xg2 = 1 - xg0 - xg1;
// write mole fractions in the fluid state
- fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
- fluidState_.setMoleFraction(gPhaseIdx, aCompIdx, xga);
- fluidState_.setMoleFraction(gPhaseIdx, cCompIdx, xgc);
+ fluidState_.setMoleFraction(gPhaseIdx, comp0Idx, xg0);
+ fluidState_.setMoleFraction(gPhaseIdx, comp2Idx, xg2);
+ fluidState_.setMoleFraction(gPhaseIdx, comp1Idx, xg1);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
@@ -331,16 +331,16 @@ public:
}
else if (phasePresence == wgPhaseOnly) {
// only water and gas phases are present
- Scalar xgc = primaryVars[switch2Idx];
- Scalar partPressH2O = fluidState_.fugacityCoefficient(wPhaseIdx, wCompIdx)*pw_;
+ Scalar xg1 = primaryVars[switch2Idx];
+ Scalar partPressH2O = fluidState_.fugacityCoefficient(wPhaseIdx, comp0Idx)*pw_;
- Scalar xgw = partPressH2O/pg_;
- Scalar xga = 1.-xgc-xgw;
+ Scalar xg0 = partPressH2O/pg_;
+ Scalar xg2 = 1.-xg1-xg0;
// write mole fractions in the fluid state
- fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
- fluidState_.setMoleFraction(gPhaseIdx, aCompIdx, xga);
- fluidState_.setMoleFraction(gPhaseIdx, cCompIdx, xgc);
+ fluidState_.setMoleFraction(gPhaseIdx, comp0Idx, xg0);
+ fluidState_.setMoleFraction(gPhaseIdx, comp2Idx, xg2);
+ fluidState_.setMoleFraction(gPhaseIdx, comp1Idx, xg1);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
@@ -380,45 +380,45 @@ public:
*/
// diffusivity coefficents
- diffusionCoefficient_[gPhaseIdx][wCompIdx] =
+ diffusionCoefficient_[gPhaseIdx][comp0Idx] =
FluidSystem::diffusionCoefficient(fluidState_,
paramCache,
gPhaseIdx,
- wCompIdx);
- diffusionCoefficient_[gPhaseIdx][cCompIdx] =
+ comp0Idx);
+ diffusionCoefficient_[gPhaseIdx][comp1Idx] =
FluidSystem::diffusionCoefficient(fluidState_,
paramCache,
gPhaseIdx,
- cCompIdx);
- diffusionCoefficient_[gPhaseIdx][aCompIdx] = 0.0; // dummy, should not be used !
+ comp1Idx);
+ diffusionCoefficient_[gPhaseIdx][comp2Idx] = 0.0; // dummy, should not be used !
- diffusionCoefficient_[wPhaseIdx][aCompIdx] =
+ diffusionCoefficient_[wPhaseIdx][comp2Idx] =
FluidSystem::diffusionCoefficient(fluidState_,
paramCache,
wPhaseIdx,
- aCompIdx);
- diffusionCoefficient_[wPhaseIdx][cCompIdx] =
+ comp2Idx);
+ diffusionCoefficient_[wPhaseIdx][comp1Idx] =
FluidSystem::diffusionCoefficient(fluidState_,
paramCache,
wPhaseIdx,
- cCompIdx);
- diffusionCoefficient_[wPhaseIdx][wCompIdx] = 0.0; // dummy, should not be used !
+ comp1Idx);
+ diffusionCoefficient_[wPhaseIdx][comp0Idx] = 0.0; // dummy, should not be used !
/* no diffusion in NAPL phase considered at the moment */
- diffusionCoefficient_[nPhaseIdx][cCompIdx] = 0.0;
- diffusionCoefficient_[nPhaseIdx][wCompIdx] = 0.0;
- diffusionCoefficient_[nPhaseIdx][aCompIdx] = 0.0;
+ diffusionCoefficient_[nPhaseIdx][comp1Idx] = 0.0;
+ diffusionCoefficient_[nPhaseIdx][comp0Idx] = 0.0;
+ diffusionCoefficient_[nPhaseIdx][comp2Idx] = 0.0;
Valgrind::CheckDefined(diffusionCoefficient_);
// porosity
- porosity_ = problem.spatialParameters().porosity(element,
+ porosity_ = problem.spatialParams().porosity(element,
elemGeom,
scvIdx);
Valgrind::CheckDefined(porosity_);
// permeability
- permeability_ = problem.spatialParameters().intrinsicPermeability(element,
+ permeability_ = problem.spatialParams().intrinsicPermeability(element,
elemGeom,
scvIdx);
Valgrind::CheckDefined(permeability_);
--
GitLab