diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
index 851ed2018e0d23d74f83b89e12abe5f014269f17..808b6d3746cd96df98909b546c39a163a03a3e69 100644
--- a/dumux/decoupled/2p2c/2p2cproperties.hh
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -150,7 +150,7 @@ SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCompressibility, true);
SET_BOOL_PROP(DecoupledTwoPTwoC, VisitFacesOnlyOnce, false);
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCapillarity, false);
//!< Restrict (no upwind) flux in transport step if direction reverses after pressure equation
-SET_BOOL_PROP(DecoupledTwoPTwoC, RestrictFluxInTransport, true);
+SET_BOOL_PROP(DecoupledTwoPTwoC, RestrictFluxInTransport, false);
SET_PROP(DecoupledTwoPTwoC, BoundaryMobility)
{
diff --git a/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh b/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh
index 9d861ccd7091be403d96ac168f0c24a62c7f7c4e..0ed1e5d62c132fa2ad996aae762f9ab5725b915c 100644
--- a/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh
@@ -208,7 +208,7 @@ public:
const Scalar capillaryPressure() const
{
if(fluidStateType_ == simple)
- return 0.;
+ return simpleFluidState_->pressure(nPhaseIdx) - simpleFluidState_->pressure(wPhaseIdx);
else
return this->fluidState_->pressure(nPhaseIdx) - this->fluidState_->pressure(wPhaseIdx);
}
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index 13c821fa03dfd71ba9bf6c67dc50bf6fd47a7055..f10da742087b630e10aa09d2ead656ef1314bf19 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -194,6 +194,17 @@ protected:
const GlobalPosition& gravity_; //!< vector including the gravity constant
static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
Dune::Timer timer_;
+
+ //! Indices of matrix and rhs entries
+ /**
+ * During the assembling of the global system of equations get-functions are called (getSource(), getFlux(), etc.), which return global matrix or right hand side entries in a vector. These can be accessed using following indices:
+ */
+ enum
+ {
+ rhs = 1,//!<index for the right hand side entry
+ matrix = 0//!<index for the global matrix entry
+
+ };
};
//! function which assembles the system of equations to be solved
@@ -208,7 +219,13 @@ protected:
*/
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
-{ // initialization: set matrix this->A_ to zero
+{
+ if(first)
+ {
+ ParentType::assemble(true);
+ return;
+ }
+ // initialization: set matrix this->A_ to zero
this->A_ = 0;
this->f_ = 0;
@@ -227,7 +244,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
else
problem().pressureModel().getSource(entries,*eIt, cellDataI, first);
- this->f_[globalIdxI] = entries[1];
+ this->f_[globalIdxI] = entries[rhs];
/***** flux term ***********/
// iterate over all faces of the cell
@@ -246,11 +263,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
//set right hand side
- this->f_[globalIdxI] -= entries[1];
+ this->f_[globalIdxI] -= entries[rhs];
// set diagonal entry
- this->A_[globalIdxI][globalIdxI] += entries[0];
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
// set off-diagonal entry
- this->A_[globalIdxI][globalIdxJ] = -entries[0];
+ this->A_[globalIdxI][globalIdxJ] = -entries[matrix];
} // end neighbor
@@ -263,9 +280,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
//set right hand side
- this->f_[globalIdxI] += entries[1];
+ this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
- this->A_[globalIdxI][globalIdxI] += entries[0];
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
}
} //end interfaces loop
// printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);
@@ -276,9 +293,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
else
problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
- this->f_[globalIdxI] += entries[1];
+ this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
- this->A_[globalIdxI][globalIdxI] += entries[0];
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
} // end grid traversal
// printmatrix(std::cout, this->A_, "global stiffness matrix after assempling", "row", 11,3);
// printvector(std::cout, this->f_, "right hand side", "row", 10);
@@ -325,7 +342,10 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
Scalar incp = 1e-2;
// numerical derivative of fluid volume with respect to pressure
- Scalar p_ = cellDataI.pressure(pressureType) + incp;
+ PhaseVector p_(incp);
+ p_[nPhaseIdx] += cellDataI.pressure(nPhaseIdx);
+ p_[wPhaseIdx] += cellDataI.pressure(wPhaseIdx);
+
Scalar sumC = (cellDataI.massConcentration(wCompIdx) + cellDataI.massConcentration(nCompIdx));
Scalar Z1 = cellDataI.massConcentration(wCompIdx) / sumC;
// initialize simple fluidstate object
@@ -356,8 +376,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
Scalar compress_term = cellDataI.dv_dp() / timestep_;
- storageEntry[0] -= compress_term*volume;
- storageEntry[1] -= cellDataI.pressure(pressureType) * compress_term * volume;
+ storageEntry[matrix] -= compress_term*volume;
+ storageEntry[rhs] -= cellDataI.pressure(pressureType) * compress_term * volume;
if (isnan(compress_term) || isinf(compress_term))
DUNE_THROW(Dune::MathError, "Compressibility term leads to NAN matrix entry at index " << globalIdxI);
@@ -366,6 +386,14 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
DUNE_THROW(Dune::NotImplemented, "Compressibility is switched off???");
}
+ // Abort error damping if there will be a possibly tiny timestep compared with last one
+ // This might be the case if the episode or simulation comes to an end.
+ if( problem().timeManager().episodeWillBeOver()
+ || problem().timeManager().willBeFinished())
+ {
+ problem().variables().cellData(globalIdxI).errorCorrection() = 0.;
+ return;
+ }
// error reduction routine: volumetric error is damped and inserted to right hand side
// if damping is not done, the solution method gets unstable!
@@ -383,12 +411,12 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxError) && (!problem().timeManager().willBeFinished()))
{
if (erri <= x_mi * maxError)
- storageEntry[1] +=
+ 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)
* cellDataI.volumeError() * volume;
else
- storageEntry[1] +=
+ storageEntry[rhs] +=
problem().variables().cellData(globalIdxI).errorCorrection() =
fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxError)
* cellDataI.volumeError() * volume;
@@ -453,9 +481,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>
int phaseIdx = std::min(cellDataI.subdomain(), cellDataJ.subdomain());
Scalar rhoMean = 0.5 * (cellDataI.density(phaseIdx) + cellDataJ.density(phaseIdx));
-
- // reset potential gradients, density
- Scalar density = 0;
+ //Scalar density = 0;
// 1p => no pC => only 1 pressure, potential
Scalar potential = (cellDataI.pressure(phaseIdx) - cellDataJ.pressure(phaseIdx)) / dist;
@@ -465,18 +491,18 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>
Scalar lambda;
if (potential >= 0.)
- {
+ {
lambda = cellDataI.mobility(phaseIdx);
- density = cellDataI.density(phaseIdx);
- }
+ //density = cellDataI.density(phaseIdx);
+ }
else
- {
+ {
lambda = cellDataJ.mobility(phaseIdx);
- density = cellDataJ.density(phaseIdx);
- }
+ //density = cellDataJ.density(phaseIdx);
+ }
entries[0] = lambda * faceArea * (permeability * unitOuterNormal) / (dist);
- entries[1] = density * lambda;
+ entries[1] = rhoMean * lambda;
entries[1] *= faceArea * (permeability * gravity_) * (unitOuterNormal * unitDistVec);
return;
}
@@ -703,8 +729,23 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
PseudoOnePTwoCFluidState<TypeTag>& pseudoFluidState
= cellData.manipulateSimpleFluidState();
- // Todo: only variables of present phase need to be updated
- Scalar press1p = this->pressure(globalIdx);
+ // prepare phase pressure for fluid state
+ // both phase pressures are necessary for the case 1p domain is assigned for
+ // the next 2p subdomain
+ PhaseVector pressure(0.);
+ Scalar pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(*eIt),
+ ((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign Sw = 1 if wPhase present, else 0
+ if(pressureType == wPhaseIdx)
+ {
+ pressure[wPhaseIdx] = this->pressure(globalIdx);
+ pressure[nPhaseIdx] = this->pressure(globalIdx)+pc;
+ }
+ else
+ {
+ pressure[wPhaseIdx] = this->pressure(globalIdx)-pc;
+ pressure[nPhaseIdx] = this->pressure(globalIdx);
+ }
+
// make shure total concentrations from solution vector are exact in fluidstate
pseudoFluidState.setMassConcentration(wCompIdx,
problem().transportModel().totalConcentration(wCompIdx,globalIdx));
@@ -716,7 +757,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
+ pseudoFluidState.massConcentration(nCompIdx);
Scalar Z1 = pseudoFluidState.massConcentration(wCompIdx)/ sumConc;
- pseudoFluidState.update(Z1, press1p, cellData.subdomain(), problem().temperatureAtPos(globalPos));
+ pseudoFluidState.update(Z1, pressure, cellData.subdomain(), problem().temperatureAtPos(globalPos));
// write stuff in fluidstate
assert(presentPhaseIdx == pseudoFluidState.presentPhaseIdx());
diff --git a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
index 04dc3957c3952772f536c45a261e63bc953272fb..207076ec77c167cabb33f1c0db7249b5a5708ec7 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
@@ -134,7 +134,9 @@ 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;
+
// resize update vector and set to zero
updateVec.resize(GET_PROP_VALUE(TypeTag, NumComponents));
updateVec[wCompIdx].resize(problem().gridView().size(0));
diff --git a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
index 5eb8be53163f9e9dae2cc0ee0cc26ff83c53a61a..abf7852cfbcb7c5ce7851549502f1f5268c1d98c 100644
--- a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
+++ b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
@@ -70,14 +70,14 @@ public:
* - Assign total concentration to the present phase
*
* \param Z1 Feed mass fraction \f$\mathrm{[-]}\f$
- * \param press1p Pressure value for present phase \f$\mathrm{[Pa]}\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$
*/
- void update(const Scalar& Z1,const Scalar& press1p,const int presentPhaseIdx, const Scalar& temperature)
+ void update(const Scalar& Z1,const Dune::FieldVector<Scalar, numPhases> phasePressure,const int presentPhaseIdx, const Scalar& temperature)
{
-
- pressure1p_=press1p;
+ pressure_[wPhaseIdx] = phasePressure[wPhaseIdx];
+ pressure_[nPhaseIdx] = phasePressure[nPhaseIdx];
temperature_ = temperature;
if (presentPhaseIdx == wPhaseIdx)
@@ -141,7 +141,7 @@ public:
* \param phaseIdx Index of the phase
*/
Scalar pressure(int phaseIdx) const
- { return pressure1p_; }
+ { return pressure_[phaseIdx]; }
Scalar density(int phaseIdx) const
{
@@ -236,8 +236,8 @@ public:
Scalar massConcentration_[numComponents];
Scalar massFractionWater_[numPhases];
Scalar moleFractionWater_[numPhases];
+ Scalar pressure_[numPhases];
Scalar Sw_;
- Scalar pressure1p_;
Scalar density_;
Scalar viscosity_;
Scalar temperature_;