From 75a5a3d135d8322768b1371134ba5349df9312b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 3 Aug 2022 08:21:18 +0200 Subject: [PATCH] Remove changes not needed for a minimal PR. --- .../blackoilintensivequantitiessimple.hh | 609 ----------------- opm/models/blackoil/blackoillocalresidual.hh | 31 +- .../common/fvbasediscretization.hh | 40 +- .../common/fvbaseintensivequantities.hh | 11 - .../common/smallelementcontext.hh | 630 ------------------ 5 files changed, 18 insertions(+), 1303 deletions(-) delete mode 100644 opm/models/blackoil/blackoilintensivequantitiessimple.hh delete mode 100644 opm/models/discretization/common/smallelementcontext.hh diff --git a/opm/models/blackoil/blackoilintensivequantitiessimple.hh b/opm/models/blackoil/blackoilintensivequantitiessimple.hh deleted file mode 100644 index da870f22f..000000000 --- a/opm/models/blackoil/blackoilintensivequantitiessimple.hh +++ /dev/null @@ -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 . - - 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 -#include -#include - -#include -#include - -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 BlackOilIntensiveQuantitiesSimple - : public GetPropType - , public GetPropType::FluxIntensiveQuantities - , public BlackOilDiffusionIntensiveQuantities() > - , public BlackOilSolventIntensiveQuantities - , public BlackOilExtboIntensiveQuantities - , public BlackOilPolymerIntensiveQuantities - , public BlackOilFoamIntensiveQuantities - , public BlackOilBrineIntensiveQuantities - , public BlackOilEnergyIntensiveQuantities - , public BlackOilMICPIntensiveQuantities -{ - using ParentType = GetPropType; - using Implementation = GetPropType; - - using Scalar = GetPropType; - using Evaluation = GetPropType; - using FluidSystem = GetPropType; - using MaterialLaw = GetPropType; - using ElementContext = GetPropType; - using PrimaryVariables = GetPropType; - using Indices = GetPropType; - using GridView = GetPropType; - using FluxModule = GetPropType; - - enum { numEq = getPropValue() }; - enum { enableSolvent = getPropValue() }; - enum { enableExtbo = getPropValue() }; - enum { enablePolymer = getPropValue() }; - enum { enableFoam = getPropValue() }; - enum { enableBrine = getPropValue() }; - enum { enableEvaporation = getPropValue() }; - enum { enableSaltPrecipitation = getPropValue() }; - enum { enableTemperature = getPropValue() }; - enum { enableEnergy = getPropValue() }; - enum { enableDiffusion = getPropValue() }; - enum { enableMICP = getPropValue() }; - enum { numPhases = getPropValue() }; - enum { numComponents = getPropValue() }; - 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; - using DimMatrix = Dune::FieldMatrix; - using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities; - using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities; - -public: - using FluidState = BlackOilFluidState; - using Problem = GetPropType; - - 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 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(*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(*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().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; - friend BlackOilExtboIntensiveQuantities; - friend BlackOilPolymerIntensiveQuantities; - friend BlackOilEnergyIntensiveQuantities; - friend BlackOilFoamIntensiveQuantities; - friend BlackOilBrineIntensiveQuantities; - friend BlackOilMICPIntensiveQuantities; - - - Implementation& asImp_() - { return *static_cast(this); } - - FluidState fluidState_; - Scalar referencePorosity_; - Evaluation porosity_; - Evaluation rockCompTransMultiplier_; - Evaluation mobility_[numPhases]; -}; - -} // namespace Opm - -#endif diff --git a/opm/models/blackoil/blackoillocalresidual.hh b/opm/models/blackoil/blackoillocalresidual.hh index 471f37621..9c28cc063 100644 --- a/opm/models/blackoil/blackoillocalresidual.hh +++ b/opm/models/blackoil/blackoillocalresidual.hh @@ -96,28 +96,29 @@ public: * \copydoc FvBaseLocalResidual::computeStorage */ template - static void computeStorage(Dune::FieldVector& storage, - const ElementContext& elemCtx, - unsigned dofIdx, - unsigned timeIdx) + void computeStorage(Dune::FieldVector& 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 - static void computeStorage(Dune::FieldVector& 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(0.0, conti0EqIdx + activeCompIdx); + else + storage[conti0EqIdx + activeCompIdx] = 0.0; + } + continue; + } + unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); LhsEval surfaceVolume = Toolbox::template decay(fs.saturation(phaseIdx)) diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh index 075e8f0f9..195020b5d 100644 --- a/opm/models/discretization/common/fvbasediscretization.hh +++ b/opm/models/discretization/common/fvbasediscretization.hh @@ -363,7 +363,7 @@ class FvBaseDiscretization using LocalLinearizer = GetPropType; using LocalResidual = GetPropType; - using Problem = GetPropType; + enum { numEq = getPropValue(), historySize = getPropValue(), @@ -402,6 +402,7 @@ class FvBaseDiscretization using DiscreteFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction; // problem restriction and prolongation operator for adaptation + using Problem = GetPropType ; 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. * diff --git a/opm/models/discretization/common/fvbaseintensivequantities.hh b/opm/models/discretization/common/fvbaseintensivequantities.hh index 4f49b320a..d8f432376 100644 --- a/opm/models/discretization/common/fvbaseintensivequantities.hh +++ b/opm/models/discretization/common/fvbaseintensivequantities.hh @@ -46,8 +46,6 @@ class FvBaseIntensiveQuantities using Implementation = GetPropType; using Scalar = GetPropType; using ElementContext = GetPropType; - using Problem = GetPropType; - using PrimaryVariables = GetPropType; 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. * diff --git a/opm/models/discretization/common/smallelementcontext.hh b/opm/models/discretization/common/smallelementcontext.hh deleted file mode 100644 index e66e127f8..000000000 --- a/opm/models/discretization/common/smallelementcontext.hh +++ /dev/null @@ -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 . - - 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 -#include - -#include - -#include - -#include - -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 SmallElementContext -{ - using Implementation = GetPropType; - - using Scalar = GetPropType; - using PrimaryVariables = GetPropType; - using IntensiveQuantities = GetPropType; - using ExtensiveQuantities = GetPropType; - - // the history size of the time discretization in number of steps - enum { timeDiscHistorySize = getPropValue() }; - - struct DofStore_ { - IntensiveQuantities intensiveQuantities[timeDiscHistorySize]; - const PrimaryVariables* priVars[timeDiscHistorySize]; - const IntensiveQuantities *thermodynamicHint[timeDiscHistorySize]; - }; - using DofVarsVector = std::vector; - using ExtensiveQuantitiesVector = std::vector; - - using Simulator = GetPropType; - using Problem = GetPropType; - using Model = GetPropType; - using Stencil = GetPropType; - using GradientCalculator = GetPropType; - using SolutionVector = GetPropType; - - using GridView = GetPropType; - using Element = typename GridView::template Codim<0>::Entity; - - static const unsigned dimWorld = GridView::dimensionworld; - static const unsigned numEq = getPropValue(); - - using CoordScalar = typename GridView::ctype; - using GlobalPosition = Dune::FieldVector; - - // 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(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(this); } - - const Implementation& asImp_() const - { return *static_cast(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 > dofVars_; - //std::vector > extensiveQuantities_; - - const Simulator *simulatorPtr_; - const Element *elemPtr_; - const GridView gridView_; - Stencil stencil_; - - int stashedDofIdx_; - int focusDofIdx_; - bool enableStorageCache_; -}; - -} // namespace Opm - -#endif