diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index c9414faba..4d6cad725 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -43,8 +43,10 @@ namespace Opm 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) }; @@ -61,65 +63,69 @@ namespace Opm AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, - const Aquancon::AquanconOutput& connection, - Simulator& ebosSimulator ) - : ebos_simulator_ (ebosSimulator), - aquct_data_ (aquct_data), - gravity_ (ebos_simulator_.problem().gravity()[2]) + const Aquancon::AquanconOutput& connection, + const Simulator& ebosSimulator) + : ebos_simulator_ (ebosSimulator) + , aquct_data_ (aquct_data) + , connection_(connection) + {} + + void initialSolutionApplied() { - initQuantities(connection); + initQuantities(connection_); } - inline void assembleAquiferEq(const SimulatorTimerInterface& timer) + void beginTimeStep() { - auto& ebosJac = ebos_simulator_.model().linearizer().matrix(); - auto& ebosResid = ebos_simulator_.model().linearizer().residual(); + 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; - size_t cellID; - for ( size_t idx = 0; idx < cell_idx_.size(); ++idx ) - { - Eval qinflow = 0.0; - cellID = cell_idx_.at(idx); - // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to - // IntensiveQuantities of that particular cell_id - const IntensiveQuantities intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellID, /*timeIdx=*/ 0)); - // This is the pressure at td + dt - updateCellPressure(pressure_current_,idx,intQuants); - updateCellDensity(idx,intQuants); - calculateInflowRate(idx, timer); - qinflow = Qai_.at(idx); - ebosResid[cellID][waterCompIdx] -= qinflow.value(); + elemCtx.updatePrimaryStencil(elem); - for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) - { - // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cellID][cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); - } + 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)); } } - inline void beforeTimeStep(const SimulatorTimerInterface& timer) + template + void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) { - auto cellID = cell_idx_.begin(); - size_t idx; - for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) - { - const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); - updateCellPressure(pressure_previous_ ,idx,intQuants); - } + 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); } - inline void afterTimeStep(const SimulatorTimerInterface& timer) + void endTimeStep() { - for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) - { - W_flux_ += (*Qai)*timer.currentStepLength(); + for (const auto& Qai: Qai_) { + W_flux_ += Qai*ebos_simulator_.timeStepSize(); } } - private: - Simulator& ebos_simulator_; + const Simulator& ebos_simulator_; // Grid variables std::vector cell_idx_; @@ -136,15 +142,17 @@ namespace Opm // Variables constants const AquiferCT::AQUCT_data aquct_data_; - Scalar mu_w_ , //water viscosity - beta_ , // Influx constant - Tc_ , // Time constant - pa0_ , // initial aquifer pressure - gravity_ ; // gravitational acceleration + Scalar mu_w_; //water viscosity + Scalar beta_; // Influx constant + Scalar Tc_; // Time constant + Scalar pa0_; // initial aquifer pressure Eval W_flux_; + Scalar gravity_() const + { return ebos_simulator_.problem().gravity()[2]; } + inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { // We use the opm-common numeric linear interpolator @@ -189,15 +197,15 @@ namespace Opm inline Scalar dpai(int idx) { - Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0) - pressure_previous_.at(idx); + Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_()*(cell_depth_.at(idx) - aquct_data_.d0) - pressure_previous_.at(idx); return dp; } // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription - inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) + inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator) { - const Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc_; - const Scalar td = timer.simulationTimeElapsed() / Tc_; + const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / Tc_; + const Scalar td = simulator.time() / Tc_; Scalar PItdprime = 0.; Scalar PItd = 0.; getInfluenceTableValues(PItd, PItdprime, td_plus_dt); @@ -206,10 +214,10 @@ namespace Opm } // This function implements Eq 5.7 of the EclipseTechnicalDescription - inline void calculateInflowRate(int idx, const SimulatorTimerInterface& timer) + inline void calculateInflowRate(int idx, const Simulator& simulator) { Scalar a, b; - calculateEqnConstants(a,b,idx,timer); + calculateEqnConstants(a,b,idx,simulator); Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) ); } @@ -256,10 +264,12 @@ namespace Opm // denom_face_areas is the sum of the areas connected to an aquifer Scalar denom_face_areas = 0.; + cellToConnectionIdx_.resize(ebos_simulator_.gridView().size(/*codim=*/0), -1); for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { + cellToConnectionIdx_[cell_idx_[idx]] = idx; + auto cellFacesRange = cell2Faces[cell_idx_.at(idx)]; - for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) { // The index of the face in the compressed grid @@ -321,11 +331,19 @@ namespace Opm 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(ebos_simulator_); + auto elemIt = 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( ebos_simulator_.model().cachedIntensiveQuantities(cell_idx_.at(0), /*timeIdx=*/ 0)->fluidState() ); + fs_aquifer.assign( iq0.fluidState() ); Eval temperature_aq, pa0_mean; temperature_aq = fs_aquifer.temperature(0); pa0_mean = pa0_; @@ -343,15 +361,26 @@ namespace Opm std::vector pw_aquifer; Scalar water_pressure_reservoir; - for (size_t idx = 0; idx < cell_idx_.size(); ++idx) - { - size_t cellIDx = cell_idx_.at(idx); - const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); - const auto& fs = intQuants.fluidState(); - + ElementContext elemCtx(ebos_simulator_); + const auto& gridView = 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 = cellToConnectionIdx_[cellIdx]; + if (idx < 0) + continue; + + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = iq0.fluidState(); + water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); - rhow_.at(idx) = fs.density(waterPhaseIdx); - pw_aquifer.push_back( (water_pressure_reservoir - rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0))*alphai_.at(idx) ); + rhow_[idx] = fs.density(waterPhaseIdx); + pw_aquifer.push_back( (water_pressure_reservoir - rhow_[idx].value()*gravity_()*(cell_depth_[idx] - aquct_data_.d0))*alphai_[idx] ); } // We take the average of the calculated equilibrium pressures. @@ -359,10 +388,11 @@ namespace Opm return aquifer_pres_avg; } - + const Aquancon::AquanconOutput connection_; + std::vector cellToConnectionIdx_; }; // class AquiferCarterTracy } // namespace Opm -#endif \ No newline at end of file +#endif diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index bc831333b..ac671ec67 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -24,6 +24,8 @@ #ifndef OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED #define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED +#include + #include #include #include @@ -34,45 +36,45 @@ namespace Opm { /// Class for handling the blackoil well model. template - class BlackoilAquiferModel { + class BlackoilAquiferModel + { + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; public: + explicit BlackoilAquiferModel(Simulator& simulator); - // --------- Types --------- - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - - typedef AquiferCarterTracy Aquifer_object; - - explicit BlackoilAquiferModel(Simulator& ebosSimulator); - - // compute the well fluxes and assemble them in to the reservoir equations as source terms - // and in the well equations. - void assemble( const SimulatorTimerInterface& timer, - const int iterationIdx ); - - // called at the end of a time step - void timeStepSucceeded(const SimulatorTimerInterface& timer); + void initialSolutionApplied(); + 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(); protected: + // --------- Types --------- + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - Simulator& ebosSimulator_; + typedef AquiferCarterTracy AquiferType; - std::vector aquifers_; + // TODO: declaring this as mutable is a hack which should be fixed in the + // long term + mutable std::vector aquifers_; + + Simulator& simulator_; // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy void init(); - void updateConnectionIntensiveQuantities() const; - - void assembleAquiferEq(const SimulatorTimerInterface& timer); - - // at the beginning of each time step (Not report step) - void prepareTimeStep(const SimulatorTimerInterface& timer); - bool aquiferActive() const; - }; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 66de30574..2c5646746 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -3,104 +3,85 @@ namespace Opm { template BlackoilAquiferModel:: - BlackoilAquiferModel(Simulator& ebosSimulator) - : ebosSimulator_(ebosSimulator) + BlackoilAquiferModel(Simulator& simulator) + : simulator_(simulator) { init(); } - - // called at the end of a time step template void - BlackoilAquiferModel:: timeStepSucceeded(const SimulatorTimerInterface& timer) + BlackoilAquiferModel::initialSolutionApplied() { - if ( !aquiferActive() ) { - return; - } - - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) - { - aquifer->afterTimeStep(timer); + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { + aquifer->initialSolutionApplied(); } } template void - BlackoilAquiferModel:: - assemble( const SimulatorTimerInterface& timer, - const int iterationIdx ) - { - if ( !aquiferActive() ) { - return; - } - - // We need to update the reservoir pressures connected to the aquifer - updateConnectionIntensiveQuantities(); - - if (iterationIdx == 0) { - // We can do the Table check and coefficients update in this function - // For now, it does nothing! - prepareTimeStep(timer); - } - - assembleAquiferEq(timer); - } - + BlackoilAquiferModel::beginEpisode() + { } template void - BlackoilAquiferModel:: updateConnectionIntensiveQuantities() const + BlackoilAquiferModel::beginTimeStep() { - ElementContext elemCtx(ebosSimulator_); - const auto& gridView = ebosSimulator_.gridView(); - const auto& elemEndIt = gridView.template end(); - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; - ++elemIt) - { - elemCtx.updatePrimaryStencil(*elemIt); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { + aquifer->beginTimeStep(); } } - // Protected function which calls the individual aquifer models template void - BlackoilAquiferModel:: assembleAquiferEq(const SimulatorTimerInterface& timer) + BlackoilAquiferModel::beginIteration() + { } + + template + template + void + BlackoilAquiferModel::addToSource(RateVector& rates, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) const { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) - { - aquifer->assembleAquiferEq(timer); + for (auto& aquifer: aquifers_) { + aquifer.addToSource(rates, context, spaceIdx, timeIdx); } } - // Protected function - // some preparation work, mostly related to group control and RESV, - // at the beginning of each time step (Not report step) template - void BlackoilAquiferModel:: prepareTimeStep(const SimulatorTimerInterface& timer) + void + BlackoilAquiferModel::endIteration() + { } + + template + void + BlackoilAquiferModel::endTimeStep() { - // Here we can ask each carter tracy aquifers to get the current previous time step's pressure - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) - { - aquifer->beforeTimeStep(timer); + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { + aquifer->endTimeStep(); } } + template + void + BlackoilAquiferModel::endEpisode() + { } + // Initialize the aquifers in the deck template void BlackoilAquiferModel:: init() { - const auto& deck = ebosSimulator_.vanguard().deck(); + const auto& deck = this->simulator_.vanguard().deck(); if ( !deck.hasKeyword("AQUCT") ) { return ; } - updateConnectionIntensiveQuantities(); - const auto& eclState = ebosSimulator_.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); @@ -115,7 +96,7 @@ namespace Opm { for (size_t i = 0; i < aquifersData.size(); ++i) { aquifers_.push_back( - AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), ebosSimulator_) + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), this->simulator_) ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index dd7af8680..5e95cb74a 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -78,6 +78,8 @@ SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume, true); SET_BOOL_PROP(EclFlowProblem, UseVolumetricResidual, false); +SET_TYPE_PROP(EclFlowProblem, EclAquiferModel, Opm::BlackoilAquiferModel); + // disable all extensions supported by black oil model. this should not really be // necessary but it makes things a bit more explicit SET_BOOL_PROP(EclFlowProblem, EnablePolymer, false); @@ -145,7 +147,6 @@ namespace Opm { BlackoilModelEbos(Simulator& ebosSimulator, const ModelParameters& param, BlackoilWellModel& well_model, - BlackoilAquiferModel& aquifer_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output) : ebosSimulator_(ebosSimulator) @@ -159,7 +160,6 @@ namespace Opm { , has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy)) , param_( param ) , well_model_ (well_model) - , aquifer_model_(aquifer_model) , terminal_output_ (terminal_output) , current_relaxation_(1.0) , dx_old_(UgGridHelpers::numCells(grid_)) @@ -342,7 +342,6 @@ namespace Opm { void afterStep(const SimulatorTimerInterface& OPM_UNUSED timer) { wellModel().timeStepSucceeded(timer.simulationTimeElapsed()); - aquiferModel().timeStepSucceeded(timer); ebosSimulator_.problem().endTimeStep(); } @@ -360,17 +359,6 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - // -------- Aquifer models ---------- - try - { - // Modify the Jacobian and residuals according to the aquifer models - aquiferModel().assemble(timer, iterationIdx); - } - catch( ... ) - { - OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models"); - } - // -------- Current time step length ---------- const double dt = timer.currentStepLength(); @@ -959,9 +947,6 @@ namespace Opm { // Well Model BlackoilWellModel& well_model_; - // Aquifer Model - BlackoilAquiferModel& aquifer_model_; - /// \brief Whether we print something to std::cout bool terminal_output_; /// \brief The number of cells of the global grid. @@ -981,9 +966,6 @@ namespace Opm { const BlackoilWellModel& wellModel() const { return well_model_; } - BlackoilAquiferModel& - aquiferModel() { return aquifer_model_; } - void beginReportStep() { ebosSimulator_.problem().beginEpisode(); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index bc930c48e..f30c6cde5 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -201,8 +201,6 @@ public: ebosSimulator_.model().addAuxiliaryModule(wellAuxMod_.get()); } - AquiferModel aquifer_model(ebosSimulator_); - // Main simulation loop. while (!timer.done()) { // Report timestep. @@ -217,7 +215,7 @@ public: wellModel.beginReportStep(timer.currentStepNum()); - auto solver = createSolver(wellModel, aquifer_model); + auto solver = createSolver(wellModel); // write the inital state at the report stage if (timer.initialStep()) { @@ -343,12 +341,11 @@ public: protected: - std::unique_ptr createSolver(WellModel& wellModel, AquiferModel& aquifer_model) + std::unique_ptr createSolver(WellModel& wellModel) { auto model = std::unique_ptr(new Model(ebosSimulator_, modelParam_, wellModel, - aquifer_model, linearSolver_, terminalOutput_));