ebos: hook up the ECL hysteresis stuff

this required to add a new callback to the problem which is invoked
when the initial solution has been applied.
This commit is contained in:
Andreas Lauser 2015-07-28 17:24:44 +02:00
parent cc420cc06d
commit 8c843dd05b

View File

@ -176,8 +176,8 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
{
typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
// Grid and world dimension
@ -199,6 +199,8 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager EclMaterialLawManager;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
@ -355,6 +357,7 @@ public:
*/
void beginTimeStep()
{ wellManager_.beginTimeStep(); }
/*!
* \brief Called by the simulator before each Newton-Raphson iteration.
*/
@ -377,6 +380,8 @@ public:
// write the summary information after each time step
summaryWriter_.write(wellManager_);
updateHysteresis_();
#ifndef NDEBUG
// in debug mode, we don't care about performance, so we check if the model does
// the right thing (i.e., the mass change inside the whole reservoir must be
@ -574,11 +579,6 @@ public:
return deck->getKeyword("PVTNUM")->getIntData()[cartesianDofIdx] - 1;
}
/*!
* \name Problem parameters
*/
//! \{
/*!
* \copydoc FvBaseProblem::name
*/
@ -597,13 +597,6 @@ public:
return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0);
}
// \}
/*!
* \name Boundary conditions
*/
//! \{
/*!
* \copydoc FvBaseProblem::boundary
*
@ -616,13 +609,6 @@ public:
int timeIdx) const
{ values.setNoFlow(); }
//! \}
/*!
* \name Volumetric terms
*/
//! \{
/*!
* \copydoc FvBaseProblem::initial
*
@ -640,6 +626,11 @@ public:
values.assignMassConservative(initialFluidStates_[globalDofIdx], matParams);
}
void initialSolutionApplied()
{
updateHysteresis_();
}
/*!
* \copydoc FvBaseProblem::source
*
@ -665,8 +656,6 @@ public:
rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
}
//! \}
private:
static bool enableEclOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
@ -1140,6 +1129,31 @@ private:
return materialLawManager_.materialLawParams(cartesianCellIdx);
}
// update the hysteresis parameters of the material laws for the whole grid
void updateHysteresis_()
{
if (!materialLawManager_.enableHysteresis())
return;
ElementContext elemCtx(this->simulator());
const auto& gridManager = this->simulator().gridManager();
auto elemIt = gridManager.gridView().template begin</*codim=*/0>();
const auto &elemEndIt = gridManager.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updateStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
int compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int cartesianDofIdx = gridManager.cartesianCellId(compressedDofIdx);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
materialLawManager_.updateHysteresis(intQuants.fluidState(), cartesianDofIdx);
}
}
std::vector<Scalar> porosity_;
std::vector<DimMatrix> intrinsicPermeability_;
EclTransmissibility<TypeTag> transmissibilities_;