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 #3039 from GitPaean/new_numerical_aquifer

Browse Source 
New numerical aquifer
This commit is contained in:
Bård Skaflestad 2021-03-16 16:30:53 +01:00 committed by GitHub
commit 0667497108
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
11 changed files with 438 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -161,6 +161,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/aquifers/AquiferInterface.hpp
opm/simulators/aquifers/AquiferCarterTracy.hpp
opm/simulators/aquifers/AquiferFetkovich.hpp
opm/simulators/aquifers/AquiferNumerical.hpp
opm/simulators/aquifers/BlackoilAquiferModel.hpp
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
opm/simulators/linalg/bda/BdaBridge.hpp

32
compareECLFiles.cmake
View File

@ -283,6 +283,22 @@ add_test_compareECLFiles(CASENAME fetkovich_2d
REL_TOL ${rel_tol}
DIR aquifer-fetkovich)
add_test_compareECLFiles(CASENAME numerical_aquifer_3d_2aqu
FILENAME 3D_2AQU_NUM
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${rel_tol}
DIR aquifer-num
TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
add_test_compareECLFiles(CASENAME numerical_aquifer_3d_1aqu
FILENAME 3D_1AQU_3CELLS
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${rel_tol}
DIR aquifer-num
TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
add_test_compareECLFiles(CASENAME spe3
FILENAME SPE3CASE1
SIMULATOR flow
@ -998,4 +1014,20 @@ if(MPI_FOUND)
REL_TOL ${rel_tol_parallel}
DIR parallel_fieldprops
TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_2aqu
FILENAME 3D_2AQU_NUM
SIMULATOR flow
ABS_TOL 0.12
REL_TOL ${coarse_rel_tol_parallel}
DIR aquifer-num
TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_1aqu
FILENAME 3D_1AQU_3CELLS
SIMULATOR flow
ABS_TOL ${abs_tol_parallel}
REL_TOL ${coarse_rel_tol_parallel}
DIR aquifer-num
TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
endif()

16
ebos/eclbasevanguard.hh
View File

@ -210,6 +210,7 @@ public:
protected:
static const int dimension = Grid::dimension;
using Element = typename GridView::template Codim<0>::Entity;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
public:
@ -637,6 +638,9 @@ public:
const std::string& caseName() const
{ return caseName_; }
const CartesianIndexMapper& cartesianMapper() const
{ return asImp_().cartesianIndexMapper(); }
/*!
* \brief Returns the number of logically Cartesian cells in each direction
*/
@ -787,10 +791,22 @@ protected:
ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
auto elemIt = this->gridView().template begin</*codim=*/0>();
const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
const auto num_aqu_cells = this->eclState_->aquifer().numericalAquifers().allAquiferCells();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellCenterDepth_[elemIdx] = cellCenterDepth(element);
if (!num_aqu_cells.empty()) {
const unsigned int global_index = cartesianIndex(elemIdx);
const auto search = num_aqu_cells.find(global_index);
if (search != num_aqu_cells.end()) {
// updating the cell depth using aquifer cell depth
cellCenterDepth_[elemIdx] = search->second->depth;
}
}
}
}

1
ebos/eclequilinitializer.hh
View File

@ -101,6 +101,7 @@ public:
EQUIL::DeckDependent::InitialStateComputer<TypeTag> initialState(materialLawManager,
eclState,
vanguard.gridView(),
vanguard.cartesianMapper(),
simulator.problem().gravity()[dimWorld - 1]);
// copy the result into the array of initial fluid states

2
ebos/ecltransmissibility.hh
View File

@ -805,7 +805,7 @@ private:
// Discard the NNC if it is between active cell and inactive cell
std::ostringstream sstr;
sstr << "NNC between active and inactive cells ("
<< low << " -> " << high << ")";
<< low << " -> " << high << ") with globalcell is (" << c1 << "->" << c2 <<")";
Opm::OpmLog::warning(sstr.str());
continue;
}

95
ebos/equil/initstateequil.hh
View File

