Merge pull request #3994 from akva2/aquifer_use_through_if

BlackOilAquiferModel: use implementations through interface
This commit is contained in:
Kai Bao 2022-08-11 11:29:10 +02:00 committed by GitHub
commit de09fbaea7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 505 additions and 553 deletions

View File

@ -247,9 +247,10 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/timestepping/SimulatorReport.hpp
opm/simulators/wells/SegmentState.hpp
opm/simulators/wells/WellContainer.hpp
opm/simulators/aquifers/AquiferInterface.hpp
opm/simulators/aquifers/AquiferAnalytical.hpp
opm/simulators/aquifers/AquiferCarterTracy.hpp
opm/simulators/aquifers/AquiferFetkovich.hpp
opm/simulators/aquifers/AquiferInterface.hpp
opm/simulators/aquifers/AquiferNumerical.hpp
opm/simulators/aquifers/BlackoilAquiferModel.hpp
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

View File

@ -0,0 +1,416 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2017 IRIS
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 OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
#define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
#include <opm/simulators/aquifers/AquiferInterface.hpp>
#include <opm/common/utility/numeric/linearInterpolation.hpp>
#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
#include <opm/output/data/Aquifer.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <numeric>
#include <unordered_map>
#include <vector>
namespace Opm
{
template <typename TypeTag>
class AquiferAnalytical : public AquiferInterface<TypeTag>
{
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
static constexpr int numEq = BlackoilIndices::numEq;
using Scalar = double;
using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
using FluidState = BlackOilFluidState<Eval,
FluidSystem,
enableTemperature,
enableEnergy,
BlackoilIndices::gasEnabled,
enableEvaporation,
enableBrine,
enableSaltPrecipitation,
BlackoilIndices::numPhases>;
// Constructor
AquiferAnalytical(int aqID,
const std::vector<Aquancon::AquancCell>& connections,
const Simulator& ebosSimulator)
: AquiferInterface<TypeTag>(aqID, ebosSimulator)
, connections_(connections)
{
}
// Destructor
virtual ~AquiferAnalytical()
{
}
void initFromRestart(const data::Aquifers& aquiferSoln) override
{
auto xaqPos = aquiferSoln.find(this->aquiferID());
if (xaqPos == aquiferSoln.end())
return;
this->assignRestartData(xaqPos->second);
this->W_flux_ = xaqPos->second.volume;
this->pa0_ = xaqPos->second.initPressure;
this->solution_set_from_restart_ = true;
}
void initialSolutionApplied() override
{
initQuantities();
}
void beginTimeStep() override
{
ElementContext elemCtx(this->ebos_simulator_);
auto elemIt = this->ebos_simulator_.gridView().template begin<0>();
const auto& elemEndIt = this->ebos_simulator_.gridView().template end<0>();
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
const int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updateIntensiveQuantities(0);
const auto& iq = elemCtx.intensiveQuantities(0, 0);
pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
}
OPM_END_PARALLEL_TRY_CATCH("AquiferAnalytical::beginTimeStep() failed: ",
this->ebos_simulator_.vanguard().grid().comm());
}
void addToSource(RateVector& rates,
const unsigned cellIdx,
const unsigned timeIdx) override
{
const auto& model = this->ebos_simulator_.model();
const int idx = this->cellToConnectionIdx_[cellIdx];
if (idx < 0)
return;
const auto* intQuantsPtr = model.cachedIntensiveQuantities(cellIdx, timeIdx);
if (intQuantsPtr == nullptr) {
throw std::logic_error("Invalid intensive quantities cache detected in AquiferAnalytical::addToSource()");
}
// This is the pressure at td + dt
this->updateCellPressure(this->pressure_current_, idx, *intQuantsPtr);
this->calculateInflowRate(idx, this->ebos_simulator_);
rates[BlackoilIndices::conti0EqIdx + compIdx_()]
+= this->Qai_[idx] / model.dofTotalVolume(cellIdx);
}
std::size_t size() const
{
return this->connections_.size();
}
protected:
virtual void assignRestartData(const data::AquiferData& xaq) = 0;
virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
virtual void calculateAquiferCondition() = 0;
virtual void calculateAquiferConstants() = 0;
virtual Scalar aquiferDepth() const = 0;
Scalar gravity_() const
{
return this->ebos_simulator_.problem().gravity()[2];
}
int compIdx_() const
{
if (this->co2store_())
return FluidSystem::oilCompIdx;
return FluidSystem::waterCompIdx;
}
void initQuantities()
{
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0
if (!this->solution_set_from_restart_) {
W_flux_ = Scalar{0};
}
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections
// function
initializeConnections();
calculateAquiferCondition();
calculateAquiferConstants();
pressure_previous_.resize(this->connections_.size(), Scalar{0});
pressure_current_.resize(this->connections_.size(), Scalar{0});
Qai_.resize(this->connections_.size(), Scalar{0});
}
void updateCellPressure(std::vector<Eval>& pressure_water,
const int idx,
const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
}
void updateCellPressure(std::vector<Scalar>& pressure_water,
const int idx,
const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
}
void initializeConnections()
{
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), Scalar{0});
// Translate the C face tag into the enum used by opm-parser's TransMult class
FaceDir::DirEnum faceDirection;
bool has_active_connection_on_proc = false;
// denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas{0};
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
for (std::size_t idx = 0; idx < this->size(); ++idx) {
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
auto elemIt = gridView.template begin</*codim=*/ 0>();
if (cell_index > 0)
std::advance(elemIt, cell_index);
//the global_index is not part of this grid
if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
continue;
has_active_connection_on_proc = true;
this->cellToConnectionIdx_[cell_index] = idx;
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get areas for all connections
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
// only deal with connections given by the aquifer
if( idx < 0)
continue;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// only deal with grid boundaries
if (!intersection.boundary())
continue;
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = FaceDir::XMinus;
break;
case 1:
faceDirection = FaceDir::XPlus;
break;
case 2:
faceDirection = FaceDir::YMinus;
break;
case 3:
faceDirection = FaceDir::YPlus;
break;
case 4:
faceDirection = FaceDir::ZMinus;
break;
case 5:
faceDirection = FaceDir::ZPlus;
break;
default:
OPM_THROW(std::logic_error,
"Internal error in initialization of aquifer.");
}
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
break;
}
}
denom_face_areas += this->faceArea_connected_.at(idx);
}
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&denom_face_areas, 1);
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (std::size_t idx = 0; idx < this->size(); ++idx) {
// Protect against division by zero NaNs.
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
? Scalar{0}
: this->faceArea_connected_.at(idx) / denom_face_areas;
}
if (this->solution_set_from_restart_) {
this->rescaleProducedVolume(has_active_connection_on_proc);
}
}
void rescaleProducedVolume(const bool has_active_connection_on_proc)
{
// Needed in parallel restart to approximate influence of aquifer
// being "owned" by a subset of the parallel processes. If the
// aquifer is fully owned by a single process--i.e., if all cells
// connecting to the aquifer are on a single process--then this_area
// is tot_area on that process and zero elsewhere.
const auto this_area = has_active_connection_on_proc
? std::accumulate(this->alphai_.begin(),
this->alphai_.end(),
Scalar{0})
: Scalar{0};
const auto tot_area = this->ebos_simulator_.vanguard()
.grid().comm().sum(this_area);
this->W_flux_ *= this_area / tot_area;
}
// This function is for calculating the aquifer properties from equilibrium state with the reservoir
Scalar calculateReservoirEquilibrium()
{
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(this->ebos_simulator_);
const auto& gridView = this->ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
const auto cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto idx = this->cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
const auto water_density = fs.density(this->phaseIdx_());
const auto gdz =
this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
pw_aquifer.push_back(this->alphai_[idx] *
(water_pressure_reservoir - water_density.value()*gdz));
}
// We take the average of the calculated equilibrium pressures.
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
Scalar vals[2];
vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
comm.sum(vals, 2);
return vals[1] / vals[0];
}
const std::vector<Aquancon::AquancCell> connections_;
// Grid variables
std::vector<Scalar> faceArea_connected_;
std::vector<int> cellToConnectionIdx_;
// Quantities at each grid id
std::vector<Scalar> cell_depth_;
std::vector<Scalar> pressure_previous_;
std::vector<Eval> pressure_current_;
std::vector<Eval> Qai_;
std::vector<Scalar> alphai_;
Scalar Tc_{}; // Time constant
Scalar pa0_{}; // initial aquifer pressure
Scalar rhow_{};
Eval W_flux_;
bool solution_set_from_restart_ {false};
};
} // namespace Opm
#endif

