From 7c81dbdaab87b7a039e8bc16b3db91c80931eac3 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 25 Oct 2018 16:47:00 +0200 Subject: [PATCH 1/6] let the aquifiers be managed by core ebos also, clean them up a bit: - do not use the intensive quantities cache directly anymore. (i.e., that code should now work if the IQ cache is disabled) - do not fiddle with the global Jacobian matrix and residual vector directly. Instead, implement the water fluxes to the reservoir as a source term like wells. one thing that did not become fully clear to me is if each aquifer ought to be assumed to be in contact with the whole reservoir or just a few cells on the boundary. The current implementation goes down the former path, while, without any deeper knowledge, I would rather suppose that the latter applies. maybe my understanding of this is just too limited, though. --- opm/autodiff/AquiferCarterTracy.hpp | 65 +++++++-------- opm/autodiff/BlackoilAquiferModel.hpp | 53 ++++++------ opm/autodiff/BlackoilAquiferModel_impl.hpp | 81 ++----------------- opm/autodiff/BlackoilModelEbos.hpp | 22 +---- .../SimulatorFullyImplicitBlackoilEbos.hpp | 7 +- 5 files changed, 69 insertions(+), 159 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index c9414faba..e6c75987a 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -45,6 +45,7 @@ namespace Opm typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; 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,8 +62,8 @@ namespace Opm AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, - const Aquancon::AquanconOutput& connection, - Simulator& ebosSimulator ) + const Aquancon::AquanconOutput& connection, + const Simulator& ebosSimulator ) : ebos_simulator_ (ebosSimulator), aquct_data_ (aquct_data), gravity_ (ebos_simulator_.problem().gravity()[2]) @@ -70,56 +71,45 @@ namespace Opm initQuantities(connection); } - inline void assembleAquiferEq(const SimulatorTimerInterface& timer) + template + void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) { - auto& ebosJac = ebos_simulator_.model().linearizer().matrix(); - auto& ebosResid = ebos_simulator_.model().linearizer().residual(); + unsigned idx = context.globalSpaceIndex(spaceIdx, timeIdx); - 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(); + // 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()); - 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); - } - } + rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] -= + Qai_.at(idx)/context.dofVolume(spaceIdx, timeIdx); } - inline void beforeTimeStep(const SimulatorTimerInterface& timer) + inline void beforeTimeStep(const Simulator& simulator) { 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)); + const auto& intQuants = *(simulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); updateCellPressure(pressure_previous_ ,idx,intQuants); } } - inline void afterTimeStep(const SimulatorTimerInterface& timer) + inline void afterTimeStep(const Simulator& simulator) { for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) { - W_flux_ += (*Qai)*timer.currentStepLength(); + W_flux_ += (*Qai)*simulator.timeStepSize(); } } private: - Simulator& ebos_simulator_; + const Simulator& ebos_simulator_; // Grid variables std::vector cell_idx_; @@ -133,6 +123,8 @@ namespace Opm std::vector rhow_; std::vector alphai_; + std::vector aquiferWaterInflux_; + // Variables constants const AquiferCT::AQUCT_data aquct_data_; @@ -164,6 +156,7 @@ namespace Opm calculateAquiferConstants(); + aquiferWaterInflux_.resize(cell_idx_.size()); pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); Qai_.resize(cell_idx_.size(), 0.0); @@ -194,10 +187,10 @@ namespace Opm } // 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 +199,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) ) ); } @@ -365,4 +358,4 @@ namespace Opm } // namespace Opm -#endif \ No newline at end of file +#endif diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index bc831333b..a0928798a 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,46 @@ namespace Opm { /// Class for handling the blackoil well model. template - class BlackoilAquiferModel { + class BlackoilAquiferModel : public Ewoms::EclBaseAquiferModel + { + typedef Ewoms::EclBaseAquiferModel ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; public: + explicit BlackoilAquiferModel(Simulator& ebosSimulator); + // at the beginning of each time step (Not report step) + void beginTimeStep(); + + // add the water rate due to aquifers to the source term. + template + void addToSource(RateVector& rates, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) const + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + aquifer->addToSource(rates, context, spaceIdx, timeIdx); + + } + + protected: // --------- 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); - - protected: - - Simulator& ebosSimulator_; - - std::vector aquifers_; + // TODO: declaring this to be mutable is a hack which should be fixed in the + // long term + mutable std::vector aquifers_; // 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..b5c21838d 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -4,87 +4,22 @@ namespace Opm { template BlackoilAquiferModel:: BlackoilAquiferModel(Simulator& ebosSimulator) - : ebosSimulator_(ebosSimulator) + : ParentType(ebosSimulator) { init(); } - - // called at the end of a time step - template - void - BlackoilAquiferModel:: timeStepSucceeded(const SimulatorTimerInterface& timer) - { - if ( !aquiferActive() ) { - return; - } - - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) - { - aquifer->afterTimeStep(timer); - } - } - - 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); - } - - - template - void - BlackoilAquiferModel:: updateConnectionIntensiveQuantities() const - { - 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); - } - } - - // Protected function which calls the individual aquifer models - template - void - BlackoilAquiferModel:: assembleAquiferEq(const SimulatorTimerInterface& timer) - { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) - { - aquifer->assembleAquiferEq(timer); - } - } - // 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::beginTimeStep() { // 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); + aquifer->beforeTimeStep(this->simulator_); } } @@ -93,14 +28,14 @@ namespace Opm { 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 +50,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_)); From 805eec95663e827c65cb479d7853e2364da51f1f Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 26 Oct 2018 12:25:02 +0200 Subject: [PATCH 2/6] make the CT aquifiers code do something --- opm/autodiff/AquiferCarterTracy.hpp | 63 +++++++++++++++++++-------- opm/autodiff/BlackoilAquiferModel.hpp | 7 +++ 2 files changed, 51 insertions(+), 19 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index e6c75987a..85de3e8d9 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -43,6 +43,7 @@ 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; @@ -63,29 +64,44 @@ namespace Opm AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, const Aquancon::AquanconOutput& connection, - const Simulator& ebosSimulator ) - : ebos_simulator_ (ebosSimulator), - aquct_data_ (aquct_data), - gravity_ (ebos_simulator_.problem().gravity()[2]) + const Simulator& ebosSimulator) + : ebos_simulator_ (ebosSimulator) + , aquct_data_ (aquct_data) + , connection_(connection) + {} + + void initialSolutionApplied() { - initQuantities(connection); + initQuantities(connection_); } template void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) { - unsigned idx = context.globalSpaceIndex(spaceIdx, timeIdx); + unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + +#warning HACK + // check if idx is in cell_idx_ + int connIdx = -1; + for (auto tmp: cell_idx_) { + if (cell_idx_[tmp] == cellIdx) { + connIdx = tmp; + break; + } + } + if (connIdx < 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()); + updateCellPressure(pressure_current_,connIdx,intQuants); + updateCellDensity(connIdx,intQuants); + calculateInflowRate(connIdx, context.simulator()); rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] -= - Qai_.at(idx)/context.dofVolume(spaceIdx, timeIdx); + Qai_.at(connIdx)/context.dofVolume(spaceIdx, timeIdx); } inline void beforeTimeStep(const Simulator& simulator) @@ -128,15 +144,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 @@ -182,7 +200,7 @@ 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; } @@ -314,11 +332,18 @@ 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_); + const auto& elem = *ebos_simulator_.gridView().template begin(); + elemCtx.updateAll(elem); + 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_; @@ -344,7 +369,7 @@ namespace Opm 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) ); + pw_aquifer.push_back( (water_pressure_reservoir - rhow_.at(idx).value()*gravity_()*(cell_depth_.at(idx) - aquct_data_.d0))*alphai_.at(idx) ); } // We take the average of the calculated equilibrium pressures. @@ -352,7 +377,7 @@ namespace Opm return aquifer_pres_avg; } - + const Aquancon::AquanconOutput connection_; }; // class AquiferCarterTracy diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index a0928798a..c1070cd9b 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -46,6 +46,13 @@ namespace Opm { public: explicit BlackoilAquiferModel(Simulator& ebosSimulator); + void initialSolutionApplied() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + aquifer->initialSolutionApplied(); + + } + // at the beginning of each time step (Not report step) void beginTimeStep(); From 339c76bcaca4ce57420a6a2c95c6ca3a521d1b4a Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 30 Oct 2018 11:21:44 +0100 Subject: [PATCH 3/6] make aquifers work again, quite a few cleanups --- opm/autodiff/AquiferCarterTracy.hpp | 137 ++++++++++++--------- opm/autodiff/BlackoilAquiferModel.hpp | 22 ++-- opm/autodiff/BlackoilAquiferModel_impl.hpp | 14 --- 3 files changed, 91 insertions(+), 82 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 85de3e8d9..6e3a19d46 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -65,7 +65,7 @@ namespace Opm AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, const Aquancon::AquanconOutput& connection, const Simulator& ebosSimulator) - : ebos_simulator_ (ebosSimulator) + : simulator_ (ebosSimulator) , aquct_data_ (aquct_data) , connection_(connection) {} @@ -75,20 +75,33 @@ namespace Opm initQuantities(connection_); } + void beginTimeStep() + { + ElementContext elemCtx(simulator_); + auto elemIt = simulator_.gridView().template begin<0>(); + const auto& elemEndIt = simulator_.gridView().template end<0>(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + + elemCtx.updatePrimaryStencil(elem); + + int cellIdx = elemCtx.globalSpaceIndex(0, 0); + int connIdx = cellToConnectionIdx_[cellIdx]; + if (connIdx < 0) + continue; + + elemCtx.updateIntensiveQuantities(0); + const auto& iq = elemCtx.intensiveQuantities(0, 0); + pressure_previous_[connIdx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); + } + } + template void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) { unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); -#warning HACK - // check if idx is in cell_idx_ - int connIdx = -1; - for (auto tmp: cell_idx_) { - if (cell_idx_[tmp] == cellIdx) { - connIdx = tmp; - break; - } - } + int connIdx = cellToConnectionIdx_[cellIdx]; if (connIdx < 0) return; @@ -100,35 +113,23 @@ namespace Opm updateCellDensity(connIdx,intQuants); calculateInflowRate(connIdx, context.simulator()); - rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] -= - Qai_.at(connIdx)/context.dofVolume(spaceIdx, timeIdx); + rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] += + Qai_[connIdx]/context.dofVolume(spaceIdx, timeIdx); } - inline void beforeTimeStep(const Simulator& simulator) + void endTimeStep() { - auto cellID = cell_idx_.begin(); - size_t idx; - for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) - { - const auto& intQuants = *(simulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); - updateCellPressure(pressure_previous_ ,idx,intQuants); - } - } - - inline void afterTimeStep(const Simulator& simulator) - { - for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) - { - W_flux_ += (*Qai)*simulator.timeStepSize(); + for (const auto& Qai: Qai_) { + totalWaterFlux_ += Qai*simulator_.timeStepSize(); } } private: - const Simulator& ebos_simulator_; + const Simulator& simulator_; // Grid variables - std::vector cell_idx_; + std::vector connectionToCellIdx_; std::vector faceArea_connected_; // Quantities at each grid id @@ -149,11 +150,11 @@ namespace Opm Scalar Tc_; // Time constant Scalar pa0_; // initial aquifer pressure - Eval W_flux_; + Eval totalWaterFlux_; Scalar gravity_() const - { return ebos_simulator_.problem().gravity()[2]; } + { return simulator_.problem().gravity()[2]; } inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { @@ -165,7 +166,7 @@ namespace Opm inline void initQuantities(const Aquancon::AquanconOutput& connection) { // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 - W_flux_ = 0.; + totalWaterFlux_ = 0.; // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function initializeConnections(connection); @@ -174,10 +175,10 @@ namespace Opm calculateAquiferConstants(); - aquiferWaterInflux_.resize(cell_idx_.size()); - pressure_previous_.resize(cell_idx_.size(), 0.); - pressure_current_.resize(cell_idx_.size(), 0.); - Qai_.resize(cell_idx_.size(), 0.0); + aquiferWaterInflux_.resize(connectionToCellIdx_.size()); + pressure_previous_.resize(connectionToCellIdx_.size(), 0.); + pressure_current_.resize(connectionToCellIdx_.size(), 0.); + Qai_.resize(connectionToCellIdx_.size(), 0.0); } inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) @@ -212,7 +213,7 @@ namespace Opm Scalar PItdprime = 0.; Scalar PItd = 0.; getInfluenceTableValues(PItd, PItdprime, td_plus_dt); - a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); + a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (totalWaterFlux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); b = beta_ / (Tc_ * ( PItd - td*PItdprime)); } @@ -241,22 +242,22 @@ namespace Opm // 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) { - const auto& eclState = ebos_simulator_.vanguard().eclState(); - const auto& ugrid = ebos_simulator_.vanguard().grid(); + const auto& eclState = simulator_.vanguard().eclState(); + const auto& ugrid = simulator_.vanguard().grid(); const auto& grid = eclState.getInputGrid(); - cell_idx_ = connection.global_index; + connectionToCellIdx_ = connection.global_index; auto globalCellIdx = ugrid.globalCell(); - assert( cell_idx_ == connection.global_index); - assert( (cell_idx_.size() == connection.influx_coeff.size()) ); + assert( connectionToCellIdx_ == connection.global_index); + assert( (connectionToCellIdx_.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 - cell_depth_.resize(cell_idx_.size(), aquct_data_.d0); - alphai_.resize(cell_idx_.size(), 1.0); - faceArea_connected_.resize(cell_idx_.size(),0.0); + cell_depth_.resize(connectionToCellIdx_.size(), aquct_data_.d0); + alphai_.resize(connectionToCellIdx_.size(), 1.0); + faceArea_connected_.resize(connectionToCellIdx_.size(),0.0); Scalar faceArea; auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); @@ -267,10 +268,12 @@ namespace Opm // denom_face_areas is the sum of the areas connected to an aquifer Scalar denom_face_areas = 0.; - for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + cellToConnectionIdx_.resize(simulator_.gridView().size(/*codim=*/0), -1); + for (size_t idx = 0; idx < connectionToCellIdx_.size(); ++idx) { - auto cellFacesRange = cell2Faces[cell_idx_.at(idx)]; - + cellToConnectionIdx_[connectionToCellIdx_[idx]] = idx; + + auto cellFacesRange = cell2Faces[connectionToCellIdx_.at(idx)]; for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) { // The index of the face in the compressed grid @@ -306,11 +309,11 @@ namespace Opm denom_face_areas += ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) ); } } - auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); + auto cellCenter = grid.getCellCenter(connectionToCellIdx_.at(idx)); cell_depth_.at(idx) = cellCenter[2]; } - for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + for (size_t idx = 0; idx < connectionToCellIdx_.size(); ++idx) { alphai_.at(idx) = ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas; } @@ -321,7 +324,7 @@ namespace Opm int pvttableIdx = aquct_data_.pvttableID - 1; - rhow_.resize(cell_idx_.size(),0.); + rhow_.resize(connectionToCellIdx_.size(),0.); if (aquct_data_.p0 < 1.0) { @@ -334,8 +337,8 @@ namespace Opm // use the thermodynamic state of the first active cell as a // reference. there might be better ways to do this... - ElementContext elemCtx(ebos_simulator_); - const auto& elem = *ebos_simulator_.gridView().template begin(); + ElementContext elemCtx(simulator_); + const auto& elem = *simulator_.gridView().template begin(); elemCtx.updateAll(elem); const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); @@ -361,15 +364,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(simulator_); + const auto& gridView = 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 connIdx = cellToConnectionIdx_[cellIdx]; + if (connIdx < 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_[connIdx] = fs.density(waterPhaseIdx); + pw_aquifer.push_back( (water_pressure_reservoir - rhow_[connIdx].value()*gravity_()*(cell_depth_[connIdx] - aquct_data_.d0))*alphai_[connIdx] ); } // We take the average of the calculated equilibrium pressures. @@ -378,6 +392,7 @@ namespace Opm } const Aquancon::AquanconOutput connection_; + std::vector cellToConnectionIdx_; }; // class AquiferCarterTracy diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index c1070cd9b..8a8bf523b 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -53,8 +53,11 @@ namespace Opm { } - // at the beginning of each time step (Not report step) - void beginTimeStep(); + void beginTimeStep() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + aquifer->beginTimeStep(); + } // add the water rate due to aquifers to the source term. template @@ -63,9 +66,14 @@ namespace Opm { unsigned spaceIdx, unsigned timeIdx) const { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) - aquifer->addToSource(rates, context, spaceIdx, timeIdx); + for (auto& aquifer: aquifers_) + aquifer.addToSource(rates, context, spaceIdx, timeIdx); + } + void endTimeStep() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + aquifer->endTimeStep(); } protected: @@ -73,11 +81,11 @@ namespace Opm { typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef AquiferCarterTracy Aquifer_object; + typedef AquiferCarterTracy AquiferType; - // TODO: declaring this to be mutable is a hack which should be fixed in the + // TODO: declaring this as mutable is a hack which should be fixed in the // long term - mutable std::vector aquifers_; + mutable std::vector aquifers_; // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy void init(); diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index b5c21838d..ae384fc1d 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -9,20 +9,6 @@ namespace Opm { init(); } - // 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::beginTimeStep() - { - // 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(this->simulator_); - } - } - // Initialize the aquifers in the deck template void From 46b52448aad45bb453102288cbe73084ce1cbf64 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 30 Oct 2018 15:50:36 +0100 Subject: [PATCH 4/6] fix a few review comments --- opm/autodiff/AquiferCarterTracy.hpp | 96 +++++++++++----------- opm/autodiff/BlackoilAquiferModel.hpp | 33 ++++++-- opm/autodiff/BlackoilAquiferModel_impl.hpp | 4 +- 3 files changed, 72 insertions(+), 61 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 6e3a19d46..8bc834faf 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -65,7 +65,7 @@ namespace Opm AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, const Aquancon::AquanconOutput& connection, const Simulator& ebosSimulator) - : simulator_ (ebosSimulator) + : ebos_simulator_ (ebosSimulator) , aquct_data_ (aquct_data) , connection_(connection) {} @@ -77,22 +77,22 @@ namespace Opm void beginTimeStep() { - ElementContext elemCtx(simulator_); - auto elemIt = simulator_.gridView().template begin<0>(); - const auto& elemEndIt = simulator_.gridView().template end<0>(); + ElementContext elemCtx(ebos_simulator_); + auto elemIt = ebos_simulator_.gridView().template begin<0>(); + const auto& elemEndIt = ebos_simulator_.gridView().template end<0>(); for (; elemIt != elemEndIt; ++elemIt) { const auto& elem = *elemIt; elemCtx.updatePrimaryStencil(elem); int cellIdx = elemCtx.globalSpaceIndex(0, 0); - int connIdx = cellToConnectionIdx_[cellIdx]; - if (connIdx < 0) + int idx = cellToConnectionIdx_[cellIdx]; + if (idx < 0) continue; elemCtx.updateIntensiveQuantities(0); const auto& iq = elemCtx.intensiveQuantities(0, 0); - pressure_previous_[connIdx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); + pressure_previous_[idx] = Opm::getValue(iq.fluidState().pressure(waterPhaseIdx)); } } @@ -101,35 +101,34 @@ namespace Opm { unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - int connIdx = cellToConnectionIdx_[cellIdx]; - if (connIdx < 0) + 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_,connIdx,intQuants); - updateCellDensity(connIdx,intQuants); - calculateInflowRate(connIdx, context.simulator()); + updateCellPressure(pressure_current_,idx,intQuants); + updateCellDensity(idx,intQuants); + calculateInflowRate(idx, context.simulator()); rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] += - Qai_[connIdx]/context.dofVolume(spaceIdx, timeIdx); + Qai_[idx]/context.dofVolume(spaceIdx, timeIdx); } void endTimeStep() { for (const auto& Qai: Qai_) { - totalWaterFlux_ += Qai*simulator_.timeStepSize(); + W_flux_ += Qai*ebos_simulator_.timeStepSize(); } } - private: - const Simulator& simulator_; + const Simulator& ebos_simulator_; // Grid variables - std::vector connectionToCellIdx_; + std::vector cell_idx_; std::vector faceArea_connected_; // Quantities at each grid id @@ -140,8 +139,6 @@ namespace Opm std::vector rhow_; std::vector alphai_; - std::vector aquiferWaterInflux_; - // Variables constants const AquiferCT::AQUCT_data aquct_data_; @@ -150,11 +147,11 @@ namespace Opm Scalar Tc_; // Time constant Scalar pa0_; // initial aquifer pressure - Eval totalWaterFlux_; + Eval W_flux_; Scalar gravity_() const - { return simulator_.problem().gravity()[2]; } + { return ebos_simulator_.problem().gravity()[2]; } inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { @@ -166,7 +163,7 @@ namespace Opm inline void initQuantities(const Aquancon::AquanconOutput& connection) { // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 - totalWaterFlux_ = 0.; + W_flux_ = 0.; // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function initializeConnections(connection); @@ -175,10 +172,9 @@ namespace Opm calculateAquiferConstants(); - aquiferWaterInflux_.resize(connectionToCellIdx_.size()); - pressure_previous_.resize(connectionToCellIdx_.size(), 0.); - pressure_current_.resize(connectionToCellIdx_.size(), 0.); - Qai_.resize(connectionToCellIdx_.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) @@ -213,7 +209,7 @@ namespace Opm Scalar PItdprime = 0.; Scalar PItd = 0.; getInfluenceTableValues(PItd, PItdprime, td_plus_dt); - a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (totalWaterFlux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); + a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); b = beta_ / (Tc_ * ( PItd - td*PItdprime)); } @@ -242,22 +238,22 @@ namespace Opm // 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) { - const auto& eclState = simulator_.vanguard().eclState(); - const auto& ugrid = simulator_.vanguard().grid(); + const auto& eclState = ebos_simulator_.vanguard().eclState(); + const auto& ugrid = ebos_simulator_.vanguard().grid(); const auto& grid = eclState.getInputGrid(); - connectionToCellIdx_ = connection.global_index; + cell_idx_ = connection.global_index; auto globalCellIdx = ugrid.globalCell(); - assert( connectionToCellIdx_ == connection.global_index); - assert( (connectionToCellIdx_.size() == connection.influx_coeff.size()) ); + assert( cell_idx_ == connection.global_index); + assert( (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 - cell_depth_.resize(connectionToCellIdx_.size(), aquct_data_.d0); - alphai_.resize(connectionToCellIdx_.size(), 1.0); - faceArea_connected_.resize(connectionToCellIdx_.size(),0.0); + cell_depth_.resize(cell_idx_.size(), aquct_data_.d0); + alphai_.resize(cell_idx_.size(), 1.0); + faceArea_connected_.resize(cell_idx_.size(),0.0); Scalar faceArea; auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); @@ -268,12 +264,12 @@ namespace Opm // denom_face_areas is the sum of the areas connected to an aquifer Scalar denom_face_areas = 0.; - cellToConnectionIdx_.resize(simulator_.gridView().size(/*codim=*/0), -1); - for (size_t idx = 0; idx < connectionToCellIdx_.size(); ++idx) + cellToConnectionIdx_.resize(ebos_simulator_.gridView().size(/*codim=*/0), -1); + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { - cellToConnectionIdx_[connectionToCellIdx_[idx]] = idx; + cellToConnectionIdx_[cell_idx_[idx]] = idx; - auto cellFacesRange = cell2Faces[connectionToCellIdx_.at(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 @@ -309,11 +305,11 @@ namespace Opm denom_face_areas += ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) ); } } - auto cellCenter = grid.getCellCenter(connectionToCellIdx_.at(idx)); + auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); cell_depth_.at(idx) = cellCenter[2]; } - for (size_t idx = 0; idx < connectionToCellIdx_.size(); ++idx) + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { alphai_.at(idx) = ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas; } @@ -324,7 +320,7 @@ namespace Opm int pvttableIdx = aquct_data_.pvttableID - 1; - rhow_.resize(connectionToCellIdx_.size(),0.); + rhow_.resize(cell_idx_.size(),0.); if (aquct_data_.p0 < 1.0) { @@ -337,8 +333,8 @@ namespace Opm // use the thermodynamic state of the first active cell as a // reference. there might be better ways to do this... - ElementContext elemCtx(simulator_); - const auto& elem = *simulator_.gridView().template begin(); + ElementContext elemCtx(ebos_simulator_); + const auto& elem = *ebos_simulator_.gridView().template begin(); elemCtx.updateAll(elem); const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); @@ -364,8 +360,8 @@ namespace Opm std::vector pw_aquifer; Scalar water_pressure_reservoir; - ElementContext elemCtx(simulator_); - const auto& gridView = simulator_.gridView(); + 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) { @@ -373,8 +369,8 @@ namespace Opm elemCtx.updatePrimaryStencil(elem); size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - int connIdx = cellToConnectionIdx_[cellIdx]; - if (connIdx < 0) + int idx = cellToConnectionIdx_[cellIdx]; + if (idx < 0) continue; elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); @@ -382,8 +378,8 @@ namespace Opm const auto& fs = iq0.fluidState(); water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); - rhow_[connIdx] = fs.density(waterPhaseIdx); - pw_aquifer.push_back( (water_pressure_reservoir - rhow_[connIdx].value()*gravity_()*(cell_depth_[connIdx] - aquct_data_.d0))*alphai_[connIdx] ); + 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. diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 8a8bf523b..c24a782b6 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -36,29 +36,34 @@ namespace Opm { /// Class for handling the blackoil well model. template - class BlackoilAquiferModel : public Ewoms::EclBaseAquiferModel + class BlackoilAquiferModel { - typedef Ewoms::EclBaseAquiferModel ParentType; - typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; public: - explicit BlackoilAquiferModel(Simulator& ebosSimulator); + explicit BlackoilAquiferModel(Simulator& simulator); void initialSolutionApplied() { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { aquifer->initialSolutionApplied(); - + } } + void beginEpisode() + { } + void beginTimeStep() { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { aquifer->beginTimeStep(); + } } + void beginIteration() + { } + // add the water rate due to aquifers to the source term. template void addToSource(RateVector& rates, @@ -66,16 +71,24 @@ namespace Opm { unsigned spaceIdx, unsigned timeIdx) const { - for (auto& aquifer: aquifers_) + for (auto& aquifer: aquifers_) { aquifer.addToSource(rates, context, spaceIdx, timeIdx); + } } + void endIteration() + { } + void endTimeStep() { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { aquifer->endTimeStep(); + } } + void endEpisode() + { } + protected: // --------- Types --------- typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; @@ -87,6 +100,8 @@ namespace Opm { // 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(); diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index ae384fc1d..3d9525c7d 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -3,8 +3,8 @@ namespace Opm { template BlackoilAquiferModel:: - BlackoilAquiferModel(Simulator& ebosSimulator) - : ParentType(ebosSimulator) + BlackoilAquiferModel(Simulator& simulator) + : simulator_(simulator) { init(); } From 6a929ca7972bba44bd1474ba2667eb1cd5ef153b Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 1 Nov 2018 14:54:30 +0100 Subject: [PATCH 5/6] try to fix subtle lifetime bug it is possible that a dune entity vanishes if its iterator gets out of scope. Whether this is a problem or not seems to be be highly depend on the used configuration... --- opm/autodiff/AquiferCarterTracy.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 8bc834faf..4d6cad725 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -334,8 +334,9 @@ namespace Opm // use the thermodynamic state of the first active cell as a // reference. there might be better ways to do this... ElementContext elemCtx(ebos_simulator_); - const auto& elem = *ebos_simulator_.gridView().template begin(); - elemCtx.updateAll(elem); + 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 From 848a65c1c8d801cc5f92a9d94f5d47d545dd35b4 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 1 Nov 2018 15:19:33 +0100 Subject: [PATCH 6/6] BlackoilAquiferModel: move the function bodies to the _impl.hpp file ... to make it more consistent with the existing code. --- opm/autodiff/BlackoilAquiferModel.hpp | 47 +++-------------- opm/autodiff/BlackoilAquiferModel_impl.hpp | 60 ++++++++++++++++++++++ 2 files changed, 68 insertions(+), 39 deletions(-) diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index c24a782b6..ac671ec67 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -44,50 +44,19 @@ namespace Opm { public: explicit BlackoilAquiferModel(Simulator& simulator); - void initialSolutionApplied() - { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->initialSolutionApplied(); - } - } - - void beginEpisode() - { } - - void beginTimeStep() - { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->beginTimeStep(); - } - } - - void beginIteration() - { } - + 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 - { - for (auto& aquifer: aquifers_) { - aquifer.addToSource(rates, context, spaceIdx, timeIdx); - } - } - - void endIteration() - { } - - void endTimeStep() - { - for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->endTimeStep(); - } - } - - void endEpisode() - { } + unsigned timeIdx) const; + void endIteration(); + void endTimeStep(); + void endEpisode(); protected: // --------- Types --------- diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 3d9525c7d..2c5646746 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -9,6 +9,66 @@ namespace Opm { init(); } + template + void + BlackoilAquiferModel::initialSolutionApplied() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { + aquifer->initialSolutionApplied(); + } + } + + template + void + BlackoilAquiferModel::beginEpisode() + { } + + template + void + BlackoilAquiferModel::beginTimeStep() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { + aquifer->beginTimeStep(); + } + } + + template + void + BlackoilAquiferModel::beginIteration() + { } + + template + template + void + BlackoilAquiferModel::addToSource(RateVector& rates, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) const + { + for (auto& aquifer: aquifers_) { + aquifer.addToSource(rates, context, spaceIdx, timeIdx); + } + } + + template + void + BlackoilAquiferModel::endIteration() + { } + + template + void + BlackoilAquiferModel::endTimeStep() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { + aquifer->endTimeStep(); + } + } + + template + void + BlackoilAquiferModel::endEpisode() + { } + // Initialize the aquifers in the deck template void