@ -48,6 +48,7 @@
#include <opm/parser/eclipse/EclipseState/Tables/PdvdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SaltvdTable.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <fmt/format.h>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
@ -62,6 +63,7 @@
#include <type_traits>
#include <utility>
#include <vector>
#include <string>
namespace Opm {
@ -1566,12 +1568,14 @@ class InitialStateComputer
using Element = typename GridView::template Codim<0>::Entity;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using Grid = GetPropType<TypeTag, Properties::Grid>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
public:
template<class MaterialLawManager>
InitialStateComputer(MaterialLawManager& materialLawManager,
const Opm::EclipseState& eclipseState,
const GridView& gridView,
const CartesianIndexMapper& cartMapper,
const double grav = Opm::unit::gravity,
const bool applySwatInit = true)
: temperature_(gridView.size(/*codim=*/0)),
@ -1581,7 +1585,8 @@ public:
sat_(FluidSystem::numPhases,
std::vector<double>(gridView.size(/*codim=*/0))),
rs_(gridView.size(/*codim=*/0)),
rv_(gridView.size(/*codim=*/0))
rv_(gridView.size(/*codim=*/0)),
cartesianIndexMapper_(cartMapper)
{
//Check for presence of kw SWATINIT
if (applySwatInit) {
@ -1590,8 +1595,10 @@ public:
}
}
//Querry cell depth, cell top-bottom.
updateCellProps_(gridView);
// Querry cell depth, cell top-bottom.
// numerical aquifer cells might be specified with different depths.
const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
updateCellProps_(gridView, num_aquifers);
// Get the equilibration records.
const std::vector<Opm::EquilRecord> rec = getEquil(eclipseState);
@ -1708,6 +1715,9 @@ public:
const auto& comm = gridView.comm();
calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
// modify the pressure and saturation for numerical aquifer cells
applyNumericalAquifers_(gridView, num_aquifers);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
}
@ -1769,12 +1779,14 @@ private:
PVec sat_;
Vec rs_;
Vec rv_;
const CartesianIndexMapper& cartesianIndexMapper_;
Vec swatInit_;
Vec cellCenterDepth_;
std::vector<std::pair<double,double>> cellZSpan_;
std::vector<std::pair<double,double>> cellZMinMax_;
void updateCellProps_(const GridView& gridView)
void updateCellProps_(const GridView& gridView,
const NumericalAquifers& aquifer)
{
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
int numElements = gridView.size(/*codim=*/0);
@ -1784,15 +1796,77 @@ private:
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
const auto num_aqu_cells = aquifer.allAquiferCells();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element);
const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
if (!num_aqu_cells.empty()) {
const auto search = num_aqu_cells.find(cartIx);
if (search != num_aqu_cells.end()) {
const auto* aqu_cell = num_aqu_cells.at(cartIx);
cellCenterDepth_[elemIdx] = aqu_cell->depth;
}
}
cellZSpan_[elemIdx] = Details::cellZSpan(element);
cellZMinMax_[elemIdx] = Details::cellZMinMax(element);
}
}
void applyNumericalAquifers_(const GridView& gridView,
const NumericalAquifers& aquifer)
{
const auto num_aqu_cells = aquifer.allAquiferCells();
if (num_aqu_cells.empty()) return;
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 Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
const auto search = num_aqu_cells.find(cartIx);
if (search != num_aqu_cells.end()) {
// numerical aquifer cells are filled with water initially
const auto watPos = FluidSystem::waterPhaseIdx;
if (FluidSystem::phaseIsActive(watPos)) {
this->sat_[watPos][elemIdx] = 1.;
} else {
throw std::logic_error { "Water phase has to be active for numerical aquifer case" };
}
const auto oilPos = FluidSystem::oilPhaseIdx;
if (FluidSystem::phaseIsActive(oilPos)) {
this->sat_[oilPos][elemIdx] = 0.;
}
const auto gasPos = FluidSystem::gasPhaseIdx;
if (FluidSystem::phaseIsActive(gasPos)) {
this->sat_[gasPos][elemIdx] = 0.;
}
const auto* aqu_cell = num_aqu_cells.at(cartIx);
const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
"AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id);
OpmLog::info(msg);
// if pressure is specified for numerical aquifers, we use these pressure values
// for numerical aquifer cells
if (aqu_cell->init_pressure) {
const double pres = *(aqu_cell->init_pressure);
this->pp_[watPos][elemIdx] = pres;
if (FluidSystem::phaseIsActive(gasPos)) {
this->pp_[gasPos][elemIdx] = pres;
}
if (FluidSystem::phaseIsActive(oilPos)) {
this->pp_[oilPos][elemIdx] = pres;
}
}
}
}
}
template<class RMap>
void setRegionPvtIdx(const Opm::EclipseState& eclState, const RMap& reg)
{
@ -1806,7 +1880,7 @@ private:
template <class RMap, class MaterialLawManager, class Comm>
void calcPressSatRsRv(const RMap& reg,
const std::vector< Opm::EquilRecord >& rec,
const std::vector<Opm::EquilRecord>& rec,
MaterialLawManager& materialLawManager,
const Comm& comm,
const double grav)
@ -1905,16 +1979,15 @@ private:
}
template <class CellRange, class PressTable, class PhaseSat>
void equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg,
const PressTable& ptable,
PhaseSat& psat)
void equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg,
const PressTable& ptable,
PhaseSat& psat)
{
using CellPos = typename PhaseSat::Position;
using CellID = std::remove_cv_t<std::remove_reference_t<
decltype(std::declval<CellPos>().cell)>>;
this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
(const CellID cell,
Details::PhaseQuantityValue& pressures,
Details::PhaseQuantityValue& saturations,

228
opm/simulators/aquifers/AquiferNumerical.hpp Normal file
View File

@ -0,0 +1,228 @@
/*
Copyright (C) 2020 Equinor ASA
Copyright (C) 2020 SINTEF Digital
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_AQUIFERNUMERICAL_HEADER_INCLUDED
#define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
#include <opm/output/data/Aquifer.hpp>
#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
namespace Opm
{
template <typename TypeTag>
class AquiferNumerical
{
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 GridView = GetPropType<TypeTag, Properties::GridView>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
enum { dimWorld = GridView::dimensionworld };
static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
static const int numEq = BlackoilIndices::numEq;
using Eval = DenseAd::Evaluation<double, numEq>;
using Toolbox = Opm::MathToolbox<Eval>;
// 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)
, flux_rate_(0.)
, cumulative_flux_(0.)
, global_cell_(global_cell)
{
this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < aquifer.numCells(); ++idx) {
const auto& cell = *(aquifer.getCellPrt(idx));
const int global_idx = cell.global_index;
const auto search = cartesian_to_compressed.find(global_idx);
// Due to parallelisation, the cell might not exist in the current process
if (search != cartesian_to_compressed.end()) {
const int cell_idx = cartesian_to_compressed.at(global_idx);
this->cell_to_aquifer_cell_idx_[cell_idx] = idx;
}
}
}
void initFromRestart([[maybe_unused]]const std::vector<data::AquiferData>& aquiferSoln)
{
// NOT handling Restart for now
}
void endTimeStep()
{
this->pressure_ = this->calculateAquiferPressure();
this->flux_rate_ = this->calculateAquiferFluxRate();
this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
}
Opm::data::AquiferData aquiferData() const
{
data::AquiferData data;
data.aquiferID = this->id_;
data.initPressure = this->init_pressure_;
data.pressure = this->pressure_;
data.fluxRate = this->flux_rate_;
data.volume = this->cumulative_flux_;
data.type = Opm::data::AquiferType::Numerical;
return data;
}
void initialSolutionApplied()
{
this->init_pressure_ = this->calculateAquiferPressure();
this->pressure_ = this->init_pressure_;
this->flux_rate_ = 0.;
this->cumulative_flux_ = 0.;
}
private:
const 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
double init_pressure_;
double pressure_; // aquifer pressure
// TODO: maybe unordered_map can also do the work to save memory?
std::vector<int> cell_to_aquifer_cell_idx_;
double calculateAquiferPressure() const
{
double sum_pressure_watervolume = 0.;
double sum_watervolume = 0.;
ElementContext elem_ctx(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;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elem_ctx.updatePrimaryStencil(elem);
const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
if (idx < 0) {
continue;
}
elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
// TODO: the porosity of the cells are still wrong for numerical aquifer cells
// Because the dofVolume still based on the grid information.
// The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
const double water_saturation = fs.saturation(waterPhaseIdx).value();
const double porosity = iq0.porosity().value();
const double volume = elem_ctx.dofTotalVolume(0, 0);
// TODO: not sure we should use water pressure here
const double water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
const double water_volume = volume * porosity * water_saturation;
sum_pressure_watervolume += water_volume * water_pressure_reservoir;
sum_watervolume += water_volume;
}
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&sum_pressure_watervolume, 1);
comm.sum(&sum_watervolume, 1);
return sum_pressure_watervolume / sum_watervolume;
}
double calculateAquiferFluxRate() const
{
double aquifer_flux = 0.;
ElementContext elem_ctx(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;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
// elem_ctx.updatePrimaryStencil(elem);
elem_ctx.updateStencil(elem);
const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
// we only need the first aquifer cell
if (idx != 0) {
continue;
}
elem_ctx.updateAllIntensiveQuantities();
elem_ctx.updateAllExtensiveQuantities();
const size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
// const auto &problem = elem_ctx.problem();
const auto &stencil = elem_ctx.stencil(0);
// const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
for (size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
const auto &face = stencil.interiorFace(face_idx);
// dof index
const size_t i = face.interiorIndex();
const size_t j = face.exteriorIndex();
// compressed index
// const size_t I = stencil.globalSpaceIndex(i);
const size_t J = stencil.globalSpaceIndex(j);
assert(stencil.globalSpaceIndex(i) == cell_index);
// we do not consider the flux within aquifer cells
// we only need the flux to the connections
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(waterPhaseIdx));
const size_t up_id = water_flux >= 0. ? i : j;
const auto &intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
const double invB = Toolbox::value(intQuantsIn.fluidState().invB(waterPhaseIdx));
const double face_area = face.area();
aquifer_flux += water_flux * invB * face_area;
}
// we only need to handle the first aquifer cell, we can exit loop here
break;
}
return aquifer_flux;
}
};
} // namespace Opm
#endif

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

@ -34,6 +34,7 @@
#include <opm/simulators/aquifers/AquiferCarterTracy.hpp>
#include <opm/simulators/aquifers/AquiferFetkovich.hpp>
#include <opm/simulators/aquifers/AquiferNumerical.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
@ -120,6 +121,7 @@ protected:
// 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;
// This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
void init();
@ -127,6 +129,7 @@ protected:
bool aquiferActive() const;
bool aquiferCarterTracyActive() const;
bool aquiferFetkovichActive() const;
bool aquiferNumericalActive() const;
};

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

@ -18,6 +18,9 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/grid/utility/cartesianToCompressed.hpp>
#include "BlackoilAquiferModel.hpp"
namespace Opm
{
@ -46,6 +49,12 @@ BlackoilAquiferModel<TypeTag>::initialSolutionApplied()
aquifer.initialSolutionApplied();
}
}
if (this->aquiferNumericalActive()) {
for (auto& aquifer : this->aquifers_numerical) {
aquifer.initialSolutionApplied();
}
}
}
template <typename TypeTag>
@ -62,6 +71,11 @@ BlackoilAquiferModel<TypeTag>::initFromRestart(const std::vector<data::AquiferDa
aquifer.initFromRestart(aquiferSoln);
}
}
if (aquiferNumericalActive()) {
for (auto& aquifer : this->aquifers_numerical) {
aquifer.initFromRestart(aquiferSoln);
}
}
}
template <typename TypeTag>
@ -132,6 +146,12 @@ BlackoilAquiferModel<TypeTag>::endTimeStep()
aquifer.endTimeStep();
}
}
if (aquiferNumericalActive()) {
for (auto& aquifer : this->aquifers_numerical) {
aquifer.endTimeStep();
this->simulator_.vanguard().grid().comm().barrier();
}
}
}
template <typename TypeTag>
void
@ -168,7 +188,7 @@ BlackoilAquiferModel<TypeTag>::init()
return;
}
// Get all the carter tracy aquifer properties data and put it in aquifers vector
const auto& connections = aquifer.connections();
for (const auto& aq : aquifer.ct()) {
aquifers_CarterTracy.emplace_back(connections[aq.aquiferID],
@ -179,12 +199,19 @@ BlackoilAquiferModel<TypeTag>::init()
aquifers_Fetkovich.emplace_back(connections[aq.aquiferID],
this->simulator_, aq);
}
}
template <typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferActive() const
{
return (aquiferCarterTracyActive() || aquiferFetkovichActive());
if (aquifer.hasNumericalAquifer()) {
const auto& num_aquifers = aquifer.numericalAquifers().aquifers();
const auto& ugrid = simulator_.vanguard().grid();
const int number_of_cells = simulator_.gridView().size(0);
const int* global_cell = Opm::UgGridHelpers::globalCell(ugrid);
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);
}
}
}
template <typename TypeTag>
bool
@ -199,6 +226,13 @@ BlackoilAquiferModel<TypeTag>::aquiferFetkovichActive() const
return !aquifers_Fetkovich.empty();
}
template<typename TypeTag>
bool
BlackoilAquiferModel<TypeTag>::aquiferNumericalActive() const
{
return !(this->aquifers_numerical.empty());
}
template<typename TypeTag>
Opm::data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const {
Opm::data::Aquifers data;
@ -215,6 +249,14 @@ Opm::data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const {
data[aqu_data.aquiferID] = aqu_data;
}
}
if (this->aquiferNumericalActive()) {
for (const auto& aqu : this->aquifers_numerical) {
Opm::data::AquiferData aqu_data = aqu.aquiferData();
data[aqu_data.aquiferID] = aqu_data;
}
}
return data;
}
} // namespace Opm

28
tests/test_equil.cc
View File

@ -475,7 +475,9 @@ BOOST_AUTO_TEST_CASE(DeckAllDead)
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 10.0);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@ -552,7 +554,9 @@ BOOST_AUTO_TEST_CASE(DeckWithCapillary)
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 10.0);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@ -590,7 +594,9 @@ BOOST_AUTO_TEST_CASE(DeckWithCapillaryOverlap)
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@ -649,7 +655,9 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveOil)
const UnstructuredGrid& grid = *(gm.c_grid());
// Initialize the fluid system
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@ -724,7 +732,9 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveGas)
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@ -802,7 +812,9 @@ BOOST_AUTO_TEST_CASE(DeckWithRSVDAndRVVD)
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@ -900,7 +912,9 @@ BOOST_AUTO_TEST_CASE(DeckWithPBVDAndPDVD)
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);

2
tests/update_reference_data.sh
View File

@ -63,6 +63,8 @@ tests[spe1_brine]="flow spe1_brine SPE1CASE1_BRINE"
tests[spe1_water]="flow_onephase spe1 SPE1CASE1_WATER"
tests[ctaquifer_2d_oilwater]="flow aquifer-oilwater 2D_OW_CTAQUIFER"
tests[fetkovich_2d]="flow aquifer-fetkovich 2D_FETKOVICHAQUIFER"
tests[numerical_aquifer_3d_1aqu]="flow aquifer-num 3D_1AQU_3CELLS"
tests[numerical_aquifer_3d_2aqu]="flow aquifer-num 3D_2AQU_NUM"
tests[msw_2d_h]="flow msw_2d_h 2D_H__"
tests[msw_3d_hfa]="flow msw_3d_hfa 3D_MSW"
tests[polymer_oilwater]="flow polymer_oilwater 2D_OILWATER_POLYMER"