View File

@ -21,7 +21,7 @@
#ifndef OPM_AQUIFERCT_HEADER_INCLUDED
#define OPM_AQUIFERCT_HEADER_INCLUDED
#include <opm/simulators/aquifers/AquiferInterface.hpp>
#include <opm/simulators/aquifers/AquiferAnalytical.hpp>
#include <opm/output/data/Aquifer.hpp>
@ -34,10 +34,10 @@ namespace Opm
{
template <typename TypeTag>
class AquiferCarterTracy : public AquiferInterface<TypeTag>
class AquiferCarterTracy : public AquiferAnalytical<TypeTag>
{
public:
typedef AquiferInterface<TypeTag> Base;
using Base = AquiferAnalytical<TypeTag>;
using typename Base::BlackoilIndices;
using typename Base::ElementContext;
@ -236,6 +236,7 @@ protected:
return this->aquct_data_.datum_depth;
}
}; // class AquiferCarterTracy
} // namespace Opm
#endif

View File

@ -21,7 +21,7 @@ along with OPM. If not, see <http://www.gnu.org/licenses/>.
#ifndef OPM_AQUIFETP_HEADER_INCLUDED
#define OPM_AQUIFETP_HEADER_INCLUDED
#include <opm/simulators/aquifers/AquiferInterface.hpp>
#include <opm/simulators/aquifers/AquiferAnalytical.hpp>
#include <opm/output/data/Aquifer.hpp>
@ -33,11 +33,11 @@ namespace Opm
{
template <typename TypeTag>
class AquiferFetkovich : public AquiferInterface<TypeTag>
class AquiferFetkovich : public AquiferAnalytical<TypeTag>
{
public:
typedef AquiferInterface<TypeTag> Base;
using Base = AquiferAnalytical<TypeTag>;
using typename Base::BlackoilIndices;
using typename Base::ElementContext;
@ -171,5 +171,7 @@ protected:
return this->aqufetp_data_.datum_depth;
}
}; // Class AquiferFetkovich
} // namespace Opm
#endif

