OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Remove changes not needed for a minimal PR.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2022-08-03 08:21:18 +02:00
parent 4001fcffcf
commit 75a5a3d135
5 changed files with 18 additions and 1303 deletions
Download Patch File Download Diff File Expand all files Collapse all files

609
opm/models/blackoil/blackoilintensivequantitiessimple.hh
View File

@@ -1,609 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::BlackOilIntensiveQuantities
*/
#ifndef EWOMS_BLACK_OIL_INTENSIVE_SIMPLE_QUANTITIES_HH
#define EWOMS_BLACK_OIL_INTENSIVE_SIMPLE_QUANTITIES_HH
#include "blackoilproperties.hh"
#include "blackoilsolventmodules.hh"
#include "blackoilextbomodules.hh"
#include "blackoilpolymermodules.hh"
#include "blackoilfoammodules.hh"
#include "blackoilbrinemodules.hh"
#include "blackoilenergymodules.hh"
#include "blackoildiffusionmodule.hh"
#include "blackoilmicpmodules.hh"
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <dune/common/fmatrix.hh>
#include <cstring>
#include <utility>
namespace Opm {
/*!
* \ingroup BlackOilModel
* \ingroup IntensiveQuantities
*
* \brief Contains the quantities which are are constant within a
* finite volume in the black-oil model.
*/
template <class TypeTag>
class BlackOilIntensiveQuantitiesSimple
: public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
, public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
, public BlackOilDiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
, public BlackOilSolventIntensiveQuantities<TypeTag>
, public BlackOilExtboIntensiveQuantities<TypeTag>
, public BlackOilPolymerIntensiveQuantities<TypeTag>
, public BlackOilFoamIntensiveQuantities<TypeTag>
, public BlackOilBrineIntensiveQuantities<TypeTag>
, public BlackOilEnergyIntensiveQuantities<TypeTag>
, public BlackOilMICPIntensiveQuantities<TypeTag>
{
using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using FluxModule = GetPropType<TypeTag, Properties::FluxModule>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { waterCompIdx = FluidSystem::waterCompIdx };
enum { oilCompIdx = FluidSystem::oilCompIdx };
enum { gasCompIdx = FluidSystem::gasCompIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { dimWorld = GridView::dimensionworld };
enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
static const bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
static const bool waterEnabled = Indices::waterEnabled;
static const bool gasEnabled = Indices::gasEnabled;
static const bool oilEnabled = Indices::oilEnabled;
using Toolbox = MathToolbox<Evaluation>;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
public:
using FluidState = BlackOilFluidState<Evaluation,
FluidSystem,
enableTemperature,
enableEnergy,
compositionSwitchEnabled,
enableEvaporation,
enableBrine,
enableSaltPrecipitation,
Indices::numPhases>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
BlackOilIntensiveQuantitiesSimple()
{
if (compositionSwitchEnabled) {
fluidState_.setRs(0.0);
fluidState_.setRv(0.0);
}
}
BlackOilIntensiveQuantitiesSimple(const BlackOilIntensiveQuantitiesSimple& other) = default;
BlackOilIntensiveQuantitiesSimple& operator=(const BlackOilIntensiveQuantitiesSimple& other) = default;
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
{
const auto& problem = elemCtx.problem();
const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
update(problem, priVars, globalSpaceIdx, timeIdx);
}
void update(const Problem& problem,const PrimaryVariables& priVars, unsigned globalSpaceIdx, unsigned timeIdx)
{
ParentType::update(problem, priVars, globalSpaceIdx, timeIdx);
const auto& linearizationType = problem.model().linearizer().getLinearizationType();
Scalar RvMax = FluidSystem::enableVaporizedOil()
? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
: 0.0;
Scalar RsMax = FluidSystem::enableDissolvedGas()
? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
: 0.0;
// asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
unsigned pvtRegionIdx = priVars.pvtRegionIndex();
fluidState_.setPvtRegionIndex(pvtRegionIdx);
//asImp_().updateSaltConcentration_(elemCtx, dofIdx, timeIdx);
// extract the water and the gas saturations for convenience
Evaluation Sw = 0.0;
if (waterEnabled) {
if (priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) {
Sw = 1.0;
} else if (priVars.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg && priVars.primaryVarsMeaning() != PrimaryVariables::Rvw_pg_Rv) {
Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
}
}
Evaluation Sg = 0.0;
if (compositionSwitchEnabled)
{
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
// -> threephase case
assert( priVars.primaryVarsMeaning() != PrimaryVariables::OnePhase_p );
Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
} else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
// -> gas-water case
Sg = 1.0 - Sw;
// deal with solvent
if (enableSolvent)
Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
} else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_po_Sg) {
// -> oil-gas case
Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
} else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
// -> gas case
Sg = 1.0;
} else
{
assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs);
// -> oil-water case
Sg = 0.0;
}
}
if (gasEnabled && waterEnabled && !oilEnabled) {
Sg = 1.0 - Sw;
}
Valgrind::CheckDefined(Sg);
Valgrind::CheckDefined(Sw);
Evaluation So = 1.0 - Sw - Sg;
// deal with solvent
if (enableSolvent)
So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
if (FluidSystem::phaseIsActive(waterPhaseIdx))
fluidState_.setSaturation(waterPhaseIdx, Sw);
if (FluidSystem::phaseIsActive(gasPhaseIdx))
fluidState_.setSaturation(gasPhaseIdx, Sg);
if (FluidSystem::phaseIsActive(oilPhaseIdx))
fluidState_.setSaturation(oilPhaseIdx, So);
//asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
// now we compute all phase pressures
Evaluation pC[numPhases];
const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
//const auto& materialParams = problem.materialLawParams(0);//NB improve speed
MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
// oil is the reference phase for pressure
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv || priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if (FluidSystem::phaseIsActive(phaseIdx))
fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
}
else {
const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
if (FluidSystem::phaseIsActive(phaseIdx))
fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
}
// calculate relative permeabilities. note that we store the result into the
// mobility_ class attribute. the division by the phase viscosity happens later.
MaterialLaw::relativePermeabilities(mobility_, materialParams, fluidState_);
Valgrind::CheckDefined(mobility_);
// update the Saturation functions for the blackoil solvent module.
//asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
// update extBO parameters
//asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
Evaluation SoMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
SoMax = max(fluidState_.saturation(oilPhaseIdx),
problem.maxOilSaturation(globalSpaceIdx));
}
// take the meaning of the switching primary variable into account for the gas
// and oil phase compositions
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
// in the threephase case, gas and oil phases are potentially present, i.e.,
// we use the compositions of the gas-saturated oil and oil-saturated gas.
if (FluidSystem::enableDissolvedGas()) {
const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRs(min(RsMax, RsSat));
}
else if (compositionSwitchEnabled)
fluidState_.setRs(0.0);
if (FluidSystem::enableVaporizedOil()) {
const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRv(min(RvMax, RvSat));
}
else if (compositionSwitchEnabled)
fluidState_.setRv(0.0);
if (FluidSystem::enableVaporizedWater()) {
const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx);
fluidState_.setRvw(RvwSat);
}
}
else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_po_Sg) {
// The switching variable is the water-gas ratio Rvw
const auto& Rvw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
fluidState_.setRvw(Rvw);
if (FluidSystem::enableVaporizedOil()) {
const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRv(min(RvMax, RvSat));
}
else if (compositionSwitchEnabled)
fluidState_.setRv(0.0);
}
else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
// The switching variable is the water-gas ratio Rvw
const auto& Rvw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
fluidState_.setRvw(Rvw);
const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRv(Rv);
if (FluidSystem::enableDissolvedGas()) {
// the oil phase is not present, but we need to compute its "composition" for
// the gravity correction anyway
const auto& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRs(min(RsMax, RsSat));
}
else {
fluidState_.setRs(0.0);
}
}
else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
// if the switching variable is the mole fraction of the gas component in the
// oil phase, we can directly set the composition of the oil phase
const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRs(min(RsMax, Rs));
if (FluidSystem::enableVaporizedOil()) {
// the gas phase is not present, but we need to compute its "composition"
// for the gravity correction anyway
const auto& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRv(min(RvMax, RvSat));
}
else
fluidState_.setRv(0.0);
if (FluidSystem::enableVaporizedWater()) {
const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx);
fluidState_.setRvw(RvwSat);
}
}
else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRv(Rv);
if (FluidSystem::enableDissolvedGas()) {
// the oil phase is not present, but we need to compute its "composition" for
// the gravity correction anyway
const auto& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRs(min(RsMax, RsSat));
} else {
fluidState_.setRs(0.0);
}
if (FluidSystem::enableVaporizedWater()) {
const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx);
fluidState_.setRvw(RvwSat);
}
} else {
assert(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p);
}
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.setRegionIndex(pvtRegionIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
paramCache.setMaxOilSat(SoMax);
}
paramCache.updateAll(fluidState_);
// compute the phase densities and transform the phase permeabilities into mobilities
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
fluidState_.setInvB(phaseIdx, b);
const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
if (enableExtbo && phaseIdx == oilPhaseIdx)
mobility_[phaseIdx] /= asImp_().oilViscosity();
else if (enableExtbo && phaseIdx == gasPhaseIdx)
mobility_[phaseIdx] /= asImp_().gasViscosity();
else
mobility_[phaseIdx] /= mu;
}
Valgrind::CheckDefined(mobility_);
// calculate the phase densities
Evaluation rho;
if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
rho = fluidState_.invB(waterPhaseIdx);
rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
fluidState_.setDensity(waterPhaseIdx, rho);
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
rho = fluidState_.invB(gasPhaseIdx);
rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableVaporizedOil()) {
rho +=
fluidState_.invB(gasPhaseIdx) *
fluidState_.Rv() *
FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
}
if (FluidSystem::enableVaporizedWater()) {
rho +=
fluidState_.invB(gasPhaseIdx) *
fluidState_.Rvw() *
FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(gasPhaseIdx, rho);
}
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
rho = fluidState_.invB(oilPhaseIdx);
rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGas()) {
rho +=
fluidState_.invB(oilPhaseIdx) *
fluidState_.Rs() *
FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
}
fluidState_.setDensity(oilPhaseIdx, rho);
}
// retrieve the porosity from the problem
referencePorosity_ = problem.porosity(globalSpaceIdx, timeIdx);
porosity_ = referencePorosity_;
// the porosity must be modified by the compressibility of the
// rock...
Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
if (rockCompressibility > 0.0) {
Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
Evaluation x;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
} else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
} else {
x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
}
porosity_ *= 1.0 + x + 0.5*x*x;
}
// deal with water induced rock compaction
porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
// the MICP processes change the porosity
if (enableMICP){
Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
porosity_ += - biofilm_ - calcite_;
}
// deal with salt-precipitation
if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::Sp) {
Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
porosity_ *= (1.0 - Sp);
}
rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
// asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
// asImp_().zPvtUpdate_();
// asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
// asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
// asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
// asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
// asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
// // update the quantities which are required by the chosen
// // velocity model
// FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
// // update the diffusion specific quantities of the intensive quantities
// DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
#ifndef NDEBUG
// some safety checks in debug mode
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
assert(isfinite(fluidState_.density(phaseIdx)));
assert(isfinite(fluidState_.saturation(phaseIdx)));
assert(isfinite(fluidState_.temperature(phaseIdx)));
assert(isfinite(fluidState_.pressure(phaseIdx)));
assert(isfinite(fluidState_.invB(phaseIdx)));
}
assert(isfinite(fluidState_.Rs()));
assert(isfinite(fluidState_.Rv()));
#endif
}
/*!
* \copydoc ImmiscibleIntensiveQuantities::fluidState
*/
const FluidState& fluidState() const
{ return fluidState_; }
/*!
* \copydoc ImmiscibleIntensiveQuantities::mobility
*/
const Evaluation& mobility(unsigned phaseIdx) const
{ return mobility_[phaseIdx]; }
/*!
* \copydoc ImmiscibleIntensiveQuantities::porosity
*/
const Evaluation& porosity() const
{ return porosity_; }
/*!
* The pressure-dependent transmissibility multiplier due to rock compressibility.
*/
const Evaluation& rockCompTransMultiplier() const
{ return rockCompTransMultiplier_; }
/*!
* \brief Returns the index of the PVT region used to calculate the thermodynamic
* quantities.
*
* This allows to specify different Pressure-Volume-Temperature (PVT) relations in
* different parts of the spatial domain. Note that this concept should be seen as a
* work-around of the fact that the black-oil model does not capture the
* thermodynamics well enough. (Because there is, err, only a single real world with
* in which all substances follow the same physical laws and hence the same
* thermodynamics.) Anyway: Since the ECL file format uses multiple PVT regions, we
* support it as well in our black-oil model. (Note that, if it is not explicitly
* specified, the PVT region index is 0.)
*/
auto pvtRegionIndex() const
-> decltype(std::declval<FluidState>().pvtRegionIndex())
{ return fluidState_.pvtRegionIndex(); }
/*!
* \copydoc ImmiscibleIntensiveQuantities::relativePermeability
*/
Evaluation relativePermeability(unsigned phaseIdx) const
{
// warning: slow
return fluidState_.viscosity(phaseIdx)*mobility(phaseIdx);
}
/*!
* \brief Returns the porosity of the rock at reference conditions.
*
* I.e., the porosity of rock which is not perturbed by pressure and temperature
* changes.
*/
Scalar referencePorosity() const
{ return referencePorosity_; }
private:
friend BlackOilSolventIntensiveQuantities<TypeTag>;
friend BlackOilExtboIntensiveQuantities<TypeTag>;
friend BlackOilPolymerIntensiveQuantities<TypeTag>;
friend BlackOilEnergyIntensiveQuantities<TypeTag>;
friend BlackOilFoamIntensiveQuantities<TypeTag>;
friend BlackOilBrineIntensiveQuantities<TypeTag>;
friend BlackOilMICPIntensiveQuantities<TypeTag>;
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
FluidState fluidState_;
Scalar referencePorosity_;
Evaluation porosity_;
Evaluation rockCompTransMultiplier_;
Evaluation mobility_[numPhases];
};
} // namespace Opm
#endif

