From 54b0c351be61201bc86dde4d1e6376dcf1aba9ba Mon Sep 17 00:00:00 2001
From: Benjamin Faigle <benjamin.faigle@posteo.de>
Date: Thu, 24 May 2012 15:19:28 +0000
Subject: [PATCH] -Naming conventions in transport module -Reformatted output
of pressure module -overhauled upwinding: Fully Upwind and central
implemented, fix for the case potential = 0. -Prepared stable modules for
merge of multiphyiscs and adaptive model: -Total Concentration (P.V.) is now
in cellData and not in FluidState -MuPhys: FluidState that is present locally
conceptually decoupled from subdomain.
Reviewed by Bernd
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8383 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/decoupled/2p2c/cellData2p2c.hh | 68 ++++--
.../2p2c/cellData2p2cmultiphysics.hh | 67 +++---
dumux/decoupled/2p2c/dec2p2cfluidstate.hh | 56 +----
dumux/decoupled/2p2c/fluxData2p2c.hh | 62 ++---
dumux/decoupled/2p2c/fvpressure2p2c.hh | 137 +++++++----
.../2p2c/fvpressure2p2cmultiphysics.hh | 199 +++++++++-------
.../decoupled/2p2c/fvpressurecompositional.hh | 92 ++++----
dumux/decoupled/2p2c/fvtransport2p2c.hh | 221 ++++++++++--------
dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh | 33 +--
9 files changed, 496 insertions(+), 439 deletions(-)
diff --git a/dumux/decoupled/2p2c/cellData2p2c.hh b/dumux/decoupled/2p2c/cellData2p2c.hh
index 21b66478c8..fb3aab44e2 100644
--- a/dumux/decoupled/2p2c/cellData2p2c.hh
+++ b/dumux/decoupled/2p2c/cellData2p2c.hh
@@ -72,6 +72,9 @@ private:
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
};
protected:
+ // primary variable (phase pressure has to be stored in fluidstate)
+ Scalar massConcentration_[numComponents];
+
Scalar mobility_[numPhases];
//numerical quantities
@@ -104,12 +107,12 @@ public:
globalIdx_ = 0;
perimeter_ = 0.;
}
-
+ //! Acess to flux data, representing information living on the intersections
FluxData& fluxData()
{
return fluxData_;
}
-
+ //! Constant acess to flux data, representing information living on the intersections
const FluxData& fluxData() const
{
return fluxData_;
@@ -134,33 +137,59 @@ public:
}
/*! Modify the phase pressure
* @param phaseIdx index of the Phase
- * @param value Value to be srored
+ * @param value Value to be stored
*/
void setPressure(int phaseIdx, Scalar value)
{
manipulateFluidState().setPressure(phaseIdx, value);
}
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::massConcentration()
+ /*!
+ * \brief Returns the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * This is equivalent to the sum of the component concentrations for all
+ * phases multiplied with the phase density.
+ *
+ * \param compIdx the index of the component
+ */
const Scalar totalConcentration(int compIdx) const
{
- return fluidState_->massConcentration(compIdx);
+ return massConcentration_[compIdx];
}
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::massConcentration()
+ //! \copydoc Dumux::CellData2P2C::totalConcentration()
const Scalar massConcentration(int compIdx) const
{
- return fluidState_->massConcentration(compIdx);
+ return massConcentration_[compIdx];
}
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::setMassConcentration()
+ /*!
+ * \brief Sets the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * @param compIdx index of the Component
+ * @param value Value to be stored
+ */
void setTotalConcentration(int compIdx, Scalar value)
{
- manipulateFluidState().setMassConcentration(compIdx, value);
+ massConcentration_[compIdx] = value;
}
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::setMassConcentration()
+ //! \copydoc Dumux::CellData2P2C::setTotalConcentration()
void setMassConcentration(int compIdx, Scalar value)
{
- manipulateFluidState().setMassConcentration(compIdx, value);
+ massConcentration_[compIdx] = value;
+ }
+ /*!
+ * \brief Calculate the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$
+ * for a given porosity (within the initialization procedure).
+ * @param porosity Porosity
+ */
+ void calculateMassConcentration(Scalar porosity)
+ {
+ massConcentration_[wCompIdx] =
+ porosity * (massFraction(wPhaseIdx,wCompIdx) * saturation(wPhaseIdx) * density(wPhaseIdx)
+ + massFraction(nPhaseIdx,wCompIdx) * saturation(nPhaseIdx) * density(nPhaseIdx));
+ massConcentration_[nCompIdx] =
+ porosity * (massFraction(wPhaseIdx,nCompIdx) * saturation(wPhaseIdx) * density(wPhaseIdx)
+ + massFraction(nPhaseIdx,nCompIdx) * saturation(nPhaseIdx) * density(nPhaseIdx));
}
//@}
@@ -282,7 +311,9 @@ public:
return fluidState_->viscosity(phaseIdx);
}
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::capillaryPressure()
+ /*!
+ * \brief Returns the capillary pressure \f$\mathrm{[Pa]}\f$
+ */
const Scalar capillaryPressure() const
{
return fluidState_->pressure(nPhaseIdx) - fluidState_->pressure(wPhaseIdx);
@@ -319,12 +350,12 @@ public:
}
//@}
-
-// void setFluidState(FluidState& fluidState)
-// DUNE_DEPRECATED
-// {
-// fluidState_ = &fluidState;
-// }
+ //! Deprecated set function, use manipulateFluidState() instead
+ void setFluidState(FluidState& fluidState)
+ DUNE_DEPRECATED
+ {
+ fluidState_ = &fluidState;
+ }
//! Returns a reference to the cells fluid state
const FluidState& fluidState() const
@@ -342,6 +373,7 @@ public:
fluidState_ = new FluidState;
return *fluidState_;
}
+
//! stores this cell datas index, only for debugging purposes!!
int& globalIdx()
{ return globalIdx_;}
diff --git a/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh b/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh
index 0ed1e5d62c..6ac0395369 100644
--- a/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh
@@ -113,30 +113,13 @@ public:
else
return this->fluidState_->pressure(phaseIdx);
}
-
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::massConcentration()
- const Scalar totalConcentration(int compIdx) const
+ //! \copydoc Dumux::CellData2P2C::setPressure()
+ void setPressure(int phaseIdx, Scalar value)
{
if(fluidStateType_ == simple)
- return simpleFluidState_->massConcentration(compIdx);
+ manipulateSimpleFluidState().setPressure(phaseIdx, value);
else
- return this->fluidState_->massConcentration(compIdx);
- }
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::massConcentration()
- const Scalar massConcentration(int compIdx) const
- {
- if(fluidStateType_ == simple)
- return simpleFluidState_->massConcentration(compIdx);
- else
- return this->fluidState_->massConcentration(compIdx);
- }
- //! \copydoc Dumux::DecoupledTwoPTwoCFluidState::setMassConcentration()
- void setTotalConcentration(int compIdx, Scalar value)
- {
- if(fluidStateType_ == simple)
- return simpleFluidState_->setMassConcentration(compIdx, value);
- else
- return this->fluidState_->setMassConcentration(compIdx, value);
+ manipulateFluidState().setPressure(phaseIdx, value);
}
//@}
@@ -159,14 +142,30 @@ public:
{
return subdomain_;
}
+ //! Specify subdomain information and fluidStateType
+ /** This function is only called if
+ */
+ void setSubdomainAndFluidStateType(int index)
+ {
+ subdomain_ = index;
+ if(index == 2)
+ fluidStateType_ = complex;
+ else
+ fluidStateType_ = simple;
+ }
/*! \name Acess to secondary variables */
//@{
//! \copydoc Dumux::DecoupledTwoPTwoCFluidState::setSaturation()
void setSaturation(int phaseIdx, Scalar value)
{
- assert(this->subdomain() == 2);
- this->fluidState_->setSaturation(phaseIdx, value);
+ if(fluidStateType_ == simple)
+ {
+ // saturation is triggered by presentPhaseIdx
+ manipulateSimpleFluidState().setPresentPhaseIdx((value==0.) ? nPhaseIdx : wPhaseIdx);
+ }
+ else
+ manipulateFluidState().setSaturation(phaseIdx, value);
}
//! \copydoc Dumux::DecoupledTwoPTwoCFluidState::saturation()
const Scalar saturation(int phaseIdx) const
@@ -185,10 +184,10 @@ public:
if(fluidStateType_ == simple)
{
assert(phaseIdx == simpleFluidState_->presentPhaseIdx());
- simpleFluidState_->setViscosity(phaseIdx, value);
+ manipulateSimpleFluidState().setViscosity(phaseIdx, value);
}
else
- this->fluidState_->setViscosity(phaseIdx, value);
+ manipulateFluidState().setViscosity(phaseIdx, value);
}
//! \copydoc Dumux::DecoupledTwoPTwoCFluidState::viscosity()
const Scalar viscosity(int phaseIdx) const
@@ -297,8 +296,13 @@ public:
//! Manipulates a simple fluidstate for a cell in the simple domain
SimpleFluidState& manipulateSimpleFluidState()
{
- assert (this->subdomain() != 2);
fluidStateType_ = simple;
+ if(this->fluidState_)
+ {
+ delete this->fluidState_;
+ this->fluidState_ = NULL;
+ }
+
if(!simpleFluidState_)
simpleFluidState_ = new SimpleFluidState;
return *simpleFluidState_;
@@ -310,13 +314,22 @@ public:
*/
FluidState& manipulateFluidState()
{
- assert(this->subdomain() == 2);
fluidStateType_ = complex;
+ if(simpleFluidState_)
+ {
+ delete simpleFluidState_;
+ simpleFluidState_ = NULL;
+ }
+
if(!this->fluidState_)
this->fluidState_ = new FluidState;
return *this->fluidState_;
}
+ //! Returns the type of the fluidState
+ const bool fluidStateType() const
+ { return fluidStateType_;}
+
};
}
#endif
diff --git a/dumux/decoupled/2p2c/dec2p2cfluidstate.hh b/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
index 46ef932ed1..a183381782 100644
--- a/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
@@ -119,7 +119,8 @@ public:
// check if there is enough of component 1 to form a phase
if (Z1 > massFraction_[nPhaseIdx][wCompIdx] && Z1 < massFraction_[wPhaseIdx][wCompIdx])
- nu_[nPhaseIdx] = -((equilRatio_[nPhaseIdx][wCompIdx]-1)*Z1 + (equilRatio_[nPhaseIdx][nCompIdx]-1)*(1-Z1)) / (equilRatio_[nPhaseIdx][wCompIdx]-1) / (equilRatio_[nPhaseIdx][nCompIdx] -1);
+ nu_[nPhaseIdx] = -((equilRatio_[nPhaseIdx][wCompIdx]-1)*Z1 + (equilRatio_[nPhaseIdx][nCompIdx]-1)*(1-Z1))
+ / (equilRatio_[nPhaseIdx][wCompIdx]-1) / (equilRatio_[nPhaseIdx][nCompIdx] -1);
else if (Z1 <= massFraction_[nPhaseIdx][wCompIdx]) // too little wComp to form a phase
{
nu_[nPhaseIdx] = 1; // only nPhase
@@ -130,9 +131,7 @@ public:
moleFraction_[nPhaseIdx][wCompIdx] /= ( massFraction_[nPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ massFraction_[nPhaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
- // corresponding phase
- massFraction_[wPhaseIdx][wCompIdx] = 1.;
- moleFraction_[wPhaseIdx][wCompIdx] = 1.;
+ // w phase is already set to equilibrium mass fraction
}
else // (Z1 >= Xw1) => no nPhase
{
@@ -144,8 +143,7 @@ public:
moleFraction_[wPhaseIdx][wCompIdx] /= ( massFraction_[wPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ massFraction_[wPhaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
- massFraction_[nPhaseIdx][wCompIdx] = 0.;
- moleFraction_[nPhaseIdx][wCompIdx] = 0.;
+ // n phase is already set to equilibrium mass fraction
}
// complete array of mass fractions
@@ -225,13 +223,6 @@ public:
// get densities with correct composition
density_[wPhaseIdx] = FluidSystem::density(*this, wPhaseIdx);
density_[nPhaseIdx] = FluidSystem::density(*this, nPhaseIdx);
-
- massConcentration_[wCompIdx] =
- poro * (massFraction_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
- + massFraction_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
- massConcentration_[nCompIdx] =
- poro * (massFraction_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
- + massFraction_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
}
//@}
/*!
@@ -273,43 +264,6 @@ public:
return moleFraction_[phaseIdx][compIdx];
}
- /*!
- * \brief Returns the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
- *
- * This is equivalent to the sum of the component concentrations for all
- * phases multiplied with the phase density.
- *
- * \param compIdx the index of the component
- */
- Scalar massConcentration(int compIdx) const
- {
- return massConcentration_[compIdx];
- };
- /*!
- * \brief Sets the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
- *
- * @param compIdx index of the Component
- * @param value Value to be stored
- */
- void setMassConcentration(int compIdx, Scalar value)
- {
- massConcentration_[compIdx] = value;
- };
- /*!
- * \brief Calculate the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$
- * for a given porosity (within the initialization procedure).
- * @param porosity Porosity
- */
- void calculateMassConcentration(Scalar porosity)
- {
- massConcentration_[wCompIdx] =
- porosity * (massFraction_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
- + massFraction_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
- massConcentration_[nCompIdx] =
- porosity * (massFraction_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
- + massFraction_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
- }
-
/*!
* \brief Returns the density of a phase \f$\mathrm{[kg/m^3]}\f$.
*
@@ -495,7 +449,7 @@ public:
DecoupledTwoPTwoCFluidState()
{ Valgrind::SetUndefined(*this); }
private:
- Scalar massConcentration_[numComponents];
+// Scalar massConcentration_[numComponents];
Scalar phasePressure_[numPhases];
Scalar temperature_;
diff --git a/dumux/decoupled/2p2c/fluxData2p2c.hh b/dumux/decoupled/2p2c/fluxData2p2c.hh
index 49fcdf339d..e36450b8c5 100644
--- a/dumux/decoupled/2p2c/fluxData2p2c.hh
+++ b/dumux/decoupled/2p2c/fluxData2p2c.hh
@@ -81,67 +81,35 @@ public:
isUpwindCell_[i] = false;
}
}
+ //! resizes the upwind vector for the case of hanging nodes
void resize(int size)
{
isUpwindCell_.resize(size);
}
//! functions returning upwind information
+ /* @param indexInInside The local inside index of the intersection
+ * @param equationIdx The equation index
+ */
const bool& isUpwindCell(int indexInInside, int equationIdx) const
{
return isUpwindCell_[indexInInside][equationIdx];
}
-
+ //! Sets the upwind information
+ /* @param indexInInside The local inside index of the intersection
+ * @param equationIdx The equation index
+ * @value value set true or false
+ */
void setUpwindCell(int indexInInside, int equationIdx, bool value)
{
isUpwindCell_[indexInInside][equationIdx] = value;
}
+ //! Console output for the FluxData
+ void outputFluxData()
+ {
+ for(int banana=0; banana<isUpwindCell_.size(); banana++)
+ printvector(std::cout, isUpwindCell_, "upwindInformation", "row", 3);
+ }
};
}
#endif
-
-
-/*** in transport module:
- * // upwind mobility
- double lambdaW, lambdaNW;
- if (potentialW >= 0.)
- {
- lambdaW = cellDataI.mobility(wPhaseIdx);
- cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, true);
- cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
- }
- else
- {
- lambdaW = cellDataJ.mobility(wPhaseIdx);
- cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, true);
- cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
- }
-
- if (potentialNW >= 0.)
- {
- lambdaNW = cellDataI.mobility(nPhaseIdx);
- cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, true);
- cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
- }
- else
- {
- lambdaW = cellDataJ.mobility(nPhaseIdx);
- cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, true);
- cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
- }
-
-
- in cellData:
-
- //! Acess to flux data
- bool& isUpwindCell(int indexInInside, int equationIdx) const
- {
- return fluxData_.isUpwindCell(indexInInside, equationIdx);
- }
-
- void setUpwindCell(int indexInInside, int equationIdx, bool value)
- {
- fluxData_.setUpwindCell(indexInInside, equationIdx, value);
- }
- */
-
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 1a331aebda..5612c0de8f 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -129,8 +129,8 @@ template<class TypeTag> class FVPressure2P2C
// convenience shortcuts for Vectors/Matrices
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
- typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
+ typedef Dune::FieldVector<Scalar, GET_PROP_VALUE(TypeTag, NumPhases)> PhaseVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
// the typenames used for the stiffness matrix and solution vector
@@ -212,9 +212,9 @@ private:
*/
template<class TypeTag>
void FVPressure2P2C<TypeTag>::getSource(Dune::FieldVector<Scalar, 2>& sourceEntry,
- const Element& elementI,
- const CellData& cellDataI,
- const bool first)
+ const Element& elementI,
+ const CellData& cellDataI,
+ const bool first)
{
sourceEntry=0.;
// cell volume & perimeter, assume linear map here
@@ -259,8 +259,8 @@ void FVPressure2P2C<TypeTag>::getSource(Dune::FieldVector<Scalar, 2>& sourceEntr
*/
template<class TypeTag>
void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEntry,
- const Element& elementI,
- const CellData& cellDataI,
+ const Element& elementI,
+ const CellData& cellDataI,
const bool first)
{
storageEntry = 0.;
@@ -366,7 +366,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
const GlobalPosition& gravity_ = problem().gravity();
// get absolute permeability
- FieldMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPtrI));
+ DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPtrI));
// get mobilities and fractional flow factors
Scalar fractionalWI=0, fractionalNWI=0;
@@ -401,11 +401,11 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
GlobalPosition unitDistVec(distVec);
unitDistVec /= dist;
- FieldMatrix permeabilityJ
+ DimMatrix permeabilityJ
= problem().spatialParams().intrinsicPermeability(*neighborPtr);
// compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
+ DimMatrix meanPermeability(0);
Dumux::harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
Dune::FieldVector<Scalar, dim> permeability(0);
@@ -479,47 +479,102 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
//do the upwinding of the mobility depending on the phase potentials
- if (potentialW >= 0.)
+ const CellData* upwindWCellData(0);
+ const CellData* upwindNCellData(0);
+ if (potentialW > 0.)
+ upwindWCellData = &cellDataI;
+ else if (potentialW < 0.)
+ upwindWCellData = &cellDataJ;
+ else
{
- dV_w = (dv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
- + dv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
- lambdaW = cellDataI.mobility(wPhaseIdx);
- gV_w = (graddv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
- + graddv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
- dV_w *= cellDataI.density(wPhaseIdx);
- gV_w *= cellDataI.density(wPhaseIdx);
+ if(cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
+ upwindWCellData = &cellDataI;
+ else if(cellDataJ.isUpwindCell(intersection.indexInOutside(), contiWEqIdx))
+ upwindWCellData = &cellDataJ;
+ //else
+ // upwinding is not done!
}
+
+ if (potentialNW > 0.)
+ upwindNCellData = &cellDataI;
+ else if (potentialNW < 0.)
+ upwindNCellData = &cellDataJ;
else
{
- dV_w = (dv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
+ if(cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
+ upwindNCellData = &cellDataI;
+ else if(cellDataJ.isUpwindCell(intersection.indexInOutside(), contiNEqIdx))
+ upwindNCellData = &cellDataJ;
+ //else
+ // upwinding is not done!
+ }
+
+ //perform upwinding if desired
+ if(!upwindWCellData)
+ {
+ // compute values for both cells
+ Scalar dV_wI = (dv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
+ + dv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
+ Scalar gV_wI = (graddv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
+ + graddv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
+ dV_wI *= cellDataI.density(wPhaseIdx);
+ gV_wI *= cellDataI.density(wPhaseIdx);
+ Scalar dV_wJ = (dv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
+ dv_dC2 * cellDataJ.massFraction(wPhaseIdx, nCompIdx));
- lambdaW = cellDataJ.mobility(wPhaseIdx);
- gV_w = (graddv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
+ Scalar gV_wJ = (graddv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
+ graddv_dC2 * cellDataJ.massFraction(wPhaseIdx, nCompIdx));
- dV_w *= cellDataJ.density(wPhaseIdx);
- gV_w *= cellDataJ.density(wPhaseIdx);
+ dV_wJ *= cellDataJ.density(wPhaseIdx);
+ gV_wJ *= cellDataJ.density(wPhaseIdx);
+
+ //compute means
+ dV_w = harmonicMean(dV_wI, dV_wJ);
+ gV_w = harmonicMean(gV_wI, gV_wJ);
+ lambdaW = harmonicMean(cellDataI.mobility(wPhaseIdx), cellDataJ.mobility(wPhaseIdx));
}
- if (potentialNW >= 0.)
+ else
+ {
+ dV_w = (dv_dC1 * upwindWCellData->massFraction(wPhaseIdx, wCompIdx)
+ + dv_dC2 * upwindWCellData->massFraction(wPhaseIdx, nCompIdx));
+ lambdaW = upwindWCellData->mobility(wPhaseIdx);
+ gV_w = (graddv_dC1 * upwindWCellData->massFraction(wPhaseIdx, wCompIdx)
+ + graddv_dC2 * upwindWCellData->massFraction(wPhaseIdx, nCompIdx));
+ dV_w *= upwindWCellData->density(wPhaseIdx);
+ gV_w *= upwindWCellData->density(wPhaseIdx);
+ }
+
+ if(!upwindNCellData)
{
- dV_n = (dv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
+ Scalar dV_nI = (dv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
+ dv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
- lambdaN = cellDataI.mobility(nPhaseIdx);
- gV_n = (graddv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
+ Scalar gV_nI = (graddv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
+ graddv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
- dV_n *= cellDataI.density(nPhaseIdx);
- gV_n *= cellDataI.density(nPhaseIdx);
+ dV_nI *= cellDataI.density(nPhaseIdx);
+ gV_nI *= cellDataI.density(nPhaseIdx);
+ Scalar dV_nJ = (dv_dC1 * cellDataJ.massFraction(nPhaseIdx, wCompIdx)
+ + dv_dC2 * cellDataJ.massFraction(nPhaseIdx, nCompIdx));
+ Scalar gV_nJ = (graddv_dC1 * cellDataJ.massFraction(nPhaseIdx, wCompIdx)
+ + graddv_dC2 * cellDataJ.massFraction(nPhaseIdx, nCompIdx));
+ dV_nJ *= cellDataJ.density(nPhaseIdx);
+ gV_nJ *= cellDataJ.density(nPhaseIdx);
+
+ //compute means
+ dV_n = harmonicMean(dV_nI, dV_nJ);
+ gV_n = harmonicMean(gV_nI, gV_nJ);
+ lambdaN = harmonicMean(cellDataI.mobility(nPhaseIdx), cellDataJ.mobility(nPhaseIdx));
}
else
{
- dV_n = (dv_dC1 * cellDataJ.massFraction(nPhaseIdx, wCompIdx)
- + dv_dC2 * cellDataJ.massFraction(nPhaseIdx, nCompIdx));
- lambdaN = cellDataJ.mobility(nPhaseIdx);
- gV_n = (graddv_dC1 * cellDataJ.massFraction(nPhaseIdx, wCompIdx)
- + graddv_dC2 * cellDataJ.massFraction(nPhaseIdx, nCompIdx));
- dV_n *= cellDataJ.density(nPhaseIdx);
- gV_n *= cellDataJ.density(nPhaseIdx);
+ dV_n = (dv_dC1 * upwindNCellData->massFraction(nPhaseIdx, wCompIdx)
+ + dv_dC2 * upwindNCellData->massFraction(nPhaseIdx, nCompIdx));
+ lambdaN = upwindNCellData->mobility(nPhaseIdx);
+ gV_n = (graddv_dC1 * upwindNCellData->massFraction(nPhaseIdx, wCompIdx)
+ + graddv_dC2 * upwindNCellData->massFraction(nPhaseIdx, nCompIdx));
+ dV_n *= upwindNCellData->density(nPhaseIdx);
+ gV_n *= upwindNCellData->density(nPhaseIdx);
}
+
+
//calculate current matrix entry
entries[matrix] = faceArea * (lambdaW * dV_w + lambdaN * dV_n)* (unitOuterNormal * unitDistVec);
if(enableVolumeIntegral)
@@ -613,7 +668,7 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
if (bcType.isDirichlet(Indices::pressureEqIdx))
{
// get absolute permeability
- FieldMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPtrI));
+ DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPtrI));
const GlobalPosition& gravity_ = problem().gravity();
//permeability vector at boundary
@@ -837,7 +892,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
CellData& cellData = problem().variables().cellData(globalIdx);
- asImp_().updateMaterialLawsInElement(*eIt);
+ problem().pressureModel().updateMaterialLawsInElement(*eIt);
maxError = std::max(maxError, fabs(cellData.volumeError()));
}
@@ -872,9 +927,9 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
// fluidState.setMassConcentration(nCompIdx,
// problem().transportModel().totalConcentration(nCompIdx,globalIdx));
// get the overall mass of component 1 Z1 = C^k / (C^1+C^2) [-]
- Scalar Z1 = fluidState.massConcentration(wCompIdx)
- / (fluidState.massConcentration(wCompIdx)
- + fluidState.massConcentration(nCompIdx));
+ Scalar Z1 = cellData.massConcentration(wCompIdx)
+ / (cellData.massConcentration(wCompIdx)
+ + cellData.massConcentration(nCompIdx));
// make shure only physical quantities enter flash calculation
#if DUNE_MINIMAL_DEBUG_LEVEL <= 3
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index e58691c1f7..74d80196a4 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -118,7 +118,7 @@ class FVPressure2P2CMultiPhysics : public FVPressure2P2C<TypeTag>
// convenience shortcuts for Vectors/Matrices
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
- typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef Dune::FieldVector<Scalar, GET_PROP_VALUE(TypeTag, NumPhases)> PhaseVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
@@ -157,6 +157,8 @@ public:
//constitutive functions are initialized and stored in the variables object
void updateMaterialLaws();
+ //updates singlephase secondary variables for one cell and stores in the variables object
+ void update1pMaterialLawsInElement(const Element&, CellData&);
//! \brief Write data files
/* \param name file name */
@@ -186,11 +188,9 @@ public:
gravity_(problem.gravity()), timer_(false)
{}
-private:
+protected:
// subdomain map
Dune::BlockVector<Dune::FieldVector<int,1> > nextSubdomain; //! vector holding next subdomain
-
-protected:
const GlobalPosition& gravity_; //!< vector including the gravity constant
static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
Dune::Timer timer_;
@@ -662,7 +662,6 @@ template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
{
//get timestep for error term
- Scalar dt = problem().timeManager().timeStepSize();
Scalar maxError = 0.;
// next subdomain map
@@ -671,14 +670,12 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
else
nextSubdomain = -1; // reduce complexity after first TS
- // iterate through leaf grid an evaluate c0 at cell center
+ // Loop A) through leaf grid
ElementIterator eItEnd = problem().gridView().template end<0> ();
for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
- const typename ElementIterator::Entity::Geometry &geo = eIt->geometry();
-
// get global coordinate of cell center
- GlobalPosition globalPos = geo.center();
+ GlobalPosition globalPos = eIt->geometry().center();
int globalIdx = problem().variables().index(*eIt);
CellData& cellData = problem().variables().cellData(globalIdx);
@@ -715,96 +712,130 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
}
timer_.stop();
// end subdomain check
- }
- else // simple
- {
- timer_.start();
- // determine which phase should be present
- int presentPhaseIdx = cellData.subdomain();
- // check if it is not yet assigned to next complex subdomain
- if(nextSubdomain[globalIdx] !=2)
- nextSubdomain[globalIdx] = presentPhaseIdx; // assign correct index
-
- // acess the simple fluid state and prepare for manipulation
- PseudoOnePTwoCFluidState<TypeTag>& pseudoFluidState
- = cellData.manipulateSimpleFluidState();
-
- // prepare phase pressure for fluid state
- // both phase pressures are necessary for the case 1p domain is assigned for
- // the next 2p subdomain
- PhaseVector pressure(0.);
- Scalar pc = MaterialLaw::pC(problem().spatialParams().materialLawParams(*eIt),
- ((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign Sw = 1 if wPhase present, else 0
- if(pressureType == wPhaseIdx)
- {
- pressure[wPhaseIdx] = this->pressure(globalIdx);
- pressure[nPhaseIdx] = this->pressure(globalIdx)+pc;
- }
- else
- {
- pressure[wPhaseIdx] = this->pressure(globalIdx)-pc;
- pressure[nPhaseIdx] = this->pressure(globalIdx);
- }
-
- // make shure total concentrations from solution vector are exact in fluidstate
- pseudoFluidState.setMassConcentration(wCompIdx,
- problem().transportModel().totalConcentration(wCompIdx,globalIdx));
- pseudoFluidState.setMassConcentration(nCompIdx,
- problem().transportModel().totalConcentration(nCompIdx,globalIdx));
-
- // get the overall mass of component 1: Z1 = C^k / (C^1+C^2) [-]
- Scalar sumConc = pseudoFluidState.massConcentration(wCompIdx)
- + pseudoFluidState.massConcentration(nCompIdx);
- Scalar Z1 = pseudoFluidState.massConcentration(wCompIdx)/ sumConc;
+ }// end complex domain
+ else if (nextSubdomain[globalIdx] != 2) //check if cell remains in simple subdomain
+ nextSubdomain[globalIdx] = cellData.subdomain();
- pseudoFluidState.update(Z1, pressure, cellData.subdomain(), problem().temperatureAtPos(globalPos));
+ } //end define complex area of next subdomain
- // write stuff in fluidstate
- assert(presentPhaseIdx == pseudoFluidState.presentPhaseIdx());
+ timer_.start();
+ // Loop B) thorugh leaf grid
+ // investigate cells that were "simple" in current TS
+ for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ int globalIdx = problem().variables().index(*eIt);
+ CellData& cellData = problem().variables().cellData(globalIdx);
- cellData.setSimpleFluidState(pseudoFluidState);
+ // store old subdomain information and assign new info
+ int oldSubdomainI = cellData.subdomain();
+ cellData.subdomain() = nextSubdomain[globalIdx];
+ //first check if simple will become complicated
+ if(oldSubdomainI != 2
+ && nextSubdomain[globalIdx] == 2)
+ {
+ // assure correct PV if subdomain used to be simple
+ if(cellData.fluidStateType() != 0); // i.e. it was simple
timer_.stop();
- // initialize viscosities
- cellData.setViscosity(presentPhaseIdx, FluidSystem::viscosity(pseudoFluidState, presentPhaseIdx));
-
- // initialize mobilities
- if(presentPhaseIdx == wPhaseIdx)
- cellData.setMobility(wPhaseIdx,
- MaterialLaw::krw(problem().spatialParams().materialLawParams(*eIt), pseudoFluidState.saturation(wPhaseIdx))
- / cellData.viscosity(wPhaseIdx));
- else
- cellData.setMobility(nPhaseIdx,
- MaterialLaw::krn(problem().spatialParams().materialLawParams(*eIt), pseudoFluidState.saturation(wPhaseIdx))
- / cellData.viscosity(nPhaseIdx));
-
- // error term handling
- Scalar vol(0.);
- vol = sumConc / pseudoFluidState.density(presentPhaseIdx);
-
- if (dt != 0)
- cellData.volumeError() = (vol - problem().spatialParams().porosity(*eIt));
-
+ this->updateMaterialLawsInElement(*eIt);
+ timer_.start();
+ }
+ else if(oldSubdomainI != 2
+ && nextSubdomain[globalIdx] != 2) // will be simple and was simple
+ {
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().center();
+// // determine which phase should be present
+// int presentPhaseIdx = cellData.subdomain(); // this is already =nextSubomainIdx
+ // perform simple update
+ this->update1pMaterialLawsInElement(*eIt, cellData);
+ timer_.stop();
}
+ //else
+ // a) will remain complex -> everything already done in loop A
+ // b) will be simple and was complex: complex FS available, so next TS
+ // will use comlex FS, next update will transform to simple
+
maxError = std::max(maxError, fabs(cellData.volumeError()));
}// end grid traversal
this->maxError_ = maxError/problem().timeManager().timeStepSize();
- // update subdomain information
- timer_.start();
- // iterate through leaf grid an copy subdomain information
- if(problem().timeManager().time() != 0.)
- {
- for (unsigned int i= 0; i < nextSubdomain.size(); i++)
- {
- problem().variables().cellData(i).subdomain() = nextSubdomain[i];
- }
- }
+
timer_.stop();
Dune::dinfo << "Subdomain routines took " << timer_.elapsed() << " seconds" << std::endl;
return;
}
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const Element& elementI, CellData& cellData)
+{
+ // get global coordinate of cell center
+ GlobalPosition globalPos = elementI.geometry().center();
+ int globalIdx = problem().variables().index(elementI);
+
+ // determine which phase should be present
+ int presentPhaseIdx = cellData.subdomain(); // this is already =nextSubomainIdx
+
+ // acess the simple fluid state and prepare for manipulation
+ PseudoOnePTwoCFluidState<TypeTag>& pseudoFluidState
+ = cellData.manipulateSimpleFluidState();
+
+ // prepare phase pressure for fluid state
+ // both phase pressures are necessary for the case 1p domain is assigned for
+ // the next 2p subdomain
+ PhaseVector pressure(0.);
+ Scalar pc = MaterialLaw::pC(problem().spatialParams().materialLawParams(elementI),
+ ((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign Sw = 1 if wPhase present, else 0
+ if(pressureType == wPhaseIdx)
+ {
+ pressure[wPhaseIdx] = this->pressure(globalIdx);
+ pressure[nPhaseIdx] = this->pressure(globalIdx)+pc;
+ }
+ else
+ {
+ pressure[wPhaseIdx] = this->pressure(globalIdx)-pc;
+ pressure[nPhaseIdx] = this->pressure(globalIdx);
+ }
+
+// // make shure total concentrations from solution vector are exact in fluidstate
+// pseudoFluidState.setMassConcentration(wCompIdx,
+// problem().transportModel().totalConcentration(wCompIdx,globalIdx));
+// pseudoFluidState.setMassConcentration(nCompIdx,
+// problem().transportModel().totalConcentration(nCompIdx,globalIdx));
+
+ // get the overall mass of component 1: Z1 = C^k / (C^1+C^2) [-]
+ Scalar sumConc = cellData.massConcentration(wCompIdx)
+ + cellData.massConcentration(nCompIdx);
+ Scalar Z1 = cellData.massConcentration(wCompIdx)/ sumConc;
+
+ pseudoFluidState.update(Z1, pressure, presentPhaseIdx, problem().temperatureAtPos(globalPos));
+
+ // write stuff in fluidstate
+ assert(presentPhaseIdx == pseudoFluidState.presentPhaseIdx());
+
+ cellData.setSimpleFluidState(pseudoFluidState);
+
+ // initialize viscosities
+ cellData.setViscosity(presentPhaseIdx, FluidSystem::viscosity(pseudoFluidState, presentPhaseIdx));
+
+ // initialize mobilities
+ if(presentPhaseIdx == wPhaseIdx)
+ cellData.setMobility(wPhaseIdx,
+ MaterialLaw::krw(problem().spatialParams().materialLawParams(elementI), pseudoFluidState.saturation(wPhaseIdx))
+ / cellData.viscosity(wPhaseIdx));
+ else
+ cellData.setMobility(nPhaseIdx,
+ MaterialLaw::krn(problem().spatialParams().materialLawParams(elementI), pseudoFluidState.saturation(wPhaseIdx))
+ / cellData.viscosity(nPhaseIdx));
+
+ // error term handling
+ Scalar vol(0.);
+ vol = sumConc / pseudoFluidState.density(presentPhaseIdx);
+
+ if (problem().timeManager().timeStepSize() != 0)
+ cellData.volumeError() = (vol - problem().spatialParams().porosity(elementI));
+ return;
+}
}//end namespace Dumux
#endif
diff --git a/dumux/decoupled/2p2c/fvpressurecompositional.hh b/dumux/decoupled/2p2c/fvpressurecompositional.hh
index c951b0ca7f..b1bd9a5f20 100644
--- a/dumux/decoupled/2p2c/fvpressurecompositional.hh
+++ b/dumux/decoupled/2p2c/fvpressurecompositional.hh
@@ -163,20 +163,25 @@ public:
ScalarSolutionType *pressureN = writer.allocateManagedBuffer(size);
ScalarSolutionType *pC = writer.allocateManagedBuffer(size);
ScalarSolutionType *saturationW = writer.allocateManagedBuffer(size);
-
ScalarSolutionType *densityWetting = writer.allocateManagedBuffer(size);
ScalarSolutionType *densityNonwetting = writer.allocateManagedBuffer(size);
- ScalarSolutionType *viscosityWetting = writer.allocateManagedBuffer(size);
- ScalarSolutionType *viscosityNonwetting = writer.allocateManagedBuffer(size);
ScalarSolutionType *mobilityW = writer.allocateManagedBuffer (size);
ScalarSolutionType *mobilityNW = writer.allocateManagedBuffer (size);
-
ScalarSolutionType *massfraction1W = writer.allocateManagedBuffer (size);
ScalarSolutionType *massfraction1NW = writer.allocateManagedBuffer (size);
-
// numerical stuff
ScalarSolutionType *volErr = writer.allocateManagedBuffer (size);
+ #if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ // add debug stuff
+ ScalarSolutionType *errorCorrPtr = writer.allocateManagedBuffer (size);
+ ScalarSolutionType *dv_dpPtr = writer.allocateManagedBuffer (size);
+ ScalarSolutionType *dV_dC1Ptr = writer.allocateManagedBuffer (size);
+ ScalarSolutionType *dV_dC2Ptr = writer.allocateManagedBuffer (size);
+ ScalarSolutionType *updEstimate1 = writer.allocateManagedBuffer (size);
+ ScalarSolutionType *updEstimate2 = writer.allocateManagedBuffer (size);
+ #endif
+
for (int i = 0; i < size; i++)
{
CellData& cellData = problem_.variables().cellData(i);
@@ -186,24 +191,27 @@ public:
(*saturationW)[i] = cellData.saturation(wPhaseIdx);
(*densityWetting)[i] = cellData.density(wPhaseIdx);
(*densityNonwetting)[i] = cellData.density(nPhaseIdx);
- (*viscosityWetting)[i] = cellData.viscosity(wPhaseIdx);
- (*viscosityNonwetting)[i] = cellData.viscosity(nPhaseIdx);
(*mobilityW)[i] = cellData.mobility(wPhaseIdx);
(*mobilityNW)[i] = cellData.mobility(nPhaseIdx);
(*massfraction1W)[i] = cellData.massFraction(wPhaseIdx,wCompIdx);
(*massfraction1NW)[i] = cellData.massFraction(nPhaseIdx,wCompIdx);
-
(*volErr)[i] = cellData.volumeError();
+
+ #if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ (*errorCorrPtr)[i] = cellData.errorCorrection();
+ (*dv_dpPtr)[i] = cellData.dv_dp();
+ (*dV_dC1Ptr)[i] = cellData.dv(wCompIdx);
+ (*dV_dC2Ptr)[i] = cellData.dv(nCompIdx);
+ (*updEstimate1)[i] = updateEstimate_[0][i];
+ (*updEstimate2)[i] = updateEstimate_[1][i];
+ #endif
}
writer.attachCellData(*pressureW, "wetting pressure");
writer.attachCellData(*pressureN, "nonwetting pressure");
writer.attachCellData(*pC, "capillary pressure");
writer.attachCellData(*saturationW, "wetting saturation");
-
writer.attachCellData(*densityWetting, "wetting density");
writer.attachCellData(*densityNonwetting, "nonwetting density");
- writer.attachCellData(*viscosityWetting, "wetting viscosity");
- writer.attachCellData(*viscosityNonwetting, "nonwetting viscosity");
writer.attachCellData(*mobilityW, "mobility w_phase");
writer.attachCellData(*mobilityNW, "mobility nw_phase");
std::ostringstream oss1, oss2;
@@ -212,44 +220,35 @@ public:
oss2 << "mass fraction " << FluidSystem::componentName(0) << " in " << FluidSystem::phaseName(1) << "-phase";
writer.attachCellData(*massfraction1NW, oss2.str());
writer.attachCellData(*volErr, "volume Error");
+ #if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ writer.attachCellData(*errorCorrPtr, "Error Correction");
+ writer.attachCellData(*dv_dpPtr, "dv_dp");
+ writer.attachCellData(*dV_dC1Ptr, "dV_dC1");
+ writer.attachCellData(*dV_dC2Ptr, "dV_dC2");
+ writer.attachCellData(*updEstimate1, "updEstimate comp 1");
+ writer.attachCellData(*updEstimate2, "updEstimate comp 2");
+ #endif
#if DUNE_MINIMAL_DEBUG_LEVEL <= 2
ScalarSolutionType *pressurePV = writer.allocateManagedBuffer(size);
ScalarSolutionType *totalConcentration1 = writer.allocateManagedBuffer (size);
ScalarSolutionType *totalConcentration2 = writer.allocateManagedBuffer (size);
-
- // add debug stuff
- ScalarSolutionType *errorCorrPtr = writer.allocateManagedBuffer (size);
- ScalarSolutionType *dv_dpPtr = writer.allocateManagedBuffer (size);
- ScalarSolutionType *dV_dC1Ptr = writer.allocateManagedBuffer (size);
- ScalarSolutionType *dV_dC2Ptr = writer.allocateManagedBuffer (size);
- ScalarSolutionType *updEstimate1 = writer.allocateManagedBuffer (size);
- ScalarSolutionType *updEstimate2 = writer.allocateManagedBuffer (size);
-
+ ScalarSolutionType *viscosityWetting = writer.allocateManagedBuffer(size);
+ ScalarSolutionType *viscosityNonwetting = writer.allocateManagedBuffer(size);
for (int i = 0; i < size; i++)
{
CellData& cellData = problem_.variables().cellData(i);
(*totalConcentration1)[i] = cellData.massConcentration(wCompIdx);
(*totalConcentration2)[i] = cellData.massConcentration(nCompIdx);
-
- (*errorCorrPtr)[i] = cellData.errorCorrection();
- (*dv_dpPtr)[i] = cellData.dv_dp();
- (*dV_dC1Ptr)[i] = cellData.dv(wCompIdx);
- (*dV_dC2Ptr)[i] = cellData.dv(nCompIdx);
- (*updEstimate1)[i] = updateEstimate_[0][i];
- (*updEstimate2)[i] = updateEstimate_[1][i];
+ (*viscosityWetting)[i] = cellData.viscosity(wPhaseIdx);
+ (*viscosityNonwetting)[i] = cellData.viscosity(nPhaseIdx);
}
*pressurePV = this->pressure();
writer.attachCellData(*pressurePV, "pressure (Primary Variable");
writer.attachCellData(*totalConcentration1, "C^w from cellData");
writer.attachCellData(*totalConcentration2, "C^n from cellData");
-
- writer.attachCellData(*errorCorrPtr, "Error Correction");
- writer.attachCellData(*dv_dpPtr, "dv_dp");
- writer.attachCellData(*dV_dC1Ptr, "dV_dC1");
- writer.attachCellData(*dV_dC2Ptr, "dV_dC2");
- writer.attachCellData(*updEstimate1, "updEstimate comp 1");
- writer.attachCellData(*updEstimate2, "updEstimate comp 2");
+ writer.attachCellData(*viscosityWetting, "wetting viscosity");
+ writer.attachCellData(*viscosityNonwetting, "nonwetting viscosity");
#endif
return;
}
@@ -323,13 +322,13 @@ protected:
Scalar ErrorTermUpperBound_; //!< Handling of error term: upper bound for error dampening
static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
private:
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this);}
-
- //! \copydoc Dumux::IMPETProblem::asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this);}
+// //! Returns the implementation of the problem (i.e. static polymorphism)
+// Implementation &asImp_()
+// { return *static_cast<Implementation *>(this);}
+//
+// //! \copydoc Dumux::IMPETProblem::asImp_()
+// const Implementation &asImp_() const
+// { return *static_cast<const Implementation *>(this);}
};
@@ -568,13 +567,12 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
= this->pressure()[globalIdx];
fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
-
- fluidState.calculateMassConcentration(problem_.spatialParams().porosity(*eIt));
-
} //end conc initial condition
+ cellData.calculateMassConcentration(problem_.spatialParams().porosity(*eIt));
+
} //end compositional
- problem_.transportModel().totalConcentration(wCompIdx,globalIdx) = fluidState.massConcentration(wCompIdx);
- problem_.transportModel().totalConcentration(nCompIdx,globalIdx) = fluidState.massConcentration(nCompIdx);
+ problem_.transportModel().totalConcentration(wCompIdx,globalIdx) = cellData.massConcentration(wCompIdx);
+ problem_.transportModel().totalConcentration(nCompIdx,globalIdx) = cellData.massConcentration(nCompIdx);
// initialize phase properties not stored in fluidstate
cellData.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, wPhaseIdx));
@@ -728,7 +726,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
DUNE_THROW(Dune::MathError, "NAN/inf of dV_dm. If that happens in first timestep, try smaller firstDt!");
}
}
- cellData.confirmVolumeDerivatives();
+ cellData.volumeDerivativesAvailable(true);
}
}//end namespace Dumux
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index 81ed17369e..9280795e03 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -98,7 +98,7 @@ class FVTransport2P2C
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix;
- typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef Dune::FieldVector<Scalar, GET_PROP_VALUE(TypeTag, NumPhases)> PhaseVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
//! Acess function for the current problem
@@ -321,7 +321,7 @@ void FVTransport2P2C<TypeTag>::updateTransportedQuantity(TransportSolutionType&
for(int compIdx = 0; compIdx < GET_PROP_VALUE(TypeTag, NumComponents); compIdx++)
{
totalConcentration_[compIdx][i] += (updateVector[compIdx][i]*=dt);
- cellDataI.setTotalConcentration(compIdx, totalConcentration_[compIdx][i]);
+ cellDataI.setMassConcentration(compIdx, totalConcentration_[compIdx][i]);
}
}
}
@@ -329,7 +329,7 @@ void FVTransport2P2C<TypeTag>::updateTransportedQuantity(TransportSolutionType&
//! Get flux at an interface between two cells
/** The flux thorugh \f$ \gamma{ij} \f$ is calculated according to the underlying pressure field,
* calculated by the pressure model.
- * \f[ - A_{\gamma_{ij}} \cdot \mathbf{K} \cdot \mathbf{u} \cdot (\mathbf{n}_{\gamma_{ij}} \cdot \mathbf{u})
+ * \f[ - A_{\gamma_{ij}} \cdot \mathbf{u} \cdot (\mathbf{n}_{\gamma_{ij}} \cdot \mathbf{u})\cdot \mathbf{K}
\sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \sum_{\kappa} X^{\kappa}_{\alpha}
\left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}\right) \f]
* Here, \f$ \mathbf{u} \f$ is the normalized vector connecting the cell centers, and \f$ \mathbf{n}_{\gamma_{ij}} \f$
@@ -384,7 +384,7 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
Dune::FieldVector<Scalar, 2> factor (0.);
Dune::FieldVector<Scalar, 2> updFactor (0.);
- Scalar potentialW(0.), potentialNW(0.);
+ PhaseVector potential(0.);
// access neighbor
ElementPointer neighborPtr = intersection.outside();
@@ -427,132 +427,151 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
{
case pw:
{
- potentialW = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
- potentialNW = (K * unitOuterNormal) * (pressI - pressJ + pcI - pcJ) / (dist);
+ potential[wPhaseIdx] = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
+ potential[nPhaseIdx] = (K * unitOuterNormal) * (pressI - pressJ + pcI - pcJ) / (dist);
break;
}
case pn:
{
- potentialW = (K * unitOuterNormal) * (pressI - pressJ - pcI + pcJ) / (dist);
- potentialNW = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
+ potential[wPhaseIdx] = (K * unitOuterNormal) * (pressI - pressJ - pcI + pcJ) / (dist);
+ potential[nPhaseIdx] = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
break;
}
}
// add gravity term
- potentialNW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;
- potentialW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
+ potential[nPhaseIdx] += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;
+ potential[wPhaseIdx] += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
- // upwind mobility
- double lambdaW(0.), lambdaNW(0.);
+ // determine upwinding direction, perform upwinding if possible
+ Dune::FieldVector<bool, GET_PROP_VALUE(TypeTag, NumPhases)> doUpwinding(true);
+ PhaseVector lambda(0.);
if(!impet_ or !restrictFluxInTransport_) // perform a simple uwpind scheme
{
- if (potentialW >= 0.)
+ if (potential[wPhaseIdx] > 0.)
{
- lambdaW = cellDataI.mobility(wPhaseIdx);
+ lambda[wPhaseIdx] = cellDataI.mobility(wPhaseIdx);
cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, true);
}
+ else if(potential[wPhaseIdx] < 0.)
+ {
+ lambda[wPhaseIdx] = cellDataJ.mobility(wPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
+ }
else
{
- lambdaW = cellDataJ.mobility(wPhaseIdx);
+ doUpwinding[wPhaseIdx] = false;
cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
}
- if (potentialNW >= 0.)
+ if (potential[nPhaseIdx] > 0.)
{
- lambdaNW = cellDataI.mobility(nPhaseIdx);
+ lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, true);
}
+ else if (potential[nPhaseIdx] < 0.)
+ {
+ lambda[nPhaseIdx] = cellDataJ.mobility(nPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
+ }
else
{
- lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ doUpwinding[nPhaseIdx] = false;
cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
}
}
else // if new potentials indicate same flow direction as in P.E., perform upwind
{
- if (potentialW >= 0. && cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
- lambdaW = cellDataI.mobility(wPhaseIdx);
- else if (potentialW < 0. && !cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
- lambdaW = cellDataJ.mobility(wPhaseIdx);
+ if (potential[wPhaseIdx] >= 0. && cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
+ lambda[wPhaseIdx] = cellDataI.mobility(wPhaseIdx);
+ else if (potential[wPhaseIdx] < 0. && !cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
+ lambda[wPhaseIdx] = cellDataJ.mobility(wPhaseIdx);
else // potential of w-phase does not coincide with that of P.E.
{
bool isUpwindCell = cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx);
//check if harmonic weithing is necessary
if (!isUpwindCell && !(cellDataI.mobility(wPhaseIdx) != 0. && cellDataJ.mobility(wPhaseIdx) == 0.)) // check if outflow induce neglected phase flux
- lambdaW = cellDataI.mobility(wPhaseIdx);
+ lambda[wPhaseIdx] = cellDataI.mobility(wPhaseIdx);
else if (isUpwindCell && !(cellDataJ.mobility(wPhaseIdx) != 0. && cellDataI.mobility(wPhaseIdx) == 0.)) // check if inflow induce neglected phase flux
- lambdaW = cellDataJ.mobility(wPhaseIdx);
+ lambda[wPhaseIdx] = cellDataJ.mobility(wPhaseIdx);
else
- {
- //a) perform harmonic averageing
- fluxEntries[wCompIdx] -= potentialW * faceArea / volume
- * harmonicMean(cellDataI.massFraction(wPhaseIdx, wCompIdx) * cellDataI.mobility(wPhaseIdx) * cellDataI.density(wPhaseIdx),
- cellDataJ.massFraction(wPhaseIdx, wCompIdx) * cellDataJ.mobility(wPhaseIdx) * cellDataJ.density(wPhaseIdx));
- fluxEntries[nCompIdx] -= potentialW * faceArea / volume
- * harmonicMean(cellDataI.massFraction(wPhaseIdx, nCompIdx) * cellDataI.mobility(wPhaseIdx) * cellDataI.density(wPhaseIdx),
- cellDataJ.massFraction(wPhaseIdx, nCompIdx) * cellDataJ.mobility(wPhaseIdx) * cellDataJ.density(wPhaseIdx));
- // b) timestep control
- // for timestep control : influx
- timestepFlux[0] += std::max(0.,
- - potentialW * faceArea / volume * harmonicMean(cellDataI.density(wPhaseIdx),cellDataJ.density(wPhaseIdx)));
- // outflux
- timestepFlux[1] += std::max(0.,
- potentialW * faceArea / volume * harmonicMean(cellDataI.density(wPhaseIdx),cellDataJ.density(wPhaseIdx)));
-
- //c) stop further standard calculations
- potentialW = 0;
-
- //d) output (only for one side)
- if(potentialW >= 0.)
- Dune::dinfo << "harmonicMean flux of phase w used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
- }
+ doUpwinding[wPhaseIdx] = false;
}
- if (potentialNW >= 0. && cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
- lambdaNW = cellDataI.mobility(nPhaseIdx);
- else if (potentialNW < 0. && !cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
- lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ if (potential[nPhaseIdx] >= 0. && cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
+ lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
+ else if (potential[nPhaseIdx] < 0. && !cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
+ lambda[nPhaseIdx] = cellDataJ.mobility(nPhaseIdx);
else // potential of n-phase does not coincide with that of P.E.
{
bool isUpwindCell = cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx);
//check if harmonic weithing is necessary
if (!isUpwindCell && !(cellDataI.mobility(nPhaseIdx) != 0. && cellDataJ.mobility(nPhaseIdx) == 0.)) // check if outflow induce neglected phase flux
- lambdaNW = cellDataI.mobility(nPhaseIdx);
+ lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
else if (isUpwindCell && !(cellDataJ.mobility(nPhaseIdx) != 0. && cellDataI.mobility(nPhaseIdx) == 0.)) // check if inflow induce neglected phase flux
- lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ lambda[nPhaseIdx] = cellDataJ.mobility(nPhaseIdx);
else
- {
- //a) perform harmonic averageing
- fluxEntries[wCompIdx] -= potentialNW * faceArea / volume
- * harmonicMean(cellDataI.massFraction(nPhaseIdx, wCompIdx) * cellDataI.mobility(nPhaseIdx) * cellDataI.density(nPhaseIdx),
- cellDataJ.massFraction(nPhaseIdx, wCompIdx) * cellDataJ.mobility(nPhaseIdx) * cellDataJ.density(nPhaseIdx));
- fluxEntries[nCompIdx] -= potentialNW * faceArea / volume
- * harmonicMean(cellDataI.massFraction(nPhaseIdx, nCompIdx) * cellDataI.mobility(nPhaseIdx) * cellDataI.density(nPhaseIdx),
- cellDataJ.massFraction(nPhaseIdx, nCompIdx) * cellDataJ.mobility(nPhaseIdx) * cellDataJ.density(nPhaseIdx));
- // b) timestep control
- // for timestep control : influx
- timestepFlux[0] += std::max(0.,
- - potentialNW * faceArea / volume * harmonicMean(cellDataI.density(nPhaseIdx),cellDataJ.density(nPhaseIdx)));
- // outflux
- timestepFlux[1] += std::max(0.,
- potentialNW * faceArea / volume * harmonicMean(cellDataI.density(nPhaseIdx),cellDataJ.density(nPhaseIdx)));
-
- //c) stop further standard calculations
- potentialNW = 0;
-
- //d) output (only for one side)
- if(potentialNW >= 0.)
- Dune::dinfo << "harmonicMean flux of phase n used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
- }
+ doUpwinding[nPhaseIdx] = false;
}
}
+ // do not perform upwinding if so desired
+ if(!doUpwinding[wPhaseIdx])
+ {
+ //a) perform harmonic averageing
+ fluxEntries[wCompIdx] -= potential[wPhaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(wPhaseIdx, wCompIdx) * cellDataI.mobility(wPhaseIdx) * cellDataI.density(wPhaseIdx),
+ cellDataJ.massFraction(wPhaseIdx, wCompIdx) * cellDataJ.mobility(wPhaseIdx) * cellDataJ.density(wPhaseIdx));
+ fluxEntries[nCompIdx] -= potential[wPhaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(wPhaseIdx, nCompIdx) * cellDataI.mobility(wPhaseIdx) * cellDataI.density(wPhaseIdx),
+ cellDataJ.massFraction(wPhaseIdx, nCompIdx) * cellDataJ.mobility(wPhaseIdx) * cellDataJ.density(wPhaseIdx));
+ // b) timestep control
+ // for timestep control : influx
+ timestepFlux[0] += std::max(0.,
+ - potential[wPhaseIdx] * faceArea / volume * harmonicMean(cellDataI.density(wPhaseIdx),cellDataJ.density(wPhaseIdx)));
+ // outflux
+ timestepFlux[1] += std::max(0.,
+ potential[wPhaseIdx] * faceArea / volume * harmonicMean(cellDataI.density(wPhaseIdx),cellDataJ.density(wPhaseIdx)));
+
+ //c) stop further standard calculations
+ potential[wPhaseIdx] = 0;
+
+ //d) output (only for one side)
+ if(potential[wPhaseIdx] >= 0.)
+ Dune::dinfo << "harmonicMean flux of phase w used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
+ }
+ if(!doUpwinding[nPhaseIdx])
+ {
+ //a) perform harmonic averageing
+ fluxEntries[wCompIdx] -= potential[nPhaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(nPhaseIdx, wCompIdx) * cellDataI.mobility(nPhaseIdx) * cellDataI.density(nPhaseIdx),
+ cellDataJ.massFraction(nPhaseIdx, wCompIdx) * cellDataJ.mobility(nPhaseIdx) * cellDataJ.density(nPhaseIdx));
+ fluxEntries[nCompIdx] -= potential[nPhaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(nPhaseIdx, nCompIdx) * cellDataI.mobility(nPhaseIdx) * cellDataI.density(nPhaseIdx),
+ cellDataJ.massFraction(nPhaseIdx, nCompIdx) * cellDataJ.mobility(nPhaseIdx) * cellDataJ.density(nPhaseIdx));
+ // b) timestep control
+ // for timestep control : influx
+ timestepFlux[0] += std::max(0.,
+ - potential[nPhaseIdx] * faceArea / volume * harmonicMean(cellDataI.density(nPhaseIdx),cellDataJ.density(nPhaseIdx)));
+ // outflux
+ timestepFlux[1] += std::max(0.,
+ potential[nPhaseIdx] * faceArea / volume * harmonicMean(cellDataI.density(nPhaseIdx),cellDataJ.density(nPhaseIdx)));
+
+ //c) stop further standard calculations
+ potential[nPhaseIdx] = 0;
+
+ //d) output (only for one side)
+ if(potential[nPhaseIdx] >= 0.)
+ Dune::dinfo << "harmonicMean flux of phase n used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
+ }
+
// calculate and standardized velocity
- double velocityJIw = std::max((-lambdaW * potentialW) * faceArea / volume, 0.0);
- double velocityIJw = std::max(( lambdaW * potentialW) * faceArea / volume, 0.0);
- double velocityJIn = std::max((-lambdaNW * potentialNW) * faceArea / volume, 0.0);
- double velocityIJn = std::max(( lambdaNW * potentialNW) * faceArea / volume, 0.0);
+ double velocityJIw = std::max((-lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea / volume, 0.0);
+ double velocityIJw = std::max(( lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea / volume, 0.0);
+ double velocityJIn = std::max((-lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea / volume, 0.0);
+ double velocityIJn = std::max(( lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea / volume, 0.0);
// for timestep control : influx
timestepFlux[0] += velocityJIw + velocityJIn;
@@ -579,7 +598,7 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
//! Get flux on Boundary
/** The flux thorugh \f$ \gamma{ij} \f$ is calculated according to the underlying pressure field,
* calculated by the pressure model.
- * \f[ - A_{\gamma_{ij}} \cdot \mathbf{K} \cdot \mathbf{u} \cdot (\mathbf{n}_{\gamma_{ij}} \cdot \mathbf{u})
+ * \f[ - A_{\gamma_{ij}} \cdot \mathbf{u} \cdot (\mathbf{n}_{\gamma_{ij}} \cdot \mathbf{u})\cdot \mathbf{K}
\sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \sum_{\kappa} X^{\kappa}_{\alpha}
\left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}\right) \f]
* Here, \f$ \mathbf{u} \f$ is the normalized vector connecting the cell centers, and \f$ \mathbf{n}_{\gamma_{ij}} \f$
@@ -631,7 +650,7 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& f
Dune::FieldVector<Scalar, 2> factor (0.);
Dune::FieldVector<Scalar, 2> updFactor (0.);
- Scalar potentialW(0.), potentialNW(0.);
+ PhaseVector potential(0.);
// center of face in global coordinates
const GlobalPosition& globalPosFace = intersection.geometry().center();
@@ -683,55 +702,55 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& f
{
case pw:
{
- potentialW = (K * unitOuterNormal) *
+ potential[wPhaseIdx] = (K * unitOuterNormal) *
(pressI - pressBound[wPhaseIdx]) / dist;
- potentialNW = (K * unitOuterNormal) *
+ potential[nPhaseIdx] = (K * unitOuterNormal) *
(pressI + pcI - pressBound[wPhaseIdx] - pcBound)
/ dist;
break;
}
case pn:
{
- potentialW = (K * unitOuterNormal) *
+ potential[wPhaseIdx] = (K * unitOuterNormal) *
(pressI - pcI - pressBound[nPhaseIdx] + pcBound)
/ dist;
- potentialNW = (K * unitOuterNormal) *
+ potential[nPhaseIdx] = (K * unitOuterNormal) *
(pressI - pressBound[nPhaseIdx]) / dist;
break;
}
}
- potentialW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;;
- potentialNW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;;
+ potential[wPhaseIdx] += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;;
+ potential[nPhaseIdx] += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;;
// do upwinding for lambdas
- double lambdaW, lambdaNW;
- if (potentialW >= 0.)
- lambdaW = cellDataI.mobility(wPhaseIdx);
+ PhaseVector lambda(0.);
+ if (potential[wPhaseIdx] >= 0.)
+ lambda[wPhaseIdx] = cellDataI.mobility(wPhaseIdx);
else
{
if(GET_PROP_VALUE(TypeTag, BoundaryMobility)==Indices::satDependent)
- lambdaW = BCfluidState.saturation(wPhaseIdx) / viscosityWBound;
+ lambda[wPhaseIdx] = BCfluidState.saturation(wPhaseIdx) / viscosityWBound;
else
- lambdaW = MaterialLaw::krw(
+ lambda[wPhaseIdx] = MaterialLaw::krw(
problem().spatialParams().materialLawParams(*elementPtrI), BCfluidState.saturation(wPhaseIdx))
/ viscosityWBound;
}
- if (potentialNW >= 0.)
- lambdaNW = cellDataI.mobility(nPhaseIdx);
+ if (potential[nPhaseIdx] >= 0.)
+ lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
else
{
if(GET_PROP_VALUE(TypeTag, BoundaryMobility)==Indices::satDependent)
- lambdaNW = BCfluidState.saturation(nPhaseIdx) / viscosityNWBound;
+ lambda[nPhaseIdx] = BCfluidState.saturation(nPhaseIdx) / viscosityNWBound;
else
- lambdaNW = MaterialLaw::krn(
+ lambda[nPhaseIdx] = MaterialLaw::krn(
problem().spatialParams().materialLawParams(*elementPtrI), BCfluidState.saturation(wPhaseIdx))
/ viscosityNWBound;
}
// calculate and standardized velocity
- double velocityJIw = std::max((-lambdaW * potentialW) * faceArea / volume, 0.0);
- double velocityIJw = std::max(( lambdaW * potentialW) * faceArea / volume, 0.0);
- double velocityJIn = std::max((-lambdaNW * potentialNW) * faceArea / volume, 0.0);
- double velocityIJn = std::max(( lambdaNW * potentialNW) * faceArea / volume, 0.0);
+ double velocityJIw = std::max((-lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea / volume, 0.0);
+ double velocityIJw = std::max(( lambda[wPhaseIdx] * potential[wPhaseIdx]) * faceArea / volume, 0.0);
+ double velocityJIn = std::max((-lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea / volume, 0.0);
+ double velocityIJn = std::max(( lambda[nPhaseIdx] * potential[nPhaseIdx]) * faceArea / volume, 0.0);
// for timestep control
timestepFlux[0] = velocityJIw + velocityJIn;
diff --git a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
index abf7852cfb..8d050b2c8d 100644
--- a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
@@ -82,7 +82,6 @@ public:
if (presentPhaseIdx == wPhaseIdx)
{
- Sw_ = 1.;
massFractionWater_[wPhaseIdx] = Z1;
massFractionWater_[nPhaseIdx] = 0.;
@@ -95,7 +94,6 @@ public:
}
else if (presentPhaseIdx == nPhaseIdx)
{
- Sw_ = 0.;
massFractionWater_[wPhaseIdx] = 0.;
massFractionWater_[nPhaseIdx] = Z1;
@@ -126,16 +124,20 @@ public:
*/
Scalar saturation(int phaseIdx) const
{
- if (phaseIdx == wPhaseIdx)
- return Sw_;
+ if (phaseIdx == presentPhaseIdx_)
+ return 1.;
else
- return Scalar(1.0) - Sw_;
+ return 0.;
};
int presentPhaseIdx() const
{
return presentPhaseIdx_;
}
+ void setPresentPhaseIdx(int index)
+ {
+ presentPhaseIdx_ = index;
+ };
/*! \brief Returns the pressure of a fluid phase \f$\mathrm{[Pa]}\f$.
* \param phaseIdx Index of the phase
@@ -197,8 +199,6 @@ public:
Scalar averageMolarMass(int phaseIdx) const
{
return aveMoMass_;
-// return moleFractionWater_[phaseIdx]*FluidSystem::molarMass(wCompIdx)
-// + (1.-moleFractionWater_[phaseIdx])*FluidSystem::molarMass(nCompIdx);
}
/*!
@@ -211,23 +211,11 @@ public:
{ return temperature_; };
/*!
- * \brief Returns the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
- *
- * This is equivalent to the sum of the component concentrations for all
- * phases multiplied with the phase density.
- *
- * \param compIdx the index of the component
- */
- Scalar massConcentration(int compIdx) const
- {
- return massConcentration_[compIdx];
- };
- /*!
- * \brief Sets the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
+ * \brief Sets the phase pressure \f$\mathrm{[Pa]}\f$.
*/
- void setMassConcentration(int compIdx, Scalar value)
+ void setPressure(int phaseIdx, Scalar value)
{
- massConcentration_[compIdx] = value;
+ pressure_[phaseIdx] = value;
};
//@}
@@ -237,7 +225,6 @@ public:
Scalar massFractionWater_[numPhases];
Scalar moleFractionWater_[numPhases];
Scalar pressure_[numPhases];
- Scalar Sw_;
Scalar density_;
Scalar viscosity_;
Scalar temperature_;
--
GitLab