View File

@ -22,119 +22,38 @@
#ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED
#define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
#include <opm/common/utility/numeric/linearInterpolation.hpp>
#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
#include <opm/output/data/Aquifer.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <numeric>
#include <unordered_map>
#include <vector>
namespace Opm
{
template <typename TypeTag>
class AquiferInterface
{
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
static const int numEq = BlackoilIndices::numEq;
typedef double Scalar;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
typedef BlackOilFluidState<Eval,
FluidSystem,
enableTemperature,
enableEnergy,
BlackoilIndices::gasEnabled,
enableEvaporation,
enableBrine,
enableSaltPrecipitation,
BlackoilIndices::numPhases>
FluidState;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
// Constructor
AquiferInterface(int aqID,
const std::vector<Aquancon::AquancCell>& connections,
const Simulator& ebosSimulator)
: aquiferID_(aqID)
, connections_(connections)
, ebos_simulator_(ebosSimulator)
{
}
// Destructor
virtual ~AquiferInterface()
{
}
virtual ~AquiferInterface() = default;
void initFromRestart(const data::Aquifers& aquiferSoln)
{
auto xaqPos = aquiferSoln.find(this->aquiferID());
if (xaqPos == aquiferSoln.end())
return;
virtual void initFromRestart(const data::Aquifers& aquiferSoln) = 0;
this->assignRestartData(xaqPos->second);
virtual void initialSolutionApplied() = 0;
this->W_flux_ = xaqPos->second.volume;
this->pa0_ = xaqPos->second.initPressure;
this->solution_set_from_restart_ = true;
}
virtual void beginTimeStep() = 0;
virtual void endTimeStep() = 0;
void initialSolutionApplied()
{
initQuantities();
}
void beginTimeStep()
{
ElementContext elemCtx(ebos_simulator_);
auto elemIt = ebos_simulator_.gridView().template begin<0>();
const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
const int idx = cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updateIntensiveQuantities(0);
const auto& iq = elemCtx.intensiveQuantities(0, 0);
pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_()));
}
OPM_END_PARALLEL_TRY_CATCH("AquiferInterface::beginTimeStep() failed: ", ebos_simulator_.vanguard().grid().comm());
}
virtual data::AquiferData aquiferData() const = 0;
template <class Context>
void addToSource(RateVector& rates,
@ -146,299 +65,30 @@ public:
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;
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();
}
virtual void addToSource(RateVector& rates,
const unsigned cellIdx,
const unsigned timeIdx) = 0;
int aquiferID() const { return this->aquiferID_; }
protected:
inline Scalar gravity_() const
{
return ebos_simulator_.problem().gravity()[2];
}
inline bool co2store_() const
bool co2store_() const
{
return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
}
inline int phaseIdx_() const
int phaseIdx_() const
{
if(co2store_())
if (co2store_())
return FluidSystem::oilPhaseIdx;
return FluidSystem::waterPhaseIdx;
}
inline int compIdx_() const
{
if(co2store_())
return FluidSystem::oilCompIdx;
return FluidSystem::waterCompIdx;
}
inline void initQuantities()
{
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0
if (!this->solution_set_from_restart_) {
W_flux_ = Scalar{0};
}
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections
// function
initializeConnections();
calculateAquiferCondition();
calculateAquiferConstants();
pressure_previous_.resize(this->connections_.size(), Scalar{0});
pressure_current_.resize(this->connections_.size(), Scalar{0});
Qai_.resize(this->connections_.size(), Scalar{0});
}
inline void
updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(phaseIdx_());
}
inline void
updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
{
const auto& fs = intQuants.fluidState();
pressure_water.at(idx) = fs.pressure(phaseIdx_()).value();
}
virtual void endTimeStep() = 0;
const int aquiferID_{};
const std::vector<Aquancon::AquancCell> connections_;
const Simulator& ebos_simulator_;
// Grid variables
std::vector<Scalar> faceArea_connected_;
std::vector<int> cellToConnectionIdx_;
// Quantities at each grid id
std::vector<Scalar> cell_depth_;
std::vector<Scalar> pressure_previous_;
std::vector<Eval> pressure_current_;
std::vector<Eval> Qai_;
std::vector<Scalar> alphai_;
Scalar Tc_{}; // Time constant
Scalar pa0_{}; // initial aquifer pressure
Scalar rhow_{};
Eval W_flux_;
bool solution_set_from_restart_ {false};
void initializeConnections()
{
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), Scalar{0});
// Translate the C face tag into the enum used by opm-parser's TransMult class
FaceDir::DirEnum faceDirection;
bool has_active_connection_on_proc = false;
// denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas{0};
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
for (std::size_t idx = 0; idx < this->size(); ++idx) {
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
auto elemIt = gridView.template begin</*codim=*/ 0>();
if (cell_index > 0)
std::advance(elemIt, cell_index);
//the global_index is not part of this grid
if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
continue;
has_active_connection_on_proc = true;
this->cellToConnectionIdx_[cell_index] = idx;
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get areas for all connections
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
// only deal with connections given by the aquifer
if( idx < 0)
continue;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// only deal with grid boundaries
if (!intersection.boundary())
continue;
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = FaceDir::XMinus;
break;
case 1:
faceDirection = FaceDir::XPlus;
break;
case 2:
faceDirection = FaceDir::YMinus;
break;
case 3:
faceDirection = FaceDir::YPlus;
break;
case 4:
faceDirection = FaceDir::ZMinus;
break;
case 5:
faceDirection = FaceDir::ZPlus;
break;
default:
OPM_THROW(std::logic_error,
"Internal error in initialization of aquifer.");
}
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
break;
}
}
denom_face_areas += this->faceArea_connected_.at(idx);
}
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&denom_face_areas, 1);
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (std::size_t idx = 0; idx < this->size(); ++idx) {
// Protect against division by zero NaNs.
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
? Scalar{0}
: this->faceArea_connected_.at(idx) / denom_face_areas;
}
if (this->solution_set_from_restart_) {
this->rescaleProducedVolume(has_active_connection_on_proc);
}
}
void rescaleProducedVolume(const bool has_active_connection_on_proc)
{
// Needed in parallel restart to approximate influence of aquifer
// being "owned" by a subset of the parallel processes. If the
// aquifer is fully owned by a single process--i.e., if all cells
// connecting to the aquifer are on a single process--then this_area
// is tot_area on that process and zero elsewhere.
const auto this_area = has_active_connection_on_proc
? std::accumulate(this->alphai_.begin(),
this->alphai_.end(),
Scalar{0})
: Scalar{0};
const auto tot_area = this->ebos_simulator_.vanguard()
.grid().comm().sum(this_area);
this->W_flux_ *= this_area / tot_area;
}
virtual void assignRestartData(const data::AquiferData& xaq) = 0;
virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
virtual void calculateAquiferCondition() = 0;
virtual void calculateAquiferConstants() = 0;
virtual Scalar aquiferDepth() const = 0;
// This function is for calculating the aquifer properties from equilibrium state with the reservoir
virtual Scalar calculateReservoirEquilibrium()
{
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(this->ebos_simulator_);
const auto& gridView = this->ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
const auto cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto idx = this->cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(phaseIdx_()).value();
const auto water_density = fs.density(phaseIdx_());
const auto gdz =
this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
pw_aquifer.push_back(this->alphai_[idx] *
(water_pressure_reservoir - water_density.value()*gdz));
}
// We take the average of the calculated equilibrium pressures.
const auto& comm = ebos_simulator_.vanguard().grid().comm();
Scalar vals[2];
vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
comm.sum(vals, 2);
return vals[1] / vals[0];
}
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
};
} // namespace Opm
#endif

