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 #1758 from akva2/janitoring

Browse Source 
convert to unix eol
This commit is contained in:
Kai Bao 2019-03-14 14:41:18 +01:00 committed by GitHub
commit a7362abb04
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 239 additions and 239 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

478
opm/autodiff/AquiferFetkovich.hpp Executable file → Normal file
View File

@ -1,240 +1,240 @@
/* /*
Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface
Copyright 2017 Statoil ASA. Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or the Free Software Foundation, either version 3 of the License, or
(at your option) any later version. (at your option) any later version.
OPM is distributed in the hope that it will be useful, OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef OPM_AQUIFETP_HEADER_INCLUDED #ifndef OPM_AQUIFETP_HEADER_INCLUDED
#define OPM_AQUIFETP_HEADER_INCLUDED #define OPM_AQUIFETP_HEADER_INCLUDED
#include <opm/autodiff/AquiferInterface.hpp> #include <opm/autodiff/AquiferInterface.hpp>
namespace Opm namespace Opm
{ {
template<typename TypeTag> template<typename TypeTag>
class AquiferFetkovich: public AquiferInterface<TypeTag> class AquiferFetkovich: public AquiferInterface<TypeTag>
{ {
public: public:
typedef AquiferInterface<TypeTag> Base; typedef AquiferInterface<TypeTag> Base;
using typename Base::Simulator; using typename Base::Simulator;
using typename Base::ElementContext; using typename Base::ElementContext;
using typename Base::FluidSystem; using typename Base::FluidSystem;
using typename Base::BlackoilIndices; using typename Base::BlackoilIndices;
using typename Base::RateVector; using typename Base::RateVector;
using typename Base::IntensiveQuantities; using typename Base::IntensiveQuantities;
using typename Base::Eval; using typename Base::Eval;
using typename Base::Scalar; using typename Base::Scalar;
using typename Base::FluidState; using typename Base::FluidState;
using Base::waterCompIdx; using Base::waterCompIdx;
using Base::waterPhaseIdx; using Base::waterPhaseIdx;
AquiferFetkovich( const Aquancon::AquanconOutput& connection, AquiferFetkovich( const Aquancon::AquanconOutput& connection,
const std::unordered_map<int, int>& cartesian_to_compressed, const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator, const Simulator& ebosSimulator,
const Aquifetp::AQUFETP_data& aqufetp_data) const Aquifetp::AQUFETP_data& aqufetp_data)
: Base(connection, cartesian_to_compressed, ebosSimulator) : Base(connection, cartesian_to_compressed, ebosSimulator)
, aqufetp_data_(aqufetp_data) , aqufetp_data_(aqufetp_data)
{} {}
void endTimeStep() void endTimeStep()
{ {
for (const auto& Qai: Base::Qai_) { for (const auto& Qai: Base::Qai_) {
Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize(); Base::W_flux_ += Qai*Base::ebos_simulator_.timeStepSize();
aquifer_pressure_ = aquiferPressure(); aquifer_pressure_ = aquiferPressure();
} }
} }
protected: protected:
// Aquifer Fetkovich Specific Variables // Aquifer Fetkovich Specific Variables
const Aquifetp::AQUFETP_data aqufetp_data_; const Aquifetp::AQUFETP_data aqufetp_data_;
Scalar aquifer_pressure_; // aquifer Scalar aquifer_pressure_; // aquifer
inline void initializeConnections(const Aquancon::AquanconOutput& connection) inline void initializeConnections(const Aquancon::AquanconOutput& connection)
{ {
const auto& eclState = Base::ebos_simulator_.vanguard().eclState(); const auto& eclState = Base::ebos_simulator_.vanguard().eclState();
const auto& ugrid = Base::ebos_simulator_.vanguard().grid(); const auto& ugrid = Base::ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid(); const auto& grid = eclState.getInputGrid();
Base::cell_idx_ = connection.global_index; Base::cell_idx_ = connection.global_index;
auto globalCellIdx = ugrid.globalCell(); auto globalCellIdx = ugrid.globalCell();
assert( Base::cell_idx_ == connection.global_index); assert( Base::cell_idx_ == connection.global_index);
assert( (Base::cell_idx_.size() == connection.influx_coeff.size()) ); assert( (Base::cell_idx_.size() == connection.influx_coeff.size()) );
assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) ); assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) );
assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.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 // 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::cell_depth_.resize(Base::cell_idx_.size(), aqufetp_data_.d0);
Base::alphai_.resize(Base::cell_idx_.size(), 1.0); Base::alphai_.resize(Base::cell_idx_.size(), 1.0);
Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0); Base::faceArea_connected_.resize(Base::cell_idx_.size(),0.0);
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid); auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
// Translate the C face tag into the enum used by opm-parser's TransMult class // Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection; Opm::FaceDir::DirEnum faceDirection;
// denom_face_areas is the sum of the areas connected to an aquifer // denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas = 0.; Scalar denom_face_areas = 0.;
Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1); Base::cellToConnectionIdx_.resize(Base::ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx)
{ {
const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]); const int cell_index = Base::cartesian_to_compressed_.at(Base::cell_idx_[idx]);
Base::cellToConnectionIdx_[cell_index] = idx; Base::cellToConnectionIdx_[cell_index] = idx;
const auto cellFacesRange = cell2Faces[cell_index]; const auto cellFacesRange = cell2Faces[cell_index];
for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter)
{ {
// The index of the face in the compressed grid // The index of the face in the compressed grid
const int faceIdx = *cellFaceIter; const int faceIdx = *cellFaceIter;
// the logically-Cartesian direction of the face // the logically-Cartesian direction of the face
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
switch(faceTag) switch(faceTag)
{ {
case 0: faceDirection = Opm::FaceDir::XMinus; case 0: faceDirection = Opm::FaceDir::XMinus;
break; break;
case 1: faceDirection = Opm::FaceDir::XPlus; case 1: faceDirection = Opm::FaceDir::XPlus;
break; break;
case 2: faceDirection = Opm::FaceDir::YMinus; case 2: faceDirection = Opm::FaceDir::YMinus;
break; break;
case 3: faceDirection = Opm::FaceDir::YPlus; case 3: faceDirection = Opm::FaceDir::YPlus;
break; break;
case 4: faceDirection = Opm::FaceDir::ZMinus; case 4: faceDirection = Opm::FaceDir::ZMinus;
break; break;
case 5: faceDirection = Opm::FaceDir::ZPlus; case 5: faceDirection = Opm::FaceDir::ZPlus;
break; break;
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer problem. Make sure faceTag is correctly defined"); default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer problem. Make sure faceTag is correctly defined");
} }
if (faceDirection == connection.reservoir_face_dir.at(idx)) if (faceDirection == connection.reservoir_face_dir.at(idx))
{ {
Base::faceArea_connected_.at(idx) = Base::getFaceArea(faceCells, ugrid, faceIdx, idx, connection); 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) ); denom_face_areas += ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) );
} }
} }
auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx)); auto cellCenter = grid.getCellCenter(Base::cell_idx_.at(idx));
Base::cell_depth_.at(idx) = cellCenter[2]; Base::cell_depth_.at(idx) = cellCenter[2];
} }
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon()); const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (size_t idx = 0; idx < Base::cell_idx_.size(); ++idx) 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 Base::alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero
0. 0.
: ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas; : ( connection.influx_multiplier.at(idx) * Base::faceArea_connected_.at(idx) )/denom_face_areas;
} }
} }
inline Scalar dpai(int idx) inline Scalar dpai(int idx)
{ {
Scalar dp = aquifer_pressure_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aqufetp_data_.d0) - Base::pressure_current_.at(idx).value() ; Scalar dp = aquifer_pressure_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aqufetp_data_.d0) - Base::pressure_current_.at(idx).value() ;
return dp; return dp;
} }
// This function implements Eq 5.12 of the EclipseTechnicalDescription // This function implements Eq 5.12 of the EclipseTechnicalDescription
inline Scalar aquiferPressure() inline Scalar aquiferPressure()
{ {
Scalar Flux = Base::W_flux_.value(); Scalar Flux = Base::W_flux_.value();
Scalar pa_ = Base::pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 ); Scalar pa_ = Base::pa0_ - Flux / ( aqufetp_data_.C_t * aqufetp_data_.V0 );
return pa_; return pa_;
} }
inline void calculateAquiferConstants() inline void calculateAquiferConstants()
{ {
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 // This function implements Eq 5.14 of the EclipseTechnicalDescription
inline void calculateInflowRate(int idx, const Simulator& simulator) inline void calculateInflowRate(int idx, const Simulator& simulator)
{ {
Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ; Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ;
Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_; Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_;
Base::Qai_.at(idx) = Base::alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_; Base::Qai_.at(idx) = Base::alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_;
} }
inline void calculateAquiferCondition() inline void calculateAquiferCondition()
{ {
int pvttableIdx = aqufetp_data_.pvttableID - 1; int pvttableIdx = aqufetp_data_.pvttableID - 1;
Base::rhow_.resize(Base::cell_idx_.size(),0.); Base::rhow_.resize(Base::cell_idx_.size(),0.);
if (!aqufetp_data_.p0) if (!aqufetp_data_.p0)
{ {
Base::pa0_ = calculateReservoirEquilibrium(); Base::pa0_ = calculateReservoirEquilibrium();
} }
else else
{ {
Base::pa0_ = *(aqufetp_data_.p0); Base::pa0_ = *(aqufetp_data_.p0);
} }
aquifer_pressure_ = Base::pa0_ ; aquifer_pressure_ = Base::pa0_ ;
// use the thermodynamic state of the first active cell as a // use the thermodynamic state of the first active cell as a
// reference. there might be better ways to do this... // reference. there might be better ways to do this...
ElementContext elemCtx(Base::ebos_simulator_); ElementContext elemCtx(Base::ebos_simulator_);
auto elemIt = Base::ebos_simulator_.gridView().template begin</*codim=*/0>(); auto elemIt = Base::ebos_simulator_.gridView().template begin</*codim=*/0>();
elemCtx.updatePrimaryStencil(*elemIt); elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
// Initialize a FluidState object first // Initialize a FluidState object first
FluidState fs_aquifer; FluidState fs_aquifer;
// We use the temperature of the first cell connected to the 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 // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs
fs_aquifer.assign( iq0.fluidState() ); fs_aquifer.assign( iq0.fluidState() );
Eval temperature_aq, pa0_mean; Eval temperature_aq, pa0_mean;
temperature_aq = fs_aquifer.temperature(0); temperature_aq = fs_aquifer.temperature(0);
pa0_mean = Base::pa0_; pa0_mean = Base::pa0_;
Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
Base::mu_w_ = mu_w_aquifer.value(); Base::mu_w_ = mu_w_aquifer.value();
} }
inline Scalar calculateReservoirEquilibrium() inline Scalar calculateReservoirEquilibrium()
{ {
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer; std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir; Scalar water_pressure_reservoir;
ElementContext elemCtx(Base::ebos_simulator_); ElementContext elemCtx(Base::ebos_simulator_);
const auto& gridView = Base::ebos_simulator_.gridView(); const auto& gridView = Base::ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>(); auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>(); const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt; const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryStencil(elem);
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int idx = Base::cellToConnectionIdx_[cellIdx]; int idx = Base::cellToConnectionIdx_[cellIdx];
if (idx < 0) if (idx < 0)
continue; continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState(); const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
Base::rhow_[idx] = fs.density(waterPhaseIdx); 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] ); 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. // 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(); Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size();
return aquifer_pres_avg; return aquifer_pres_avg;
} }
}; //Class AquiferFetkovich }; //Class AquiferFetkovich
} // namespace Opm } // namespace Opm
#endif #endif

0
opm/autodiff/FlowMainEbos.hpp Executable file → Normal file
View File