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

Merge pull request #3988 from atgeirr/faster_assembly_minimal

Browse Source 
Faster assembly minimal
This commit is contained in:
Atgeirr Flø Rasmussen 2022-08-10 11:20:59 +02:00 committed by GitHub
commit a242dadb6f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
13 changed files with 445 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -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})

196
ebos/eclfluxmodule.hh
View File

@ -100,6 +100,7 @@ class EclTransExtensiveQuantities
{
using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@ -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,35 +237,88 @@ 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;
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;
}
const bool upwindIsInterior = (static_cast<unsigned>(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 = up.rockCompTransMultiplier();
if (upwindIsInterior)
volumeFlux[phaseIdx] =
pressureDifferences[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea);
else
volumeFlux[phaseIdx] =
pressureDifferences[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea));
}
}
template<class EvalType>
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_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
pressureDifference_[phaseIdx] = 0.0;
volumeFlux_[phaseIdx] = 0.0;
continue;
upIdx = interiorDofIdx;
dnIdx = exteriorDofIdx;
pressureDifference = 0.0;
return;
}
// do the gravity correction: compute the hydrostatic pressure for the
@ -274,48 +331,46 @@ protected:
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);
pressureExterior += Toolbox::value(rhoAvg)*(distZg);
else
pressureExterior += rhoAvg*(distZ*g);
pressureExterior += rhoAvg*(distZg);
pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
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_[phaseIdx] > 0.0) {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
if (pressureDifference > 0.0) {
upIdx = exteriorDofIdx;
dnIdx = interiorDofIdx;
}
else if (pressureDifference_[phaseIdx] < 0.0) {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
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
Scalar Vin = elemCtx.dofVolume(interiorDofIdx_, /*timeIdx=*/0);
Scalar Vex = elemCtx.dofVolume(exteriorDofIdx_, /*timeIdx=*/0);
if (Vin > Vex) {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
upIdx = interiorDofIdx;
dnIdx = exteriorDofIdx;
}
else if (Vin < Vex) {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
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 (I < J) {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
if (globalIndexIn < globalIndexEx) {
upIdx = interiorDofIdx;
dnIdx = exteriorDofIdx;
}
else {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
upIdx = exteriorDofIdx;
dnIdx = interiorDofIdx;
}
}
}
@ -324,36 +379,26 @@ protected:
// 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;
if (std::abs(Toolbox::value(pressureDifference)) > thpres) {
if (pressureDifference < 0.0)
pressureDifference += thpres;
else
pressureDifference_[phaseIdx] -= thpres;
pressureDifference -= thpres;
}
else {
pressureDifference_[phaseIdx] = 0.0;
volumeFlux_[phaseIdx] = 0.0;
continue;
pressureDifference = 0.0;
}
}
// 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);
protected:
/*!
* \brief Update the required gradients for interior faces
*/
void calculateGradients_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{
Valgrind::SetUndefined(*this);
// TODO: should the rock compaction transmissibility multiplier be upstreamed
// or averaged? all fluids should see the same compaction?!
const Evaluation& transMult =
problem.template rockCompTransMultiplier<Evaluation>(up, stencil.globalSpaceIndex(upstreamIdx));
if (upstreamIdx == interiorDofIdx_)
volumeFlux_[phaseIdx] =
pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea);
else
volumeFlux_[phaseIdx] =
pressureDifference_[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea));
}
volumeAndPhasePressureDifferences(upIdx_ , dnIdx_, volumeFlux_, pressureDifference_, elemCtx, scvfIdx, timeIdx);
}
/*!
@ -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<double>(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];
};

88
ebos/eclproblem.hh
View File

@ -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;

25
flow/flow_blackoil_tpfa.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <flow/flow_ebos_blackoil_tpfa.hpp>
int main(int argc, char** argv)
{
return Opm::flowEbosBlackoilTpfaMainStandalone(argc, argv);
}

77
flow/flow_ebos_blackoil_tpfa.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <flow/flow_ebos_blackoil_tpfa.hpp>
#include <opm/material/common/ResetLocale.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/simulators/flow/Main.hpp>
#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
#include <opm/models/discretization/common/tpfalinearizer.hh>
namespace Opm {
namespace Properties {
namespace TTag {
struct EclFlowProblemTPFA {
using InheritsFrom = std::tuple<EclFlowProblem>;
};
}
}
}
namespace Opm {
namespace Properties {
template<class TypeTag>
struct Linearizer<TypeTag, TTag::EclFlowProblemTPFA> { using type = TpfaLinearizer<TypeTag>; };
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::EclFlowProblemTPFA> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
template<class TypeTag>
struct EnableDiffusion<TypeTag, TTag::EclFlowProblemTPFA> { 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<Properties::TTag::EclFlowProblemTPFA>
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<TypeTag>();
}
} // namespace Opm

30
flow/flow_ebos_blackoil_tpfa.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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

29
opm/simulators/aquifers/AquiferInterface.hpp
View File

@ -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);
// This is the pressure at td + dt
this->updateCellPressure(this->pressure_current_, idx, intQuants);
this->calculateInflowRate(idx, context.simulator());
rates[BlackoilIndices::conti0EqIdx + compIdx_()]
+= this->Qai_[idx] / context.dofVolume(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, *intQuantsPtr);
this->calculateInflowRate(idx, ebos_simulator_);
rates[BlackoilIndices::conti0EqIdx + compIdx_()]
+= this->Qai_[idx] / model.dofTotalVolume(cellIdx);
}
std::size_t size() const {
return this->connections_.size();
}

19
opm/simulators/aquifers/AquiferNumerical.hpp
View File

@ -49,12 +49,14 @@ public:
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = FluidSystem::numPhases };
static const int numEq = BlackoilIndices::numEq;
using Eval = DenseAd::Evaluation<double, numEq>;
using Toolbox = MathToolbox<Eval>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
// Constructor
AquiferNumerical(const SingleNumericalAquifer& aquifer,
const std::unordered_map<int, int>& cartesian_to_compressed,
@ -250,6 +252,14 @@ private:
return sum_pressure_watervolume / sum_watervolume;
}
template <class ElemCtx>
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_()));

1
opm/simulators/aquifers/BlackoilAquiferModel.hpp
View File

@ -95,6 +95,7 @@ public:
// add the water rate due to aquifers to the source term.
template <class Context>
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();

18
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
View File

@ -126,6 +126,24 @@ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
}
}
template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::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 <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::endIteration()

3
opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -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();

5
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -128,8 +128,6 @@ namespace Opm {
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::BlockVector<VectorBlockType> BVector;
// typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef BlackOilPolymerModule<TypeTag> PolymerModule;
typedef BlackOilMICPModule<TypeTag> MICPModule;
@ -211,6 +209,9 @@ namespace Opm {
endReportStep();
}
void computeTotalRatesForDof(RateVector& rate,
unsigned globalIdx) const;
template <class Context>
void computeTotalRatesForDof(RateVector& rate,
const Context& context,

16
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -500,6 +500,22 @@ namespace Opm {
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
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<typename TypeTag>
template <class Context>
void