View File

@ -25,6 +25,7 @@
#include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
#include <opm/simulators/aquifers/AquiferInterface.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <algorithm>
@ -37,33 +38,33 @@
namespace Opm
{
template <typename TypeTag>
class AquiferNumerical
class AquiferNumerical : public AquiferInterface<TypeTag>
{
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = FluidSystem::numPhases };
static const int numEq = BlackoilIndices::numEq;
static constexpr 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>;
using typename AquiferInterface<TypeTag>::RateVector;
// Constructor
AquiferNumerical(const SingleNumericalAquifer& aquifer,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebos_simulator,
const int* global_cell)
: id_ (aquifer.id())
, ebos_simulator_ (ebos_simulator)
: AquiferInterface<TypeTag>(aquifer.id(), ebos_simulator)
, flux_rate_ (0.0)
, cumulative_flux_(0.0)
, global_cell_ (global_cell)
@ -88,7 +89,7 @@ public:
}
}
void initFromRestart(const data::Aquifers& aquiferSoln)
void initFromRestart(const data::Aquifers& aquiferSoln) override
{
auto xaqPos = aquiferSoln.find(this->aquiferID());
if (xaqPos == aquiferSoln.end())
@ -107,14 +108,17 @@ public:
this->solution_set_from_restart_ = true;
}
void endTimeStep()
void beginTimeStep() override {}
void addToSource(RateVector&, const unsigned, const unsigned) override {}
void endTimeStep() override
{
this->pressure_ = this->calculateAquiferPressure();
this->flux_rate_ = this->calculateAquiferFluxRate();
this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
}
data::AquiferData aquiferData() const
data::AquiferData aquiferData() const override
{
data::AquiferData data;
data.aquiferID = this->aquiferID();
@ -128,7 +132,7 @@ public:
return data;
}
void initialSolutionApplied()
void initialSolutionApplied() override
{
if (this->solution_set_from_restart_) {
return;
@ -139,38 +143,7 @@ public:
this->cumulative_flux_ = 0.;
}
int aquiferID() const
{
return static_cast<int>(this->id_);
}
private:
const std::size_t id_;
const Simulator& ebos_simulator_;
double flux_rate_; // aquifer influx rate
double cumulative_flux_; // cumulative aquifer influx
const int* global_cell_; // mapping to global index
std::vector<double> init_pressure_{};
double pressure_; // aquifer pressure
bool solution_set_from_restart_ {false};
bool connects_to_reservoir_ {false};
// TODO: maybe unordered_map can also do the work to save memory?
std::vector<int> cell_to_aquifer_cell_idx_;
inline bool co2store_() const
{
return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
}
inline int phaseIdx_() const
{
if(co2store_())
return FluidSystem::oilPhaseIdx;
return FluidSystem::waterPhaseIdx;
}
void checkConnectsToReservoir()
{
ElementContext elem_ctx(this->ebos_simulator_);
@ -325,6 +298,19 @@ private:
return aquifer_flux;
}
double flux_rate_; // aquifer influx rate
double cumulative_flux_; // cumulative aquifer influx
const int* global_cell_; // mapping to global index
std::vector<double> init_pressure_{};
double pressure_; // aquifer pressure
bool solution_set_from_restart_ {false};
bool connects_to_reservoir_ {false};
// TODO: maybe unordered_map can also do the work to save memory?
std::vector<int> cell_to_aquifer_cell_idx_;
};
} // namespace Opm
#endif

