diff --git a/dumux/decoupled/2p2c/2p2cfluidstate.hh b/dumux/decoupled/2p2c/2p2cfluidstate.hh
index 125641d02b9309c682ad2ad2f710e1c2ed5c3beb..60bef96508b3232136c2ee5365b7805e646fdc64 100644
--- a/dumux/decoupled/2p2c/2p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/2p2cfluidstate.hh
@@ -63,168 +63,6 @@ public:
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
public:
- /*!
- * \name flash calculation routines
- * Routines to determine the phase composition after the transport step.
- */
- //@{
-
- //! the update routine equals the 2p2c - concentration flash for constant p & t.
- /*!
- * Routine goes as follows:
- * - determination of the equilibrium constants from the fluid system
- * - determination of maximum solubilities (mole fractions) according to phase pressures
- * - comparison with Z1 to determine phase presence => phase mass fractions
- * - round off fluid properties
- * \param Z1 Feed mass fraction: Mass of comp1 per total mass \f$\mathrm{[-]}\f$
- * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$
- * \param poro Porosity \f$\mathrm{[-]}\f$
- * \param temperature Temperature \f$\mathrm{[K]}\f$
- */
- void update(Scalar Z1, Dune::FieldVector<Scalar, numPhases> phasePressure, Scalar poro, Scalar temperature)
- {
- if (pressureType == Indices::pressureGlobal)
- {
- DUNE_THROW(Dune::NotImplemented, "Pressure type not supported in fluidState!");
- }
- else
- phasePressure_[wPhaseIdx] = phasePressure[wPhaseIdx];
- phasePressure_[nPhaseIdx] = phasePressure[nPhaseIdx];
- temperature_=temperature;
-
-
- //mole equilibrium ratios K for in case wPhase is reference phase
- double k1 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, wCompIdx); // = p^wComp_vap / p
- double k2 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, nCompIdx); // = H^nComp_w / p
-
- // get mole fraction from equilibrium konstants
- moleFraction_[wPhaseIdx][wCompIdx] = (1. - k2) / (k1 -k2);
- moleFraction_[nPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * k1;
-
-
- // transform mole to mass fractions
- massFraction_[wPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
- / ( moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[wPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
- massFraction_[nPhaseIdx][wCompIdx] = moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
- / ( moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[nPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
-
-
- //mass equilibrium ratios
- equilRatio_[nPhaseIdx][wCompIdx] = massFraction_[nPhaseIdx][wCompIdx] / massFraction_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
- equilRatio_[nPhaseIdx][nCompIdx] = (1.-massFraction_[nPhaseIdx][wCompIdx])/ (1.-massFraction_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
- equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
-
- // phase fraction of nPhase [mass/totalmass]
- nu_[nPhaseIdx] = 0;
-
- // 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);
- else if (Z1 <= massFraction_[nPhaseIdx][wCompIdx]) // too little wComp to form a phase
- {
- nu_[nPhaseIdx] = 1; // only nPhase
- massFraction_[nPhaseIdx][wCompIdx] = Z1; // hence, assign complete mass soluted into nPhase
- massFraction_[nPhaseIdx][nCompIdx] = 1. - massFraction_[nPhaseIdx][wCompIdx];
- // store as moleFractions
- moleFraction_[nPhaseIdx][wCompIdx] = ( massFraction_[nPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx) ); // = moles of compIdx
- moleFraction_[nPhaseIdx][wCompIdx] /= ( massFraction_[nPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
- + massFraction_[nPhaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
-
- // w phase is already set to equilibrium mass fraction
- }
- else // (Z1 >= Xw1) => no nPhase
- {
- nu_[nPhaseIdx] = 0; // no second phase
- massFraction_[wPhaseIdx][wCompIdx] = Z1;
- massFraction_[wPhaseIdx][nCompIdx] = 1. - massFraction_[wPhaseIdx][wCompIdx];
- // store as moleFractions
- moleFraction_[wPhaseIdx][wCompIdx] = ( massFraction_[wPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx) ); // = moles of compIdx
- moleFraction_[wPhaseIdx][wCompIdx] /= ( massFraction_[wPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
- + massFraction_[wPhaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
-
- // n phase is already set to equilibrium mass fraction
- }
-
- // complete array of mass fractions
- massFraction_[wPhaseIdx][nCompIdx] = 1. - massFraction_[wPhaseIdx][wCompIdx];
- massFraction_[nPhaseIdx][nCompIdx] = 1. - massFraction_[nPhaseIdx][wCompIdx];
- // complete array of mole fractions
- moleFraction_[wPhaseIdx][nCompIdx] = 1. - moleFraction_[wPhaseIdx][wCompIdx];
- moleFraction_[nPhaseIdx][nCompIdx] = 1. - moleFraction_[nPhaseIdx][wCompIdx];
-
- // complete phase mass fractions
- nu_[wPhaseIdx] = 1. - nu_[nPhaseIdx];
-
- // get densities with correct composition
- density_[wPhaseIdx] = FluidSystem::density(*this, wPhaseIdx);
- density_[nPhaseIdx] = FluidSystem::density(*this, nPhaseIdx);
-
- Sw_ = (nu_[wPhaseIdx]) / density_[wPhaseIdx];
- Sw_ /= (nu_[wPhaseIdx]/density_[wPhaseIdx] + nu_[nPhaseIdx]/density_[nPhaseIdx]);
- }
-
- //! a flash routine for 2p2c if the saturation instead of total concentration is known.
- /*!
- * Routine goes as follows:
- * - determination of the equilibrium constants from the fluid system
- * - determination of maximum solubilities (mole fractions) according to phase pressures
- * - round off fluid properties
- * \param sat Saturation of phase 1 \f$\mathrm{[-]}\f$
- * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$
- * \param poro Porosity \f$\mathrm{[-]}\f$
- * \param temperature Temperature \f$\mathrm{[K]}\f$
- */
- void satFlash(Scalar sat, Dune::FieldVector<Scalar, numPhases> phasePressure, Scalar poro, Scalar temperature)
- {
- if (pressureType == Indices::pressureGlobal)
- {
- DUNE_THROW(Dune::NotImplemented, "Pressure type not supported in fluidState!");
- }
- else if (sat <= 0. || sat >= 1.)
- Dune::dinfo << "saturation initial and boundary conditions set to zero or one!"
- << " assuming fully saturated compositional conditions" << std::endl;
-
- // assign values
- Sw_ = sat;
- phasePressure_[wPhaseIdx] = phasePressure[wPhaseIdx];
- phasePressure_[nPhaseIdx] = phasePressure[nPhaseIdx];
- temperature_=temperature;
- nu_[nPhaseIdx] = nu_[wPhaseIdx] = NAN; //in contrast to the standard update() method, satflash() does not calculate nu.
-
-
- //mole equilibrium ratios K for in case wPhase is reference phase
- double k1 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, wCompIdx); // = p^wComp_vap / p
- double k2 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, nCompIdx); // = H^nComp_w / p
-
- // get mole fraction from equilibrium konstants
- moleFraction_[wPhaseIdx][wCompIdx] = (1. - k2) / (k1 -k2);
- moleFraction_[nPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * k1;
-
-
- // transform mole to mass fractions
- massFraction_[wPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
- / ( moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[wPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
- massFraction_[nPhaseIdx][wCompIdx] = moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
- / ( moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[nPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
-
- // complete array of mass fractions
- massFraction_[wPhaseIdx][nCompIdx] = 1. - massFraction_[wPhaseIdx][wCompIdx];
- massFraction_[nPhaseIdx][nCompIdx] = 1. - massFraction_[nPhaseIdx][wCompIdx];
- // complete array of mole fractions
- moleFraction_[wPhaseIdx][nCompIdx] = 1. - moleFraction_[wPhaseIdx][wCompIdx];
- moleFraction_[nPhaseIdx][nCompIdx] = 1. - moleFraction_[nPhaseIdx][wCompIdx];
-
- //mass equilibrium ratios
- equilRatio_[nPhaseIdx][wCompIdx] = massFraction_[nPhaseIdx][wCompIdx] / massFraction_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
- equilRatio_[nPhaseIdx][nCompIdx] = (1.-massFraction_[nPhaseIdx][wCompIdx])/ (1.-massFraction_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
- equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
-
- // get densities with correct composition
- density_[wPhaseIdx] = FluidSystem::density(*this, wPhaseIdx);
- density_[nPhaseIdx] = FluidSystem::density(*this, nPhaseIdx);
- }
- //@}
/*!
* \name acess functions
*/
@@ -296,6 +134,7 @@ public:
return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, wCompIdx);
else
DUNE_THROW(Dune::NotImplemented, "component not found in fluidState!");
+ return 0.;
}
/*!
@@ -355,6 +194,16 @@ public:
else
return nu_[phaseIdx];
}
+ /*!
+ * \brief Returns the phase mass fraction \f$ \nu \f$:
+ * phase mass per total mass \f$\mathrm{[kg/kg]}\f$.
+ *
+ * \param phaseIdx the index of the phase
+ */
+ Scalar& nu(int phaseIdx) const
+ {
+ return phaseMassFraction(phaseIdx);
+ }
/*!
* \name Functions to set Data
diff --git a/dumux/decoupled/2p2c/2p2cproblem.hh b/dumux/decoupled/2p2c/2p2cproblem.hh
index f20e1582ae6a894924f3caf3da892d1f110496aa..fc24a59f527c06058edf560227ff887f9ea4138d 100644
--- a/dumux/decoupled/2p2c/2p2cproblem.hh
+++ b/dumux/decoupled/2p2c/2p2cproblem.hh
@@ -112,7 +112,9 @@ public:
if (this->adaptiveGrid)
{
this->gridAdapt().adaptGrid();
- asImp_().pressureModel().updateMaterialLaws(false);
+
+ if(this->gridAdapt().wasAdapted())
+ asImp_().pressureModel().updateMaterialLaws(false);
}
}
/*!
@@ -203,9 +205,9 @@ protected:
}
else if (equation == -1)
{
- values[Indices::pressureEqIdx] = 0.;
- values[Indices::contiNEqIdx] =0.;
- values[Indices::contiWEqIdx] =0.;
+ // set everything to zero
+ for (int i = 0; i < values.size(); i++)
+ values[i] = 0.;
}
else
DUNE_THROW(Dune::InvalidStateException, "vector of primary variables can not be set properly");
diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
index 14dcd5e3a27cd8f06a65a34eed806f99d8f364f3..afbeb76cd3bc94acd504341ea1819672ab096949 100644
--- a/dumux/decoupled/2p2c/2p2cproperties.hh
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -89,9 +89,10 @@ NEW_PROP_TAG( FluidSystem ); //!< The fluid system
NEW_PROP_TAG( FluidState ); //!< The fluid state
NEW_PROP_TAG( ImpetEnableVolumeIntegral ); //!< Enables volume integral in the pressure equation (volume balance formulation)
NEW_PROP_TAG( EnableVolumeIntegral ); //!< DEPRECATED Enables volume integral in the pressure equation (volume balance formulation)
-NEW_PROP_TAG( EnableMultiPointFluxApproximationOnAdaptiveGrids ); //!< Two-point flux approximation (false) or mpfa (true)
+NEW_PROP_TAG( GridAdaptEnableMultiPointFluxApproximation); //!< HangingNode: Two-point flux approximation (false) or mpfa (true)
+NEW_PROP_TAG( EnableMultiPointFluxApproximationOnAdaptiveGrids ); //Deprecated
NEW_PROP_TAG( EnableSecondHalfEdge ); //!< Uses second interaction volume for second half-edge in 2D
-
+NEW_PROP_TAG( GridAdaptEnableSecondHalfEdge ); //!< Uses second interaction volume for second half-edge in 2D
}}
//DUMUX includes
@@ -155,8 +156,8 @@ SET_PROP(DecoupledTwoPTwoC, TransportSolutionType)
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCompressibility, true); //!< Compositional models are very likely compressible
SET_BOOL_PROP(DecoupledTwoPTwoC, VisitFacesOnlyOnce, false); //!< Faces are regarded from both sides
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCapillarity, false); //!< Capillarity is enabled
-SET_BOOL_PROP(DecoupledTwoPTwoC, ImpetRestrictFluxInTransport, GET_PROP_VALUE(TypeTag, RestrictFluxInTransport)); //!< Restrict (no upwind) flux in transport step if direction reverses after pressure equation
-SET_BOOL_PROP(DecoupledTwoPTwoC, RestrictFluxInTransport, false); //!< DEPRECATED Restrict (no upwind) flux in transport step if direction reverses after pressure equation
+SET_INT_PROP(DecoupledTwoPTwoC, ImpetRestrictFluxInTransport, GET_PROP_VALUE(TypeTag, RestrictFluxInTransport)); //!< Restrict (no upwind) flux in transport step if direction reverses after pressure equation
+SET_INT_PROP(DecoupledTwoPTwoC, RestrictFluxInTransport, 0); //!< DEPRECATED Restrict (no upwind) flux in transport step if direction reverses after pressure equation
SET_PROP(DecoupledTwoPTwoC, BoundaryMobility) //!< Saturation scales flux on Dirichlet B.C.
{ static const int value = DecoupledTwoPTwoCIndices<TypeTag>::satDependent;};
@@ -168,7 +169,9 @@ SET_TYPE_PROP(DecoupledTwoPTwoC, FluidState, TwoPTwoCFluidState<TypeTag>);
SET_TYPE_PROP(DecoupledTwoPTwoC, SpatialParameters, typename GET_PROP_TYPE(TypeTag, SpatialParams));//!< DEPRECATED SpatialParameters property
-SET_BOOL_PROP(DecoupledTwoPTwoC, EnableMultiPointFluxApproximationOnAdaptiveGrids, false); //!< MPFA disabled on adaptive grids
+SET_BOOL_PROP(DecoupledTwoPTwoC, GridAdaptEnableMultiPointFluxApproximation,
+ GET_PROP_VALUE(TypeTag,EnableMultiPointFluxApproximationOnAdaptiveGrids)); //!< MPFA disabled on adaptive grids
+SET_BOOL_PROP(DecoupledTwoPTwoC, EnableMultiPointFluxApproximationOnAdaptiveGrids, false); //!< DEPRECATED MPFA disabled on adaptive grids
SET_BOOL_PROP(DecoupledTwoPTwoC, ImpetEnableVolumeIntegral, GET_PROP_VALUE(TypeTag,EnableVolumeIntegral)); //!< Regard volume integral in pressure equation
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableVolumeIntegral, true); //!< DEPRECATED Regard volume integral in pressure equation
diff --git a/dumux/decoupled/2p2c/cellData2p2c.hh b/dumux/decoupled/2p2c/cellData2p2c.hh
index 28be1a2f56727ed1d6b7e1970711a81b05bd27ea..d248437103bbfa290137417fb60166bfb1b78bf2 100644
--- a/dumux/decoupled/2p2c/cellData2p2c.hh
+++ b/dumux/decoupled/2p2c/cellData2p2c.hh
@@ -83,6 +83,7 @@ protected:
Scalar dv_dp_;
Scalar dv_[numComponents];
bool volumeDerivativesAvailable_;
+ bool wasRefined_;
int globalIdx_;
Scalar perimeter_;
@@ -104,6 +105,7 @@ public:
errorCorrection_= 0.;
dv_dp_ = 0.;
volumeDerivativesAvailable_ = false;
+ wasRefined_ = false;
globalIdx_ = 0;
perimeter_ = 0.;
}
@@ -118,7 +120,6 @@ public:
return fluxData_;
}
-
/*! \name Acess to primary variables */
//@{
/*! Acess to the phase pressure
@@ -351,8 +352,7 @@ public:
//@}
//! Deprecated set function, use manipulateFluidState() instead
- void setFluidState(FluidState& fluidState)
- DUNE_DEPRECATED
+ void setFluidState(FluidState& fluidState) DUNE_DEPRECATED
{
fluidState_ = &fluidState;
}
@@ -389,9 +389,15 @@ public:
//! Resets the cell data after a timestep was completed: No volume derivatives yet available
void reset()
{
+ wasRefined_ = false;
volumeDerivativesAvailable_ = false;
}
-
+ //! Indicates if current cell was refined at this time step
+ bool& wasRefined()
+ { return wasRefined_;}
+ //!\copydoc wasRefined()
+ const bool& wasRefined() const
+ { return wasRefined_;}
//! Indicates if current cell is the upwind cell for a given interface
/* @param indexInInside Local face index seen from current cell
* @param phaseIdx The index of the phase
diff --git a/dumux/decoupled/2p2c/fluxData2p2c.hh b/dumux/decoupled/2p2c/fluxData2p2c.hh
index e36450b8c541379171a0c77050bf27b2a9bf8627..cd0346ece0bc7e4d11989d8a9b771fc56361acee 100644
--- a/dumux/decoupled/2p2c/fluxData2p2c.hh
+++ b/dumux/decoupled/2p2c/fluxData2p2c.hh
@@ -86,6 +86,11 @@ public:
{
isUpwindCell_.resize(size);
}
+ //! returns the size of the upwind vector which equals number of faces
+ int size()
+ {
+ return isUpwindCell_.size();
+ }
//! functions returning upwind information
/* @param indexInInside The local inside index of the intersection
* @param equationIdx The equation index
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 8c03597f90c8dfff63e6b7afe4f905c53adbd0f4..dc84219b050dc2ad3b68027857ed809c37835b5c 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -300,7 +300,6 @@ void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEn
// error reduction routine: volumetric error is damped and inserted to right hand side
// if damping is not done, the solution method gets unstable!
problem().variables().cellData(globalIdxI).volumeError() /= timestep_;
- Scalar maxError = maxError_;
Scalar erri = fabs(cellDataI.volumeError());
Scalar x_lo = ErrorTermLowerBound_;
Scalar x_mi = ErrorTermUpperBound_;
@@ -310,17 +309,17 @@ void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEn
Scalar lofac = 0.;
Scalar hifac = 0.;
- if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxError))
+ if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxError_))
{
- if (erri <= x_mi * maxError)
+ if (erri <= x_mi * maxError_)
storageEntry[rhs] +=
problem().variables().cellData(globalIdxI).errorCorrection() =
- fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxError)
+ fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxError_)
* cellDataI.volumeError() * volume;
else
storageEntry[rhs] +=
problem().variables().cellData(globalIdxI).errorCorrection() =
- fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxError)
+ fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxError_)
* cellDataI.volumeError() * volume;
}
else
@@ -510,28 +509,29 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
}
//perform upwinding if desired
- if(!upwindWCellData)
+ if(!upwindWCellData or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father()))
{
- // 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));
- Scalar gV_wJ = (graddv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
- + graddv_dC2 * cellDataJ.massFraction(wPhaseIdx, nCompIdx));
- dV_wJ *= cellDataJ.density(wPhaseIdx);
- gV_wJ *= cellDataJ.density(wPhaseIdx);
+ if (cellDataI.wasRefined() && cellDataJ.wasRefined())
+ {
+ problem().variables().cellData(problem().variables().index(*elementPtrI)).setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
+ }
+
+ Scalar averagedMassFraction[2];
+ averagedMassFraction[wCompIdx]
+ = harmonicMean(cellDataI.massFraction(wPhaseIdx, wCompIdx), cellDataJ.massFraction(wPhaseIdx, wCompIdx));
+ averagedMassFraction[nCompIdx]
+ = harmonicMean(cellDataI.massFraction(wPhaseIdx, nCompIdx), cellDataJ.massFraction(wPhaseIdx, nCompIdx));
+ Scalar averageDensity = harmonicMean(cellDataI.density(wPhaseIdx), cellDataJ.density(wPhaseIdx));
//compute means
- dV_w = harmonicMean(dV_wI, dV_wJ);
- gV_w = harmonicMean(gV_wI, gV_wJ);
+ dV_w = dv_dC1 * averagedMassFraction[wCompIdx] + dv_dC2 * averagedMassFraction[nCompIdx];
+ dV_w *= averageDensity;
+ gV_w = graddv_dC1 * averagedMassFraction[wCompIdx] + graddv_dC2 * averagedMassFraction[nCompIdx];
+ gV_w *= averageDensity;
lambdaW = harmonicMean(cellDataI.mobility(wPhaseIdx), cellDataJ.mobility(wPhaseIdx));
}
- else
+ else //perform upwinding
{
dV_w = (dv_dC1 * upwindWCellData->massFraction(wPhaseIdx, wCompIdx)
+ dv_dC2 * upwindWCellData->massFraction(wPhaseIdx, nCompIdx));
@@ -542,28 +542,29 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
gV_w *= upwindWCellData->density(wPhaseIdx);
}
- if(!upwindNCellData)
+ if(!upwindNCellData or (cellDataI.wasRefined() && cellDataJ.wasRefined()))
{
- Scalar dV_nI = (dv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
- + dv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
- Scalar gV_nI = (graddv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
- + graddv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
- 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);
+ if (cellDataI.wasRefined() && cellDataJ.wasRefined())
+ {
+ problem().variables().cellData(problem().variables().index(*elementPtrI)).setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
+ }
+ Scalar averagedMassFraction[2];
+ averagedMassFraction[wCompIdx]
+ = harmonicMean(cellDataI.massFraction(nPhaseIdx, wCompIdx), cellDataJ.massFraction(nPhaseIdx, wCompIdx));
+ averagedMassFraction[nCompIdx]
+ = harmonicMean(cellDataI.massFraction(nPhaseIdx, nCompIdx), cellDataJ.massFraction(nPhaseIdx, nCompIdx));
+ Scalar averageDensity = harmonicMean(cellDataI.density(nPhaseIdx), cellDataJ.density(nPhaseIdx));
//compute means
- dV_n = harmonicMean(dV_nI, dV_nJ);
- gV_n = harmonicMean(gV_nI, gV_nJ);
+ dV_n = dv_dC1 * averagedMassFraction[wCompIdx] + dv_dC2 * averagedMassFraction[nCompIdx];
+ dV_n *= averageDensity;
+ gV_n = graddv_dC1 * averagedMassFraction[wCompIdx] + graddv_dC2 * averagedMassFraction[nCompIdx];
+ gV_n *= averageDensity;
lambdaN = harmonicMean(cellDataI.mobility(nPhaseIdx), cellDataJ.mobility(nPhaseIdx));
}
else
- {
+ {
dV_n = (dv_dC1 * upwindNCellData->massFraction(nPhaseIdx, wCompIdx)
+ dv_dC2 * upwindNCellData->massFraction(nPhaseIdx, nCompIdx));
lambdaN = upwindNCellData->mobility(nPhaseIdx);
@@ -573,8 +574,6 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
gV_n *= upwindNCellData->density(nPhaseIdx);
}
-
-
//calculate current matrix entry
entries[matrix] = faceArea * (lambdaW * dV_w + lambdaN * dV_n)* (unitOuterNormal * unitDistVec);
if(enableVolumeIntegral)
@@ -983,7 +982,8 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
}
//complete fluid state
- fluidState.update(Z1, pressure, problem().spatialParams().porosity(elementI), temperature_);
+ CompositionalFlash<TypeTag> flashSolver;
+ flashSolver.concentrationFlash2p2c(fluidState, Z1, pressure, problem().spatialParams().porosity(elementI), temperature_);
// iterations part in case of enabled capillary pressure
Scalar pc(0.), oldPc(0.);
@@ -1015,8 +1015,8 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
//store old pc
oldPc = pc;
//update with better pressures
- fluidState.update(Z1, pressure, problem().spatialParams().porosity(elementI),
- problem().temperatureAtPos(globalPos));
+ flashSolver.concentrationFlash2p2c(fluidState, Z1, pressure,
+ problem().spatialParams().porosity(elementI), problem().temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem().spatialParams().materialLawParams(elementI),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index d70d7cd5719b658946a2089cb73f464ffde30653..77ce4fb3adfc947c5920ecb5080083f0ffb95483 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -208,7 +208,7 @@ protected:
};
//! function which assembles the system of equations to be solved
-/** for first == true, this function assembles the matrix and right hand side for
+/*! for first == true, this function assembles the matrix and right hand side for
* the solution of the pressure field in the same way as in the class FVPressure2P.
* for first == false, the approach is changed to \f[-\frac{\partial V}{\partial p}
* \frac{\partial p}{\partial t}+\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}\nabla\cdot
@@ -350,10 +350,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
Scalar Z1 = cellDataI.massConcentration(wCompIdx) / sumC;
// initialize simple fluidstate object
PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
- pseudoFluidState.update(Z1, p_, cellDataI.subdomain(),
+ CompositionalFlash<TypeTag> flashSolver;
+ flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
problem().temperatureAtPos(elementI.geometry().center()));
- Scalar v_ = 1. / FluidSystem::density(pseudoFluidState, presentPhaseIdx);
+ Scalar v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
cellDataI.dv_dp() = (sumC * ( v_ - (1. /cellDataI.density(presentPhaseIdx)))) /incp;
if (cellDataI.dv_dp()>0)
@@ -362,9 +363,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
Dune::dinfo << "dv_dp larger 0 at Idx " << globalIdxI << " , try and invert secant"<< std::endl;
p_ -= 2*incp;
- pseudoFluidState.update(Z1, p_, cellDataI.subdomain(),
+ flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
problem().temperatureAtPos(elementI.geometry().center()));
- v_ = 1. / FluidSystem::density(pseudoFluidState, presentPhaseIdx);
+ v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
cellDataI.dv_dp() = (sumC * ( v_ - (1. /cellDataI.density(presentPhaseIdx)))) /incp;
// dV_dp > 0 is unphysical: Try inverse increment for secant
if (cellDataI.dv_dp()>0)
@@ -479,6 +480,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>
// due to "safety cell" around subdomain, both cells I and J
// have single-phase conditions, although one is in 2p domain.
int phaseIdx = std::min(cellDataI.subdomain(), cellDataJ.subdomain());
+ // determine respective equation Idx
+ int contiEqIdx = (phaseIdx == wPhaseIdx) ? contiWEqIdx : contiNEqIdx;
Scalar rhoMean = 0.5 * (cellDataI.density(phaseIdx) + cellDataJ.density(phaseIdx));
//Scalar density = 0;
@@ -490,14 +493,22 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>
Scalar lambda;
- if (potential >= 0.)
+ if (potential > 0.)
{
lambda = cellDataI.mobility(phaseIdx);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiEqIdx, false); // store in cellJ since cellI is const
//density = cellDataI.density(phaseIdx);
}
- else
+ else if (potential < 0.)
{
lambda = cellDataJ.mobility(phaseIdx);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiEqIdx, true);
+ //density = cellDataJ.density(phaseIdx);
+ }
+ else
+ {
+ lambda = harmonicMean(cellDataI.mobility(phaseIdx) , cellDataJ.mobility(phaseIdx));
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiEqIdx, false);
//density = cellDataJ.density(phaseIdx);
}
@@ -732,23 +743,16 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
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();
- }
+ // use complex update of the fluidstate
+ timer_.stop();
this->updateMaterialLawsInElement(*eIt, postTimeStep);
timer_.start();
}
else if(oldSubdomainI != 2
&& nextSubdomain[globalIdx] != 2) // will be simple and was simple
{
-// // determine which phase should be present
-// int presentPhaseIdx = cellData.subdomain(); // this is already =nextSubomainIdx
-
- // perform simple update
+ // perform simple update
this->update1pMaterialLawsInElement(*eIt, cellData, postTimeStep);
- timer_.stop();
}
//else
// a) will remain complex -> everything already done in loop A
@@ -760,7 +764,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
this->maxError_ = maxError/problem().timeManager().timeStepSize();
timer_.stop();
- Dune::dinfo << "Subdomain routines took " << timer_.elapsed() << " seconds" << std::endl;
+
+ if(problem().timeManager().willBeFinished() or problem().timeManager().episodeWillBeOver())
+ Dune::dinfo << "Subdomain routines took " << timer_.elapsed() << " seconds" << std::endl;
return;
}
@@ -799,18 +805,13 @@ void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const El
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));
+ CompositionalFlash<TypeTag> flashSolver;
+ flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, pressure, presentPhaseIdx, problem().temperatureAtPos(globalPos));
// write stuff in fluidstate
assert(presentPhaseIdx == pseudoFluidState.presentPhaseIdx());
diff --git a/dumux/decoupled/2p2c/fvpressurecompositional.hh b/dumux/decoupled/2p2c/fvpressurecompositional.hh
index bb0826849ad6cf052b2bca333ad690d3027886b7..24c11c4082c83a12122e6163378f89d2b99b5f1b 100644
--- a/dumux/decoupled/2p2c/fvpressurecompositional.hh
+++ b/dumux/decoupled/2p2c/fvpressurecompositional.hh
@@ -27,6 +27,7 @@
// dumux environment
#include "dumux/common/math.hh"
#include <dumux/decoupled/common/fv/fvpressure.hh>
+#include <dumux/material/constraintsolvers/compositionalflash.hh>
#include <dumux/decoupled/2p2c/2p2cproperties.hh>
#include <dumux/io/vtkmultiwriter.hh>
@@ -235,6 +236,10 @@ public:
ScalarSolutionType *totalConcentration2 = writer.allocateManagedBuffer (size);
ScalarSolutionType *viscosityWetting = writer.allocateManagedBuffer(size);
ScalarSolutionType *viscosityNonwetting = writer.allocateManagedBuffer(size);
+// ScalarSolutionType *nun = writer.allocateManagedBuffer(size);
+// ScalarSolutionType *nuw = writer.allocateManagedBuffer(size);
+ ScalarSolutionType *faceUpwindW = writer.allocateManagedBuffer(size);
+ ScalarSolutionType *faceUpwindN = writer.allocateManagedBuffer(size);
for (int i = 0; i < size; i++)
{
CellData& cellData = problem_.variables().cellData(i);
@@ -242,8 +247,25 @@ public:
(*totalConcentration2)[i] = cellData.massConcentration(nCompIdx);
(*viscosityWetting)[i] = cellData.viscosity(wPhaseIdx);
(*viscosityNonwetting)[i] = cellData.viscosity(nPhaseIdx);
+// (*nun)[i] = cellData.phaseMassFraction(nPhaseIdx);
+// (*nuw)[i] = cellData.phaseMassFraction(wPhaseIdx);
+ (*faceUpwindW)[i] = 0;
+ (*faceUpwindN)[i] = 0;
+ // run thorugh all local face idx and collect upwind information
+ for(int faceIdx = 0; faceIdx<cellData.fluxData().size(); faceIdx++)
+ {
+ if(cellData.isUpwindCell(faceIdx, contiWEqIdx))
+ (*faceUpwindW)[i] += pow(10,static_cast<double>(1-faceIdx));
+ if(cellData.isUpwindCell(faceIdx, contiNEqIdx))
+ (*faceUpwindN)[i] += pow(10,static_cast<double>(1-faceIdx));
+ }
}
+// writer.attachCellData(*nun, "phase mass fraction n-phase");
+// writer.attachCellData(*nuw, "phase mass fraction w-phase");
*pressurePV = this->pressure();
+
+ writer.attachCellData(*faceUpwindW, "isUpwind w-phase");
+ writer.attachCellData(*faceUpwindN, "isUpwind n-phase");
writer.attachCellData(*pressurePV, "pressure (Primary Variable");
writer.attachCellData(*totalConcentration1, "C^w from cellData");
writer.attachCellData(*totalConcentration2, "C^n from cellData");
@@ -454,6 +476,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
CellData& cellData = problem_.variables().cellData(globalIdx);
// acess the fluid state and prepare for manipulation
FluidState& fluidState = cellData.manipulateFluidState();
+ CompositionalFlash<TypeTag> flashSolver;
// initial conditions
PhaseVector pressure(0.);
@@ -470,12 +493,13 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
if (icFormulation == Indices::saturation) // saturation initial condition
{
sat_0 = problem_.initSat(*eIt);
- fluidState.satFlash(sat_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
+ flashSolver.saturationFlash2p2c(fluidState, sat_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
else if (icFormulation == Indices::concentration) // concentration initial condition
{
Scalar Z1_0 = problem_.initConcentration(*eIt);
- fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
+ flashSolver.concentrationFlash2p2c(fluidState, Z1_0, pressure,
+ problem_.spatialParams().porosity(*eIt), temperature_);
}
}
else if(compositional) //means we regard compositional effects since we know an estimate pressure field
@@ -509,8 +533,8 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
break;
}
}
-
- fluidState.satFlash(sat_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
+ flashSolver.saturationFlash2p2c(fluidState, sat_0, pressure,
+ problem_.spatialParams().porosity(*eIt), temperature_);
}
else if (icFormulation == Indices::concentration) // concentration initial condition
{
@@ -548,8 +572,8 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
//store old pc
Scalar oldPc = pc;
//update with better pressures
- fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt),
- problem_.temperatureAtPos(globalPos));
+ flashSolver.concentrationFlash2p2c(fluidState, Z1_0, pressure,
+ problem_.spatialParams().porosity(*eIt), problem_.temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
@@ -565,7 +589,8 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
{
pressure[wPhaseIdx] = pressure[nPhaseIdx]
= this->pressure()[globalIdx];
- fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
+ flashSolver.concentrationFlash2p2c(fluidState, Z1_0,
+ pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
} //end conc initial condition
cellData.calculateMassConcentration(problem_.spatialParams().porosity(*eIt));
@@ -629,8 +654,9 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
// get cell temperature
Scalar temperature_ = problem_.temperatureAtPos(globalPos);
- // initialize an Fluid state for the update
+ // initialize an Fluid state and a flash solver
FluidState updFluidState;
+ CompositionalFlash<TypeTag> flashSolver;
/**********************************
* a) get necessary variables
@@ -667,7 +693,6 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
}
Scalar incp = 1e-2;
-
/**********************************
* c) Secant method for derivatives
**********************************/
@@ -676,7 +701,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
PhaseVector p_(incp);
p_ += pressure;
Scalar Z1 = mass[0] / mass.one_norm();
- updFluidState.update(Z1,
+ flashSolver.concentrationFlash2p2c(updFluidState, Z1,
p_, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
@@ -691,7 +716,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
Dune::dinfo << "dv_dp larger 0 at Idx " << globalIdx << " , try and invert secant"<< std::endl;
p_ -= 2*incp;
- updFluidState.update(Z1,
+ flashSolver.concentrationFlash2p2c(updFluidState, Z1,
p_, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
@@ -711,7 +736,9 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
{
mass[compIdx] += massIncrement[compIdx];
Z1 = mass[0] / mass.one_norm();
- updFluidState.update(Z1, pressure, problem_.spatialParams().porosity(element), temperature_);
+
+ flashSolver.concentrationFlash2p2c(updFluidState, Z1,
+ pressure, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
for(int phaseIdx = 0; phaseIdx< numPhases; phaseIdx++)
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index b97ed8d9b63a34b8d46d57cf248200f5dac7f497..2e3158b60bd27bf5602176d5f1227d5c099d97eb 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -73,7 +73,7 @@ class FVTransport2P2C
enum
{
- dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld,
};
enum
{
@@ -86,7 +86,8 @@ class FVTransport2P2C
{
wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
- contiWEqIdx=Indices::contiWEqIdx, contiNEqIdx=Indices::contiNEqIdx
+ contiWEqIdx=Indices::contiWEqIdx, contiNEqIdx=Indices::contiNEqIdx,
+ NumPhases = GET_PROP_VALUE(TypeTag, NumPhases)
};
typedef typename GridView::Traits::template Codim<0>::Entity Element;
@@ -98,7 +99,7 @@ class FVTransport2P2C
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix;
- typedef Dune::FieldVector<Scalar, GET_PROP_VALUE(TypeTag, NumPhases)> PhaseVector;
+ typedef Dune::FieldVector<Scalar, NumPhases> PhaseVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
//! Acess function for the current problem
@@ -173,6 +174,11 @@ public:
//! \copydoc transportedQuantity()
void getTransportedQuantity(TransportSolutionType& transportedQuantity)
{
+ // resize update vector and set to zero
+ transportedQuantity.resize(GET_PROP_VALUE(TypeTag, NumComponents));
+ transportedQuantity[wCompIdx].resize(problem_.gridView().size(0));
+ transportedQuantity[nCompIdx].resize(problem_.gridView().size(0));
+
transportedQuantity = totalConcentration_;
}
//! \copydoc transportedQuantity()
@@ -195,7 +201,7 @@ public:
totalConcentration_[wCompIdx].resize(problem.gridView().size(0));
totalConcentration_[nCompIdx].resize(problem.gridView().size(0));
- restrictFluxInTransport_ = GET_PARAM_FROM_GROUP(TypeTag,bool, Impet, RestrictFluxInTransport);
+ restrictFluxInTransport_ = GET_PARAM_FROM_GROUP(TypeTag,int, Impet, RestrictFluxInTransport);
}
virtual ~FVTransport2P2C()
@@ -206,9 +212,10 @@ protected:
TransportSolutionType totalConcentration_;
Problem& problem_;
bool impet_;
+ int averagedFaces_;
static const 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$)
- bool restrictFluxInTransport_;
+ int restrictFluxInTransport_;
bool switchNormals;
};
@@ -238,6 +245,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt,
dt = 1E100;
// store if we do update Estimate for flux functions
impet_ = impet;
+ averagedFaces_ = 0;
// resize update vector and set to zero
updateVec.resize(GET_PROP_VALUE(TypeTag, NumComponents));
@@ -302,7 +310,12 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt,
}
} // end grid traversal
if(impet)
- Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell << " leads to dt = "<<dt * GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, CFLFactor)<< std::endl;
+ {
+ Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell << " leads to dt = "
+ <<dt * GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, CFLFactor)<< std::endl;
+ if(averagedFaces_ != 0)
+ Dune::dinfo << " Averageing done for " << averagedFaces_ << " faces. "<< std::endl;
+ }
return;
}
/* Updates the transported quantity once an update is calculated.
@@ -363,10 +376,11 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
Scalar pcI = cellDataI.capillaryPressure();
DimMatrix K_I(problem().spatialParams().intrinsicPermeability(*elementPtrI));
- Scalar SwmobI = std::max((cellDataI.saturation(wPhaseIdx)
+ PhaseVector SmobI(0.);
+ SmobI[wPhaseIdx] = std::max((cellDataI.saturation(wPhaseIdx)
- problem().spatialParams().materialLawParams(*elementPtrI).Swr())
, 1e-2);
- Scalar SnmobI = std::max((cellDataI.saturation(nPhaseIdx)
+ SmobI[nPhaseIdx] = std::max((cellDataI.saturation(nPhaseIdx)
- problem().spatialParams().materialLawParams(*elementPtrI).Snr())
, 1e-2);
@@ -380,6 +394,10 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
unitOuterNormal *= -1.0;
Scalar faceArea = intersection.geometry().volume();
+// // local interface index
+// const int indexInInside = intersection.indexInInside();
+// const int indexInOutside = intersection.indexInOutside();
+
// create vector for timestep and for update
Dune::FieldVector<Scalar, 2> factor (0.);
Dune::FieldVector<Scalar, 2> updFactor (0.);
@@ -443,128 +461,101 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
potential[wPhaseIdx] += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
// determine upwinding direction, perform upwinding if possible
- Dune::FieldVector<bool, GET_PROP_VALUE(TypeTag, NumPhases)> doUpwinding(true);
+ Dune::FieldVector<bool, NumPhases> doUpwinding(true);
PhaseVector lambda(0.);
- if(!impet_ or !restrictFluxInTransport_) // perform a simple uwpind scheme
+ for(int phaseIdx = 0; phaseIdx < NumPhases; phaseIdx++)
{
- if (potential[wPhaseIdx] > 0.)
- {
- 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);
- }
+ int contiEqIdx = 0;
+ if(phaseIdx == wPhaseIdx)
+ contiEqIdx = contiWEqIdx;
else
- {
- doUpwinding[wPhaseIdx] = false;
- cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
- cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
- }
+ contiEqIdx = contiNEqIdx;
- if (potential[nPhaseIdx] > 0.)
- {
- 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
+ if(!impet_ or restrictFluxInTransport_==0) // perform a strict uwpind scheme
{
- doUpwinding[nPhaseIdx] = false;
- cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
- cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
+ if (potential[phaseIdx] > 0.)
+ {
+ lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiEqIdx, true);
+
+ }
+ else if(potential[phaseIdx] < 0.)
+ {
+ lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiEqIdx, false);
+ }
+ else
+ {
+ doUpwinding[phaseIdx] = false;
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiEqIdx, false);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiEqIdx, false);
+ }
}
- }
- else // if new potentials indicate same flow direction as in P.E., perform upwind
- {
- 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.
+ else // Transport after PE with check on flow direction
{
- 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
- lambda[wPhaseIdx] = cellDataI.mobility(wPhaseIdx);
- else if (isUpwindCell && !(cellDataJ.mobility(wPhaseIdx) != 0. && cellDataI.mobility(wPhaseIdx) == 0.)) // check if inflow induce neglected phase flux
- lambda[wPhaseIdx] = cellDataJ.mobility(wPhaseIdx);
- else
- doUpwinding[wPhaseIdx] = false;
+ bool wasUpwindCell = cellDataI.isUpwindCell(intersection.indexInInside(), contiEqIdx);
+
+ if (potential[phaseIdx] > 0. && wasUpwindCell)
+ lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
+ else if (potential[phaseIdx] < 0. && !wasUpwindCell)
+ lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
+ // potential direction does not coincide with that of P.E.
+ else if(restrictFluxInTransport_ == 2) // perform central averageing for all direction changes
+ doUpwinding[phaseIdx] = false;
+ else // i.e. restrictFluxInTransport == 1
+ {
+ //check if harmonic weithing is necessary
+ if (!wasUpwindCell && (cellDataJ.mobility(phaseIdx) != 0. // check if outflow induce neglected (i.e. mob=0) phase flux
+ or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father())))
+ lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
+ else if (wasUpwindCell && (cellDataI.mobility(phaseIdx) != 0. // check if inflow induce neglected phase flux
+ or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father())))
+ lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
+ else
+ doUpwinding[phaseIdx] = false;
+ }
+
}
- 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.
+ // do not perform upwinding if so desired
+ if(!doUpwinding[phaseIdx])
{
- 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
- lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
- else if (isUpwindCell && !(cellDataJ.mobility(nPhaseIdx) != 0. && cellDataI.mobility(nPhaseIdx) == 0.)) // check if inflow induce neglected phase flux
- lambda[nPhaseIdx] = cellDataJ.mobility(nPhaseIdx);
- else
- doUpwinding[nPhaseIdx] = false;
- }
- }
+ //a) no flux if there wont be any flux regardless how to average/upwind
+ if(cellDataI.mobility(phaseIdx)+cellDataJ.mobility(phaseIdx)==0.)
+ {
+ potential[phaseIdx] = 0;
+ continue;
+ }
- // 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(globalIdxI > globalIdxJ)
- 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(globalIdxI > globalIdxJ)
- Dune::dinfo << "harmonicMean flux of phase n used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
+ //b) perform harmonic averageing
+ fluxEntries[wCompIdx] -= potential[phaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(phaseIdx, wCompIdx) * cellDataI.mobility(phaseIdx) * cellDataI.density(phaseIdx),
+ cellDataJ.massFraction(phaseIdx, wCompIdx) * cellDataJ.mobility(phaseIdx) * cellDataJ.density(phaseIdx));
+ fluxEntries[nCompIdx] -= potential[phaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(phaseIdx, nCompIdx) * cellDataI.mobility(phaseIdx) * cellDataI.density(phaseIdx),
+ cellDataJ.massFraction(phaseIdx, nCompIdx) * cellDataJ.mobility(phaseIdx) * cellDataJ.density(phaseIdx));
+ // c) timestep control
+ // for timestep control : influx
+ timestepFlux[0] += std::max(0.,
+ - potential[phaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.mobility(phaseIdx),cellDataJ.mobility(phaseIdx)));
+ // outflux
+ timestepFlux[1] += std::max(0.,
+ potential[phaseIdx] * faceArea / volume
+ * harmonicMean(cellDataI.mobility(phaseIdx),cellDataJ.mobility(phaseIdx))/SmobI[phaseIdx]);
+
+ //d) output (only for one side)
+ averagedFaces_++;
+ #if DUNE_MINIMAL_DEBUG_LEVEL < 3
+ // verbose (only for one side)
+ if(globalIdxI > globalIdxJ)
+ Dune::dinfo << "harmonicMean flux of phase" << phaseIdx <<" used from cell" << globalIdxI<< " into " << globalIdxJ
+ << " ; TE upwind I = "<< cellDataI.isUpwindCell(intersection.indexInInside(), contiEqIdx) << " but pot = "<< potential[phaseIdx] << " \n";
+ #endif
+
+ //e) stop further standard calculations
+ potential[phaseIdx] = 0;
+ }
}
// calculate and standardized velocity
@@ -576,8 +567,8 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
// for timestep control : influx
timestepFlux[0] += velocityJIw + velocityJIn;
- double foutw = velocityIJw/SwmobI;
- double foutn = velocityIJn/SnmobI;
+ double foutw = velocityIJw/SmobI[wPhaseIdx];
+ double foutn = velocityIJn/SmobI[nPhaseIdx];
if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
timestepFlux[1] += foutw + foutn;
@@ -825,6 +816,9 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
FluidState &BCfluidState,
PhaseVector &pressBound)
{
+ // prepare a flash solver
+ CompositionalFlash<TypeTag> flashSolver;
+
const ElementPointer eIt= intersection.inside();
// read boundary values
PrimaryVariables primaryVariablesOnBoundary(0.);
@@ -859,18 +853,16 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
else // capillarity neglected
pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
-
- BCfluidState.satFlash(satBound, pressBound, problem().spatialParams().porosity(*eIt),
- problem().temperatureAtPos(globalPosFace));
-
+ flashSolver.saturationFlash2p2c(BCfluidState, satBound, pressBound,
+ problem().spatialParams().porosity(*eIt), problem().temperatureAtPos(globalPosFace));
}
else if (bcType == Indices::concentration)
{
// saturation and hence pc and hence corresponding pressure unknown
pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
Scalar Z1Bound = primaryVariablesOnBoundary[contiWEqIdx];
- BCfluidState.update(Z1Bound, pressBound, problem().spatialParams().porosity(*eIt),
- problem().temperatureAtPos(globalPosFace));
+ flashSolver.concentrationFlash2p2c(BCfluidState, Z1Bound, pressBound,
+ problem().spatialParams().porosity(*eIt), problem().temperatureAtPos(globalPosFace));
if(GET_PROP_VALUE(TypeTag, EnableCapillarity))
{
@@ -902,8 +894,8 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
//store old pc
Scalar oldPc = pcBound;
//update with better pressures
- BCfluidState.update(Z1Bound, pressBound, problem().spatialParams().porosity(*eIt),
- problem().temperatureAtPos(globalPosFace));
+ flashSolver.concentrationFlash2p2c(BCfluidState, Z1Bound, pressBound,
+ problem().spatialParams().porosity(*eIt), problem().temperatureAtPos(globalPosFace));
pcBound = MaterialLaw::pC(problem().spatialParams().materialLawParams(*eIt),
BCfluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
diff --git a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
index 83f1eab9be22061c4b2ad356d43e96e0fd6bbc0f..9248c32150e42093d54ec10169629fe3fb07acf5 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
@@ -135,7 +135,8 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
// initialize dt very large
dt = 1E100;
// store if we do update Estimate for flux functions
- this->impet_ = impet;
+ this->impet_ = impet;
+ this->averagedFaces_ = 0.;
// resize update vector and set to zero
updateVec.resize(GET_PROP_VALUE(TypeTag, NumComponents));
@@ -203,7 +204,11 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
}
} // end grid traversal
if(impet)
+ {
Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell << " leads to dt = "<<dt * GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, CFLFactor)<< std::endl;
+ if(this->averagedFaces_ != 0)
+ Dune::dinfo << " Averageing done for " << this->averagedFaces_ << " faces. "<< std::endl;
+ }
return;
}
}
diff --git a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
index 8d050b2c8d2e360eccfb53efad481418b23f6e5f..7652d6e5af7203d607f0ef167cd11498d5da7ede 100644
--- a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
@@ -74,49 +74,6 @@ public:
* \param presentPhaseIdx Subdomain Index = Indication which phase is present
* \param temperature Temperature \f$\mathrm{[K]}\f$
*/
- void update(const Scalar& Z1,const Dune::FieldVector<Scalar, numPhases> phasePressure,const int presentPhaseIdx, const Scalar& temperature)
- {
- pressure_[wPhaseIdx] = phasePressure[wPhaseIdx];
- pressure_[nPhaseIdx] = phasePressure[nPhaseIdx];
- temperature_ = temperature;
-
- if (presentPhaseIdx == wPhaseIdx)
- {
- massFractionWater_[wPhaseIdx] = Z1;
- massFractionWater_[nPhaseIdx] = 0.;
-
- moleFractionWater_[wPhaseIdx] = Z1 / FluidSystem::molarMass(0);
- moleFractionWater_[wPhaseIdx] /= ( Z1 / FluidSystem::molarMass(0)
- + (1.-Z1) / FluidSystem::molarMass(1));
- moleFractionWater_[nPhaseIdx] = 0.;
-
- presentPhaseIdx_ = wPhaseIdx;
- }
- else if (presentPhaseIdx == nPhaseIdx)
- {
- massFractionWater_[wPhaseIdx] = 0.;
- massFractionWater_[nPhaseIdx] = Z1;
-
- // interested in nComp => 1-X1
- moleFractionWater_[nPhaseIdx] = ( Z1 / FluidSystem::molarMass(0) ); // = moles of compIdx
- moleFractionWater_[nPhaseIdx] /= (Z1/ FluidSystem::molarMass(0)
- + (1.-Z1) / FluidSystem::molarMass(1) ); // /= total moles in phase
- moleFractionWater_[nPhaseIdx] = 0.;
-
- presentPhaseIdx_ = nPhaseIdx;
- }
- else
- Dune::dgrave << __FILE__ <<": Twophase conditions in single-phase flash!"
- << " Z1 is "<<Z1<< std::endl;
-
- density_ = FluidSystem::density(*this, presentPhaseIdx);
-
- aveMoMass_ = moleFractionWater_[presentPhaseIdx]*FluidSystem::molarMass(wCompIdx)
- + (1.-moleFractionWater_[presentPhaseIdx])*FluidSystem::molarMass(nCompIdx);
-
- return;
- }
- //@}
/*! \name Acess functions */
//@{
/*! \brief Returns the saturation of a phase.
@@ -134,10 +91,6 @@ public:
{
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
@@ -150,7 +103,7 @@ public:
if(phaseIdx == presentPhaseIdx_)
return density_;
else
- return FluidSystem::density(*this, phaseIdx);
+ return 0.; // only required for output if phase is not present anyhow
}
/*!
@@ -160,6 +113,10 @@ public:
*/
Scalar massFraction(int phaseIdx, int compIdx) const
{
+ if(phaseIdx != presentPhaseIdx_)
+ return 0.;
+
+
if (compIdx == wPhaseIdx)
return massFractionWater_[phaseIdx];
else
@@ -174,6 +131,8 @@ public:
*/
Scalar moleFraction(int phaseIdx, int compIdx) const
{
+ if(phaseIdx != presentPhaseIdx_)
+ return 0.;
if (compIdx == wPhaseIdx)
return moleFractionWater_[phaseIdx];
else
@@ -190,11 +149,6 @@ public:
assert(phaseIdx == presentPhaseIdx_);
return viscosity_;
}
- void setViscosity(int phaseIdx, Scalar value)
- {
- assert(phaseIdx == presentPhaseIdx_);
- viscosity_ = value;
- }
Scalar averageMolarMass(int phaseIdx) const
{
@@ -209,7 +163,86 @@ public:
*/
Scalar temperature(int phaseIdx) const
{ return temperature_; };
+ //@}
+ /*!
+ * \name Functions to set Data
+ */
+ //@{
+ /*!
+ * \brief Sets the viscosity of a phase \f$\mathrm{[Pa*s]}\f$.
+ *
+ * \param phaseIdx the index of the phase
+ * @param value Value to be stored
+ */
+ void setViscosity(int phaseIdx, Scalar value)
+ {
+ assert(phaseIdx == presentPhaseIdx_);
+ viscosity_ = value;
+ }
+ /*!
+ * \brief Sets the mass fraction of a component in a phase.
+ *
+ * \param phaseIdx the index of the phase
+ * \param compIdx the index of the component
+ * @param value Value to be stored
+ */
+ void setMassFraction(int phaseIdx, int compIdx, Scalar value)
+ {
+ if (compIdx == wCompIdx)
+ massFractionWater_[phaseIdx] = value;
+ else
+ massFractionWater_[phaseIdx] = 1- value;
+ }
+
+ /*!
+ * \brief Sets the molar fraction of a component in a fluid phase.
+ *
+ * \param phaseIdx the index of the phase
+ * \param compIdx the index of the component
+ * @param value Value to be stored
+ */
+ void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
+ {
+ if (compIdx == wCompIdx)
+ moleFractionWater_[phaseIdx] = value;
+ else
+ moleFractionWater_[phaseIdx] = 1-value;
+ }
+ /*!
+ * \brief Sets the density of a phase \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param phaseIdx the index of the phase
+ * @param value Value to be stored
+ */
+ void setDensity(int phaseIdx, Scalar value)
+ {
+ assert(phaseIdx == presentPhaseIdx_);
+ density_ = value;
+ }
+ /*!
+ * \brief Sets the phase Index that is present in this fluidState.
+ * @param phaseIdx the index of the phase
+ */
+ void setPresentPhaseIdx(int phaseIdx)
+ {
+ presentPhaseIdx_ = phaseIdx;
+ }
+
+ /*!
+ * \brief Sets the temperature
+ *
+ * @param value Value to be stored
+ */
+ void setTemperature(Scalar value)
+ {
+ temperature_ = value;
+ }
+ //TODO: doc me
+ void setAverageMolarMass(int phaseIdx, Scalar value)
+ {
+ aveMoMass_ = value;
+ }
/*!
* \brief Sets the phase pressure \f$\mathrm{[Pa]}\f$.
*/
diff --git a/dumux/material/constraintsolvers/compositionalflash.hh b/dumux/material/constraintsolvers/compositionalflash.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dc7e6eb43b630a5253920d22476dcf0463a15cea
--- /dev/null
+++ b/dumux/material/constraintsolvers/compositionalflash.hh
@@ -0,0 +1,322 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Benjamin Faigle *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Determines the pressures and saturations of all fluid phases
+ * given the total mass of all components.
+ */
+#ifndef DUMUX_COMPOSITIONAL_FLASH_HH
+#define DUMUX_COMPOSITIONAL_FLASH_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/exceptions.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+#include <dumux/decoupled/2p2c/2p2cfluidstate.hh>
+#include <dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh>
+
+namespace Dumux
+{
+/*!
+ * \brief Flash calculation routines for compositional decoupled models
+ *
+ * Routines for isothermal and isobaric 2p2c and 1p2c flash.
+ * \tparam TypeTag The property Type Tag
+ */
+template <class TypeTag>
+class CompositionalFlash
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef PseudoOnePTwoCFluidState<TypeTag> FluidState1p2c;
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
+
+ enum{
+ wPhaseIdx = FluidSystem::wPhaseIdx,
+ nPhaseIdx = FluidSystem::nPhaseIdx,
+ wCompIdx = FluidSystem::wCompIdx,
+ nCompIdx = FluidSystem::nCompIdx
+ };
+
+public:
+ typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
+ typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
+/*!
+ * \name Concentration flash for a given feed fraction
+ */
+//@{
+ //! 2p2c Flash for constant p & t if concentrations (feed mass fraction) is given.
+ /*!
+ * Routine goes as follows:
+ * - determination of the equilibrium constants from the fluid system
+ * - determination of maximum solubilities (mole fractions) according to phase pressures
+ * - comparison with Z1 to determine phase presence => phase mass fractions
+ * - round off fluid properties
+ * \param fluidState The decoupled fluid State
+ * \param Z1 Feed mass fraction: Mass of comp1 per total mass \f$\mathrm{[-]}\f$
+ * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$
+ * \param poro Porosity \f$\mathrm{[-]}\f$
+ * \param temperature Temperature \f$\mathrm{[K]}\f$
+ */
+ static void concentrationFlash2p2c(FluidState &fluidState,
+ const Scalar &Z1,
+ const PhaseVector &phasePressure,
+ const Scalar &porosity,
+ const Scalar &temperature)
+ {
+ // set the temperature, pressure
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(wPhaseIdx, phasePressure[wPhaseIdx]);
+ fluidState.setPressure(nPhaseIdx, phasePressure[nPhaseIdx]);
+
+ //mole equilibrium ratios K for in case wPhase is reference phase
+ double k1 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, wCompIdx); // = p^wComp_vap
+ double k2 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, nCompIdx); // = H^nComp_w
+
+ // get mole fraction from equilibrium konstants
+ fluidState.setMoleFraction(wPhaseIdx,wCompIdx, ((1. - k2) / (k1 -k2)));
+ fluidState.setMoleFraction(nPhaseIdx,wCompIdx, (fluidState.moleFraction(wPhaseIdx,wCompIdx) * k1));
+
+ // transform mole to mass fractions
+ fluidState.setMassFraction(wPhaseIdx, wCompIdx,
+ (fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ / ( fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ + (1.-fluidState.moleFraction(wPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
+ fluidState.setMassFraction(nPhaseIdx,wCompIdx,
+ (fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ / ( fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ + (1.-fluidState.moleFraction(nPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
+
+ //mass equilibrium ratios
+ Scalar equilRatio_[numPhases][numComponents];
+ equilRatio_[nPhaseIdx][wCompIdx] = fluidState.massFraction(nPhaseIdx,wCompIdx)
+ / fluidState.massFraction(wPhaseIdx, wCompIdx); // = Xn1 / Xw1 = K1
+ equilRatio_[nPhaseIdx][nCompIdx] = (1.-fluidState.massFraction(nPhaseIdx, wCompIdx))
+ / (1.-fluidState.massFraction(wPhaseIdx, wCompIdx)); // =(1.-Xn1) / (1.-Xw1) = K2
+ equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
+
+ // phase fraction of nPhase [mass/totalmass]
+ fluidState.setNu(nPhaseIdx, 0.);
+
+ // check if there is enough of component 1 to form a phase
+ if (Z1 > fluidState.massFraction(nPhaseIdx,wCompIdx)
+ && Z1 < fluidState.massFraction(wPhaseIdx,wCompIdx))
+ fluidState.setNu(nPhaseIdx, -((equilRatio_[nPhaseIdx][wCompIdx]-1)*Z1 + (equilRatio_[nPhaseIdx][nCompIdx]-1)*(1-Z1))
+ / (equilRatio_[nPhaseIdx][wCompIdx]-1) / (equilRatio_[nPhaseIdx][nCompIdx] -1));
+ else if (Z1 <= fluidState.massFraction(nPhaseIdx,wCompIdx)) // too little wComp to form a phase
+ {
+ fluidState.setNu(nPhaseIdx, 1.); // only nPhase
+ fluidState.setMassFraction(nPhaseIdx,wCompIdx, Z1); // hence, assign complete mass dissolved into nPhase
+
+ // store as moleFractions
+ Scalar xw_n = Z1 /*=Xw_n*/ / FluidSystem::molarMass(wCompIdx); // = moles of compIdx
+ xw_n /= ( Z1 /*=Xw_n*/ / FluidSystem::molarMass(wCompIdx)
+ +(1- Z1 /*=Xn_n*/) / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
+
+ fluidState.setMoleFraction(nPhaseIdx,wCompIdx, xw_n);
+
+// // opposing non-present phase is already set to equilibrium mass fraction
+// fluidState.setMassFraction(wPhaseIdx,wCompIdx, 1.); // non present phase is set to be pure
+// fluidState.setMoleFraction(wPhaseIdx,wCompIdx, 1.); // non present phase is set to be pure
+ }
+ else // (Z1 >= Xw1) => no nPhase
+ {
+ fluidState.setNu(nPhaseIdx, 0.); // no second phase
+ fluidState.setMassFraction(wPhaseIdx, wCompIdx, Z1);
+
+ // store as moleFractions
+ Scalar xw_w = Z1 /*=Xw_w*/ / FluidSystem::molarMass(wCompIdx); // = moles of compIdx
+ xw_w /= ( Z1 /*=Xw_w*/ / FluidSystem::molarMass(wCompIdx)
+ +(1- Z1 /*=Xn_w*/) / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
+ fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xw_w);
+
+// // opposing non-present phase is already set to equilibrium mass fraction
+// fluidState.setMassFraction(nPhaseIdx,wCompIdx, 0.); // non present phase is set to be pure
+// fluidState.setMoleFraction(nPhaseIdx,wCompIdx, 0.); // non present phase is set to be pure
+ }
+
+ // complete array of mass fractions
+ fluidState.setMassFraction(wPhaseIdx, nCompIdx, 1. - fluidState.massFraction(wPhaseIdx,wCompIdx));
+ fluidState.setMassFraction(nPhaseIdx, nCompIdx, 1. - fluidState.massFraction(nPhaseIdx,wCompIdx));
+ // complete array of mole fractions
+ fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(wPhaseIdx,wCompIdx));
+ fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(nPhaseIdx,wCompIdx));
+
+ // complete phase mass fractions
+ fluidState.setNu(wPhaseIdx, 1. - fluidState.phaseMassFraction(nPhaseIdx));
+
+ // get densities with correct composition
+ fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, wPhaseIdx));
+ fluidState.setDensity(nPhaseIdx, FluidSystem::density(fluidState, nPhaseIdx));
+
+ Scalar Sw = fluidState.phaseMassFraction(wPhaseIdx) / fluidState.density(wPhaseIdx);
+ Sw /= (fluidState.phaseMassFraction(wPhaseIdx)/fluidState.density(wPhaseIdx)
+ + fluidState.phaseMassFraction(nPhaseIdx)/fluidState.density(nPhaseIdx));
+ fluidState.setSaturation(wPhaseIdx, Sw);
+ };
+
+ //! The simplest possible update routine for 1p2c "flash" calculations
+ /*!
+ * Routine goes as follows:
+ * - Check if we are in single phase condition
+ * - Assign total concentration to the present phase
+ *
+ * \param Z1 Feed mass fraction \f$\mathrm{[-]}\f$
+ * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$
+ * \param presentPhaseIdx Subdomain Index = Indication which phase is present
+ * \param temperature Temperature \f$\mathrm{[K]}\f$
+ */
+ static void concentrationFlash1p2c(FluidState1p2c& fluidState, const Scalar& Z1,const Dune::FieldVector<Scalar, numPhases> phasePressure,const int presentPhaseIdx, const Scalar& temperature)
+ {
+ // set the temperature, pressure
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(wPhaseIdx, phasePressure[wPhaseIdx]);
+ fluidState.setPressure(nPhaseIdx, phasePressure[nPhaseIdx]);
+
+ fluidState.setPresentPhaseIdx(presentPhaseIdx);
+ fluidState.setMassFraction(presentPhaseIdx,wCompIdx, Z1);
+
+ // calculate mole fraction and average molar mass
+ Scalar xw_alpha= Z1 / FluidSystem::molarMass(wCompIdx);
+ xw_alpha /= ( Z1 / FluidSystem::molarMass(wCompIdx)
+ + (1.-Z1) / FluidSystem::molarMass(nCompIdx));
+ fluidState.setMoleFraction(presentPhaseIdx, wCompIdx, xw_alpha);
+
+// if (presentPhaseIdx == wPhaseIdx)
+// {
+//
+//// fluidState.setMassFraction(wPhaseIdx,wCompIdx, 0.;
+//
+//
+//
+//
+//
+//// fluidState.moleFractionWater_[nPhaseIdx] = 0.;
+//
+// fluidState.setPresentPhaseIdx(presentPhaseIdx);
+// }
+// else if (presentPhaseIdx == nPhaseIdx)
+// {
+// fluidState.setMassFraction[wPhaseIdx] = 0.;
+// fluidState.setMassFraction[nPhaseIdx] = Z1;
+//
+// // interested in nComp => 1-X1
+// fluidState.moleFractionWater_[nPhaseIdx] = ( Z1 / FluidSystem::molarMass(0) ); // = moles of compIdx
+// fluidState.moleFractionWater_[nPhaseIdx] /= (Z1/ FluidSystem::molarMass(0)
+// + (1.-Z1) / FluidSystem::molarMass(1) ); // /= total moles in phase
+// fluidState.moleFractionWater_[nPhaseIdx] = 0.;
+//
+// fluidState.presentPhaseIdx_ = nPhaseIdx;
+// }
+// else
+// Dune::dgrave << __FILE__ <<": Twophase conditions in single-phase flash!"
+// << " Z1 is "<<Z1<< std::endl;
+
+ fluidState.setAverageMolarMass(presentPhaseIdx,
+ fluidState.massFraction(presentPhaseIdx, wCompIdx)*FluidSystem::molarMass(wCompIdx)
+ +fluidState.massFraction(presentPhaseIdx, nCompIdx)*FluidSystem::molarMass(nCompIdx));
+
+ fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx));
+
+ return;
+ }
+//@}
+
+/*!
+ * \name Saturation flash for a given saturation (e.g. at boundary)
+ */
+//@{
+ //! a flash routine for 2p2c systems if the saturation instead of total concentration is known.
+ /*!
+ * Routine goes as follows:
+ * - determination of the equilibrium constants from the fluid system
+ * - determination of maximum solubilities (mole fractions) according to phase pressures
+ * - round off fluid properties
+ * \param fluidState The decoupled fluid state
+ * \param sat Saturation of phase 1 \f$\mathrm{[-]}\f$
+ * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$
+ * \param poro Porosity \f$\mathrm{[-]}\f$
+ * \param temperature Temperature \f$\mathrm{[K]}\f$
+ */
+ static void saturationFlash2p2c(FluidState &fluidState,
+ const Scalar &saturation,
+ const PhaseVector &phasePressure,
+ const Scalar &porosity,
+ const Scalar &temperature)
+ {
+ if (saturation == 0. || saturation == 1.)
+ Dune::dinfo << "saturation initial and boundary conditions set to zero or one!"
+ << " assuming fully saturated compositional conditions" << std::endl;
+
+ // set the temperature, pressure
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(wPhaseIdx, phasePressure[wPhaseIdx]);
+ fluidState.setPressure(nPhaseIdx, phasePressure[nPhaseIdx]);
+
+ //in contrast to the standard update() method, satflash() does not calculate nu.
+ fluidState.setNu(wPhaseIdx, NAN);
+ fluidState.setNu(nPhaseIdx, NAN);
+
+ //mole equilibrium ratios K for in case wPhase is reference phase
+ double k1 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, wCompIdx); // = p^wComp_vap
+ double k2 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, nCompIdx); // = H^nComp_w
+
+ // get mole fraction from equilibrium konstants
+ fluidState.setMoleFraction(wPhaseIdx,wCompIdx, ((1. - k2) / (k1 -k2)));
+ fluidState.setMoleFraction(nPhaseIdx,wCompIdx, (fluidState.moleFraction(wPhaseIdx,wCompIdx) * k1));
+
+ // transform mole to mass fractions
+ fluidState.setMassFraction(wPhaseIdx, wCompIdx,
+ (fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ / ( fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ + (1.-fluidState.moleFraction(wPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
+ fluidState.setMassFraction(nPhaseIdx,wCompIdx,
+ (fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ / ( fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ + (1.-fluidState.moleFraction(nPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
+
+ // complete array of mass fractions
+ fluidState.setMassFraction(wPhaseIdx, nCompIdx, 1. - fluidState.massFraction(wPhaseIdx,wCompIdx));
+ fluidState.setMassFraction(nPhaseIdx, nCompIdx, 1. - fluidState.massFraction(nPhaseIdx,wCompIdx));
+ // complete array of mole fractions
+ fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(wPhaseIdx,wCompIdx));
+ fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(nPhaseIdx,wCompIdx));
+
+ // get densities with correct composition
+ fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, wPhaseIdx));
+ fluidState.setDensity(nPhaseIdx, FluidSystem::density(fluidState, nPhaseIdx));
+
+ // set saturation
+ fluidState.setSaturation(wPhaseIdx, saturation);
+ }
+//@}
+};
+
+} // end namespace Dumux
+
+#endif