From 3fb96deb368b5af3ce264324d88e1b14779d3be9 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 20 Dec 2019 14:30:13 +0100 Subject: [PATCH 1/2] re-formatting aquifer files with clang-format --- .../aquifers/AquiferCarterTracy.hpp | 421 +++++++++--------- opm/simulators/aquifers/AquiferFetkovich.hpp | 284 ++++++------ opm/simulators/aquifers/AquiferInterface.hpp | 223 +++++----- .../aquifers/BlackoilAquiferModel.hpp | 87 ++-- .../aquifers/BlackoilAquiferModel_impl.hpp | 349 +++++++-------- 5 files changed, 679 insertions(+), 685 deletions(-) diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index 471de6071..5e2a56e49 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -31,249 +31,250 @@ namespace Opm { - template - class AquiferCarterTracy: public AquiferInterface +template +class AquiferCarterTracy : public AquiferInterface +{ +public: + typedef AquiferInterface Base; + + using typename Base::BlackoilIndices; + using typename Base::ElementContext; + using typename Base::Eval; + using typename Base::FluidState; + using typename Base::FluidSystem; + using typename Base::IntensiveQuantities; + using typename Base::RateVector; + using typename Base::Scalar; + using typename Base::Simulator; + + using Base::waterCompIdx; + using Base::waterPhaseIdx; + AquiferCarterTracy(const Aquancon::AquanconOutput& connection, + const std::unordered_map& cartesian_to_compressed, + const Simulator& ebosSimulator, + const AquiferCT::AQUCT_data& aquct_data) + : Base(connection, cartesian_to_compressed, ebosSimulator) + , aquct_data_(aquct_data) { - public: - typedef AquiferInterface Base; + } - using typename Base::Simulator; - using typename Base::ElementContext; - using typename Base::FluidSystem; - using typename Base::BlackoilIndices; - using typename Base::RateVector; - using typename Base::IntensiveQuantities; - using typename Base::Eval; - using typename Base::Scalar; - using typename Base::FluidState; + void endTimeStep() override + { + for (const auto& Qai : Base::Qai_) { + Base::W_flux_ += Qai * Base::ebos_simulator_.timeStepSize(); + } + } - using Base::waterCompIdx; - using Base::waterPhaseIdx; - AquiferCarterTracy( const Aquancon::AquanconOutput& connection, - const std::unordered_map& cartesian_to_compressed, - const Simulator& ebosSimulator, - const AquiferCT::AQUCT_data& aquct_data) - : Base(connection, cartesian_to_compressed, ebosSimulator) - , aquct_data_(aquct_data) - {} +protected: + // Variables constants + const AquiferCT::AQUCT_data aquct_data_; + Scalar beta_; // Influx constant + // TODO: it is possible it should be a AD variable + Scalar mu_w_; // water viscosity - void endTimeStep() override - { - for (const auto& Qai: Base::Qai_) { - Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize(); - } - } + // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer + inline void initializeConnections(const Aquancon::AquanconOutput& connection) override + { + const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); + const auto& ugrid = Base::ebos_simulator_.vanguard().grid(); + const auto& grid = eclState.getInputGrid(); - protected: - // Variables constants - const AquiferCT::AQUCT_data aquct_data_; - Scalar beta_; // Influx constant - // TODO: it is possible it should be a AD variable - Scalar mu_w_; // water viscosity + Base::cell_idx_ = connection.global_index; + auto globalCellIdx = ugrid.globalCell(); - // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer - inline void initializeConnections(const Aquancon::AquanconOutput& connection) override - { - const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); - const auto& ugrid = Base::ebos_simulator_.vanguard().grid(); - const auto& grid = eclState.getInputGrid(); + assert(Base::cell_idx_ == connection.global_index); + assert((Base::cell_idx_.size() <= connection.influx_coeff.size())); + assert((connection.influx_coeff.size() == connection.influx_multiplier.size())); + assert((connection.influx_multiplier.size() == connection.reservoir_face_dir.size())); - Base::cell_idx_ = connection.global_index; - auto globalCellIdx = ugrid.globalCell(); + // We hack the cell depth values for now. We can actually get it from elementcontext pos + Base::cell_depth_.resize(Base::cell_idx_.size(), aquct_data_.d0); + Base::alphai_.resize(Base::cell_idx_.size(), 1.0); + Base::faceArea_connected_.resize(Base::cell_idx_.size(), 0.0); - assert( Base::cell_idx_ == connection.global_index); - assert( (Base::cell_idx_.size() <= connection.influx_coeff.size()) ); - assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) ); - assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); + auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); + auto faceCells = Opm::UgGridHelpers::faceCells(ugrid); - // We hack the cell depth values for now. We can actually get it from elementcontext pos - Base::cell_depth_.resize(Base::cell_idx_.size(), aquct_data_.d0); - Base::alphai_.resize(Base::cell_idx_.size(), 1.0); - Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0); + // Translate the C face tag into the enum used by opm-parser's TransMult class + Opm::FaceDir::DirEnum faceDirection; - auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); - auto faceCells = Opm::UgGridHelpers::faceCells(ugrid); + // denom_face_areas is the sum of the areas connected to an aquifer + Scalar denom_face_areas = 0.; + Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1); + for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) { + const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]); + Base::cellToConnectionIdx_[cell_index] = idx; - // Translate the C face tag into the enum used by opm-parser's TransMult class - Opm::FaceDir::DirEnum faceDirection; + const auto cellFacesRange = cell2Faces[cell_index]; + for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) { + // The index of the face in the compressed grid + const int faceIdx = *cellFaceIter; - // denom_face_areas is the sum of the areas connected to an aquifer - Scalar denom_face_areas = 0.; - Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1); - for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) - { - const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]); - Base::cellToConnectionIdx_[cell_index] = idx; + // the logically-Cartesian direction of the face + const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); - const auto cellFacesRange = cell2Faces[cell_index]; - for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) - { - // The index of the face in the compressed grid - const int faceIdx = *cellFaceIter; - - // the logically-Cartesian direction of the face - const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); - - switch(faceTag) - { - case 0: faceDirection = Opm::FaceDir::XMinus; - break; - case 1: faceDirection = Opm::FaceDir::XPlus; - break; - case 2: faceDirection = Opm::FaceDir::YMinus; - break; - case 3: faceDirection = Opm::FaceDir::YPlus; - break; - case 4: faceDirection = Opm::FaceDir::ZMinus; - break; - case 5: faceDirection = Opm::FaceDir::ZPlus; - break; - default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined"); - } - - if (faceDirection == connection.reservoir_face_dir.at(idx)) - { - Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection); - denom_face_areas += ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) ); - } - } - auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx)); - Base::cell_depth_.at(idx) = cellCenter[2]; + switch (faceTag) { + case 0: + faceDirection = Opm::FaceDir::XMinus; + break; + case 1: + faceDirection = Opm::FaceDir::XPlus; + break; + case 2: + faceDirection = Opm::FaceDir::YMinus; + break; + case 3: + faceDirection = Opm::FaceDir::YPlus; + break; + case 4: + faceDirection = Opm::FaceDir::ZMinus; + break; + case 5: + faceDirection = Opm::FaceDir::ZPlus; + break; + default: + OPM_THROW(Opm::NumericalIssue, + "Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined"); } - const double eps_sqrt = std::sqrt(std::numeric_limits::epsilon()); - for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) - { - Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero - 0. - : ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas; + if (faceDirection == connection.reservoir_face_dir.at(idx)) { + Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection); + denom_face_areas += (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)); } } + auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx)); + Base::cell_depth_.at(idx) = cellCenter[2]; + } - void assignRestartData(const data::AquiferData& /* xaq */) override - { - throw std::runtime_error { - "Restart-based initialization not currently supported " - "for Carter-Tracey analytic aquifers" - }; - } + const double eps_sqrt = std::sqrt(std::numeric_limits::epsilon()); + for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) { + Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt) + ? // Prevent no connection NaNs due to division by zero + 0. + : (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)) / denom_face_areas; + } + } - inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) - { - // We use the opm-common numeric linear interpolator - pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td); - pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td); - } + void assignRestartData(const data::AquiferData& /* xaq */) override + { + throw std::runtime_error {"Restart-based initialization not currently supported " + "for Carter-Tracey analytic aquifers"}; + } - inline Scalar dpai(int idx) - { - Scalar dp = Base::pa0_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aquct_data_.d0) - Base::pressure_previous_.at(idx); - return dp; - } + inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) + { + // We use the opm-common numeric linear interpolator + pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td); + pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td); + } - // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription - inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator) - { - const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Base::Tc_; - const Scalar td = simulator.time() / Base::Tc_; - Scalar PItdprime = 0.; - Scalar PItd = 0.; - getInfluenceTableValues(PItd, PItdprime, td_plus_dt); - a = 1.0/Base::Tc_ * ( (beta_ * dpai(idx)) - (Base::W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); - b = beta_ / (Base::Tc_ * ( PItd - td*PItdprime)); - } + inline Scalar dpai(int idx) + { + Scalar dp = Base::pa0_ + + Base::rhow_.at(idx).value() * Base::gravity_() * (Base::cell_depth_.at(idx) - aquct_data_.d0) + - Base::pressure_previous_.at(idx); + return dp; + } - // This function implements Eq 5.7 of the EclipseTechnicalDescription - inline void calculateInflowRate(int idx, const Simulator& simulator) override - { - Scalar a, b; - calculateEqnConstants(a,b,idx,simulator); - Base::Qai_.at(idx) = Base::alphai_.at(idx)*( a - b * ( Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx) ) ); - } + // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription + inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator) + { + const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Base::Tc_; + const Scalar td = simulator.time() / Base::Tc_; + Scalar PItdprime = 0.; + Scalar PItd = 0.; + getInfluenceTableValues(PItd, PItdprime, td_plus_dt); + a = 1.0 / Base::Tc_ * ((beta_ * dpai(idx)) - (Base::W_flux_.value() * PItdprime)) / (PItd - td * PItdprime); + b = beta_ / (Base::Tc_ * (PItd - td * PItdprime)); + } - inline void calculateAquiferConstants() override - { - // We calculate the influx constant - beta_ = aquct_data_.c2 * aquct_data_.h - * aquct_data_.theta * aquct_data_.phi_aq - * aquct_data_.C_t - * aquct_data_.r_o * aquct_data_.r_o; - // We calculate the time constant - Base::Tc_ = mu_w_ * aquct_data_.phi_aq - * aquct_data_.C_t - * aquct_data_.r_o * aquct_data_.r_o - / ( aquct_data_.k_a * aquct_data_.c1 ); - } + // This function implements Eq 5.7 of the EclipseTechnicalDescription + inline void calculateInflowRate(int idx, const Simulator& simulator) override + { + Scalar a, b; + calculateEqnConstants(a, b, idx, simulator); + Base::Qai_.at(idx) + = Base::alphai_.at(idx) * (a - b * (Base::pressure_current_.at(idx) - Base::pressure_previous_.at(idx))); + } - inline void calculateAquiferCondition() override - { + inline void calculateAquiferConstants() override + { + // We calculate the influx constant + beta_ = aquct_data_.c2 * aquct_data_.h * aquct_data_.theta * aquct_data_.phi_aq * aquct_data_.C_t + * aquct_data_.r_o * aquct_data_.r_o; + // We calculate the time constant + Base::Tc_ = mu_w_ * aquct_data_.phi_aq * aquct_data_.C_t * aquct_data_.r_o * aquct_data_.r_o + / (aquct_data_.k_a * aquct_data_.c1); + } - int pvttableIdx = aquct_data_.pvttableID - 1; - Base::rhow_.resize(Base::cell_idx_.size(),0.); - if (!aquct_data_.p0) - { - Base::pa0_ = calculateReservoirEquilibrium(); - } - else - { - Base::pa0_ = *(aquct_data_.p0); - } + inline void calculateAquiferCondition() override + { - // use the thermodynamic state of the first active cell as a - // reference. there might be better ways to do this... - ElementContext elemCtx(Base::ebos_simulator_); - auto elemIt = Base::ebos_simulator_.gridView().template begin(); - elemCtx.updatePrimaryStencil(*elemIt); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - // Initialize a FluidState object first - FluidState fs_aquifer; - // We use the temperature of the first cell connected to the aquifer - // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs - fs_aquifer.assign( iq0.fluidState() ); - Eval temperature_aq, pa0_mean; - temperature_aq = fs_aquifer.temperature(0); - pa0_mean = Base::pa0_; - Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); - mu_w_ = mu_w_aquifer.value(); + int pvttableIdx = aquct_data_.pvttableID - 1; + Base::rhow_.resize(Base::cell_idx_.size(), 0.); + if (!aquct_data_.p0) { + Base::pa0_ = calculateReservoirEquilibrium(); + } else { + Base::pa0_ = *(aquct_data_.p0); + } - } + // use the thermodynamic state of the first active cell as a + // reference. there might be better ways to do this... + ElementContext elemCtx(Base::ebos_simulator_); + auto elemIt = Base::ebos_simulator_.gridView().template begin(); + elemCtx.updatePrimaryStencil(*elemIt); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + // Initialize a FluidState object first + FluidState fs_aquifer; + // We use the temperature of the first cell connected to the aquifer + // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs + fs_aquifer.assign(iq0.fluidState()); + Eval temperature_aq, pa0_mean; + temperature_aq = fs_aquifer.temperature(0); + pa0_mean = Base::pa0_; + Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); + mu_w_ = mu_w_aquifer.value(); + } - // This function is for calculating the aquifer properties from equilibrium state with the reservoir - // TODO: this function can be moved to the Inteface class, since it is the same for both Aquifer models - inline Scalar calculateReservoirEquilibrium() override - { - // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices - std::vector pw_aquifer; - Scalar water_pressure_reservoir; + // This function is for calculating the aquifer properties from equilibrium state with the reservoir + // TODO: this function can be moved to the Inteface class, since it is the same for both Aquifer models + inline Scalar calculateReservoirEquilibrium() override + { + // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices + std::vector pw_aquifer; + Scalar water_pressure_reservoir; - ElementContext elemCtx(Base::ebos_simulator_); - const auto& gridView = Base::ebos_simulator_.gridView(); - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - elemCtx.updatePrimaryStencil(elem); + ElementContext elemCtx(Base::ebos_simulator_); + const auto& gridView = Base::ebos_simulator_.gridView(); + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); - size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - int idx = Base::cellToConnectionIdx_[cellIdx]; - if (idx < 0) - continue; + size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + int idx = Base::cellToConnectionIdx_[cellIdx]; + if (idx < 0) + continue; - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = iq0.fluidState(); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = iq0.fluidState(); - water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); - Base::rhow_[idx] = fs.density(waterPhaseIdx); - pw_aquifer.push_back( (water_pressure_reservoir - Base::rhow_[idx].value()*Base::gravity_()*(Base::cell_depth_[idx] - aquct_data_.d0))*Base::alphai_[idx] ); - } + water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); + Base::rhow_[idx] = fs.density(waterPhaseIdx); + pw_aquifer.push_back( + (water_pressure_reservoir + - Base::rhow_[idx].value() * Base::gravity_() * (Base::cell_depth_[idx] - aquct_data_.d0)) + * Base::alphai_[idx]); + } - // We take the average of the calculated equilibrium pressures. - Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size(); - return aquifer_pres_avg; - } - }; // class AquiferCarterTracy + // We take the average of the calculated equilibrium pressures. + Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / pw_aquifer.size(); + return aquifer_pres_avg; + } +}; // class AquiferCarterTracy } // namespace Opm #endif diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index 90843fec7..558a9d626 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -31,43 +31,44 @@ along with OPM. If not, see . namespace Opm { - template - class AquiferFetkovich: public AquiferInterface - { +template +class AquiferFetkovich : public AquiferInterface +{ - public: +public: typedef AquiferInterface Base; - using typename Base::Simulator; - using typename Base::ElementContext; - using typename Base::FluidSystem; using typename Base::BlackoilIndices; - using typename Base::RateVector; - using typename Base::IntensiveQuantities; + using typename Base::ElementContext; using typename Base::Eval; - using typename Base::Scalar; using typename Base::FluidState; + using typename Base::FluidSystem; + using typename Base::IntensiveQuantities; + using typename Base::RateVector; + using typename Base::Scalar; + using typename Base::Simulator; using Base::waterCompIdx; using Base::waterPhaseIdx; - AquiferFetkovich( const Aquancon::AquanconOutput& connection, - const std::unordered_map& cartesian_to_compressed, - const Simulator& ebosSimulator, - const Aquifetp::AQUFETP_data& aqufetp_data) - : Base(connection, cartesian_to_compressed, ebosSimulator) - , aqufetp_data_(aqufetp_data) - {} + AquiferFetkovich(const Aquancon::AquanconOutput& connection, + const std::unordered_map& cartesian_to_compressed, + const Simulator& ebosSimulator, + const Aquifetp::AQUFETP_data& aqufetp_data) + : Base(connection, cartesian_to_compressed, ebosSimulator) + , aqufetp_data_(aqufetp_data) + { + } void endTimeStep() override { - for (const auto& Qai: Base::Qai_) { - Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize(); - aquifer_pressure_ = aquiferPressure(); - } + for (const auto& Qai : Base::Qai_) { + Base::W_flux_ += Qai * Base::ebos_simulator_.timeStepSize(); + aquifer_pressure_ = aquiferPressure(); + } } - protected: +protected: // Aquifer Fetkovich Specific Variables // TODO: using const reference here will cause segmentation fault, which is very strange const Aquifetp::AQUFETP_data aqufetp_data_; @@ -75,173 +76,174 @@ namespace Opm inline void initializeConnections(const Aquancon::AquanconOutput& connection) override { - const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); - const auto& ugrid = Base::ebos_simulator_.vanguard().grid(); - const auto& grid = eclState.getInputGrid(); + const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); + const auto& ugrid = Base::ebos_simulator_.vanguard().grid(); + const auto& grid = eclState.getInputGrid(); - Base::cell_idx_ = connection.global_index; - auto globalCellIdx = ugrid.globalCell(); + Base::cell_idx_ = connection.global_index; + auto globalCellIdx = ugrid.globalCell(); - assert( Base::cell_idx_ == connection.global_index); - assert( (Base::cell_idx_.size() == connection.influx_coeff.size()) ); - assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) ); - assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); + assert(Base::cell_idx_ == connection.global_index); + assert((Base::cell_idx_.size() == connection.influx_coeff.size())); + assert((connection.influx_coeff.size() == connection.influx_multiplier.size())); + assert((connection.influx_multiplier.size() == connection.reservoir_face_dir.size())); - // We hack the cell depth values for now. We can actually get it from elementcontext pos - Base::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0); - Base::alphai_.resize(Base::cell_idx_.size(), 1.0); - Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0); + // We hack the cell depth values for now. We can actually get it from elementcontext pos + Base::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0); + Base::alphai_.resize(Base::cell_idx_.size(), 1.0); + Base::faceArea_connected_.resize(Base::cell_idx_.size(), 0.0); - auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); - auto faceCells = Opm::UgGridHelpers::faceCells(ugrid); + auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); + auto faceCells = Opm::UgGridHelpers::faceCells(ugrid); - // Translate the C face tag into the enum used by opm-parser's TransMult class - Opm::FaceDir::DirEnum faceDirection; + // Translate the C face tag into the enum used by opm-parser's TransMult class + Opm::FaceDir::DirEnum faceDirection; - // denom_face_areas is the sum of the areas connected to an aquifer - Scalar denom_face_areas = 0.; - Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1); - for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) - { - const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]); - Base::cellToConnectionIdx_[cell_index] = idx; + // denom_face_areas is the sum of the areas connected to an aquifer + Scalar denom_face_areas = 0.; + Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1); + for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) { + const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]); + Base::cellToConnectionIdx_[cell_index] = idx; - const auto cellFacesRange = cell2Faces[cell_index]; - for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) - { - // The index of the face in the compressed grid - const int faceIdx = *cellFaceIter; + const auto cellFacesRange = cell2Faces[cell_index]; + for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) { + // The index of the face in the compressed grid + const int faceIdx = *cellFaceIter; - // the logically-Cartesian direction of the face - const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); + // the logically-Cartesian direction of the face + const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); - switch(faceTag) - { - case 0: faceDirection = Opm::FaceDir::XMinus; - break; - case 1: faceDirection = Opm::FaceDir::XPlus; - break; - case 2: faceDirection = Opm::FaceDir::YMinus; - break; - case 3: faceDirection = Opm::FaceDir::YPlus; - break; - case 4: faceDirection = Opm::FaceDir::ZMinus; - break; - case 5: faceDirection = Opm::FaceDir::ZPlus; - break; - default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer problem. Make sure faceTag is correctly defined"); - } + switch (faceTag) { + case 0: + faceDirection = Opm::FaceDir::XMinus; + break; + case 1: + faceDirection = Opm::FaceDir::XPlus; + break; + case 2: + faceDirection = Opm::FaceDir::YMinus; + break; + case 3: + faceDirection = Opm::FaceDir::YPlus; + break; + case 4: + faceDirection = Opm::FaceDir::ZMinus; + break; + case 5: + faceDirection = Opm::FaceDir::ZPlus; + break; + default: + OPM_THROW(Opm::NumericalIssue, + "Initialization of Aquifer problem. Make sure faceTag is correctly defined"); + } - if (faceDirection == connection.reservoir_face_dir.at(idx)) - { - Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection); - denom_face_areas += ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) ); - } + if (faceDirection == connection.reservoir_face_dir.at(idx)) { + Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection); + denom_face_areas += (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)); + } + } + auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx)); + Base::cell_depth_.at(idx) = cellCenter[2]; } - auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx)); - Base::cell_depth_.at(idx) = cellCenter[2]; - } - const double eps_sqrt = std::sqrt(std::numeric_limits::epsilon()); - for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) - { - Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero - 0. - : ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas; - } + const double eps_sqrt = std::sqrt(std::numeric_limits::epsilon()); + for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) { + Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt) + ? // Prevent no connection NaNs due to division by zero + 0. + : (connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx)) / denom_face_areas; + } } void assignRestartData(const data::AquiferData& xaq) override { - if (xaq.type != data::AquiferType::Fetkovich) - { - throw std::invalid_argument { - "Analytic aquifer data for unexpected aquifer type " - "passed to Fetkovich aquifer" - }; - } + if (xaq.type != data::AquiferType::Fetkovich) { + throw std::invalid_argument {"Analytic aquifer data for unexpected aquifer type " + "passed to Fetkovich aquifer"}; + } - this->aquifer_pressure_ = xaq.pressure; + this->aquifer_pressure_ = xaq.pressure; } inline Eval dpai(int idx) { - const Eval dp = aquifer_pressure_ - Base::pressure_current_.at(idx) - + Base::rhow_[idx] * Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0); - return dp; + const Eval dp = aquifer_pressure_ - Base::pressure_current_.at(idx) + + Base::rhow_[idx] * Base::gravity_() * (Base::cell_depth_[idx] - aqufetp_data_.d0); + return dp; } // This function implements Eq 5.12 of the EclipseTechnicalDescription inline Scalar aquiferPressure() { - Scalar Flux = Base::W_flux_.value(); - Scalar pa_ = Base::pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 ); - return pa_; + Scalar Flux = Base::W_flux_.value(); + Scalar pa_ = Base::pa0_ - Flux / (aqufetp_data_.C_t * aqufetp_data_.V0); + return pa_; } inline void calculateAquiferConstants() override { - Base::Tc_ = ( aqufetp_data_.C_t * aqufetp_data_.V0 ) / aqufetp_data_.J ; + Base::Tc_ = (aqufetp_data_.C_t * aqufetp_data_.V0) / aqufetp_data_.J; } // This function implements Eq 5.14 of the EclipseTechnicalDescription inline void calculateInflowRate(int idx, const Simulator& simulator) override { - const Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ; - const Scalar coef = (1 - exp(-td_Tc_)) / td_Tc_; - Base::Qai_.at(idx) = Base::alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef; + const Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_; + const Scalar coef = (1 - exp(-td_Tc_)) / td_Tc_; + Base::Qai_.at(idx) = Base::alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef; } inline void calculateAquiferCondition() override { - Base::rhow_.resize(Base::cell_idx_.size(),0.); + Base::rhow_.resize(Base::cell_idx_.size(), 0.); - if (this->solution_set_from_restart_) { - return; - } + if (this->solution_set_from_restart_) { + return; + } - if (!aqufetp_data_.p0) - { - Base::pa0_ = calculateReservoirEquilibrium(); - } - else - { - Base::pa0_ = *(aqufetp_data_.p0); - } - aquifer_pressure_ = Base::pa0_ ; + if (!aqufetp_data_.p0) { + Base::pa0_ = calculateReservoirEquilibrium(); + } else { + Base::pa0_ = *(aqufetp_data_.p0); + } + aquifer_pressure_ = Base::pa0_; } inline Scalar calculateReservoirEquilibrium() override { - // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices - std::vector pw_aquifer; - Scalar water_pressure_reservoir; + // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices + std::vector pw_aquifer; + Scalar water_pressure_reservoir; - ElementContext elemCtx(Base::ebos_simulator_); - const auto& gridView = Base::ebos_simulator_.gridView(); - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - elemCtx.updatePrimaryStencil(elem); - size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - int idx = Base::cellToConnectionIdx_[cellIdx]; - if (idx < 0) - continue; + ElementContext elemCtx(Base::ebos_simulator_); + const auto& gridView = Base::ebos_simulator_.gridView(); + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + int idx = Base::cellToConnectionIdx_[cellIdx]; + if (idx < 0) + continue; - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = iq0.fluidState(); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = iq0.fluidState(); - water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); - Base::rhow_[idx] = fs.density(waterPhaseIdx); - pw_aquifer.push_back( (water_pressure_reservoir - Base::rhow_[idx].value()*Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0))*Base::alphai_[idx] ); - } + water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); + Base::rhow_[idx] = fs.density(waterPhaseIdx); + pw_aquifer.push_back( + (water_pressure_reservoir + - Base::rhow_[idx].value() * Base::gravity_() * (Base::cell_depth_[idx] - aqufetp_data_.d0)) + * Base::alphai_[idx]); + } - // We take the average of the calculated equilibrium pressures. - const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.); - const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha; - return aquifer_pres_avg; + // We take the average of the calculated equilibrium pressures. + const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.); + const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha; + return aquifer_pres_avg; } - }; //Class AquiferFetkovich +}; // Class AquiferFetkovich } // namespace Opm #endif diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 6ea145308..ee18165a9 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -22,28 +22,28 @@ #ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED #define OPM_AQUIFERINTERFACE_HEADER_INCLUDED +#include +#include #include #include -#include -#include #include #include -#include #include +#include #include -#include #include #include +#include namespace Opm { - template - class AquiferInterface - { - public: +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; @@ -59,154 +59,165 @@ namespace Opm typedef DenseAd::Evaluation Eval; - typedef Opm::BlackOilFluidState FluidState; + 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) - {} + 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() {} + virtual ~AquiferInterface() + { + } void initFromRestart(const std::vector& aquiferSoln) { - auto xaqPos = std::find_if(aquiferSoln.begin(), aquiferSoln.end(), - [this](const data::AquiferData& xaq) -> bool - { - return xaq.aquiferID == this->connection_.aquiferID; - }); + auto xaqPos + = std::find_if(aquiferSoln.begin(), aquiferSoln.end(), [this](const data::AquiferData& xaq) -> bool { + return xaq.aquiferID == this->connection_.aquiferID; + }); - if (xaqPos == aquiferSoln.end()) { - // No restart value applies to this aquifer. Nothing to do. - return; - } + if (xaqPos == aquiferSoln.end()) { + // No restart value applies to this aquifer. Nothing to do. + return; + } - this->assignRestartData(*xaqPos); + this->assignRestartData(*xaqPos); - this->W_flux_ = xaqPos->volume; - this->pa0_ = xaqPos->initPressure; - this->solution_set_from_restart_ = true; + this->W_flux_ = xaqPos->volume; + this->pa0_ = xaqPos->initPressure; + this->solution_set_from_restart_ = true; } void initialSolutionApplied() { - initQuantities(connection_); + 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; + 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); + elemCtx.updatePrimaryStencil(elem); - int cellIdx = elemCtx.globalSpaceIndex(0, 0); - int idx = cellToConnectionIdx_[cellIdx]; - if (idx < 0) - continue; + 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)); - } + 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); + unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - int idx = cellToConnectionIdx_[cellIdx]; - if (idx < 0) - return; + 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); + // 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: +protected: inline Scalar gravity_() const { - return ebos_simulator_.problem().gravity()[2]; + 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 - if (!this->solution_set_from_restart_) - { - W_flux_ = 0.; - } + // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 + if (!this->solution_set_from_restart_) { + W_flux_ = 0.; + } - // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function - initializeConnections(connection); - calculateAquiferCondition(); - calculateAquiferConstants(); + // 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); + 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) + 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); + 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) + 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(); + 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); + 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, + 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; + // 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; @@ -232,9 +243,9 @@ namespace Opm Eval W_flux_; - bool solution_set_from_restart_{false}; + bool solution_set_from_restart_ {false}; - virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0; + virtual void initializeConnections(const Aquancon::AquanconOutput& connection) = 0; virtual void assignRestartData(const data::AquiferData& xaq) = 0; @@ -244,8 +255,8 @@ namespace Opm 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 - }; + 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.hpp b/opm/simulators/aquifers/BlackoilAquiferModel.hpp index 19158f9c8..13bd0d9ca 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel.hpp @@ -26,9 +26,9 @@ #include +#include #include #include -#include #include @@ -39,62 +39,59 @@ #include -namespace Opm { +namespace Opm +{ - /// Class for handling the blackoil well model. - template - class BlackoilAquiferModel - { - typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; - typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; +/// Class for handling the blackoil well model. +template +class BlackoilAquiferModel +{ + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; - public: - explicit BlackoilAquiferModel(Simulator& simulator); +public: + explicit BlackoilAquiferModel(Simulator& simulator); - void initialSolutionApplied(); - void initFromRestart(const std::vector& aquiferSoln); + void initialSolutionApplied(); + void initFromRestart(const std::vector& aquiferSoln); - void beginEpisode(); - void beginTimeStep(); - void beginIteration(); - // add the water rate due to aquifers to the source term. - template - void addToSource(RateVector& rates, - const Context& context, - unsigned spaceIdx, - unsigned timeIdx) const; - void endIteration(); - void endTimeStep(); - void endEpisode(); + void beginEpisode(); + void beginTimeStep(); + void beginIteration(); + // add the water rate due to aquifers to the source term. + template + void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) const; + void endIteration(); + void endTimeStep(); + void endEpisode(); - template - void serialize(Restarter& res); + template + void serialize(Restarter& res); - template - void deserialize(Restarter& res); + template + void deserialize(Restarter& res); - protected: - // --------- Types --------- - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +protected: + // --------- Types --------- + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef AquiferCarterTracy AquiferCarterTracy_object; - typedef AquiferFetkovich AquiferFetkovich_object; + typedef AquiferCarterTracy AquiferCarterTracy_object; + typedef AquiferFetkovich AquiferFetkovich_object; - Simulator& simulator_; + Simulator& simulator_; - std::unordered_map cartesian_to_compressed_; - mutable std::vector aquifers_CarterTracy; - mutable std::vector aquifers_Fetkovich; + std::unordered_map cartesian_to_compressed_; + mutable std::vector aquifers_CarterTracy; + mutable std::vector aquifers_Fetkovich; - // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy - void init(); + // 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 aquiferActive() const; + bool aquiferCarterTracyActive() const; + bool aquiferFetkovichActive() const; +}; } // namespace Opm diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index 594bf42d9..0f1513b05 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -1,225 +1,208 @@ #include -namespace Opm { +namespace Opm +{ - template - BlackoilAquiferModel:: - BlackoilAquiferModel(Simulator& simulator) - : simulator_(simulator) - { +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(); - } +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(); - } + if (aquiferFetkovichActive()) { + for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) { + aquifer->initialSolutionApplied(); + } } - } +} - template - void - BlackoilAquiferModel::initFromRestart(const std::vector& aquiferSoln) - { - if(aquiferCarterTracyActive()) - { - for (auto& aquifer : aquifers_CarterTracy) - { - aquifer.initFromRestart(aquiferSoln); - } +template +void +BlackoilAquiferModel::initFromRestart(const std::vector& aquiferSoln) +{ + if (aquiferCarterTracyActive()) { + for (auto& aquifer : aquifers_CarterTracy) { + aquifer.initFromRestart(aquiferSoln); + } } - if(aquiferFetkovichActive()) - { - for (auto& aquifer : aquifers_Fetkovich) - { - aquifer.initFromRestart(aquiferSoln); - } + if (aquiferFetkovichActive()) { + for (auto& aquifer : aquifers_Fetkovich) { + aquifer.initFromRestart(aquiferSoln); + } } - } +} - template - void - BlackoilAquiferModel::beginEpisode() - { } +template +void +BlackoilAquiferModel::beginEpisode() +{ +} - template - void - BlackoilAquiferModel::beginIteration() - { } +template +void +BlackoilAquiferModel::beginIteration() +{ +} - template - void BlackoilAquiferModel:: beginTimeStep() - { - if(aquiferCarterTracyActive()) - { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) - { - aquifer->beginTimeStep(); - } +template +void +BlackoilAquiferModel::beginTimeStep() +{ + 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(); - } + 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); - } +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); - } + if (aquiferFetkovichActive()) { + for (auto& aquifer : aquifers_Fetkovich) { + aquifer.addToSource(rates, context, spaceIdx, timeIdx); + } } - } +} - template - void - BlackoilAquiferModel::endIteration() - { } +template +void +BlackoilAquiferModel::endIteration() +{ +} - template - void BlackoilAquiferModel:: endTimeStep() - { - if(aquiferCarterTracyActive()) - { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) - { - aquifer->endTimeStep(); - } +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(); - } + if (aquiferFetkovichActive()) { + for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) { + aquifer->endTimeStep(); + } } - } - template - void - BlackoilAquiferModel::endEpisode() - { } +} +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::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"); - } +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() - { +// 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(); + // 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); + // 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(); + 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); + 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)); + 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)) - ); - } + 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(); + 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); + // 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(); + 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); + 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)); + 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)) - ); - } + 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 - { +} +template +bool +BlackoilAquiferModel::aquiferActive() const +{ return (aquiferCarterTracyActive() || aquiferFetkovichActive()); - } - template - bool - BlackoilAquiferModel:: aquiferCarterTracyActive() const - { +} +template +bool +BlackoilAquiferModel::aquiferCarterTracyActive() const +{ return !aquifers_CarterTracy.empty(); - } - template - bool - BlackoilAquiferModel:: aquiferFetkovichActive() const - { +} +template +bool +BlackoilAquiferModel::aquiferFetkovichActive() const +{ return !aquifers_Fetkovich.empty(); - } +} } // namespace Opm From 7fc81ae22844394b36f54996ebe2a67282c3d121 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 20 Dec 2019 15:10:36 +0100 Subject: [PATCH 2/2] using range for a few occasions in BlackoilAquiferModel_impl.hpp --- .../aquifers/BlackoilAquiferModel_impl.hpp | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index 0f1513b05..0c4b3ac4f 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -14,13 +14,13 @@ void BlackoilAquiferModel::initialSolutionApplied() { if (aquiferCarterTracyActive()) { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) { - aquifer->initialSolutionApplied(); + for (auto& aquifer : aquifers_CarterTracy) { + aquifer.initialSolutionApplied(); } } if (aquiferFetkovichActive()) { - for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) { - aquifer->initialSolutionApplied(); + for (auto& aquifer : aquifers_Fetkovich) { + aquifer.initialSolutionApplied(); } } } @@ -58,13 +58,13 @@ void BlackoilAquiferModel::beginTimeStep() { if (aquiferCarterTracyActive()) { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) { - aquifer->beginTimeStep(); + for (auto& aquifer : aquifers_CarterTracy) { + aquifer.beginTimeStep(); } } if (aquiferFetkovichActive()) { - for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) { - aquifer->beginTimeStep(); + for (auto& aquifer : aquifers_Fetkovich) { + aquifer.beginTimeStep(); } } } @@ -100,13 +100,13 @@ void BlackoilAquiferModel::endTimeStep() { if (aquiferCarterTracyActive()) { - for (auto aquifer = aquifers_CarterTracy.begin(); aquifer != aquifers_CarterTracy.end(); ++aquifer) { - aquifer->endTimeStep(); + for (auto& aquifer : aquifers_CarterTracy) { + aquifer.endTimeStep(); } } if (aquiferFetkovichActive()) { - for (auto aquifer = aquifers_Fetkovich.begin(); aquifer != aquifers_Fetkovich.end(); ++aquifer) { - aquifer->endTimeStep(); + for (auto& aquifer : aquifers_Fetkovich) { + aquifer.endTimeStep(); } } }