View File

@ -113,28 +113,17 @@ protected:
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
typedef AquiferCarterTracy<TypeTag> AquiferCarterTracy_object;
typedef AquiferFetkovich<TypeTag> AquiferFetkovich_object;
Simulator& simulator_;
// TODO: probably we can use one variable to store both types of aquifers, because
// they share the base class
mutable std::vector<AquiferCarterTracy_object> aquifers_CarterTracy;
mutable std::vector<AquiferFetkovich_object> aquifers_Fetkovich;
std::vector<AquiferNumerical<TypeTag>> aquifers_numerical;
std::vector<std::unique_ptr<AquiferInterface<TypeTag>>> aquifers;
// This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
void init();
bool aquiferActive() const;
bool aquiferCarterTracyActive() const;
bool aquiferFetkovichActive() const;
bool aquiferNumericalActive() const;
};
} // namespace Opm
#include "BlackoilAquiferModel_impl.hpp"
#endif

View File

@ -19,7 +19,6 @@
*/
#include <opm/grid/utility/cartesianToCompressed.hpp>
#include "BlackoilAquiferModel.hpp"
namespace Opm
{
@ -39,45 +38,16 @@ template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::initialSolutionApplied()
{
if (aquiferCarterTracyActive()) {
for (auto& aquifer : aquifers_CarterTracy) {
aquifer.initialSolutionApplied();
}
}
if (aquiferFetkovichActive()) {
for (auto& aquifer : aquifers_Fetkovich) {
aquifer.initialSolutionApplied();
}
}
if (this->aquiferNumericalActive()) {
for (auto& aquifer : this->aquifers_numerical) {
aquifer.initialSolutionApplied();
}
}
for (auto& aquifer : aquifers)
aquifer->initialSolutionApplied();
}
template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln)
{
if (this->aquiferCarterTracyActive()) {
for (auto& aquifer : this->aquifers_CarterTracy) {
aquifer.initFromRestart(aquiferSoln);
}
}
if (this->aquiferFetkovichActive()) {
for (auto& aquifer : this->aquifers_Fetkovich) {
aquifer.initFromRestart(aquiferSoln);
}
}
if (this->aquiferNumericalActive()) {
for (auto& aquifer : this->aquifers_numerical) {
aquifer.initFromRestart(aquiferSoln);
}
}
for (auto& aquifer : this->aquifers)
aquifer->initFromRestart(aquiferSoln);
}
template <typename TypeTag>
@ -94,16 +64,8 @@ template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::beginTimeStep()
{
if (aquiferCarterTracyActive()) {
for (auto& aquifer : aquifers_CarterTracy) {
aquifer.beginTimeStep();
}
}
if (aquiferFetkovichActive()) {
for (auto& aquifer : aquifers_Fetkovich) {
aquifer.beginTimeStep();
}
}
for (auto& aquifer : aquifers)
aquifer->beginTimeStep();
}
template <typename TypeTag>
@ -114,16 +76,8 @@ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
unsigned spaceIdx,
unsigned timeIdx) const
{
if (aquiferCarterTracyActive()) {
for (auto& aquifer : aquifers_CarterTracy) {
aquifer.addToSource(rates, context, spaceIdx, timeIdx);
}
}
if (aquiferFetkovichActive()) {
for (auto& aquifer : aquifers_Fetkovich) {
aquifer.addToSource(rates, context, spaceIdx, timeIdx);
}
}
for (auto& aquifer : aquifers)
aquifer->addToSource(rates, context, spaceIdx, timeIdx);
}
template <typename TypeTag>
@ -132,16 +86,8 @@ 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);
}
}
for (auto& aquifer : aquifers)
aquifer->addToSource(rates, globalSpaceIdx, timeIdx);
}
template <typename TypeTag>
@ -153,21 +99,12 @@ template <typename TypeTag>
void
BlackoilAquiferModel<TypeTag>::endTimeStep()
{
if (aquiferCarterTracyActive()) {
for (auto& aquifer : aquifers_CarterTracy) {
aquifer.endTimeStep();
}
}
if (aquiferFetkovichActive()) {
for (auto& aquifer : aquifers_Fetkovich) {
aquifer.endTimeStep();
}
}
if (aquiferNumericalActive()) {
for (auto& aquifer : this->aquifers_numerical) {
aquifer.endTimeStep();
for (auto& aquifer : aquifers) {
aquifer->endTimeStep();
using NumAq = AquiferNumerical<TypeTag>;
NumAq* num = dynamic_cast<NumAq*>(aquifer.get());
if (num)
this->simulator_.vanguard().grid().comm().barrier();
}
}
}
@ -207,13 +144,15 @@ BlackoilAquiferModel<TypeTag>::init()
const auto& connections = aquifer.connections();
for (const auto& aq : aquifer.ct()) {
aquifers_CarterTracy.emplace_back(connections[aq.aquiferID],
this->simulator_, aq);
auto aqf = std::make_unique<AquiferCarterTracy<TypeTag>>(connections[aq.aquiferID],
this->simulator_, aq);
aquifers.push_back(std::move(aqf));
}
for (const auto& aq : aquifer.fetp()) {
aquifers_Fetkovich.emplace_back(connections[aq.aquiferID],
this->simulator_, aq);
auto aqf = std::make_unique<AquiferFetkovich<TypeTag>>(connections[aq.aquiferID],
this->simulator_, aq);
aquifers.push_back(std::move(aqf));
}
if (aquifer.hasNumericalAquifer()) {
@ -223,55 +162,23 @@ BlackoilAquiferModel<TypeTag>::init()
const std::unordered_map<int, int> cartesian_to_compressed = cartesianToCompressed(number_of_cells,
global_cell);
for ([[maybe_unused]]const auto& [id, aqu] : num_aquifers) {
this->aquifers_numerical.emplace_back(aqu,
cartesian_to_compressed, this->simulator_, global_cell);
auto aqf = std::make_unique<AquiferNumerical<TypeTag>>(aqu,
cartesian_to_compressed,
this->simulator_,
global_cell);
aquifers.push_back(std::move(aqf));
}
}
}
template <typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferCarterTracyActive() const
{
return !aquifers_CarterTracy.empty();
}
template <typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferFetkovichActive() const
{
return !aquifers_Fetkovich.empty();
}
template<typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferNumericalActive() const
{
return !(this->aquifers_numerical.empty());
}
template<typename TypeTag>
data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const
{
data::Aquifers data;
if (this->aquiferCarterTracyActive()) {
for (const auto& aqu : this->aquifers_CarterTracy) {
data.insert_or_assign(aqu.aquiferID(), aqu.aquiferData());
}
}
if (this->aquiferFetkovichActive()) {
for (const auto& aqu : this->aquifers_Fetkovich) {
data.insert_or_assign(aqu.aquiferID(), aqu.aquiferData());
}
}
if (this->aquiferNumericalActive()) {
for (const auto& aqu : this->aquifers_numerical) {
data.insert_or_assign(aqu.aquiferID(), aqu.aquiferData());
}
}
for (const auto& aqu : this->aquifers)
data.insert_or_assign(aqu->aquiferID(), aqu->aquiferData());
return data;
}
} // namespace Opm