Commit d38a4f4b authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

improve central weighting

	-consistent storage of flux directions in fluxdata -> also applicable for adaptive grids.
	- central in PE for the rare case of potential = 0
	- improve future scalability for 3p3c model
- Introduced flash solver (as had been intended already for dumux 2.1) and removed methods from fluidstates
- properties put in groups
discussed with Philipp

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8893 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 1ba4ef46
......@@ -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
......
......@@ -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");
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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) [-]