diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 91c192834..5a8e8cf94 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -1,225 +1,225 @@ -/* - 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 . -*/ - -#ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED -#define OPM_AQUIFERINTERFACE_HEADER_INCLUDED - -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - -namespace Opm -{ - template - class AquiferInterface - { - public: - typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; - typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; - - enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; - enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; - - static const int numEq = BlackoilIndices::numEq; - typedef double Scalar; - - typedef DenseAd::Evaluation Eval; - - typedef Opm::BlackOilFluidState FluidState; - - static const auto waterCompIdx = FluidSystem::waterCompIdx; - static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; - - // Constructor - AquiferInterface( const Aquancon::AquanconOutput& connection, - const std::unordered_map& cartesian_to_compressed, - const Simulator& ebosSimulator) - : connection_(connection) - , ebos_simulator_(ebosSimulator) - , cartesian_to_compressed_(cartesian_to_compressed) - {} - - // Deconstructor - virtual ~AquiferInterface() {} - - void initialSolutionApplied() - { - initQuantities(connection_); - } - - void beginTimeStep() - { - ElementContext elemCtx(ebos_simulator_); - auto elemIt = ebos_simulator_.gridView().template begin<0>(); - const auto& elemEndIt = ebos_simulator_.gridView().template end<0>(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - - elemCtx.updatePrimaryStencil(elem); - - int cellIdx = elemCtx.globalSpaceIndex(0, 0); - int idx = cellToConnectionIdx_[cellIdx]; - if (idx < 0) - continue; - - elemCtx.updateIntensiveQuantities(0); - const auto& iq = elemCtx.intensiveQuantities(0, 0); - pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); - } - } - - template - void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) - { - unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - - int idx = 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 IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx); - // This is the pressure at td + dt - updateCellPressure(pressure_current_,idx,intQuants); - updateCellDensity(idx,intQuants); - calculateInflowRate(idx, context.simulator()); - rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] += - Qai_[idx]/context.dofVolume(spaceIdx, timeIdx); - } - - protected: - inline Scalar gravity_() const - { - return ebos_simulator_.problem().gravity()[2]; - } - - inline void initQuantities(const Aquancon::AquanconOutput& connection) - { - // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 - W_flux_ = 0.; - - // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function - initializeConnections(connection); - calculateAquiferCondition(); - calculateAquiferConstants(); - - pressure_previous_.resize(cell_idx_.size(), 0.); - pressure_current_.resize(cell_idx_.size(), 0.); - Qai_.resize(cell_idx_.size(), 0.0); - } - - inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) - { - const auto& fs = intQuants.fluidState(); - pressure_water.at(idx) = fs.pressure(waterPhaseIdx); - } - - inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) - { - const auto& fs = intQuants.fluidState(); - pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); - } - - inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants) - { - const auto& fs = intQuants.fluidState(); - rhow_.at(idx) = fs.density(waterPhaseIdx); - } - - template - inline double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid, - const int faceIdx, const int idx, - const Aquancon::AquanconOutput& connection) const - { - // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell - // Do not make the connection if the product of the two cellIdx > 0. This is because the - // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) - double faceArea = 0.; - const auto cellNeighbour0 = faceCells(faceIdx,0); - const auto cellNeighbour1 = faceCells(faceIdx,1); - const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); - const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))? - defaultFaceArea : - *(connection.influx_coeff.at(idx)); - faceArea = (cellNeighbour0 * cellNeighbour1 > 0)? 0. : calculatedFaceArea; - if (cellNeighbour1 == 0){ - faceArea = (cellNeighbour0 < 0)? faceArea : 0.; - } - else if (cellNeighbour0 == 0){ - faceArea = (cellNeighbour1 < 0)? faceArea : 0.; - } - return faceArea; - } - - virtual void endTimeStep() = 0; - - const Aquancon::AquanconOutput connection_; - const Simulator& ebos_simulator_; - const std::unordered_map cartesian_to_compressed_; - - // Grid variables - std::vector cell_idx_; - std::vector faceArea_connected_; - std::vector cellToConnectionIdx_; - // Quantities at each grid id - std::vector cell_depth_; - std::vector pressure_previous_; - std::vector pressure_current_; - std::vector Qai_; - std::vector rhow_; - std::vector alphai_; - - Scalar mu_w_; //water viscosity - Scalar Tc_; // Time constant - Scalar pa0_; // initial aquifer pressure - - Eval W_flux_; - - virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0; - - virtual Scalar dpai(int idx) = 0; - - virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0; - - virtual void calculateAquiferCondition() = 0; - - virtual void calculateAquiferConstants() = 0; - - virtual Scalar calculateReservoirEquilibrium() =0; - // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer - }; -} // namespace Opm -#endif +/* + 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 . +*/ + +#ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED +#define OPM_AQUIFERINTERFACE_HEADER_INCLUDED + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +namespace Opm +{ + template + class AquiferInterface + { + public: + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + + static const int numEq = BlackoilIndices::numEq; + typedef double Scalar; + + typedef DenseAd::Evaluation Eval; + + typedef Opm::BlackOilFluidState FluidState; + + static const auto waterCompIdx = FluidSystem::waterCompIdx; + static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; + + // Constructor + AquiferInterface( const Aquancon::AquanconOutput& connection, + const std::unordered_map& cartesian_to_compressed, + const Simulator& ebosSimulator) + : connection_(connection) + , ebos_simulator_(ebosSimulator) + , cartesian_to_compressed_(cartesian_to_compressed) + {} + + // Deconstructor + virtual ~AquiferInterface() {} + + void initialSolutionApplied() + { + initQuantities(connection_); + } + + void beginTimeStep() + { + ElementContext elemCtx(ebos_simulator_); + auto elemIt = ebos_simulator_.gridView().template begin<0>(); + const auto& elemEndIt = ebos_simulator_.gridView().template end<0>(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + + elemCtx.updatePrimaryStencil(elem); + + int cellIdx = elemCtx.globalSpaceIndex(0, 0); + int idx = cellToConnectionIdx_[cellIdx]; + if (idx < 0) + continue; + + elemCtx.updateIntensiveQuantities(0); + const auto& iq = elemCtx.intensiveQuantities(0, 0); + pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); + } + } + + template + void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) + { + unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + + int idx = 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 IntensiveQuantities intQuants = context.intensiveQuantities(spaceIdx, timeIdx); + // This is the pressure at td + dt + updateCellPressure(pressure_current_,idx,intQuants); + updateCellDensity(idx,intQuants); + calculateInflowRate(idx, context.simulator()); + rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] += + Qai_[idx]/context.dofVolume(spaceIdx, timeIdx); + } + + protected: + inline Scalar gravity_() const + { + return ebos_simulator_.problem().gravity()[2]; + } + + inline void initQuantities(const Aquancon::AquanconOutput& connection) + { + // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 + W_flux_ = 0.; + + // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function + initializeConnections(connection); + calculateAquiferCondition(); + calculateAquiferConstants(); + + pressure_previous_.resize(cell_idx_.size(), 0.); + pressure_current_.resize(cell_idx_.size(), 0.); + Qai_.resize(cell_idx_.size(), 0.0); + } + + inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + pressure_water.at(idx) = fs.pressure(waterPhaseIdx); + } + + inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); + } + + inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + rhow_.at(idx) = fs.density(waterPhaseIdx); + } + + template + inline double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid, + const int faceIdx, const int idx, + const Aquancon::AquanconOutput& connection) const + { + // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell + // Do not make the connection if the product of the two cellIdx > 0. This is because the + // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) + double faceArea = 0.; + const auto cellNeighbour0 = faceCells(faceIdx,0); + const auto cellNeighbour1 = faceCells(faceIdx,1); + const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))? + defaultFaceArea : + *(connection.influx_coeff.at(idx)); + faceArea = (cellNeighbour0 * cellNeighbour1 > 0)? 0. : calculatedFaceArea; + if (cellNeighbour1 == 0){ + faceArea = (cellNeighbour0 < 0)? faceArea : 0.; + } + else if (cellNeighbour0 == 0){ + faceArea = (cellNeighbour1 < 0)? faceArea : 0.; + } + return faceArea; + } + + virtual void endTimeStep() = 0; + + const Aquancon::AquanconOutput connection_; + const Simulator& ebos_simulator_; + const std::unordered_map cartesian_to_compressed_; + + // Grid variables + std::vector cell_idx_; + std::vector faceArea_connected_; + std::vector cellToConnectionIdx_; + // Quantities at each grid id + std::vector cell_depth_; + std::vector pressure_previous_; + std::vector pressure_current_; + std::vector Qai_; + std::vector rhow_; + std::vector alphai_; + + Scalar mu_w_; //water viscosity + Scalar Tc_; // Time constant + Scalar pa0_; // initial aquifer pressure + + Eval W_flux_; + + virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0; + + virtual Scalar dpai(int idx) = 0; + + virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0; + + virtual void calculateAquiferCondition() = 0; + + virtual void calculateAquiferConstants() = 0; + + virtual Scalar calculateReservoirEquilibrium() =0; + // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer + }; +} // namespace Opm +#endif diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index 90f5dee89..b81e7ca64 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -1,173 +1,146 @@ -#include -namespace Opm { +#include +namespace Opm { - template - BlackoilAquiferModel:: - BlackoilAquiferModel(Simulator& simulator) - : simulator_(simulator) - { - init(); - } - - template - void - BlackoilAquiferModel::initialSolutionApplied() - { - if(aquiferCarterTracyActive()) - { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) - { - aquifer->initialSolutionApplied(); - } - } - if(aquiferFetkovichActive()) - { - for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) - { - aquifer->initialSolutionApplied(); - } - } - } - - template - void - BlackoilAquiferModel::beginEpisode() - { } - - template - void - BlackoilAquiferModel::beginIteration() - { } - template - void BlackoilAquiferModel:: beginTimeStep() + BlackoilAquiferModel:: + BlackoilAquiferModel(Simulator& simulator) + : simulator_(simulator) + { + init(); + } + + template + void + BlackoilAquiferModel::initialSolutionApplied() { if(aquiferCarterTracyActive()) - { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) - { - aquifer->beginTimeStep(); - } - } - if(aquiferFetkovichActive()) - { - for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) - { - aquifer->beginTimeStep(); - } - } - } - - template - template - void BlackoilAquiferModel:: addToSource(RateVector& rates, const Context& context, 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); - } - } - } - - template - void - BlackoilAquiferModel::endIteration() - { } - - template - void BlackoilAquiferModel:: endTimeStep() - { - if(aquiferCarterTracyActive()) - { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) - { - aquifer->endTimeStep(); - } - } - if(aquiferFetkovichActive()) - { - for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) - { - aquifer->endTimeStep(); - } - } - } - template - void - BlackoilAquiferModel::endEpisode() - { } - - template - template - void - BlackoilAquiferModel::serialize(Restarter& /* res */) - { - // TODO (?) - throw std::logic_error("BlackoilAquiferModel::serialize() is not yet implemented"); - } - - template - template - void - BlackoilAquiferModel::deserialize(Restarter& /* res */) - { - // TODO (?) - throw std::logic_error("BlackoilAquiferModel::deserialize() is not yet implemented"); - } - - // Initialize the aquifers in the deck - template - void - BlackoilAquiferModel:: init() + { + for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) + { + aquifer->initialSolutionApplied(); + } + } + if(aquiferFetkovichActive()) + { + for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) + { + aquifer->initialSolutionApplied(); + } + } + } + + template + void + BlackoilAquiferModel::beginEpisode() + { } + + template + void + BlackoilAquiferModel::beginIteration() + { } + + template + void BlackoilAquiferModel:: beginTimeStep() { - const auto& deck = this->simulator_.vanguard().deck(); + if(aquiferCarterTracyActive()) + { + for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) + { + aquifer->beginTimeStep(); + } + } + if(aquiferFetkovichActive()) + { + for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) + { + aquifer->beginTimeStep(); + } + } + } + + template + template + void BlackoilAquiferModel:: addToSource(RateVector& rates, const Context& context, 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); + } + } + } + + template + void + BlackoilAquiferModel::endIteration() + { } + + template + void BlackoilAquiferModel:: endTimeStep() + { + if(aquiferCarterTracyActive()) + { + for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) + { + aquifer->endTimeStep(); + } + } + if(aquiferFetkovichActive()) + { + for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) + { + aquifer->endTimeStep(); + } + } + } + template + void + BlackoilAquiferModel::endEpisode() + { } + + template + template + void + BlackoilAquiferModel::serialize(Restarter& /* res */) + { + // TODO (?) + throw std::logic_error("BlackoilAquiferModel::serialize() is not yet implemented"); + } + + template + template + void + BlackoilAquiferModel::deserialize(Restarter& /* res */) + { + // TODO (?) + throw std::logic_error("BlackoilAquiferModel::deserialize() is not yet implemented"); + } + + // Initialize the aquifers in the deck + template + void + BlackoilAquiferModel:: init() + { + const auto& deck = this->simulator_.vanguard().deck(); if (deck.hasKeyword("AQUCT")) { - //updateConnectionIntensiveQuantities(); - const auto& eclState = this->simulator_.vanguard().eclState(); - - // Get all the carter tracy aquifer properties data and put it in aquifers vector - const AquiferCT aquiferct = AquiferCT(eclState,deck); - const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); - - std::vector aquifersData = aquiferct.getAquifers(); - std::vector aquifer_connection = aquifer_connect.getAquOutput(); - - assert( aquifersData.size() == aquifer_connection.size() ); - const auto& ugrid = simulator_.vanguard().grid(); - const auto& gridView = simulator_.gridView(); - const int number_of_cells = gridView.size(0); - - cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, - Opm::UgGridHelpers::globalCell(ugrid)); - - for (size_t i = 0; i < aquifersData.size(); ++i) - { - aquifers_CarterTracy.push_back( - AquiferCarterTracy (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i)) - ); - } - } - if(deck.hasKeyword("AQUFETP")) - { - //updateConnectionIntensiveQuantities(); - const auto& eclState = this->simulator_.vanguard().eclState(); - - // Get all the carter tracy aquifer properties data and put it in aquifers vector - const Aquifetp aquifetp = Aquifetp(deck); - const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); - - std::vector aquifersData = aquifetp.getAquifers(); - std::vector aquifer_connection = aquifer_connect.getAquOutput(); - + //updateConnectionIntensiveQuantities(); + const auto& eclState = this->simulator_.vanguard().eclState(); + + // Get all the carter tracy aquifer properties data and put it in aquifers vector + const AquiferCT aquiferct = AquiferCT(eclState,deck); + const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); + + std::vector aquifersData = aquiferct.getAquifers(); + std::vector aquifer_connection = aquifer_connect.getAquOutput(); + assert( aquifersData.size() == aquifer_connection.size() ); const auto& ugrid = simulator_.vanguard().grid(); const auto& gridView = simulator_.gridView(); @@ -176,30 +149,57 @@ namespace Opm { cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, Opm::UgGridHelpers::globalCell(ugrid)); - for (size_t i = 0; i < aquifersData.size(); ++i) - { - aquifers_Fetkovich.push_back( - AquiferFetkovich (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i)) - ); - } - } - } - template - bool - BlackoilAquiferModel:: aquiferActive() const + for (size_t i = 0; i < aquifersData.size(); ++i) + { + aquifers_CarterTracy.push_back( + AquiferCarterTracy (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i)) + ); + } + } + if(deck.hasKeyword("AQUFETP")) + { + //updateConnectionIntensiveQuantities(); + const auto& eclState = this->simulator_.vanguard().eclState(); + + // Get all the carter tracy aquifer properties data and put it in aquifers vector + const Aquifetp aquifetp = Aquifetp(deck); + const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); + + std::vector aquifersData = aquifetp.getAquifers(); + std::vector aquifer_connection = aquifer_connect.getAquOutput(); + + assert( aquifersData.size() == aquifer_connection.size() ); + const auto& ugrid = simulator_.vanguard().grid(); + const auto& gridView = simulator_.gridView(); + const int number_of_cells = gridView.size(0); + + cartesian_to_compressed_ = cartesianToCompressed(number_of_cells, + Opm::UgGridHelpers::globalCell(ugrid)); + + for (size_t i = 0; i < aquifersData.size(); ++i) + { + aquifers_Fetkovich.push_back( + AquiferFetkovich (aquifer_connection.at(i), cartesian_to_compressed_, this->simulator_ , aquifersData.at(i)) + ); + } + } + } + template + bool + BlackoilAquiferModel:: aquiferActive() const { return (aquiferCarterTracyActive() || aquiferFetkovichActive()); - } - template - bool - BlackoilAquiferModel:: aquiferCarterTracyActive() const - { - return !aquifers_CarterTracy.empty(); - } - template - bool - BlackoilAquiferModel:: aquiferFetkovichActive() const - { - return !aquifers_Fetkovich.empty(); - } -} // namespace Opm + } + template + bool + BlackoilAquiferModel:: aquiferCarterTracyActive() const + { + return !aquifers_CarterTracy.empty(); + } + template + bool + BlackoilAquiferModel:: aquiferFetkovichActive() const + { + return !aquifers_Fetkovich.empty(); + } +} // namespace Opm