diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d4fcea02..391579832 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -423,7 +423,7 @@ set(FLOW_MODELS blackoil brine energy extbo foam gasoil gaswater oilwater oilwater_brine gaswater_brine oilwater_polymer oilwater_polymer_injectivity micp polymer solvent gasoil_energy brine_saltprecipitation gaswater_saltprec_vapwat brine_precsalt_vapwat) -set(FLOW_VARIANT_MODELS brine_energy onephase onephase_energy) +set(FLOW_VARIANT_MODELS brine_energy onephase onephase_energy blackoil_tpfa) set(FLOW_TGTS) foreach(OBJ ${COMMON_MODELS} ${FLOW_MODELS} ${FLOW_VARIANT_MODELS}) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index f8f7b2f13..9b63bac89 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -100,6 +100,7 @@ class EclTransExtensiveQuantities { using Implementation = GetPropType; + using IntensiveQuantities = GetPropType; using FluidSystem = GetPropType; using ElementContext = GetPropType; using Scalar = GetPropType; @@ -206,25 +207,28 @@ protected: void updatePolymer(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); } - /*! - * \brief Update the required gradients for interior faces - */ - void calculateGradients_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) - { - Valgrind::SetUndefined(*this); +public: + static void volumeAndPhasePressureDifferences(short (&upIdx)[numPhases], + short (&dnIdx)[numPhases], + Evaluation (&volumeFlux)[numPhases], + Evaluation (&pressureDifferences)[numPhases], + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) + { const auto& problem = elemCtx.problem(); const auto& stencil = elemCtx.stencil(timeIdx); const auto& scvf = stencil.interiorFace(scvfIdx); + unsigned interiorDofIdx = scvf.interiorIndex(); + unsigned exteriorDofIdx = scvf.exteriorIndex(); - interiorDofIdx_ = scvf.interiorIndex(); - exteriorDofIdx_ = scvf.exteriorIndex(); - assert(interiorDofIdx_ != exteriorDofIdx_); + assert(interiorDofIdx != exteriorDofIdx); - unsigned I = stencil.globalSpaceIndex(interiorDofIdx_); - unsigned J = stencil.globalSpaceIndex(exteriorDofIdx_); + unsigned I = stencil.globalSpaceIndex(interiorDofIdx); + unsigned J = stencil.globalSpaceIndex(exteriorDofIdx); - Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx_, exteriorDofIdx_); + Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx); Scalar faceArea = scvf.area(); Scalar thpres = problem.thresholdPressure(I, J); @@ -233,129 +237,170 @@ protected: // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) Scalar g = elemCtx.problem().gravity()[dimWorld - 1]; - const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx); - const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx_, timeIdx); + const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx); + const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx); // this is quite hacky because the dune grid interface does not provide a // cellCenterDepth() method (so we ask the problem to provide it). The "good" // solution would be to take the Z coordinate of the element centroids, but since // ECL seems to like to be inconsistent on that front, it needs to be done like // here... - Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx_, timeIdx); - Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx_, timeIdx); + Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx); + Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx); // the distances from the DOF's depths. (i.e., the additional depth of the // exterior DOF) Scalar distZ = zIn - zEx; + Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0); + Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0); + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { if (!FluidSystem::phaseIsActive(phaseIdx)) continue; - - // check shortcut: if the mobility of the phase is zero in the interior as - // well as the exterior DOF, we can skip looking at the phase. - if (intQuantsIn.mobility(phaseIdx) <= 0.0 && - intQuantsEx.mobility(phaseIdx) <= 0.0) - { - upIdx_[phaseIdx] = interiorDofIdx_; - dnIdx_[phaseIdx] = exteriorDofIdx_; - pressureDifference_[phaseIdx] = 0.0; - volumeFlux_[phaseIdx] = 0.0; + calculatePhasePressureDiff_(upIdx[phaseIdx], + dnIdx[phaseIdx], + pressureDifferences[phaseIdx], + intQuantsIn, + intQuantsEx, + phaseIdx,//input + interiorDofIdx,//input + exteriorDofIdx,//intput + Vin, + Vex, + I, + J, + distZ*g, + thpres); + if (pressureDifferences[phaseIdx] == 0) { + volumeFlux[phaseIdx] = 0.0; continue; } - // do the gravity correction: compute the hydrostatic pressure for the - // external at the depth of the internal one - const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx); - Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx)); - Evaluation rhoAvg = (rhoIn + rhoEx)/2; - - const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx); - Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx)); - if (enableExtbo) // added stability; particulary useful for solvent migrating in pure water - // where the solvent fraction displays a 0/1 behaviour ... - pressureExterior += Toolbox::value(rhoAvg)*(distZ*g); - else - pressureExterior += rhoAvg*(distZ*g); - - pressureDifference_[phaseIdx] = pressureExterior - pressureInterior; - - // decide the upstream index for the phase. for this we make sure that the - // degree of freedom which is regarded upstream if both pressures are equal - // is always the same: if the pressure is equal, the DOF with the lower - // global index is regarded to be the upstream one. - if (pressureDifference_[phaseIdx] > 0.0) { - upIdx_[phaseIdx] = exteriorDofIdx_; - dnIdx_[phaseIdx] = interiorDofIdx_; - } - else if (pressureDifference_[phaseIdx] < 0.0) { - upIdx_[phaseIdx] = interiorDofIdx_; - dnIdx_[phaseIdx] = exteriorDofIdx_; - } - else { - // if the pressure difference is zero, we chose the DOF which has the - // larger volume associated to it as upstream DOF - Scalar Vin = elemCtx.dofVolume(interiorDofIdx_, /*timeIdx=*/0); - Scalar Vex = elemCtx.dofVolume(exteriorDofIdx_, /*timeIdx=*/0); - if (Vin > Vex) { - upIdx_[phaseIdx] = interiorDofIdx_; - dnIdx_[phaseIdx] = exteriorDofIdx_; - } - else if (Vin < Vex) { - upIdx_[phaseIdx] = exteriorDofIdx_; - dnIdx_[phaseIdx] = interiorDofIdx_; - } - else { - assert(Vin == Vex); - // if the volumes are also equal, we pick the DOF which exhibits the - // smaller global index - if (I < J) { - upIdx_[phaseIdx] = interiorDofIdx_; - dnIdx_[phaseIdx] = exteriorDofIdx_; - } - else { - upIdx_[phaseIdx] = exteriorDofIdx_; - dnIdx_[phaseIdx] = interiorDofIdx_; - } - } - } - - // apply the threshold pressure for the intersection. note that the concept - // of threshold pressure is a quite big hack that only makes sense for ECL - // datasets. (and even there, its physical justification is quite - // questionable IMO.) - if (std::abs(Toolbox::value(pressureDifference_[phaseIdx])) > thpres) { - if (pressureDifference_[phaseIdx] < 0.0) - pressureDifference_[phaseIdx] += thpres; - else - pressureDifference_[phaseIdx] -= thpres; - } - else { - pressureDifference_[phaseIdx] = 0.0; - volumeFlux_[phaseIdx] = 0.0; - continue; - } - - // this is slightly hacky because in the automatic differentiation case, it - // only works for the element centered finite volume method. for ebos this - // does not matter, though. - unsigned upstreamIdx = upstreamIndex_(phaseIdx); - const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); - + const bool upwindIsInterior = (static_cast(upIdx[phaseIdx]) == interiorDofIdx); + const IntensiveQuantities& up = upwindIsInterior ? intQuantsIn : intQuantsEx; // TODO: should the rock compaction transmissibility multiplier be upstreamed // or averaged? all fluids should see the same compaction?! - const Evaluation& transMult = - problem.template rockCompTransMultiplier(up, stencil.globalSpaceIndex(upstreamIdx)); + const Evaluation& transMult = up.rockCompTransMultiplier(); - if (upstreamIdx == interiorDofIdx_) - volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea); + if (upwindIsInterior) + volumeFlux[phaseIdx] = + pressureDifferences[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea); else - volumeFlux_[phaseIdx] = - pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea)); + volumeFlux[phaseIdx] = + pressureDifferences[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea)); } } + template + static void calculatePhasePressureDiff_(short& upIdx, + short& dnIdx, + EvalType& pressureDifference, + const IntensiveQuantities& intQuantsIn, + const IntensiveQuantities& intQuantsEx, + const unsigned phaseIdx, + const unsigned interiorDofIdx, + const unsigned exteriorDofIdx, + const Scalar Vin, + const Scalar Vex, + const unsigned globalIndexIn, + const unsigned globalIndexEx, + const Scalar distZg, + const Scalar thpres + ) + { + + // check shortcut: if the mobility of the phase is zero in the interior as + // well as the exterior DOF, we can skip looking at the phase. + if (intQuantsIn.mobility(phaseIdx) <= 0.0 && + intQuantsEx.mobility(phaseIdx) <= 0.0) + { + upIdx = interiorDofIdx; + dnIdx = exteriorDofIdx; + pressureDifference = 0.0; + return; + } + + // do the gravity correction: compute the hydrostatic pressure for the + // external at the depth of the internal one + const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx); + Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx)); + Evaluation rhoAvg = (rhoIn + rhoEx)/2; + + const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx); + Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx)); + if (enableExtbo) // added stability; particulary useful for solvent migrating in pure water + // where the solvent fraction displays a 0/1 behaviour ... + pressureExterior += Toolbox::value(rhoAvg)*(distZg); + else + pressureExterior += rhoAvg*(distZg); + + pressureDifference = pressureExterior - pressureInterior; + + // decide the upstream index for the phase. for this we make sure that the + // degree of freedom which is regarded upstream if both pressures are equal + // is always the same: if the pressure is equal, the DOF with the lower + // global index is regarded to be the upstream one. + if (pressureDifference > 0.0) { + upIdx = exteriorDofIdx; + dnIdx = interiorDofIdx; + } + else if (pressureDifference < 0.0) { + upIdx = interiorDofIdx; + dnIdx = exteriorDofIdx; + } + else { + // if the pressure difference is zero, we chose the DOF which has the + // larger volume associated to it as upstream DOF + if (Vin > Vex) { + upIdx = interiorDofIdx; + dnIdx = exteriorDofIdx; + } + else if (Vin < Vex) { + upIdx = exteriorDofIdx; + dnIdx = interiorDofIdx; + } + else { + assert(Vin == Vex); + // if the volumes are also equal, we pick the DOF which exhibits the + // smaller global index + if (globalIndexIn < globalIndexEx) { + upIdx = interiorDofIdx; + dnIdx = exteriorDofIdx; + } + else { + upIdx = exteriorDofIdx; + dnIdx = interiorDofIdx; + } + } + } + + // apply the threshold pressure for the intersection. note that the concept + // of threshold pressure is a quite big hack that only makes sense for ECL + // datasets. (and even there, its physical justification is quite + // questionable IMO.) + if (std::abs(Toolbox::value(pressureDifference)) > thpres) { + if (pressureDifference < 0.0) + pressureDifference += thpres; + else + pressureDifference -= thpres; + } + else { + pressureDifference = 0.0; + } + } + +protected: + /*! + * \brief Update the required gradients for interior faces + */ + void calculateGradients_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) + { + Valgrind::SetUndefined(*this); + + volumeAndPhasePressureDifferences(upIdx_ , dnIdx_, volumeFlux_, pressureDifference_, elemCtx, scvfIdx, timeIdx); + } + /*! * \brief Update the required gradients for boundary faces */ @@ -374,7 +419,7 @@ protected: const auto& stencil = elemCtx.stencil(timeIdx); const auto& scvf = stencil.boundaryFace(scvfIdx); - interiorDofIdx_ = scvf.interiorIndex(); + unsigned interiorDofIdx = scvf.interiorIndex(); Scalar trans = problem.transmissibilityBoundary(elemCtx, scvfIdx); Scalar faceArea = scvf.area(); @@ -384,14 +429,14 @@ protected: // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) Scalar g = elemCtx.problem().gravity()[dimWorld - 1]; - const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx); + const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx); // this is quite hacky because the dune grid interface does not provide a // cellCenterDepth() method (so we ask the problem to provide it). The "good" // solution would be to take the Z coordinate of the element centroids, but since // ECL seems to like to be inconsistent on that front, it needs to be done like // here... - Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx_, timeIdx); + Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx); Scalar zEx = scvf.integrationPos()[dimWorld - 1]; // the distances from the DOF's depths. (i.e., the additional depth of the @@ -420,17 +465,17 @@ protected: // global index is regarded to be the upstream one. if (pressureDifference_[phaseIdx] > 0.0) { upIdx_[phaseIdx] = -1; - dnIdx_[phaseIdx] = interiorDofIdx_; + dnIdx_[phaseIdx] = interiorDofIdx; } else { - upIdx_[phaseIdx] = interiorDofIdx_; + upIdx_[phaseIdx] = interiorDofIdx; dnIdx_[phaseIdx] = -1; } Evaluation transModified = trans; - short upstreamIdx = upstreamIndex_(phaseIdx); - if (upstreamIdx == interiorDofIdx_) { + unsigned upstreamIdx = upstreamIndex_(phaseIdx); + if (upstreamIdx == interiorDofIdx) { // this is slightly hacky because in the automatic differentiation case, it // only works for the element centered finite volume method. for ebos this @@ -438,7 +483,8 @@ protected: const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); // deal with water induced rock compaction - transModified *= problem.template rockCompTransMultiplier(up, stencil.globalSpaceIndex(upstreamIdx)); + const double transMult = Toolbox::value(up.rockCompTransMultiplier()); + transModified *= transMult; volumeFlux_[phaseIdx] = pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); @@ -451,7 +497,7 @@ protected: // interior element. TODO: this could probably be done more efficiently const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, - interiorDofIdx_, + interiorDofIdx, /*timeIdx=*/0); typename FluidState::Scalar kr[numPhases]; MaterialLaw::relativePermeabilities(kr, matParams, exFluidState); @@ -491,8 +537,6 @@ private: Evaluation pressureDifference_[numPhases]; // the local indices of the interior and exterior degrees of freedom - unsigned short interiorDofIdx_; - unsigned short exteriorDofIdx_; short upIdx_[numPhases]; short dnIdx_[numPhases]; }; diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 84e27c445..a94d2eabb 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -1438,6 +1438,14 @@ public: return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility; } + /*! + * \brief Direct access to the transmissibility between two elements. + */ + Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const + { + return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx); + } + /*! * \copydoc EclTransmissiblity::diffusivity */ @@ -1534,6 +1542,19 @@ public: Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const { unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + return this->porosity(globalSpaceIdx, timeIdx); + } + + /*! + * \brief Direct indexed access to the porosity. + * + * For the EclProblem, this method is identical to referencePorosity(). The intensive + * quantities object may apply various multipliers (e.g. ones which model rock + * compressibility and water induced rock compaction) to it which depend on the + * current physical conditions. + */ + Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const + { return this->referencePorosity_[timeIdx][globalSpaceIdx]; } @@ -1547,11 +1568,11 @@ public: Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const { unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - return dofCenterDepth(globalSpaceIdx); + return this->dofCenterDepth(globalSpaceIdx); } /*! - * \brief Returns the depth of an degree of freedom [m] + * \brief Direct indexed acces to the depth of an degree of freedom [m] * * For ECL problems this is defined as the average of the depth of an element and is * thus slightly different from the depth of an element's centroid. @@ -1561,7 +1582,6 @@ public: return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx); } - /*! * \copydoc BlackoilProblem::rockCompressibility */ @@ -1580,6 +1600,25 @@ public: return this->rockParams_[tableIdx].compressibility; } + /*! + * Direct access to rock compressibility. + * + * While the above overload could be implemented in terms of this method, + * that would require always looking up the global space index, which + * is not always needed. + */ + Scalar rockCompressibility(unsigned globalSpaceIdx) const + { + if (this->rockParams_.empty()) + return 0.0; + + unsigned tableIdx = 0; + if (!this->rockTableIdx_.empty()) { + tableIdx = this->rockTableIdx_[globalSpaceIdx]; + } + return this->rockParams_[tableIdx].compressibility; + } + /*! * \copydoc BlackoilProblem::rockReferencePressure */ @@ -1598,6 +1637,25 @@ public: return this->rockParams_[tableIdx].referencePressure; } + /*! + * Direct access to rock reference pressure. + * + * While the above overload could be implemented in terms of this method, + * that would require always looking up the global space index, which + * is not always needed. + */ + Scalar rockReferencePressure(unsigned globalSpaceIdx) const + { + if (this->rockParams_.empty()) + return 1e5; + + unsigned tableIdx = 0; + if (!this->rockTableIdx_.empty()) { + tableIdx = this->rockTableIdx_[globalSpaceIdx]; + } + return this->rockParams_[tableIdx].referencePressure; + } + /*! * \copydoc FvBaseMultiPhaseProblem::materialLawParams */ @@ -1606,11 +1664,13 @@ public: unsigned spaceIdx, unsigned timeIdx) const { unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - return materialLawParams(globalSpaceIdx); + return this->materialLawParams(globalSpaceIdx); } const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const - { return materialLawManager_->materialLawParams(globalDofIdx); } + { + return materialLawManager_->materialLawParams(globalDofIdx); + } /*! * \brief Return the parameters for the energy storage law of the rock @@ -1953,14 +2013,21 @@ public: const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + source(rate, globalDofIdx, timeIdx); + } + + void source(RateVector& rate, + unsigned globalDofIdx, + unsigned timeIdx) const { rate = 0.0; - wellModel_.computeTotalRatesForDof(rate, context, spaceIdx, timeIdx); + wellModel_.computeTotalRatesForDof(rate, globalDofIdx); // convert the source term from the total mass rate of the // cell to the one per unit of volume as used by the model. - const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx); for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) { rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx); @@ -1969,12 +2036,11 @@ public: } if (enableAquifers_) - aquiferModel_.addToSource(rate, context, spaceIdx, timeIdx); + aquiferModel_.addToSource(rate, globalDofIdx, timeIdx); // if requested, compensate systematic mass loss for cells which were "well // behaved" in the last time step if (enableDriftCompensation_) { - const auto& intQuants = context.intensiveQuantities(spaceIdx, timeIdx); const auto& simulator = this->simulator(); const auto& model = this->model(); @@ -1982,11 +2048,11 @@ public: // current time step might be shorter than the last one Scalar maxCompensation = 10.0*model.newtonMethod().tolerance(); - Scalar poro = intQuants.referencePorosity(); + Scalar poro = this->porosity(globalDofIdx, timeIdx); Scalar dt = simulator.timeStepSize(); EqVector dofDriftRate = drift_[globalDofIdx]; - dofDriftRate /= dt*context.dofTotalVolume(spaceIdx, timeIdx); + dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx); // compute the weighted total drift rate Scalar totalDriftRate = 0.0; diff --git a/flow/flow_blackoil_tpfa.cpp b/flow/flow_blackoil_tpfa.cpp new file mode 100644 index 000000000..e644171bc --- /dev/null +++ b/flow/flow_blackoil_tpfa.cpp @@ -0,0 +1,25 @@ +/* + 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 3 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 . +*/ + +#include "config.h" + +#include + +int main(int argc, char** argv) +{ + return Opm::flowEbosBlackoilTpfaMainStandalone(argc, argv); +} diff --git a/flow/flow_ebos_blackoil_tpfa.cpp b/flow/flow_ebos_blackoil_tpfa.cpp new file mode 100644 index 000000000..5227d56f8 --- /dev/null +++ b/flow/flow_ebos_blackoil_tpfa.cpp @@ -0,0 +1,77 @@ +/* + 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 3 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 . +*/ +#include "config.h" + +#include + +#include +#include +#include +#include + +#include +#include + +namespace Opm { + namespace Properties { + namespace TTag { + struct EclFlowProblemTPFA { + using InheritsFrom = std::tuple; + }; + } + } +} + +namespace Opm { + namespace Properties { + + template + struct Linearizer { using type = TpfaLinearizer; }; + + template + struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; + + template + struct EnableDiffusion { static constexpr bool value = false; }; + + } +} + + +namespace Opm +{ + +// ----------------- Main program ----------------- +int flowEbosBlackoilTpfaMain(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + resetLocale(); + + FlowMainEbos + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +int flowEbosBlackoilTpfaMainStandalone(int argc, char** argv) +{ + using TypeTag = Properties::TTag::EclFlowProblemTPFA; + auto mainObject = Opm::Main(argc, argv); + return mainObject.runStatic(); +} + +} // namespace Opm diff --git a/flow/flow_ebos_blackoil_tpfa.hpp b/flow/flow_ebos_blackoil_tpfa.hpp new file mode 100644 index 000000000..b262be98d --- /dev/null +++ b/flow/flow_ebos_blackoil_tpfa.hpp @@ -0,0 +1,30 @@ +/* + 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 3 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 . +*/ +#ifndef FLOW_EBOS_BLACKOIL_TPFA_HPP +#define FLOW_EBOS_BLACKOIL_TPFA_HPP + +namespace Opm { + +//! \brief Main function used in flow binary. +int flowEbosBlackoilTpfaMain(int argc, char** argv, bool outputCout, bool outputFiles); + +//! \brief Main function used in flow_brine binary. +int flowEbosBlackoilTpfaMainStandalone(int argc, char** argv); + +} + +#endif // FLOW_EBOS_BLACKOIL_TPFA_HPP diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index fcc997f8c..74180a2d9 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -143,24 +143,31 @@ public: const unsigned timeIdx) { const unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + addToSource(rates, cellIdx, timeIdx); + } + + void addToSource(RateVector& rates, + const unsigned cellIdx, + const unsigned timeIdx) + { + const auto& model = ebos_simulator_.model(); const int idx = this->cellToConnectionIdx_[cellIdx]; if (idx < 0) return; - // We are dereferencing the value of IntensiveQuantities because - // cachedIntensiveQuantities return a const pointer to - // IntensiveQuantities of that particular cell_id - const auto& intQuants = context.intensiveQuantities(spaceIdx, timeIdx); + const auto* intQuantsPtr = model.cachedIntensiveQuantities(cellIdx, timeIdx); + if (intQuantsPtr == nullptr) { + throw std::logic_error("Invalid intensive quantities cache detected in AquiferInterface::addToSource()"); + } // This is the pressure at td + dt - this->updateCellPressure(this->pressure_current_, idx, intQuants); - this->calculateInflowRate(idx, context.simulator()); + this->updateCellPressure(this->pressure_current_, idx, *intQuantsPtr); + this->calculateInflowRate(idx, ebos_simulator_); rates[BlackoilIndices::conti0EqIdx + compIdx_()] - += this->Qai_[idx] / context.dofVolume(spaceIdx, timeIdx); + += this->Qai_[idx] / model.dofTotalVolume(cellIdx); } - std::size_t size() const { return this->connections_.size(); } diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp index 7f17701ff..f06ffe75c 100644 --- a/opm/simulators/aquifers/AquiferNumerical.hpp +++ b/opm/simulators/aquifers/AquiferNumerical.hpp @@ -49,12 +49,14 @@ public: using MaterialLaw = GetPropType; enum { dimWorld = GridView::dimensionworld }; - + enum { numPhases = FluidSystem::numPhases }; static const int numEq = BlackoilIndices::numEq; using Eval = DenseAd::Evaluation; using Toolbox = MathToolbox; + using IntensiveQuantities = GetPropType; + using ExtensiveQuantities = GetPropType; // Constructor AquiferNumerical(const SingleNumericalAquifer& aquifer, const std::unordered_map& cartesian_to_compressed, @@ -250,6 +252,14 @@ private: return sum_pressure_watervolume / sum_watervolume; } + template + double getWaterFlux(const ElemCtx& elem_ctx, unsigned face_idx) const + { + const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0); + const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_())); + return water_flux; + } + double calculateAquiferFluxRate() const { double aquifer_flux = 0.0; @@ -276,8 +286,6 @@ private: if (idx != 0) { continue; } - elem_ctx.updateAllIntensiveQuantities(); - elem_ctx.updateAllExtensiveQuantities(); const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0); // const auto &problem = elem_ctx.problem(); @@ -300,9 +308,10 @@ private: if (this->cell_to_aquifer_cell_idx_[J] > 0) { continue; } - const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0); - const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_())); + elem_ctx.updateAllIntensiveQuantities(); + elem_ctx.updateAllExtensiveQuantities(); + const double water_flux = getWaterFlux(elem_ctx,face_idx); const std::size_t up_id = water_flux >= 0.0 ? i : j; const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0); const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_())); diff --git a/opm/simulators/aquifers/BlackoilAquiferModel.hpp b/opm/simulators/aquifers/BlackoilAquiferModel.hpp index 1e2ae5468..342f614af 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel.hpp @@ -95,6 +95,7 @@ public: // add the water rate due to aquifers to the source term. template void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) const; + void addToSource(RateVector& rates, unsigned globalSpaceIdx, unsigned timeIdx) const; void endIteration(); void endTimeStep(); void endEpisode(); diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index 3ad985fc2..829c5835c 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -126,6 +126,24 @@ BlackoilAquiferModel::addToSource(RateVector& rates, } } +template +void +BlackoilAquiferModel::addToSource(RateVector& rates, + unsigned globalSpaceIdx, + unsigned timeIdx) const +{ + if (aquiferCarterTracyActive()) { + for (auto& aquifer : aquifers_CarterTracy) { + aquifer.addToSource(rates, globalSpaceIdx, timeIdx); + } + } + if (aquiferFetkovichActive()) { + for (auto& aquifer : aquifers_Fetkovich) { + aquifer.addToSource(rates, globalSpaceIdx, timeIdx); + } + } +} + template void BlackoilAquiferModel::endIteration() diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp index de628562c..6d3947730 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp @@ -218,7 +218,8 @@ public: ebosSimulator_.setEpisodeIndex(-1); ebosSimulator_.setEpisodeLength(0.0); ebosSimulator_.setTimeStepSize(0.0); - + // Make cache up to date. No need for updating it in elementCtx. + ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); wellModel_().beginReportStep(timer.currentStepNum()); ebosSimulator_.problem().writeOutput(); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 13cd0d42a..ec36549cd 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -128,8 +128,6 @@ namespace Opm { typedef Dune::FieldVector VectorBlockType; typedef Dune::BlockVector BVector; - // typedef Dune::FieldMatrix MatrixBlockType; - typedef BlackOilPolymerModule PolymerModule; typedef BlackOilMICPModule MICPModule; @@ -211,6 +209,9 @@ namespace Opm { endReportStep(); } + void computeTotalRatesForDof(RateVector& rate, + unsigned globalIdx) const; + template void computeTotalRatesForDof(RateVector& rate, const Context& context, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index ebde65f63..066da1432 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -500,6 +500,22 @@ namespace Opm { } + template + void + BlackoilWellModel:: + computeTotalRatesForDof(RateVector& rate, + unsigned elemIdx) const + { + rate = 0; + + if (!is_cell_perforated_[elemIdx]) + return; + + for (const auto& well : well_container_) + well->addCellRates(rate, elemIdx); + } + + template template void