31
opm/models/blackoil/blackoillocalresidual.hh
View File

@@ -96,28 +96,29 @@ public:
* \copydoc FvBaseLocalResidual::computeStorage
*/
template <class LhsEval>
static void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx) const
{
// retrieve the intensive quantities for the SCV at the specified point in time
const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
computeStorage(storage,
intQuants);
}
const auto& fs = intQuants.fluidState();
/*!
* Compute storage term without using the element context.
*/
template <class LhsEval>
static void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const IntensiveQuantities& intQuants)
{
storage = 0.0;
const auto& fs = intQuants.fluidState();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
if (Indices::numPhases == 3) { // add trivial equation for the pseudo phase
unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
if (timeIdx == 0)
storage[conti0EqIdx + activeCompIdx] = variable<LhsEval>(0.0, conti0EqIdx + activeCompIdx);
else
storage[conti0EqIdx + activeCompIdx] = 0.0;
}
continue;
}
unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
LhsEval surfaceVolume =
Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))

40
opm/models/discretization/common/fvbasediscretization.hh
View File

@@ -363,7 +363,7 @@ class FvBaseDiscretization
using LocalLinearizer = GetPropType<TypeTag, Properties::LocalLinearizer>;
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
enum {
numEq = getPropValue<TypeTag, Properties::NumEq>(),
historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
@@ -402,6 +402,7 @@ class FvBaseDiscretization
using DiscreteFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables>;
// problem restriction and prolongation operator for adaptation
using Problem = GetPropType<TypeTag, Properties::Problem> ;
using ProblemRestrictProlongOperator = typename Problem :: RestrictProlongOperator ;
// discrete function restriction and prolongation operator for adaptation
@@ -783,43 +784,6 @@ public:
}
}
void invalidateAndUpdateIntensiveSingleQuantitiesSimple(const Problem& problem,
const PrimaryVariables& primaryVar,
unsigned dofIdx,
unsigned timeIdx) const
{
//invalidateIntensiveQuantitiesCache(timeIdx);
auto& intquant = intensiveQuantityCache_[timeIdx][dofIdx];
intquant.update(problem, primaryVar, dofIdx, timeIdx);
intensiveQuantityCacheUpToDate_[timeIdx][dofIdx] = true;
}
void invalidateAndUpdateIntensiveQuantitiesSimple(const Problem& problem,
const SolutionVector& primaryVars,
unsigned timeIdx) const
{
//invalidateIntensiveQuantitiesCache(timeIdx);
size_t numGridDof = primaryVars.size();
// omp_sched_t kind;
// int chunk;
//omp_get_schedule(&kind, &chunk);
//printf("%d %d\n", kind, chunk);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
const auto& primaryVar = primaryVars[dofIdx];
auto& intquant = intensiveQuantityCache_[timeIdx][dofIdx];
intquant.update(problem, primaryVar, dofIdx, timeIdx);
}
std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
intensiveQuantityCacheUpToDate_[timeIdx].end(),
/*value=*/true);
// loop over all elements...
}
/*!
* \brief Move the intensive quantities for a given time index to the back.
*

11
opm/models/discretization/common/fvbaseintensivequantities.hh
View File

@@ -46,8 +46,6 @@ class FvBaseIntensiveQuantities
using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
public:
// default constructor
@@ -71,15 +69,6 @@ public:
unsigned timeIdx)
{ extrusionFactor_ = elemCtx.problem().extrusionFactor(elemCtx, dofIdx, timeIdx); }
/*!
* \brief Update all quantities for a given control volume.
*/
void update(const Problem& problem,
const PrimaryVariables& /* primaryVars */,
unsigned /* globalSpaceIdx */,
unsigned /* timeIdx */)
{ extrusionFactor_ = problem.extrusionFactor(); }
/*!
* \brief Return how much a given sub-control volume is extruded.
*

630
opm/models/discretization/common/smallelementcontext.hh
View File

@@ -1,630 +0,0 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::FvBaseElementContext
*/
#ifndef EWOMS_SMALL_ELEMENT_CONTEXT_HH
#define EWOMS_SMALL_ELEMENT_CONTEXT_HH
#include "fvbaseproperties.hh"
#include <opm/models/discretization/common/linearizationtype.hh>
#include <opm/models/utils/alignedallocator.hh>
#include <opm/material/common/Unused.hpp>
#include <dune/common/fvector.hh>
#include <vector>
namespace Opm {
/*!
* \ingroup FiniteVolumeDiscretizations
*
* \brief This class stores an array of IntensiveQuantities objects, one
* intensive quantities object for each of the element's vertices
*/
template<class TypeTag>
class SmallElementContext
{
using Implementation = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
// the history size of the time discretization in number of steps
enum { timeDiscHistorySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>() };
struct DofStore_ {
IntensiveQuantities intensiveQuantities[timeDiscHistorySize];
const PrimaryVariables* priVars[timeDiscHistorySize];
const IntensiveQuantities *thermodynamicHint[timeDiscHistorySize];
};
using DofVarsVector = std::vector<DofStore_>;
using ExtensiveQuantitiesVector = std::vector<ExtensiveQuantities>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using Model = GetPropType<TypeTag, Properties::Model>;
using Stencil = GetPropType<TypeTag, Properties::Stencil>;
using GradientCalculator = GetPropType<TypeTag, Properties::GradientCalculator>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
static const unsigned dimWorld = GridView::dimensionworld;
static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
using CoordScalar = typename GridView::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
// we don't allow copies of element contexts!
SmallElementContext(const SmallElementContext& ) = delete;
public:
/*!
* \brief The constructor.
*/
explicit SmallElementContext(const Simulator& simulator)
: gridView_(simulator.gridView())
, stencil_(gridView_, simulator.model().dofMapper() )
{
// remember the simulator object
simulatorPtr_ = &simulator;
enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
stashedDofIdx_ = -1;
focusDofIdx_ = -1;
}
static void *operator new(size_t size)
{ return aligned_alloc(alignof(SmallElementContext), size); }
static void operator delete(void *ptr)
{ aligned_free(ptr); }
/*!
* \brief Construct all volume and extensive quantities of an element
* from scratch.
*
* \param elem The DUNE Codim<0> entity for which the volume
* variables ought to be calculated
*/
void updateAll(const Element& elem)
{
throw std::logic_error("Only use update stencil");
// asImp_().updateStencil(elem);
// asImp_().updateAllIntensiveQuantities();
// asImp_().updateAllExtensiveQuantities();
}
/*!
* \brief Compute the finite volume geometry for an element.
*
* \param elem The grid element for which the finite volume geometry ought to be
* computed.
*/
void updateStencil(const Element& elem)
{
// remember the current element
elemPtr_ = &elem;
// update the stencil. the center gradients are quite expensive to calculate and
// most models don't need them, so that we only do this if the model explicitly
// enables them
stencil_.update(elem);
// resize the arrays containing the flux and the volume variables
//dofVars_.resize(stencil_.numDof());
//extensiveQuantities_.resize(stencil_.numInteriorFaces());
}
/*!
* \brief Update the primary topological part of the stencil, but nothing else.
*
* \param elem The grid element for which the finite volume geometry ought to be
* computed.
*/
void updatePrimaryStencil(const Element& elem)
{
// remember the current element
elemPtr_ = &elem;
// update the finite element geometry
stencil_.updatePrimaryTopology(elem);
dofVars_.resize(stencil_.numPrimaryDof());
}
/*!
* \brief Update the topological part of the stencil, but nothing else.
*
* \param elem The grid element for which the finite volume geometry ought to be
* computed.
*/
void updateStencilTopology(const Element& elem)
{
// remember the current element
elemPtr_ = &elem;
// update the finite element geometry
stencil_.updateTopology(elem);
}
/*!
* \brief Compute the intensive quantities of all sub-control volumes of the current
* element for all time indices.
*/
void updateAllIntensiveQuantities()
{
// if (!enableStorageCache_) {
// // if the storage cache is disabled, we need to calculate the storage term
// // from scratch, i.e. we need the intensive quantities of all of the history.
// for (unsigned timeIdx = 0; timeIdx < timeDiscHistorySize; ++ timeIdx)
// asImp_().updateIntensiveQuantities(timeIdx);
// }
// else
// // if the storage cache is enabled, we only need to recalculate the storage
// // term for the most recent point of history (i.e., for the current iterative
// // solution)
// asImp_().updateIntensiveQuantities(/*timeIdx=*/0);
}
/*!
* \brief Compute the intensive quantities of all sub-control volumes of the current
* element for a single time index.
*
* \param timeIdx The index of the solution vector used by the time discretization.
*/
void updateIntensiveQuantities(unsigned timeIdx)
{ //updateIntensiveQuantities_(timeIdx, numDof(timeIdx));
}
/*!
* \brief Compute the intensive quantities of all sub-control volumes of the current
* element for a single time index.
*
* \param timeIdx The index of the solution vector used by the time discretization.
*/
void updatePrimaryIntensiveQuantities(unsigned timeIdx)
{
updateIntensiveQuantities_(timeIdx, numPrimaryDof(timeIdx));
}
/*!
* \brief Compute the intensive quantities of a single sub-control volume of the
* current element for a single time index.
*
* \param priVars The PrimaryVariables which should be used to calculate the
* intensive quantities.
* \param dofIdx The local index in the current element of the sub-control volume
* which should be updated.
* \param timeIdx The index of the solution vector used by the time discretization.
*/
void updateIntensiveQuantities(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx)
{ asImp_().updateSingleIntQuants_(priVars, dofIdx, timeIdx); }
/*!
* \brief Compute the extensive quantities of all sub-control volume
* faces of the current element for all time indices.
*/
void updateAllExtensiveQuantities()
{ //asImp_().updateExtensiveQuantities(/*timeIdx=*/0);
}
/*!
* \brief Compute the extensive quantities of all sub-control volume
* faces of the current element for a single time index.
*
* \param timeIdx The index of the solution vector used by the
* time discretization.
*/
void updateExtensiveQuantities(unsigned timeIdx)
{
throw std::logic_error("Extensive quantities should not be used");
// gradientCalculator_.prepare(/*context=*/asImp_(), timeIdx);
// for (unsigned fluxIdx = 0; fluxIdx < numInteriorFaces(timeIdx); fluxIdx++) {
// extensiveQuantities_[fluxIdx].update(/*context=*/asImp_(),
// /*localIndex=*/fluxIdx,
// timeIdx);
// }
}
/*!
* \brief Sets the degree of freedom on which the simulator is currently "focused" on
*
* I.e., in the case of automatic differentiation, all derivatives are with regard to
* the primary variables of that degree of freedom. Only "primary" DOFs can be
* focused on.
*/
void setFocusDofIndex(unsigned dofIdx)
{ focusDofIdx_ = dofIdx; }
/*!
* \brief Returns the degree of freedom on which the simulator is currently "focused" on
*
* \copydetails setFocusDof()
*/
unsigned focusDofIndex() const
{ return focusDofIdx_; }
/*!
* \brief Returns the linearization type.
*
* \copydetails setLinearizationType()
*/
LinearizationType linearizationType() const
{ return this->model().linearizer().getLinearizationType(); }
/*!
* \brief Return a reference to the simulator.
*/
const Simulator& simulator() const
{ return *simulatorPtr_; }
/*!
* \brief Return a reference to the problem.
*/
const Problem& problem() const
{ return simulatorPtr_->problem(); }
/*!
* \brief Return a reference to the model.
*/
const Model& model() const
{ return simulatorPtr_->model(); }
/*!
* \brief Return a reference to the grid view.
*/
const GridView& gridView() const
{ return gridView_; }
/*!
* \brief Return the current element.
*/
const Element& element() const
{ return *elemPtr_; }
/*!
* \brief Return the number of sub-control volumes of the current element.
*/
size_t numDof(unsigned timeIdx) const
{ return stencil(timeIdx).numDof(); }
/*!
* \brief Return the number of primary degrees of freedom of the current element.
*/
size_t numPrimaryDof(unsigned timeIdx) const
{ return stencil(timeIdx).numPrimaryDof(); }
/*!
* \brief Return the number of non-boundary faces which need to be
* considered for the flux apporixmation.
*/
size_t numInteriorFaces(unsigned timeIdx) const
{ return stencil(timeIdx).numInteriorFaces(); }
/*!
* \brief Return the number of boundary faces which need to be
* considered for the flux apporixmation.
*/
size_t numBoundaryFaces(unsigned timeIdx) const
{ return stencil(timeIdx).numBoundaryFaces(); }
/*!
* \brief Return the current finite element geometry.
*
* \param timeIdx The index of the solution vector used by the
* time discretization.
*/
const Stencil& stencil(unsigned) const
{ return stencil_; }
/*!
* \brief Return the position of a local entities in global coordinates
*
* \param dofIdx The local index of the degree of freedom
* in the current element.
* \param timeIdx The index of the solution vector used by the
* time discretization.
*/
const GlobalPosition& pos(unsigned dofIdx, unsigned) const
{ return stencil_.subControlVolume(dofIdx).globalPos(); }
/*!
* \brief Return the global spatial index for a sub-control volume
*
* \param dofIdx The local index of the degree of freedom
* in the current element.
* \param timeIdx The index of the solution vector used by the
* time discretization.
*/
unsigned globalSpaceIndex(unsigned dofIdx, unsigned timeIdx) const
{ return stencil(timeIdx).globalSpaceIndex(dofIdx); }
/*!
* \brief Return the element-local volume associated with a degree of freedom
*
* In the case of the vertex-centered finite volume method, this is different from
* the total volume because a finite volume usually spans multiple elements...
*
* \param dofIdx The local index of the degree of freedom in the current element.
* \param timeIdx The index of the solution vector used by the time discretization.
*/
Scalar dofVolume(unsigned dofIdx, unsigned timeIdx) const
{ return stencil(timeIdx).subControlVolume(dofIdx).volume(); }
/*!
* \brief Return the total volume associated with a degree of freedom
*
* "Total" means the volume controlled by a degree of freedom disregarding the
* element. (For example in the vertex-centered finite volume method, a control
* volume typically encompasses parts of multiple elements.)
*
* \param dofIdx The local index of the degree of freedom in the current element.
* \param timeIdx The index of the solution vector used by the time discretization.
*/
Scalar dofTotalVolume(unsigned dofIdx, unsigned timeIdx) const
{ return model().dofTotalVolume(globalSpaceIndex(dofIdx, timeIdx)); }
/*!
* \brief Returns whether the current element is on the domain's
* boundary.
*/
bool onBoundary() const
{ return element().hasBoundaryIntersections(); }
/*!
* \brief Return a reference to the intensive quantities of a
* sub-control volume at a given time.
*
* If the time step index is not given, return the volume
* variables for the current time.
*
* \param dofIdx The local index of the degree of freedom in the current element.
* \param timeIdx The index of the solution vector used by the time discretization.
*/
const IntensiveQuantities& intensiveQuantities(unsigned dofIdx, unsigned timeIdx) const
{
#ifndef NDEBUG
assert(dofIdx < numDof(timeIdx));
if (enableStorageCache_ && timeIdx != 0 && problem().recycleFirstIterationStorage())
throw std::logic_error("If caching of the storage term is enabled, only the intensive quantities "
"for the most-recent substep (i.e. time index 0) are available!");
#endif
unsigned globalIdx = globalSpaceIndex(dofIdx, timeIdx);
const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
if(cachedIntQuants){
return *cachedIntQuants;//dofVars_[dofIdx].intensiveQuantities[timeIdx];
}else{
throw std::logic_error("Cached quantities need to be calcualted before for this element context");
}
}
/*!
* \brief Return the thermodynamic hint for a given local index.
*
* \sa Discretization::thermodynamicHint(int, int)
*
* \param dofIdx The local index of the degree of freedom in the current element.
* \param timeIdx The index of the solution vector used by the time discretization.
*/
// const IntensiveQuantities *thermodynamicHint(unsigned dofIdx, unsigned timeIdx) const
// {
// assert(dofIdx < numDof(timeIdx));
// return dofVars_[dofIdx].thermodynamicHint[timeIdx];
// }
/*!
* \copydoc intensiveQuantities()
*/
// IntensiveQuantities& intensiveQuantities(unsigned dofIdx, unsigned timeIdx)
// {
// assert(dofIdx < numDof(timeIdx));
// return dofVars_[dofIdx].intensiveQuantities[timeIdx];
// }
/*!
* \brief Return the primary variables for a given local index.
*
* \param dofIdx The local index of the degree of freedom
* in the current element.
* \param timeIdx The index of the solution vector used by the
* time discretization.
*/
const PrimaryVariables& primaryVars(unsigned dofIdx, unsigned timeIdx) const
{
//need for compilation with vtk
//throw std::logic_error("Cached quantities need to be calcualted before for this element context");
assert(dofIdx < numDof(timeIdx));
return *dofVars_[dofIdx].priVars[timeIdx];
}
/*!
* \brief Returns true if no intensive quanties are stashed
*
* In most cases quantities are stashed only if a partial derivative is to be
* calculated via finite difference methods.
*/
bool haveStashedIntensiveQuantities() const
{ return stashedDofIdx_ != -1; }
/*!
* \brief Return the (local) index of the DOF for which the primary variables were
* stashed
*
* If none, then this returns -1.
*/
int stashedDofIdx() const
{ return stashedDofIdx_; }
/*!
* \brief Stash the intensive quantities for a degree of freedom on internal memory.
*
* \param dofIdx The local index of the degree of freedom in the current element.
*/
// void stashIntensiveQuantities(unsigned dofIdx)
// {
// assert(dofIdx < numDof(/*timeIdx=*/0));
// intensiveQuantitiesStashed_ = dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0];
// priVarsStashed_ = *dofVars_[dofIdx].priVars[/*timeIdx=*/0];
// stashedDofIdx_ = static_cast<int>(dofIdx);
// }
/*!
* \brief Restores the intensive quantities for a degree of freedom from internal memory.
*
* \param dofIdx The local index of the degree of freedom in the current element.
*/
// void restoreIntensiveQuantities(unsigned dofIdx)
// {
// dofVars_[dofIdx].priVars[/*timeIdx=*/0] = &priVarsStashed_;
// dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0] = intensiveQuantitiesStashed_;
// stashedDofIdx_ = -1;
// }
/*!
* \brief Return a reference to the gradient calculation class of
* the chosen spatial discretization.
*/
const GradientCalculator& gradientCalculator() const
{
throw std::logic_error("Gradient calculator should not be used");
// return gradientCalculator_;
}
/*!
* \brief Return a reference to the extensive quantities of a
* sub-control volume face.
*
* \param fluxIdx The local index of the sub-control volume face for which the
* extensive quantities are requested
* \param timeIdx The index of the solution vector used by the time discretization.
*/
const ExtensiveQuantities& extensiveQuantities(unsigned fluxIdx, unsigned) const
{
throw std::logic_error("Extensive quantiteis shoud not be used");
// return extensiveQuantities_[fluxIdx];
}
/*!
* \brief Returns true iff the cache for the storage term ought to be used for this context.
*
* If it is used, intensive quantities can only be accessed for the most recent time
* index. (time index 0.)
*/
bool enableStorageCache() const
{ return enableStorageCache_; }
/*!
* \brief Specifies if the cache for the storage term ought to be used for this context.
*/
void setEnableStorageCache(bool yesno)
{ enableStorageCache_ = yesno; }
private:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
protected:
/*!
* \brief Update the first 'n' intensive quantities objects from the primary variables.
*
* This method considers the intensive quantities cache.
*/
void updateIntensiveQuantities_(unsigned timeIdx, size_t numDof)
{
//throw std::invalid_argument("Intenisive quantiteis should not be updated for this element context");
// update the intensive quantities for the whole history
const SolutionVector& globalSol = model().solution(timeIdx);
// update the non-gradient quantities
for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
unsigned globalIdx = globalSpaceIndex(dofIdx, timeIdx);
const PrimaryVariables& dofSol = globalSol[globalIdx];
dofVars_[dofIdx].priVars[timeIdx] = &dofSol;
dofVars_[dofIdx].thermodynamicHint[timeIdx] =
model().thermodynamicHint(globalIdx, timeIdx);
const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
if (cachedIntQuants) {
dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
}
else {
updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
model().updateCachedIntensiveQuantities(dofVars_[dofIdx].intensiveQuantities[timeIdx],
globalIdx,
timeIdx);
}
}
}
void updateSingleIntQuants_(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx)
{
//throw std::invalid_argument("Intenisive quantiteis should not be updated for this element context");
#ifndef NDEBUG
if (enableStorageCache_ && timeIdx != 0 && problem().recycleFirstIterationStorage())
throw std::logic_error("If caching of the storage term is enabled, only the intensive quantities "
"for the most-recent substep (i.e. time index 0) are available!");
#endif
dofVars_[dofIdx].priVars[timeIdx] = &priVars;
dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx);
}
//IntensiveQuantities intensiveQuantitiesStashed_;
//PrimaryVariables priVarsStashed_;
//GradientCalculator gradientCalculator_;
std::vector<DofStore_, aligned_allocator<DofStore_, alignof(DofStore_)> > dofVars_;
//std::vector<ExtensiveQuantities, aligned_allocator<ExtensiveQuantities, alignof(ExtensiveQuantities)> > extensiveQuantities_;
const Simulator *simulatorPtr_;
const Element *elemPtr_;
const GridView gridView_;
Stencil stencil_;
int stashedDofIdx_;
int focusDofIdx_;
bool enableStorageCache_;
};
} // namespace Opm
#endif