From 72be136529651cc7cf8fb15fd4f3cc850b8fd6ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 12 Oct 2021 08:40:31 +0200 Subject: [PATCH 1/5] Stop using EclWellManager in test_equil.cc. --- tests/test_equil.cc | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/tests/test_equil.cc b/tests/test_equil.cc index 9148e99365..e626bbff5c 100644 --- a/tests/test_equil.cc +++ b/tests/test_equil.cc @@ -26,9 +26,11 @@ #include #include -#include #include +#include +#include + #include #include @@ -64,16 +66,23 @@ namespace Opm::Properties { + +template +struct EnableTerminalOutput { + static constexpr bool value = true; +}; + namespace TTag { + struct TestEquilTypeTag { - using InheritsFrom = std::tuple; + using InheritsFrom = std::tuple; }; } template struct EclWellModel { - using type = EclWellManager; + using type = BlackoilWellModel; }; } // namespace Opm::Properties @@ -199,15 +208,20 @@ namespace { struct EquilFixture { EquilFixture() { - int argc = boost::unit_test::framework::master_test_suite().argc; - char** argv = boost::unit_test::framework::master_test_suite().argv; + int argc = boost::unit_test::framework::master_test_suite().argc; + char** argv = boost::unit_test::framework::master_test_suite().argv; #if HAVE_DUNE_FEM - Dune::Fem::MPIManager::initialize(argc, argv); + Dune::Fem::MPIManager::initialize(argc, argv); #else - Dune::MPIHelper::instance(argc, argv); + Dune::MPIHelper::instance(argc, argv); #endif Opm::EclGenericVanguard::setCommunication(std::make_unique()); using TypeTag = Opm::Properties::TTag::TestEquilTypeTag; + Opm::BlackoilModelParametersEbos::registerParameters(); + Opm::Parameters::registerParam("EnableTerminalOutput", + "EnableTerminalOutput", + Opm::getPropValue(), + "Dummy added for the well model to compile."); Opm::registerAllParameters_(); } }; From 2da4165b3ecde1a57e92e6d92b4d5605ff83ee87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 13 Oct 2021 10:55:05 +0200 Subject: [PATCH 2/5] Stop using EclWellManager in test_ecl_output.cc. --- tests/test_ecl_output.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_ecl_output.cc b/tests/test_ecl_output.cc index 92c6fa1bf8..d175b760fc 100644 --- a/tests/test_ecl_output.cc +++ b/tests/test_ecl_output.cc @@ -26,7 +26,6 @@ #include #include -#include #include #include @@ -41,6 +40,7 @@ #include #include #include +#include #if HAVE_DUNE_FEM #include @@ -86,7 +86,7 @@ struct EnableAsyncEclOutput { template struct EclWellModel { - using type = EclWellManager; + using type = BlackoilWellModel; }; } // namespace Opm::Properties From e6e8eab6d2ea14276f3e1e9d564610c306c415a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 13 Oct 2021 10:56:33 +0200 Subject: [PATCH 3/5] Remove unused EclWellManager class. --- ebos/eclwellmanager.hh | 890 ----------------------------------------- 1 file changed, 890 deletions(-) delete mode 100644 ebos/eclwellmanager.hh diff --git a/ebos/eclwellmanager.hh b/ebos/eclwellmanager.hh deleted file mode 100644 index f067a535ec..0000000000 --- a/ebos/eclwellmanager.hh +++ /dev/null @@ -1,890 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/* - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 2 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . - - Consult the COPYING file in the top-level source directory of this - module for the precise wording of the license and the list of - copyright holders. -*/ -/** - * \file - * - * \copydoc Opm::EclWellManager - */ -#ifndef EWOMS_ECL_WELL_MANAGER_HH -#define EWOMS_ECL_WELL_MANAGER_HH - -#include "eclpeacemanwell.hh" - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -#include - -#include -#include -#include -#include - -namespace Opm { - -/*! - * \ingroup EclBlackOilSimulator - * - * \brief A class which handles well controls as specified by an - * Eclipse deck - */ -template -class EclWellManager -{ - using Simulator = GetPropType; - using GridView = GetPropType; - using Grid = GetPropType; - using Scalar = GetPropType; - using Evaluation = GetPropType; - using FluidSystem = GetPropType; - using ElementContext = GetPropType; - using RateVector = GetPropType; - - enum { numEq = getPropValue() }; - enum { numPhases = FluidSystem::numPhases }; - enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; - enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; - enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; - - typedef typename GridView::template Codim<0>::Entity Element; - - using Well = EclPeacemanWell; - - typedef std::map > > WellConnectionsMap; - - typedef Dune::FieldVector EvalEqVector; - -public: - EclWellManager(Simulator& simulator) - : simulator_(simulator) - { } - - /*! - * \brief This sets up the basic properties of all wells. - * - * I.e., well positions, names etc... - */ - void init() - { - const Schedule& deckSchedule = simulator_.vanguard().schedule(); - const auto& summaryState = simulator_.vanguard().summaryState(); - // create the wells which intersect with the current process' grid - const auto wellsatEnd = deckSchedule.getWellsatEnd(); - for (size_t deckWellIdx = 0; deckWellIdx < deckSchedule.numWells(); ++deckWellIdx) - { - const auto& deckWell = wellsatEnd[deckWellIdx]; - const std::string& wellName = deckWell.name(); - Scalar wellTemperature = 273.15 + 15.56; // [K] - if (deckWell.isInjector()) - wellTemperature = deckWell.injectionControls(summaryState).temperature; - - // set the name of the well but not much else. (i.e., if it is not completed, - // the well primarily serves as a placeholder.) The big rest of the well is - // specified by the updateWellCompletions_() method - auto well = std::make_shared(simulator_); - well->setName(wellName); - well->setWellStatus(Well::Shut); - well->setTemperature(wellTemperature); - - wells_.push_back(well); - wellNameToIndex_[well->name()] = wells_.size() - 1; - } - } - - /*! - * \brief This should be called the problem before each simulation - * episode to adapt the well controls. - */ - void beginEpisode(bool wasRestarted=false) - { - const EclipseState& eclState = simulator_.vanguard().eclState(); - const Schedule& deckSchedule = simulator_.vanguard().schedule(); - const auto& summaryState = simulator_.vanguard().summaryState(); - unsigned episodeIdx = simulator_.episodeIndex(); - WellConnectionsMap wellCompMap; - computeWellConnectionsMap_(episodeIdx, wellCompMap); - - if (wasRestarted || wellTopologyChanged_(eclState, deckSchedule, episodeIdx)) - updateWellTopology_(episodeIdx, wellCompMap, gridDofIsPenetrated_); - - // set those parameters of the wells which do not change the topology of the - // linearized system of equations - updateWellParameters_(episodeIdx, wellCompMap); - - const auto& deckWells = deckSchedule.getWells(episodeIdx); - // set the injection data for the respective wells. - for (const auto& deckWell : deckWells) { - if (!hasWell(deckWell.name())) - continue; - - auto well = this->well(deckWell.name()); - - if (deckWell.isInjector( )) - well->setTemperature(deckWell.injectionControls(summaryState).temperature); - - auto deckWellStatus = deckWell.getStatus( ); - switch (deckWellStatus) { - case ::Opm::Well::Status::AUTO: - // TODO: for now, auto means open... - case ::Opm::Well::Status::OPEN: - well->setWellStatus(Well::Open); - break; - case ::Opm::Well::Status::STOP: - well->setWellStatus(Well::Closed); - break; - case ::Opm::Well::Status::SHUT: - well->setWellStatus(Well::Shut); - break; - } - - // make sure that the well is either an injector or a - // producer for the current episode. (it is not allowed to - // be neither or to be both...) - assert((deckWell.isInjector( )?1:0) + - (deckWell.isProducer( )?1:0) == 1); - - if (deckWell.isInjector( )) { - well->setWellType(Well::Injector); - const auto controls = deckWell.injectionControls(summaryState); - switch (controls.injector_type) { - case InjectorType::WATER: - well->setInjectedPhaseIndex(FluidSystem::waterPhaseIdx); - break; - case InjectorType::GAS: - well->setInjectedPhaseIndex(FluidSystem::gasPhaseIdx); - break; - case InjectorType::OIL: - well->setInjectedPhaseIndex(FluidSystem::oilPhaseIdx); - break; - case InjectorType::MULTI: - throw std::runtime_error("Not implemented: Multi-phase injector wells"); - } - - using InjectorCMode = ::Opm::Well::InjectorCMode; - switch (controls.cmode) { - case InjectorCMode::RATE: - well->setControlMode(Well::ControlMode::VolumetricSurfaceRate); - break; - - case InjectorCMode::RESV: - well->setControlMode(Well::ControlMode::VolumetricReservoirRate); - break; - - case InjectorCMode::BHP: - well->setControlMode(Well::ControlMode::BottomHolePressure); - break; - - case InjectorCMode::THP: - well->setControlMode(Well::ControlMode::TubingHeadPressure); - break; - - case InjectorCMode::GRUP: - throw std::runtime_error("Not implemented: Well groups"); - - case InjectorCMode::CMODE_UNDEFINED: - std::cout << "Warning: Control mode of injection well " << well->name() - << " is undefined. Assuming well to be shut.\n"; - well->setWellStatus(Well::WellStatus::Shut); - continue; - } - - switch (controls.injector_type) { - case InjectorType::WATER: - well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/0.0, /*water=*/1.0); - break; - - case InjectorType::OIL: - well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/0.0, /*water=*/0.0); - break; - - case InjectorType::GAS: - well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/1.0, /*water=*/0.0); - break; - - case InjectorType::MULTI: - throw std::runtime_error("Not implemented: Multi-phase injection wells"); - } - - well->setMaximumSurfaceRate(controls.surface_rate); - well->setMaximumReservoirRate(controls.reservoir_rate); - well->setTargetBottomHolePressure(controls.bhp_limit); - - // TODO - well->setTargetTubingHeadPressure(1e30); - //well->setTargetTubingHeadPressure(controls.thp_limit); - } - - if (deckWell.isProducer( )) { - well->setWellType(Well::Producer); - const auto controls = deckWell.productionControls(summaryState); - - using ProducerCMode = ::Opm::Well::ProducerCMode; - switch (controls.cmode) { - case ProducerCMode::ORAT: - well->setControlMode(Well::ControlMode::VolumetricSurfaceRate); - well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/0.0, /*water=*/0.0); - well->setMaximumSurfaceRate(controls.oil_rate); - break; - - case ProducerCMode::GRAT: - well->setControlMode(Well::ControlMode::VolumetricSurfaceRate); - well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/1.0, /*water=*/0.0); - well->setMaximumSurfaceRate(controls.gas_rate); - break; - - case ProducerCMode::WRAT: - well->setControlMode(Well::ControlMode::VolumetricSurfaceRate); - well->setVolumetricPhaseWeights(/*oil=*/0.0, /*gas=*/0.0, /*water=*/1.0); - well->setMaximumSurfaceRate(controls.water_rate); - break; - - case ProducerCMode::LRAT: - well->setControlMode(Well::ControlMode::VolumetricSurfaceRate); - well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/0.0, /*water=*/1.0); - well->setMaximumSurfaceRate(controls.liquid_rate); - break; - - case ProducerCMode::CRAT: - throw std::runtime_error("Not implemented: Linearly combined rates"); - - case ProducerCMode::RESV: - well->setControlMode(Well::ControlMode::VolumetricReservoirRate); - well->setVolumetricPhaseWeights(/*oil=*/1.0, /*gas=*/1.0, /*water=*/1.0); - well->setMaximumSurfaceRate(controls.resv_rate); - break; - - case ProducerCMode::BHP: - well->setControlMode(Well::ControlMode::BottomHolePressure); - break; - - case ProducerCMode::THP: - well->setControlMode(Well::ControlMode::TubingHeadPressure); - break; - - case ProducerCMode::GRUP: - throw std::runtime_error("Not implemented: Well groups"); - - case ProducerCMode::NONE: - // fall-through - case ProducerCMode::CMODE_UNDEFINED: - std::cout << "Warning: Control mode of production well " << well->name() - << " is undefined. Assuming well to be shut."; - well->setWellStatus(Well::WellStatus::Shut); - continue; - } - - well->setTargetBottomHolePressure(controls.bhp_limit); - - // TODO - well->setTargetTubingHeadPressure(-1e30); - //well->setTargetTubingHeadPressure(controls.thp_limit); - } - } - } - - /*! - * \brief Return the number of wells considered by the EclWellManager. - */ - unsigned numWells() const - { return wells_.size(); } - - /*! - * \brief Return if a given well name is known to the wells manager - */ - bool hasWell(const std::string& wellName) const - { - return wellNameToIndex_.find( wellName ) != wellNameToIndex_.end(); - } - - /*! - * \brief Returns true iff a given degree of freedom is currently penetrated by any well. - */ - bool gridDofIsPenetrated(unsigned globalDofIdx) const - { return gridDofIsPenetrated_[globalDofIdx]; } - - /*! - * \brief Given a well name, return the corresponding index. - * - * A std::runtime_error will be thrown if the well name is unknown. - */ - unsigned wellIndex(const std::string& wellName) const - { - assert( hasWell( wellName ) ); - const auto& it = wellNameToIndex_.find(wellName); - if (it == wellNameToIndex_.end()) - throw std::runtime_error("No well called '"+wellName+"'found"); - return it->second; - } - - /*! - * \brief Given a well name, return the corresponding well. - * - * A std::runtime_error will be thrown if the well name is unknown. - */ - std::shared_ptr well(const std::string& wellName) const - { return wells_[wellIndex(wellName)]; } - - /*! - * \brief Given a well name, return the corresponding well. - * - * A std::runtime_error will be thrown if the well name is unknown. - */ - std::shared_ptr well(const std::string& wellName) - { return wells_[wellIndex(wellName)]; } - - /*! - * \brief Given a well index, return the corresponding well. - */ - std::shared_ptr well(size_t wellIdx) const - { return wells_[wellIdx]; } - - /*! - * \brief Given a well index, return the corresponding well. - */ - std::shared_ptr well(size_t wellIdx) - { return wells_[wellIdx]; } - - /*! - * \brief Informs the well manager that a time step has just begun. - */ - void beginTimeStep() - { - // iterate over all wells and notify them individually - for (size_t wellIdx = 0; wellIdx < wells_.size(); ++wellIdx) - wells_[wellIdx]->beginTimeStep(); - } - - /*! - * \brief Informs the well that an iteration has just begun. - * - * In this method, the well calculates the bottom hole and tubing head pressures, the - * actual unconstraint production and injection rates, etc. - */ - void beginIteration() - { - // call the preprocessing routines - const size_t wellSize = wells_.size(); - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) - wells_[wellIdx]->beginIterationPreProcess(); - - // call the accumulation routines - ElementContext elemCtx(simulator_); - const auto gridView = simulator_.vanguard().gridView(); - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - OPM_BEGIN_PARALLEL_TRY_CATCH(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) - continue; - - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) - wells_[wellIdx]->beginIterationAccumulate(elemCtx, /*timeIdx=*/0); - } - OPM_END_PARALLEL_TRY_CATCH("EclWellManager::beginIteration() failed: ", simulator_.vanguard().grid().comm()); - - // call the postprocessing routines - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) - wells_[wellIdx]->beginIterationPostProcess(); - } - - /*! - * \brief Informs the well manager that an iteration has just been finished. - */ - void endIteration() - { - // iterate over all wells and notify them individually - const size_t wellSize = wells_.size(); - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) - wells_[wellIdx]->endIteration(); - } - - /*! - * \brief Informs the well manager that a time step has just been finished. - */ - void endTimeStep() - { - Scalar dt = simulator_.timeStepSize(); - - // iterate over all wells and notify them individually. also, update the - // production/injection totals for the active wells. - const size_t wellSize = wells_.size(); - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) { - auto well = wells_[wellIdx]; - well->endTimeStep(); - - // update the surface volumes of the produced/injected fluids - std::array* injectedVolume; - if (wellTotalInjectedVolume_.count(well->name()) == 0) { - injectedVolume = &wellTotalInjectedVolume_[well->name()]; - std::fill(injectedVolume->begin(), injectedVolume->end(), 0.0); - } - else - injectedVolume = &wellTotalInjectedVolume_[well->name()]; - - std::array* producedVolume; - if (wellTotalProducedVolume_.count(well->name()) == 0) { - producedVolume = &wellTotalProducedVolume_[well->name()]; - std::fill(producedVolume->begin(), producedVolume->end(), 0.0); - } - else - producedVolume = &wellTotalProducedVolume_[well->name()]; - - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - // this assumes that the implicit Euler method is used for time - // integration. TODO: Once the time discretization becomes pluggable, - // this integration needs to be done by the time discretization code! - Scalar vol = dt * well->surfaceRate(phaseIdx); - - if (vol < 0) - (*producedVolume)[phaseIdx] += -vol; - else - (*injectedVolume)[phaseIdx] += vol; - } - } - } - - /*! - * \brief Informs the well manager that a simulation episode has just been finished. - */ - void endEpisode() - { } - - /*! - * \brief Returns the surface volume of a fluid phase produced by a well. - */ - Scalar totalProducedVolume(const std::string& wellName, unsigned phaseIdx) const - { - if (wellTotalProducedVolume_.count(wellName) == 0) - return 0.0; // well not yet seen - return wellTotalProducedVolume_.at(wellName)[phaseIdx]; - } - - /*! - * \brief Returns the surface volume of a fluid phase injected by a well. - */ - Scalar totalInjectedVolume(const std::string& wellName, unsigned phaseIdx) const - { - if (wellTotalInjectedVolume_.count(wellName) == 0) - return 0.0; // well not yet seen - return wellTotalInjectedVolume_.at(wellName)[phaseIdx]; - } - - /*! - * \brief Computes the source term due to wells for a degree of - * freedom. - */ - template - void computeTotalRatesForDof(EvalEqVector& q, - const Context& context, - unsigned dofIdx, - unsigned timeIdx) const - { - q = 0.0; - - if (!gridDofIsPenetrated(context.globalSpaceIndex(dofIdx, timeIdx))) - return; - - RateVector wellRate; - - // iterate over all wells and add up their individual rates - const size_t wellSize = wells_.size(); - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) { - wellRate = 0.0; - wells_[wellIdx]->computeTotalRatesForDof(wellRate, context, dofIdx, timeIdx); - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) - q[eqIdx] += wellRate[eqIdx]; - } - } - - data::Wells wellData() const - { - data::Wells wellDat; - - using rt = data::Rates::opt; - for (unsigned wellIdx = 0; wellIdx < numWells(); ++wellIdx) { - const auto& ebosWell = well(wellIdx); - auto& wellOut = wellDat[ebosWell->name()]; - - wellOut.bhp = ebosWell->bottomHolePressure(); - wellOut.thp = ebosWell->tubingHeadPressure(); - wellOut.temperature = 0; - wellOut.rates.set( rt::wat, ebosWell->surfaceRate(waterPhaseIdx) ); - wellOut.rates.set( rt::oil, ebosWell->surfaceRate(oilPhaseIdx) ); - wellOut.rates.set( rt::gas, ebosWell->surfaceRate(gasPhaseIdx) ); - - const int numConnections = ebosWell->numConnections(); - wellOut.connections.resize(numConnections); - - for( int i = 0; i < numConnections; ++i ) { - auto& connection = wellOut.connections[ i ]; - connection.index = 0; - connection.pressure = 0.0; - connection.reservoir_rate = 0.0; - connection.rates.set( rt::wat, 0.0 ); - connection.rates.set( rt::oil, 0.0 ); - connection.rates.set( rt::gas, 0.0 ); - } - } - - return wellDat; - } - - data::GroupAndNetworkValues - groupAndNetworkData(const int /* reportStepIdx */) - { - return {}; - } - - /*! - * \brief This method writes the complete state of all wells - * to the hard disk. - */ - template - void serialize(Restarter&) - { - /* do nothing: Everything which we need here is provided by the deck->.. */ - } - - /*! - * \brief This method restores the complete state of the all wells - * from disk. - * - * It is the inverse of the serialize() method. - */ - template - void deserialize(Restarter&) - { - // initialize the wells for the current episode - beginEpisode(/*wasRestarted=*/true); - } - - /*! - * \brief Returns true if something in a well changed compared to the previous report - * step. - * - * "Something" can either be the well topology (i.e., which grid blocks are contained - * in which well) or it can be a well parameter like the bottom hole pressure... - */ - bool wellsChanged(const EclipseState& eclState, - const Schedule& schedule, - unsigned reportStepIdx) const - { - if (wellTopologyChanged_(eclState, reportStepIdx)) - return true; - - if ((schedule.size() - 1) <= (unsigned) reportStepIdx) - // for the "until the universe dies" episode, the wells don't change - return false; - - const Events& events = schedule[reportStepIdx].events(); - return events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE | - ScheduleEvents::INJECTION_UPDATE | - ScheduleEvents::WELL_STATUS_CHANGE ); - } - - void initFromRestartFile(const RestartValue&) { - throw std::logic_error("initFromRestartFile() method not implemented for class eclwellmanager"); - } - - const WellState& wellState() const - { - throw std::logic_error("wellState() method not implemented for class eclwellmanager"); - } - - WellState& wellState() - { - throw std::logic_error("wellState() method not implemented for class eclwellmanager"); - } - - void commitWGState() - { - throw std::logic_error("commitWellState() method not implemented for class eclwellmanager"); - } - - void commitWGState(WGState) - { - throw std::logic_error("commitWGState() method not implemented for class eclwellmanager"); - } - - WellTestState& wellTestState() { - throw std::logic_error("wellTestState() method not implemented for class eclwellmanager"); - } - - - - void resetWGState() - { - throw std::logic_error("resetWGState() method not implemented for class eclwellmanager"); - } - - void updateNupcolWGState() - { - throw std::logic_error("updateNupcolWGState() method not implemented for class eclwellmanager"); - } - - void - updateEclWell(int, int) - { - throw std::logic_error("wellPI() method not implemented for class eclwellmanager"); - } - - - void - updateEclWells(int, const std::unordered_set&) { - throw std::logic_error("wellPI() method not implemented for class eclwellmanager"); - } - - - double - wellPI(int ) const - { - throw std::logic_error("wellPI() method not implemented for class eclwellmanager"); - } - - double - wellPI(const std::string& ) const - { - throw std::logic_error("wellPI() method not implemented for class eclwellmanager"); - } - - - - - -protected: - bool wellTopologyChanged_(const EclipseState&, - const Schedule& schedule, - unsigned reportStepIdx) const - { - if (reportStepIdx == 0) { - // the well topology has always been changed relative to before the - // simulation is started... - return true; - } - - if ((schedule.size() - 1) <= (unsigned) reportStepIdx) - // for the "until the universe dies" episode, the wells don't change - return false; - - const Events& events = schedule[reportStepIdx].events(); - return events.hasEvent(ScheduleEvents::NEW_WELL | - ScheduleEvents::COMPLETION_CHANGE); - } - - void updateWellTopology_(unsigned, - const WellConnectionsMap& wellConnections, - std::vector& gridDofIsPenetrated) const - { - auto& model = simulator_.model(); - const auto& vanguard = simulator_.vanguard(); - - // first, remove all wells from the reservoir - model.clearAuxiliaryModules(); - auto wellIt = wells_.begin(); - const auto& wellEndIt = wells_.end(); - for (; wellIt != wellEndIt; ++wellIt) { - (*wellIt)->clear(); - (*wellIt)->beginSpec(); - } - - ////// - // tell the active wells which DOFs they contain - const auto gridView = simulator_.vanguard().gridView(); - - gridDofIsPenetrated.resize(model.numGridDof()); - std::fill(gridDofIsPenetrated.begin(), gridDofIsPenetrated.end(), false); - - ElementContext elemCtx(simulator_); - auto elemIt = gridView.template begin(); - const auto elemEndIt = gridView.template end(); - std::set > wells; - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) - continue; // non-local entities need to be skipped - - elemCtx.updateStencil(elem); - for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) { - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - unsigned cartesianDofIdx = vanguard.cartesianIndex(globalDofIdx); - - if (wellConnections.count(cartesianDofIdx) == 0) - // the current DOF is not contained in any well, so we must skip - // it... - continue; - - gridDofIsPenetrated[globalDofIdx] = true; - - auto eclWell = wellConnections.at(cartesianDofIdx).second; - eclWell->addDof(elemCtx, dofIdx); - - wells.insert(eclWell); - } - ////// - } - - // register all wells at the model as auxiliary equations - wellIt = wells_.begin(); - for (; wellIt != wellEndIt; ++wellIt) { - (*wellIt)->endSpec(); - model.addAuxiliaryModule(wellIt->get()); - } - } - - void computeWellConnectionsMap_(unsigned reportStepIdx, - WellConnectionsMap& cartesianIdxToConnectionMap) - { - const auto& deckSchedule = simulator_.vanguard().schedule(); - -#ifndef NDEBUG - const auto& eclState = simulator_.vanguard().eclState(); - const auto& eclGrid = eclState.getInputGrid(); - assert( int(eclGrid.getNX()) == simulator_.vanguard().cartesianDimensions()[ 0 ] ); - assert( int(eclGrid.getNY()) == simulator_.vanguard().cartesianDimensions()[ 1 ] ); - assert( int(eclGrid.getNZ()) == simulator_.vanguard().cartesianDimensions()[ 2 ] ); -#endif - - // compute the mapping from logically Cartesian indices to the well the - // respective connection. - const auto deckWells = deckSchedule.getWells(reportStepIdx); - for (const auto& deckWell : deckWells) { - const std::string& wellName = deckWell.name(); - - if (!hasWell(wellName)) - { -#ifndef NDEBUG - if( simulator_.vanguard().grid().comm().size() == 1 ) - { - std::cout << "Well '" << wellName << "' suddenly appears in the connection " - << "for the report step, but has not been previously specified. " - << "Ignoring.\n"; - } -#endif - continue; - } - - std::array cartesianCoordinate; - // set the well parameters defined by the current set of connections - const auto& connectionSet = deckWell.getConnections(); - for (size_t connIdx = 0; connIdx < connectionSet.size(); connIdx ++) { - const auto& connection = connectionSet.get(connIdx); - cartesianCoordinate[ 0 ] = connection.getI(); - cartesianCoordinate[ 1 ] = connection.getJ(); - cartesianCoordinate[ 2 ] = connection.getK(); - unsigned cartIdx = simulator_.vanguard().cartesianIndex( cartesianCoordinate ); - - // in this code we only support each cell to be part of at most a single - // well. TODO (?) change this? - assert(cartesianIdxToConnectionMap.count(cartIdx) == 0); - - auto eclWell = wells_[wellIndex(wellName)]; - cartesianIdxToConnectionMap[cartIdx] = std::make_pair(&connection, eclWell); - } - } - } - - - void updateWellParameters_(unsigned reportStepIdx, const WellConnectionsMap& wellConnections) - { - const auto& deckSchedule = simulator_.vanguard().schedule(); - const auto deckWells = deckSchedule.getWells(reportStepIdx); - - // set the reference depth for all wells - for (const auto& deckWell : deckWells) { - const std::string& wellName = deckWell.name(); - - if( hasWell( wellName ) ) - { - wells_[wellIndex(wellName)]->clear(); - wells_[wellIndex(wellName)]->setReferenceDepth(deckWell.getRefDepth()); - } - } - - // associate the well connections with grid cells and register them in the - // Peaceman well object - const auto& vanguard = simulator_.vanguard(); - const GridView gridView = vanguard.gridView(); - - ElementContext elemCtx(simulator_); - auto elemIt = gridView.template begin(); - const auto elemEndIt = gridView.template end(); - - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) - continue; // non-local entities need to be skipped - - elemCtx.updateStencil(elem); - for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) - { - assert( elemCtx.numPrimaryDof(/*timeIdx=*/0) == 1 ); - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - unsigned cartesianDofIdx = vanguard.cartesianIndex(globalDofIdx); - - if (wellConnections.count(cartesianDofIdx) == 0) - // the current DOF is not contained in any well, so we must skip - // it... - continue; - - const auto& connInfo = wellConnections.at(cartesianDofIdx); - const Connection* connection = connInfo.first; - std::shared_ptr eclWell = connInfo.second; - eclWell->addDof(elemCtx, dofIdx); - eclWell->setConnectionTransmissibilityFactor(elemCtx, dofIdx, connection->CF()); - eclWell->setRadius(elemCtx, dofIdx, connection->rw()); - //eclWell->setEffectivePermeability(elemCtx, dofIdx, connection->Kh()); - } - } - } - - Simulator& simulator_; - - std::vector > wells_; - std::vector gridDofIsPenetrated_; - std::map wellNameToIndex_; - std::map > wellTotalInjectedVolume_; - std::map > wellTotalProducedVolume_; -}; -} // namespace Opm - -#endif From 432f4d8247c1b1f1bdd9ae9588e8255b54f7905c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 13 Oct 2021 11:00:08 +0200 Subject: [PATCH 4/5] Remove unused class EclPeacemanWell --- ebos/eclpeacemanwell.hh | 1782 --------------------------------------- 1 file changed, 1782 deletions(-) delete mode 100644 ebos/eclpeacemanwell.hh diff --git a/ebos/eclpeacemanwell.hh b/ebos/eclpeacemanwell.hh deleted file mode 100644 index 92efdecfa5..0000000000 --- a/ebos/eclpeacemanwell.hh +++ /dev/null @@ -1,1782 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/* - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 2 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . - - Consult the COPYING file in the top-level source directory of this - module for the precise wording of the license and the list of - copyright holders. -*/ -/** - * \file - * - * \copydoc Opm::EclPeacemanWell - */ -#ifndef EWOMS_ECL_PEACEMAN_WELL_HH -#define EWOMS_ECL_PEACEMAN_WELL_HH - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include - -namespace Opm { - -template -class EcfvDiscretization; - -/*! - * \ingroup EclBlackOilSimulator - * - * \brief The well model of Peaceman. - * - * This class is tailored for the element centered finite volume - * discretization, assumes a vertical borehole and is intended to be - * used by the EclWellManager. - * - * See: - * - * Z. Chen, G. Huan, Y. Ma: Computational Methods for Multiphase - * Flows in Porous Media, 1st edition, SIAM, 2006, pp. 445-446 - * - * and - * - * D. W. Peaceman: Interpretation of well-block pressures in numerical - * reservoir simulation, The 52nd Annual SPE Fall Technical Conference - * and Exhibition, Denver, CO., 1977 - */ -template -class EclPeacemanWell : public BaseAuxiliaryModule -{ - using AuxModule = BaseAuxiliaryModule; - - using NeighborSet = typename AuxModule::NeighborSet; - using SparseMatrixAdapter = GetPropType; - using SolutionVector = GetPropType; - using GlobalEqVector = GetPropType; - - using Scalar = GetPropType; - using Evaluation = GetPropType; - using Discretization = GetPropType; - using FluidSystem = GetPropType; - using Simulator = GetPropType; - using ElementContext = GetPropType; - using IntensiveQuantities = GetPropType; - using RateVector = GetPropType; - using GridView = GetPropType; - - typedef MathToolbox Toolbox; - - typedef typename GridView::template Codim<0>::Entity Element; - typedef Element ElementStorage; - - // the dimension of the simulator's world - static const int dimWorld = GridView::dimensionworld; - - // convenient access to the number of phases and the number of - // components - static const unsigned numComponents = getPropValue(); - static const unsigned numPhases = getPropValue(); - - // convenient access to the phase and component indices. If the compiler bails out - // here, you're probably using an incompatible fluid system. This class has only been - // tested with Opm::FluidSystems::BlackOil... - static const unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx; - static const unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx; - static const unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; - - static const unsigned oilCompIdx = FluidSystem::oilCompIdx; - static const unsigned waterCompIdx = FluidSystem::waterCompIdx; - static const unsigned gasCompIdx = FluidSystem::gasCompIdx; - - static const unsigned numModelEq = getPropValue(); - static const unsigned conti0EqIdx = GetPropType::conti0EqIdx; - static const unsigned contiEnergyEqIdx = GetPropType::contiEnergyEqIdx; - - static constexpr unsigned historySize = getPropValue(); - - static constexpr bool enableEnergy = getPropValue(); - - typedef CompositionalFluidState FluidState; - typedef Dune::FieldMatrix DimMatrix; - - // all quantities that need to be stored per degree of freedom that intersects the - // well. - struct DofVariables { - DofVariables() = default; - DofVariables(const DofVariables&) = default; - - // retrieve the solution dependent quantities that are only updated at the - // beginning of a time step from the IntensiveQuantities of the model - void updateBeginTimestep(const IntensiveQuantities&) - {} - - // retrieve the solution dependent quantities from the IntensiveQuantities of the - // model - void update(const IntensiveQuantities& intQuants) - { - const auto& fs = intQuants.fluidState(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - pressure[phaseIdx] = fs.pressure(phaseIdx); - density[phaseIdx] = fs.density(phaseIdx); - mobility[phaseIdx] = intQuants.mobility(phaseIdx); - } - - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - oilMassFraction[compIdx] = fs.massFraction(oilPhaseIdx, compIdx); - gasMassFraction[compIdx] = fs.massFraction(gasPhaseIdx, compIdx); - } - } - - // the depth of the centroid of the DOF - Scalar depth; - - // the volume in m^3 of the DOF - Scalar totalVolume; - - // the effective size of an element in each direction. This is defined as the - // distance of the face centers along the respective axis. - std::array effectiveSize; - - // the intrinsic permeability matrix for the degree of freedom - DimMatrix permeability; - - // the effective permeability of the connection. usually that's the geometric - // mean of the X and Y permeabilities of the DOF times the DOF's height - Scalar effectivePermeability; - - // The connection transmissibility factor to be used for a given DOF. this is - // usually computed from the values above but it can be explicitly specified by - // the user... - Scalar connectionTransmissibilityFactor; - - // the radius of the well for the given degree of freedom - Scalar boreholeRadius; - - // The skin factor of the well at the given degree of freedom - Scalar skinFactor; - - ////////////// - // the following quantities depend on the considered solution and are thus updated - // at the beginning of each Newton-Raphson iteration. - ////////////// - - // the phase pressures inside a DOF - std::array pressure; - - // the phase densities at the DOF - std::array density; - - // the phase mobilities of the DOF - std::array mobility; - - // the composition of the oil phase at the DOF - std::array oilMassFraction; - - // the composition of the gas phase at the DOF - std::array gasMassFraction; - - ElementStorage element; - unsigned pvtRegionIdx; - unsigned localDofIdx; - }; - - // some safety checks/caveats - static_assert(std::is_same >::value, - "The Peaceman well model is only implemented for the " - "element-centered finite volume discretization!"); - static_assert(dimWorld == 3, - "The Peaceman well model is only implemented for 3D grids!"); - -public: - enum ControlMode { - BottomHolePressure, - TubingHeadPressure, - VolumetricSurfaceRate, - VolumetricReservoirRate - }; - - enum WellType { - Undefined, - Injector, - Producer - }; - - enum WellStatus { - // production/injection is ongoing - Open, - - // no production/injection, but well is only closed above the reservoir, so cross - // flow is possible - Closed, - - // well is completely separated from the reservoir, e.g. by filling it with - // concrete. - Shut - }; - - EclPeacemanWell(const Simulator& simulator) - : simulator_(simulator) - { - // set the initial status of the well - wellType_ = Undefined; - wellStatus_ = Shut; - controlMode_ = BottomHolePressure; - - wellTotalVolume_ = 0.0; - - bhpLimit_ = 0.0; - thpLimit_ = 0.0; - - targetBottomHolePressure_ = 0.0; - actualBottomHolePressure_ = 0.0; - maximumSurfaceRate_ = 0.0; - maximumReservoirRate_ = 0.0; - - actualWeightedSurfaceRate_ = 0.0; - actualWeightedResvRate_ = 0.0; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - actualSurfaceRates_[phaseIdx] = 0.0; - actualResvRates_[phaseIdx] = 0.0; - - volumetricWeight_[phaseIdx] = 0.0; - } - - refDepth_ = 0.0; - - // set the composition of the injected fluids based. If - // somebody is stupid enough to inject oil, we assume he wants - // to loose his fortune on dry oil... - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) - injectionFluidState_.setMoleFraction(phaseIdx, compIdx, 0.0); - injectionFluidState_.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0); - injectionFluidState_.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0); - injectionFluidState_.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0); - - // set the temperature to 25 deg C, just so that it is set - injectionFluidState_.setTemperature(273.15 + 25); - - injectedPhaseIdx_ = oilPhaseIdx; - } - - /*! - * \copydoc Opm::BaseAuxiliaryModule::numDofs() - */ - unsigned numDofs() const override - { return 1; } - - /*! - * \copydoc Opm::BaseAuxiliaryModule::addNeighbors() - */ - void addNeighbors(std::vector& neighbors) const override - { - int wellGlobalDof = AuxModule::localToGlobalDof(/*localDofIdx=*/0); - - // the well's bottom hole pressure always affects itself... - neighbors[wellGlobalDof].insert(wellGlobalDof); - - // add the grid DOFs which are influenced by the well, and add the well dof to - // the ones neighboring the grid ones - auto wellDofIt = dofVariables_.begin(); - const auto& wellDofEndIt = dofVariables_.end(); - for (; wellDofIt != wellDofEndIt; ++ wellDofIt) { - neighbors[wellGlobalDof].insert(wellDofIt->first); - neighbors[wellDofIt->first].insert(wellGlobalDof); - } - } - - /*! - * \copydoc Opm::BaseAuxiliaryModule::addNeighbors() - */ - void applyInitial() override - { - auto& sol = const_cast(simulator_.model().solution(/*timeIdx=*/0)); - - int wellGlobalDof = AuxModule::localToGlobalDof(/*localDofIdx=*/0); - sol[wellGlobalDof] = 0.0; - - // make valgrind shut up about the DOFs for the well even if the PrimaryVariables - // class contains some "holes" due to alignment - Valgrind::SetDefined(sol[wellGlobalDof]); - - // also apply the initial solution of the well to the "old" time steps - for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx) { - auto& oldSol = const_cast(simulator_.model().solution(timeIdx)); - - oldSol[wellGlobalDof] = sol[wellGlobalDof]; - } - } - - /*! - * \copydoc Opm::BaseAuxiliaryModule::linearize() - */ - void linearize(SparseMatrixAdapter& matrix, GlobalEqVector& residual) override - { - const SolutionVector& curSol = simulator_.model().solution(/*timeIdx=*/0); - - typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlock; - - unsigned wellGlobalDofIdx = AuxModule::localToGlobalDof(/*localDofIdx=*/0); - residual[wellGlobalDofIdx] = 0.0; - - MatrixBlock diagBlock(0.0); - for (unsigned i = 0; i < numModelEq; ++ i) - diagBlock[i][i] = 1.0; - - MatrixBlock block(0.0); - - if (wellStatus() == Shut) { - // if the well is shut, make the auxiliary DOFs a trivial equation in the - // matrix: the main diagonal is already set to the identity matrix, the - // off-diagonal matrix entries must be set to 0. - auto wellDofIt = dofVariables_.begin(); - const auto& wellDofEndIt = dofVariables_.end(); - for (; wellDofIt != wellDofEndIt; ++ wellDofIt) { - matrix.setBlock(wellGlobalDofIdx, wellDofIt->first, block); - matrix.setBlock(wellDofIt->first, wellGlobalDofIdx, block); - } - matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock); - residual[wellGlobalDofIdx] = 0.0; - return; - } - else if (dofVariables_.empty()) { - // the well does not feature any perforations on the local process - matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock); - residual[wellGlobalDofIdx] = 0.0; - return; - } - - Scalar wellResid = wellResidual_(actualBottomHolePressure_); - residual[wellGlobalDofIdx][0] = wellResid; - - // account for the effect of the grid DOFs which are influenced by the well on - // the well equation and the effect of the well on the grid DOFs - auto wellDofIt = dofVariables_.begin(); - const auto& wellDofEndIt = dofVariables_.end(); - - ElementContext elemCtx(simulator_); - for (; wellDofIt != wellDofEndIt; ++ wellDofIt) { - unsigned gridDofIdx = wellDofIt->first; - const auto& dofVars = *dofVariables_[gridDofIdx]; - DofVariables tmpDofVars(dofVars); - auto priVars(curSol[gridDofIdx]); - - ///////////// - // influence of grid on well - elemCtx.updateStencil( dofVars.element ); - // reset block from previous values - block = 0.0; - for (unsigned priVarIdx = 0; priVarIdx < numModelEq; ++priVarIdx) { - // calculate the derivative of the well equation w.r.t. the current - // primary variable using forward differences - Scalar eps = - 1e3 - *std::numeric_limits::epsilon() - *std::max(1.0, priVars[priVarIdx]); - priVars[priVarIdx] += eps; - - elemCtx.updateIntensiveQuantities(priVars, dofVars.localDofIdx, /*timeIdx=*/0); - tmpDofVars.update(elemCtx.intensiveQuantities(dofVars.localDofIdx, /*timeIdx=*/0)); - - Scalar dWellEq_dPV = - (wellResidual_(actualBottomHolePressure_, &tmpDofVars, gridDofIdx) - wellResid) - / eps; - block[0][priVarIdx] = dWellEq_dPV; - - // go back to the original primary variables - priVars[priVarIdx] -= eps; - } - matrix.setBlock(wellGlobalDofIdx, gridDofIdx, block); - - // - ///////////// - - ///////////// - // influence of well on grid: - RateVector q(0.0); - RateVector modelRate; - std::array resvRates; - - elemCtx.updateIntensiveQuantities(priVars, dofVars.localDofIdx, /*timeIdx=*/0); - - const auto& fluidState = elemCtx.intensiveQuantities(dofVars.localDofIdx, /*timeIdx=*/0).fluidState(); - - // first, we need the source term of the grid for the slightly disturbed well. - Scalar eps = - 1e3 - *std::numeric_limits::epsilon() - *std::max(1e5, actualBottomHolePressure_); - computeVolumetricDofRates_(resvRates, actualBottomHolePressure_ + eps, *dofVariables_[gridDofIdx]); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - modelRate.setVolumetricRate(fluidState, phaseIdx, resvRates[phaseIdx]); - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) - q[compIdx] += modelRate[compIdx]; - } - - // then, we subtract the source rates for a undisturbed well. - computeVolumetricDofRates_(resvRates, actualBottomHolePressure_, *dofVariables_[gridDofIdx]); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - modelRate.setVolumetricRate(fluidState, phaseIdx, resvRates[phaseIdx]); - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) - q[compIdx] -= modelRate[compIdx]; - } - - // and finally, we divide by the epsilon to get the derivative - for (unsigned eqIdx = 0; eqIdx < numModelEq; ++eqIdx) - q[eqIdx] /= eps; - - // now we put this derivative into the right place in the Jacobian - // matrix. This is a bit hacky because it assumes that the model uses a mass - // rate for each component as its first conservation equation, but we require - // the black-oil model for now anyway, so this should not be too much of a - // problem... - Valgrind::CheckDefined(q); - block = 0.0; - for (unsigned eqIdx = 0; eqIdx < numModelEq; ++ eqIdx) - block[eqIdx][0] = - getValue(q[eqIdx])/dofVars.totalVolume; - - matrix.setBlock(gridDofIdx, wellGlobalDofIdx, block); - - // - ///////////// - } - - // effect of changing the well's bottom hole pressure on the well equation - Scalar eps = - 1e3 - *std::numeric_limits::epsilon() - *std::max(1e7, targetBottomHolePressure_); - Scalar wellResidStar = wellResidual_(actualBottomHolePressure_ + eps); - diagBlock[0][0] = (wellResidStar - wellResid)/eps; - - matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock); - } - - Scalar volumetricSurfaceRateForConnection(int globalDofIdx, int phaseIdx) const - { - const DofVariables& dofVars = *dofVariables_.at(globalDofIdx); - std::array volumetricReservoirRates; - computeVolumetricDofRates_(volumetricReservoirRates, actualBottomHolePressure_, dofVars); - std::array volumetricSurfaceRates; - computeSurfaceRates_(volumetricSurfaceRates, volumetricReservoirRates, dofVars); - return volumetricSurfaceRates[phaseIdx]; - } - - - // reset the well to the initial state, i.e. remove all degrees of freedom... - void clear() - { - dofVarsStore_.clear(); - dofVariables_.clear(); - } - - /*! - * \brief Begin the specification of the well. - * - * The specification process is the following: - * - * beginSpec() - * setName("FOO"); - * // add degrees of freedom to the well - * for (dof in wellDofs) - * addDof(dof); - * endSpec() - * - * // set the radius of the well at the dof [m]. - * // optional, if not specified, it is assumed to be 0.1524m - * setRadius(dof, someRadius); - * - * // set the skin factor of the well. - * // optional, if not specified, it is assumed to be 0 - * setSkinFactor(dof, someSkinFactor); - * - * // specify the phase which is supposed to be injected. (Optional, - * // if unspecified, the well will throw an - * // exception if it would inject something.) - * setInjectedPhaseIndex(phaseIdx); - * - * // set maximum production rate at reservoir conditions - * // (kg/s, optional, if not specified, the well is assumed to be - * // shut for production) - * setMaximumReservoirRate(someMassRate); - * - * // set maximum injection rate at reservoir conditions - * // (kg/s, optional, if not specified, the well is assumed to be - * // shut for injection) - * setMinmumReservoirRate(someMassRate); - * - * // set the relative weight of the mass rate of a fluid phase. - * // (Optional, if unspecified each phase exhibits a weight of 1) - * setPhaseWeight(phaseIdx, someWeight); - * - * // set maximum production rate at surface conditions - * // (kg/s, optional, if not specified, the well is assumed to be - * // not limited by the surface rate) - * setMaximumSurfaceRate(someMassRate); - * - * // set maximum production rate at surface conditions - * // (kg/s, optional, if not specified, the well is assumed to be - * // not limited by the surface rate) - * setMinimumSurfaceRate(someMassRate); - * - * // set the minimum pressure at the bottom of the well (Pa, - * // optional, if not specified, the well is assumes it estimates - * // the bottom hole pressure based on the tubing head pressure - * // assuming hydrostatic conditions.) - * setMinimumBottomHolePressure(somePressure); - * - * // set the pressure at the top of the well (Pa, - * // optional, if not specified, the tubing head pressure is - * // assumed to be 1 bar) - * setTubingHeadPressure(somePressure); - * - * // set the control mode of the well [m]. - * // optional, if not specified, it is assumed to be "BottomHolePressure" - * setControlMode(Well::TubingHeadPressure); - * - * // set the tubing head pressure of the well [Pa] - * // only require if the control mode is "TubingHeadPressure" - * setTubingHeadPressure(1e5); - */ - void beginSpec() - { - // this is going to be set to a real value by any realistic grid. Shall we bet? - refDepth_ = 1e100; - - // By default, take the bottom hole pressure as a given - controlMode_ = ControlMode::BottomHolePressure; - - // use one bar for the default bottom hole and tubing head - // pressures. For the bottom hole pressure, this is probably - // off by at least one magnitude... - bhpLimit_ = 1e5; - thpLimit_ = 1e5; - - // reset the actually observed bottom hole pressure - actualBottomHolePressure_ = 0.0; - - // By default, all fluids exhibit the weight 1.0 - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - volumetricWeight_[phaseIdx] = 1.0; - - wellType_ = Undefined; - - wellTotalVolume_ = 0.0; - } - - /*! - * \brief Set the relative weight of the volumetric phase rates. - */ - void setVolumetricPhaseWeights(Scalar oilWeight, Scalar gasWeight, Scalar waterWeight) - { - volumetricWeight_[oilPhaseIdx] = oilWeight; - volumetricWeight_[gasPhaseIdx] = gasWeight; - volumetricWeight_[waterPhaseIdx] = waterWeight; - } - - /*! - * \brief Return the human-readable name of the well - * - * Well, let's say "readable by some humans". - */ - const std::string& name() const - { return name_; } - - /*! - * \brief Set the human-readable name of the well - */ - void setName(const std::string& newName) - { name_ = newName; } - - /*! - * \brief Add a degree of freedom to the well. - */ - template - void addDof(const Context& context, unsigned dofIdx) - { - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - if (applies(globalDofIdx)) - // we already have this DOF in the well! - return; - - const auto& dofPos = context.pos(dofIdx, /*timeIdx=*/0); - - dofVarsStore_.push_back(DofVariables()); - dofVariables_[globalDofIdx] = &dofVarsStore_.back(); - DofVariables& dofVars = *dofVariables_[globalDofIdx]; - wellTotalVolume_ += context.model().dofTotalVolume(globalDofIdx); - - dofVars.element = context.element(); - - dofVars.localDofIdx = dofIdx; - dofVars.pvtRegionIdx = context.problem().pvtRegionIndex(context, dofIdx, /*timeIdx=*/0); - assert(dofVars.pvtRegionIdx == 0); - - // determine the size of the element - dofVars.effectiveSize.fill(0.0); - - // we assume all elements to be hexahedrons! - assert(context.element().subEntities(/*codim=*/dimWorld) == 8); - - const auto& refElem = Dune::ReferenceElements::cube(); - - // determine the current element's effective size - const auto& elem = context.element(); - unsigned faceIdx = 0; - unsigned numFaces = refElem.size(/*codim=*/1); - for (; faceIdx < numFaces; ++faceIdx) { - const auto& faceCenterLocal = refElem.position(faceIdx, /*codim=*/1); - const auto& faceCenter = elem.geometry().global(faceCenterLocal); - - switch (faceIdx) { - case 0: - dofVars.effectiveSize[0] -= faceCenter[0]; - break; - case 1: - dofVars.effectiveSize[0] += faceCenter[0]; - break; - case 2: - dofVars.effectiveSize[1] -= faceCenter[1]; - break; - case 3: - dofVars.effectiveSize[1] += faceCenter[1]; - break; - case 4: - dofVars.depth += faceCenter[2]; - dofVars.effectiveSize[2] -= faceCenter[2]; - break; - case 5: - dofVars.depth += faceCenter[2]; - dofVars.effectiveSize[2] += faceCenter[2]; - break; - } - } - - // the volume associated with the DOF - dofVars.totalVolume = context.model().dofTotalVolume(globalDofIdx); - - // the depth of the degree of freedom - dofVars.depth /= 2; - - // default borehole radius: 1/2 foot - dofVars.boreholeRadius = 0.3048/2; - - // default skin factor: 0 - dofVars.skinFactor = 0; - - // the intrinsic permeability tensor of the DOF - const auto& K = context.problem().intrinsicPermeability(context, dofIdx, /*timeIdx=*/0); - dofVars.permeability = K; - - // default the effective permeability: Geometric mean of the x and y components - // of the intrinsic permeability of DOF times the DOF's height. - assert(K[0][0] > 0); - assert(K[1][1] > 0); - dofVars.effectivePermeability = - std::sqrt(K[0][0]*K[1][1])*dofVars.effectiveSize[2]; - - // from that, compute the default connection transmissibility factor - computeConnectionTransmissibilityFactor_(globalDofIdx); - - // we assume that the z-coordinate represents depth (and not - // height) here... - if (dofPos[2] < refDepth_) - refDepth_ = dofPos[2]; - } - - int numConnections() const - { return dofVariables_.size(); } - - /*! - * \brief Finalize the specification of the borehole. - */ - void endSpec() - { - const auto& comm = simulator_.gridView().comm(); - - int nTotal = dofVariables_.size(); - nTotal = comm.sum(nTotal); - if (nTotal == 0) { - // well does not penetrate any active cell on any process. notify the - // user about this. - std::cout << "Well " << name() << " does not penetrate any active cell." - << " Assuming it to be shut!\n"; - setWellStatus(WellStatus::Shut); - } - - // determine the maximum depth of the well over all processes - refDepth_ = comm.min(refDepth_); - - // the total volume of the well must also be summed over all processes - wellTotalVolume_ = comm.sum(wellTotalVolume_); - } - - /*! - * \brief Set the control mode of the well. - * - * This specifies which quantities are assumed to be externally - * given and which must be calculated based on those. - */ - void setControlMode(ControlMode controlMode) - { controlMode_ = controlMode; } - - /*! - * \brief Set the temperature of the injected fluids [K] - */ - void setTemperature(Scalar value) - { wellTemperature_ = value; } - - /*! - * \brief Set the connection transmissibility factor for a given degree of freedom. - */ - template - void setConnectionTransmissibilityFactor(const Context& context, unsigned dofIdx, Scalar value) - { - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - dofVariables_[globalDofIdx]->connectionTransmissibilityFactor = value; - } - - /*! - * \brief Set the effective permeability Kh to be used for a given degree of freedom. - * - * By default, Kh is sqrt(K_xx * K_yy) * h, where K_xx and K_yy is the permeability - * for the DOF in X and Y directions and h is the height associated with the degree - * of freedom. - * - * Note: The connection transmissibility factor is updated after calling this method, - * so if setConnectionTransmissibilityFactor() is to have any effect, it should - * be called after setEffectivePermeability()! - */ - template - void setEffectivePermeability(const Context& context, unsigned dofIdx, Scalar value) - { - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - dofVariables_[globalDofIdx].effectivePermeability = value; - - computeConnectionTransmissibilityFactor_(globalDofIdx); - } - - /*! - * \brief Set the type of the well (i.e., injector or producer). - */ - void setWellType(WellType wellType) - { wellType_ = wellType; } - - /*! - * \brief Returns the type of the well (i.e., injector or producer). - */ - WellType wellType() const - { return wellType_; } - - /*! - * \brief Set the index of fluid phase to be injected. - * - * This is only relevant if the well type is an injector. - */ - void setInjectedPhaseIndex(unsigned injPhaseIdx) - { injectedPhaseIdx_ = injPhaseIdx; } - - /*! - * \brief Sets the reference depth for the bottom hole pressure [m] - */ - void setReferenceDepth(Scalar value) - { refDepth_ = value; } - - /*! - * \brief The reference depth for the bottom hole pressure [m] - */ - Scalar referenceDepth() const - { return refDepth_; } - - /*! - * \brief Set whether the well is open,closed or shut - */ - void setWellStatus(WellStatus status) - { wellStatus_ = status; } - - /*! - * \brief Return whether the well is open,closed or shut - */ - WellStatus wellStatus() const - { return wellStatus_; } - - /*! - * \brief Return true iff a degree of freedom is directly affected - * by the well - */ - bool applies(unsigned globalDofIdx) const - { return dofVariables_.count(globalDofIdx) > 0; } - - /*! - * \brief Set the maximum/minimum bottom hole pressure [Pa] of the well. - */ - void setTargetBottomHolePressure(Scalar val) - { bhpLimit_ = val; } - - /*! - * \brief Return the maximum/minimum bottom hole pressure [Pa] of the well. - * - * For injectors, this is the maximum, for producers it's the minimum. - */ - Scalar targetBottomHolePressure() const - { return bhpLimit_; } - - /*! - * \brief Return the maximum/minimum bottom hole pressure [Pa] of the well. - */ - Scalar bottomHolePressure() const - { return actualBottomHolePressure_; } - - /*! - * \brief Set the tubing head pressure [Pa] of the well. - */ - void setTargetTubingHeadPressure(Scalar val) - { thpLimit_ = val; } - - /*! - * \brief Return the maximum/minimum tubing head pressure [Pa] of the well. - * - * For injectors, this is the maximum, for producers it's the minimum. - */ - Scalar targetTubingHeadPressure() const - { return thpLimit_; } - - /*! - * \brief Return the maximum/minimum tubing head pressure [Pa] of the well. - */ - Scalar tubingHeadPressure() const - { - // warning: this is a bit hacky... - Scalar rho = 650; // kg/m^3 - Scalar g = 9.81; // m/s^2 - return actualBottomHolePressure_ + rho*refDepth_*g; - } - - /*! - * \brief Set the maximum combined rate of the fluids at the surface. - */ - void setMaximumSurfaceRate(Scalar value) - { maximumSurfaceRate_ = value; } - - /*! - * \brief Return the weighted maximum surface rate [m^3/s] of the well. - */ - Scalar maximumSurfaceRate() const - { return maximumSurfaceRate_; } - - /*! - * \brief Set the maximum combined rate of the fluids at the surface. - */ - void setMaximumReservoirRate(Scalar value) - { maximumReservoirRate_ = value; } - - /*! - * \brief Return the weighted maximum reservoir rate [m^3/s] of the well. - */ - Scalar maximumReservoirRate() const - { return maximumReservoirRate_; } - - /*! - * \brief Return the reservoir rate [m^3/s] actually seen by the well in the current time - * step. - */ - Scalar reservoirRate() const - { return actualWeightedResvRate_; } - - /*! - * \brief Return the weighted surface rate [m^3/s] actually seen by the well in the current time - * step. - */ - Scalar surfaceRate() const - { return actualWeightedSurfaceRate_; } - - /*! - * \brief Return the reservoir rate [m^3/s] of a given fluid which is actually seen - * by the well in the current time step. - */ - Scalar reservoirRate(unsigned phaseIdx) const - { return actualResvRates_[phaseIdx]; } - - /*! - * \brief Return the weighted surface rate [m^3/s] of a given fluid which is actually - * seen by the well in the current time step. - */ - Scalar surfaceRate(unsigned phaseIdx) const - { return actualSurfaceRates_[phaseIdx]; } - - /*! - * \brief Set the skin factor of the well - * - * Note: The connection transmissibility factor is updated after calling this method, - * so if setConnectionTransmissibilityFactor() is to have any effect, it should - * be called after setSkinFactor()! - */ - template - void setSkinFactor(const Context& context, unsigned dofIdx, Scalar value) - { - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - dofVariables_[globalDofIdx].skinFactor = value; - - computeConnectionTransmissibilityFactor_(globalDofIdx); - } - - /*! - * \brief Return the well's skin factor at a DOF [-]. - */ - Scalar skinFactor(unsigned gridDofIdx) const - { return dofVariables_.at(gridDofIdx).skinFactor_; } - - /*! - * \brief Set the borehole radius of the well - * - * Note: The connection transmissibility factor is updated after calling this method, - * so if setConnectionTransmissibilityFactor() is to have any effect, it should - * be called after setRadius()! - */ - template - void setRadius(const Context& context, unsigned dofIdx, Scalar value) - { - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - dofVariables_[globalDofIdx]->boreholeRadius = value; - - computeConnectionTransmissibilityFactor_(globalDofIdx); - } - - /*! - * \brief Return the well's radius at a cell [m]. - */ - Scalar radius(unsigned gridDofIdx) const - { return dofVariables_.at(gridDofIdx)->radius_; } - - /*! - * \brief Informs the well that a time step has just begun. - */ - void beginTimeStep() - { - if (wellStatus() == Shut) - return; - - // calculate the bottom hole pressure to be actually used - if (controlMode_ == ControlMode::TubingHeadPressure) { - // assume a density of 650 kg/m^3 for the bottom hole pressure - // calculation - Scalar rho = 650.0; - targetBottomHolePressure_ = thpLimit_ + rho*refDepth_; - } - else if (controlMode_ == ControlMode::BottomHolePressure) - targetBottomHolePressure_ = bhpLimit_; - else - // TODO: also take the tubing head pressure limit into account... - targetBottomHolePressure_ = bhpLimit_; - - // make it very likely that we screw up if we control for {surface,reservoir} - // rate, but depend on the {reservoir,surface} rate somewhere... - if (controlMode_ == ControlMode::VolumetricSurfaceRate) - maximumReservoirRate_ = 1e100; - else if (controlMode_ == ControlMode::VolumetricReservoirRate) - maximumSurfaceRate_ = 1e100; - - // reset the iteration index - iterationIdx_ = 0; - } - - /*! - * \brief Informs the well that an iteration has just begun. - * - * The beginIteration*() methods, the well calculates the bottom - * and tubing head pressures, the actual unconstraint production and - * injection rates, etc. The callback is split into three parts as - * this arrangement avoids iterating over the whole grid and to - * re-calculate the volume variables for each well. - * - * This is supposed to prepare the well object to do the - * computations which are required to do the DOF specific - * things. - */ - void beginIterationPreProcess() - { } - - /*! - * \brief Do the DOF specific part at the beginning of each iteration - */ - template - void beginIterationAccumulate(Context& context, unsigned timeIdx) - { - if (wellStatus() == Shut) - return; - - for (unsigned dofIdx = 0; dofIdx < context.numPrimaryDof(timeIdx); ++dofIdx) { - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, timeIdx); - if (!applies(globalDofIdx)) - continue; - - DofVariables& dofVars = *dofVariables_.at(globalDofIdx); - const auto& intQuants = context.intensiveQuantities(dofIdx, timeIdx); - - if (iterationIdx_ == 0) - dofVars.updateBeginTimestep(intQuants); - - dofVars.update(intQuants); - } - } - - /*! - * \brief Informs the well that an iteration has just begun. - * - * This is the post-processing part which uses the results of the - * accumulation callback. - */ - void beginIterationPostProcess() - { - if (wellStatus() == Shut) - return; - - auto& sol = const_cast(simulator_.model().solution(/*timeIdx=*/0)); - int wellGlobalDof = AuxModule::localToGlobalDof(/*localDofIdx=*/0); - - if (!dofVariables_.empty()) { - // retrieve the bottom hole pressure from the global system of equations - actualBottomHolePressure_ = Toolbox::value(dofVariables_.begin()->second->pressure[0]); - actualBottomHolePressure_ = computeRateEquivalentBhp_(); - } - else - // start with 300 bars if we don't have anything better - actualBottomHolePressure_ = 300 * 1e5; - - sol[wellGlobalDof][0] = actualBottomHolePressure_; - - computeOverallRates_(actualBottomHolePressure_, - actualResvRates_, - actualSurfaceRates_); - - actualWeightedResvRate_ = computeWeightedRate_(actualResvRates_); - actualWeightedSurfaceRate_ = computeWeightedRate_(actualSurfaceRates_); - } - - /*! - * \brief Called by the simulator after each Newton-Raphson iteration. - */ - void endIteration() - { ++ iterationIdx_; } - - /*! - * \brief Called by the simulator after each time step. - */ - void endTimeStep() - { - if (wellStatus() == Shut) - return; - - // we use a condition that is always false here to prevent the code below from - // bitrotting. (i.e., at least it stays compileable) - if (false && simulator_.gridView().comm().rank() == 0) { - std::cout << "Well '" << name() << "':\n"; - std::cout << " Control mode: " << controlMode_ << "\n"; - std::cout << " BHP limit: " << bhpLimit_/1e5 << " bar\n"; - std::cout << " Observed BHP: " << actualBottomHolePressure_/1e5 << " bar\n"; - std::cout << " Weighted surface rate limit: " << maximumSurfaceRate_ << "\n"; - std::cout << " Weighted surface rate: " << std::abs(actualWeightedSurfaceRate_) << " (=" - << 100*std::abs(actualWeightedSurfaceRate_)/maximumSurfaceRate_ << "%)\n"; - - std::cout << " Surface rates:\n"; - std::cout << " oil: " - << actualSurfaceRates_[oilPhaseIdx] << " m^3/s = " - << actualSurfaceRates_[oilPhaseIdx]*(24*60*60) << " m^3/day = " - << actualSurfaceRates_[oilPhaseIdx]*(24*60*60)/0.15898729 << " STB/day = " - << actualSurfaceRates_[oilPhaseIdx]*(24*60*60) - *FluidSystem::referenceDensity(oilPhaseIdx, /*pvtRegionIdx=*/0) << " kg/day" - << "\n"; - std::cout << " gas: " - << actualSurfaceRates_[gasPhaseIdx] << " m^3/s = " - << actualSurfaceRates_[gasPhaseIdx]*(24*60*60) << " m^3/day = " - << actualSurfaceRates_[gasPhaseIdx]*(24*60*60)/28.316847 << " MCF/day = " - << actualSurfaceRates_[gasPhaseIdx]*(24*60*60) - *FluidSystem::referenceDensity(gasPhaseIdx, /*pvtRegionIdx=*/0) << " kg/day" - << "\n"; - std::cout << " water: " - << actualSurfaceRates_[waterPhaseIdx] << " m^3/s = " - << actualSurfaceRates_[waterPhaseIdx]*(24*60*60) << " m^3/day = " - << actualSurfaceRates_[waterPhaseIdx]*(24*60*60)/0.15898729 << " STB/day = " - << actualSurfaceRates_[waterPhaseIdx]*(24*60*60) - *FluidSystem::referenceDensity(waterPhaseIdx, /*pvtRegionIdx=*/0) << " kg/day" - << "\n"; - } - } - - /*! - * \brief Computes the source term for a degree of freedom. - */ - template - void computeTotalRatesForDof(RateVector& q, - const Context& context, - unsigned dofIdx, - unsigned timeIdx) const - { - q = 0.0; - - unsigned globalDofIdx = context.globalSpaceIndex(dofIdx, timeIdx); - if (wellStatus() == Shut || !applies(globalDofIdx)) - return; - - // create a DofVariables object for the current evaluation point - DofVariables tmp(*dofVariables_.at(globalDofIdx)); - - tmp.update(context.intensiveQuantities(dofIdx, timeIdx)); - - std::array volumetricRates; - computeVolumetricDofRates_(volumetricRates, actualBottomHolePressure_, tmp); - - // convert to mass rates - RateVector modelRate(0.0); - const auto& intQuants = context.intensiveQuantities(dofIdx, timeIdx); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - // energy is disabled or we have production for the given phase, i.e., we - // can use the intensive quantities' fluid state - modelRate.setVolumetricRate(intQuants.fluidState(), phaseIdx, volumetricRates[phaseIdx]); - - if (enableEnergy) { - if (volumetricRates[phaseIdx] < 0.0) { - // producer - const auto& fs = intQuants.fluidState(); - modelRate[contiEnergyEqIdx] += volumetricRates[phaseIdx]*fs.density(phaseIdx)*fs.enthalpy(phaseIdx); - } - else if (volumetricRates[phaseIdx] > 0.0 - && injectedPhaseIdx_ == phaseIdx) - { - // injector for the right phase. we need to use the thermodynamic - // quantities from the borehole as upstream - // - // TODO: This is not implemented in a very efficient way, the - // required quantities could be precomputed at initialization! - auto fs = injectionFluidState_; - - // TODO: maybe we need to use a depth dependent pressure here. the - // difference is probably not very large, and for wells that span - // multiple perforations it is unclear what "well temperature" means - // anyway. - fs.setPressure(phaseIdx, actualBottomHolePressure_); - - fs.setTemperature(wellTemperature_); - - typename FluidSystem::template ParameterCache paramCache; - unsigned globalSpaceIdx = context.globalSpaceIndex(dofIdx, timeIdx); - unsigned pvtRegionIdx = context.primaryVars(dofIdx, timeIdx).pvtRegionIndex(); - paramCache.setRegionIndex(pvtRegionIdx); - paramCache.setMaxOilSat(context.problem().maxOilSaturation(globalSpaceIdx)); - paramCache.updatePhase(fs, phaseIdx); - - const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx); - fs.setDensity(phaseIdx, rho); - - const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx); - fs.setEnthalpy(phaseIdx, h); - - modelRate[contiEnergyEqIdx] += volumetricRates[phaseIdx]*fs.density(phaseIdx)*fs.enthalpy(phaseIdx); - } - } - - for (unsigned eqIdx = 0; eqIdx < modelRate.size(); ++eqIdx) - q[conti0EqIdx + eqIdx] += modelRate[conti0EqIdx + eqIdx]; - } - - Valgrind::CheckDefined(q); - } - -protected: - // compute the connection transmissibility factor based on the effective permeability - // of a connection, the radius of the borehole and the skin factor. - void computeConnectionTransmissibilityFactor_(unsigned globalDofIdx) - { - auto& dofVars = *dofVariables_[globalDofIdx]; - - const auto& D = dofVars.effectiveSize; - const auto& K = dofVars.permeability; - Scalar Kh = dofVars.effectivePermeability; - Scalar S = dofVars.skinFactor; - Scalar rWell = dofVars.boreholeRadius; - - // compute the "equivalence radius" r_0 of the connection - assert(K[0][0] > 0.0); - assert(K[1][1] > 0.0); - Scalar tmp1 = std::sqrt(K[1][1]/K[0][0]); - Scalar tmp2 = 1.0 / tmp1; - Scalar r0 = std::sqrt(D[0]*D[0]*tmp1 + D[1]*D[1]*tmp2); - r0 /= std::sqrt(tmp1) + std::sqrt(tmp2); - r0 *= 0.28; - - // we assume the well borehole in the center of the dof and that it is vertical, - // i.e., the area which is exposed to the flow is 2*pi*r0*h. (for non-vertical - // wells this would need to be multiplied with the cosine of the angle and the - // height must be adapted...) - const Scalar exposureFactor = 2*M_PI; - - dofVars.connectionTransmissibilityFactor = exposureFactor*Kh/(std::log(r0 / rWell) + S); - } - - template - void computeVolumetricDofRates_(std::array& volRates, - const BhpEval& bottomHolePressure, - const DofVariables& dofVars) const - { - typedef MathToolbox DofVarsToolbox; - typedef typename std::conditional::value, - ResultEval, - Scalar>::type DofEval; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - volRates[phaseIdx] = 0.0; - - // connection transmissibility factor for the current DOF. - Scalar Twj = dofVars.connectionTransmissibilityFactor; - - // bottom hole pressure and depth of the degree of freedom - ResultEval pbh = bottomHolePressure; - Scalar depth = dofVars.depth; - - // gravity constant - Scalar g = simulator_.problem().gravity()[dimWorld - 1]; - - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - // well model due to Peaceman; see Chen et al., p. 449 - - // phase pressure in grid cell - const DofEval& p = DofVarsToolbox::template decay(dofVars.pressure[phaseIdx]); - - // density and mobility of fluid phase - const DofEval& rho = DofVarsToolbox::template decay(dofVars.density[phaseIdx]); - DofEval lambda; - if (wellType_ == Producer) { - //assert(p < pbh); - lambda = DofVarsToolbox::template decay(dofVars.mobility[phaseIdx]); - } - else if (wellType_ == Injector) { - //assert(p > pbh); - if (phaseIdx != injectedPhaseIdx_) - continue; - - // use the total mobility, i.e. the sum of all phase mobilities at the - // injector cell. this seems a bit weird: at the wall of the borehole, - // there should only be injected phase present, so its mobility should be - // 1/viscosity... - lambda = 0.0; - for (unsigned phase2Idx = 0; phase2Idx < numPhases; ++phase2Idx) { - if (!FluidSystem::phaseIsActive(phase2Idx)) - continue; - - lambda += DofVarsToolbox::template decay(dofVars.mobility[phase2Idx]); - } - } - else - throw std::logic_error("Type of well \""+name()+"\" is undefined"); - - Valgrind::CheckDefined(pbh); - Valgrind::CheckDefined(p); - Valgrind::CheckDefined(g); - Valgrind::CheckDefined(rho); - Valgrind::CheckDefined(lambda); - Valgrind::CheckDefined(depth); - Valgrind::CheckDefined(refDepth_); - - // pressure in the borehole ("hole pressure") at the given location - ResultEval ph = pbh + rho*g*(depth - refDepth_); - - // volumetric reservoir rate for the phase - volRates[phaseIdx] = Twj*lambda*(ph - p); - - Valgrind::CheckDefined(g); - Valgrind::CheckDefined(ph); - Valgrind::CheckDefined(volRates[phaseIdx]); - } - } - - /*! - * \brief Given the volumetric rates for all phases, return the - * corresponding weighted rate - * - * The weights are user-specified and can be set using - * setVolumetricPhaseWeights() - */ - template - Eval computeWeightedRate_(const std::array& volRates) const - { - Eval result = 0; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - result += volRates[phaseIdx]*volumetricWeight_[phaseIdx]; - } - return result; - } - - /*! - * \brief Convert volumetric reservoir rates into volumetric volume rates. - * - * This requires the density and composition of the phases and - * thus the applicable fluid state. - */ - template - void computeSurfaceRates_(std::array& surfaceRates, - const std::array& reservoirRate, - const DofVariables& dofVars) const - { - // the array for the surface rates and the one for the reservoir rates must not - // be the same! - assert(&surfaceRates != &reservoirRate); - - int regionIdx = dofVars.pvtRegionIdx; - - // If your compiler bails out here, you have not chosen the correct fluid - // system. Currently, only Opm::FluidSystems::BlackOil is supported, sorry... - Scalar rhoOilSurface = FluidSystem::referenceDensity(oilPhaseIdx, regionIdx); - Scalar rhoGasSurface = FluidSystem::referenceDensity(gasPhaseIdx, regionIdx); - Scalar rhoWaterSurface = FluidSystem::referenceDensity(waterPhaseIdx, regionIdx); - - // oil - if (FluidSystem::phaseIsActive(oilPhaseIdx)) - surfaceRates[oilPhaseIdx] = - // oil in gas phase - reservoirRate[gasPhaseIdx] - * Toolbox::value(dofVars.density[gasPhaseIdx]) - * Toolbox::value(dofVars.gasMassFraction[oilCompIdx]) - / rhoOilSurface - + - // oil in oil phase - reservoirRate[oilPhaseIdx] - * Toolbox::value(dofVars.density[oilPhaseIdx]) - * Toolbox::value(dofVars.oilMassFraction[oilCompIdx]) - / rhoOilSurface; - - // gas - if (FluidSystem::phaseIsActive(gasPhaseIdx)) - surfaceRates[gasPhaseIdx] = - // gas in gas phase - reservoirRate[gasPhaseIdx] - * Toolbox::value(dofVars.density[gasPhaseIdx]) - * Toolbox::value(dofVars.gasMassFraction[gasCompIdx]) - / rhoGasSurface - + - // gas in oil phase - reservoirRate[oilPhaseIdx] - * Toolbox::value(dofVars.density[oilPhaseIdx]) - * Toolbox::value(dofVars.oilMassFraction[gasCompIdx]) - / rhoGasSurface; - - // water - if (FluidSystem::phaseIsActive(waterPhaseIdx)) - surfaceRates[waterPhaseIdx] = - reservoirRate[waterPhaseIdx] - * Toolbox::value(dofVars.density[waterPhaseIdx]) - / rhoWaterSurface; - } - - const WellState& wellState() const - { - throw std::logic_error("wellState() method not implemented for class eclpeacemanwell"); - } - - WellState& wellState() - { - throw std::logic_error("wellState() method not implemented for class eclpeacemanwell"); - } - - void commitWGState() - { - throw std::logic_error("commitWellState() method not implemented for class eclpeacemanwell"); - } - - void commitWGState(WGState) - { - throw std::logic_error("commitWellState() method not implemented for class eclpeacemanwell"); - } - - void resetWGState() - { - throw std::logic_error("resetWellState() method not implemented for class eclpeacemanwell"); - } - - void updateNupcolWGState() - { - throw std::logic_error("updateNupcolWellState() method not implemented for class eclpeacemanwell"); - } - - void - updateEclWell(int, int) - { - throw std::logic_error("updateEclWell() method not implemented for class eclpeacemanwell"); - } - - - void - updateEclWells(int, const std::unordered_set&) { - throw std::logic_error("updateEclWells() method not implemented for class eclpeacemanwell"); - } - - - double - wellPI(int) const - { - throw std::logic_error("wellPI() method not implemented for class eclpeacemanwell"); - } - - double - wellPI(const std::string& ) const - { - throw std::logic_error("wellPI() method not implemented for class eclpeacemanwell"); - } - - - - /*! - * \brief Compute the volumetric phase rate of the complete well given a bottom hole - * pressure. - * - * A single degree of freedom may be different from the evaluation point. - */ - void computeOverallRates_(Scalar bottomHolePressure, - std::array& overallResvRates, - std::array& overallSurfaceRates, - const DofVariables *evalDofVars = 0, - int globalEvalDofIdx = -1) const - - { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - overallResvRates[phaseIdx] = 0.0; - overallSurfaceRates[phaseIdx] = 0.0; - } - - auto dofVarsIt = dofVariables_.begin(); - const auto& dofVarsEndIt = dofVariables_.end(); - for (; dofVarsIt != dofVarsEndIt; ++ dofVarsIt) { - std::array volumetricReservoirRates; - const DofVariables *tmp; - if (dofVarsIt->first == globalEvalDofIdx) - tmp = evalDofVars; - else - tmp = dofVarsIt->second; - - computeVolumetricDofRates_(volumetricReservoirRates, bottomHolePressure, *tmp); - - std::array volumetricSurfaceRates; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - volumetricSurfaceRates[ phaseIdx ] = 0; - } - computeSurfaceRates_(volumetricSurfaceRates, volumetricReservoirRates, *tmp); - - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - overallResvRates[phaseIdx] += volumetricReservoirRates[phaseIdx]; - overallSurfaceRates[phaseIdx] += volumetricSurfaceRates[phaseIdx]; - } - } - } - - /*! - * \brief Compute the weighted volumetric rate of the complete well given a bottom - * hole pressure. - * - * A single degree of freedom may be different from the evaluation point. - */ - Scalar computeOverallWeightedSurfaceRate_(Scalar bottomHolePressure, - std::array& overallSurfaceRates, - const DofVariables& evalDofVars, - int globalEvalDofIdx) const - - { - static std::array resvRatesDummy; - computeOverallRates_(bottomHolePressure, - overallSurfaceRates, - resvRatesDummy, - evalDofVars, - globalEvalDofIdx); - return computeWeightedRate_(overallSurfaceRates); - } - - // this is a more convenient version of the method above if all degrees of freedom - // are supposed to be at their evaluation points. - Scalar computeOverallWeightedSurfaceRate_(Scalar bottomHolePressure, - std::array& overallSurfaceRates) const - { - // create a dummy DofVariables object and call the method above using an index - // that is guaranteed to never be part of a well... - static DofVariables dummyDofVars; - return computeOverallWeightedSurfaceRate_(bottomHolePressure, - overallSurfaceRates, - dummyDofVars, - /*globalEvalDofIdx=*/-1); - } - - /*! - * \brief Compute the "rate-equivalent bottom hole pressure" - * - * I.e. The bottom hole pressure where the well rate is exactly the one which is - * targeted. This is zero of the "rate-equivalent bottom hole pressure" would be - * smaller than 1 bar. - */ - Scalar computeRateEquivalentBhp_() const - { - if (wellStatus() == Shut) - // there is no flow happening in the well, so we return 0... - return 0.0; - - // initialize the bottom hole pressure which we would like to calculate - Scalar bhpScalar = actualBottomHolePressure_; - if (bhpScalar > 1e8) - bhpScalar = 1e8; - if (bhpScalar < 1e5) - bhpScalar = 1e5; - - // if the BHP goes below 1 bar for the first time, we reset it to 10 bars and - // are "on bail", i.e. if it goes below 1 bar again, we give up because the - // final pressure would be below 1 bar... - bool onBail = false; - - // Newton-Raphson method - typedef DenseAd::Evaluation BhpEval; - - BhpEval bhpEval(bhpScalar); - bhpEval.setDerivative(0, 1.0); - const Scalar tolerance = 1e3*std::numeric_limits::epsilon(); - const int maxIter = 20; - for (int iterNum = 0; iterNum < maxIter; ++iterNum) { - const auto& f = wellResidual_(bhpEval); - - if (std::abs(f.derivative(0)) < 1e-20) - throw NumericalIssue("Cannot determine the bottom hole pressure for well "+name() - +": Derivative of the well residual is too small"); - Scalar delta = f.value()/f.derivative(0); - - bhpEval.setValue(bhpEval.value() - delta); - if (bhpEval < 1e5) { - bhpEval.setValue(1e5); - if (onBail) - return bhpEval.value(); - else - onBail = true; - } - else - onBail = false; - - if (std::abs(delta/bhpEval.value()) < tolerance) - return bhpEval.value(); - } - - throw NumericalIssue("Could not determine the bottom hole pressure of well '"+name() - +"' within " + std::to_string(maxIter) + " iterations."); - } - - template - BhpEval wellResidual_(const BhpEval& bhp, - const DofVariables *replacementDofVars = 0, - int replacedGridIdx = -1) const - { - typedef MathToolbox BhpEvalToolbox; - - // compute the volumetric reservoir and surface rates for the complete well - BhpEval resvRate = 0.0; - - std::array totalSurfaceRates; - std::fill(totalSurfaceRates.begin(), totalSurfaceRates.end(), 0.0); - - auto dofVarsIt = dofVariables_.begin(); - const auto& dofVarsEndIt = dofVariables_.end(); - for (; dofVarsIt != dofVarsEndIt; ++ dofVarsIt) { - std::array resvRates; - const DofVariables *dofVars = dofVarsIt->second; - if (replacedGridIdx == dofVarsIt->first) - dofVars = replacementDofVars; - computeVolumetricDofRates_(resvRates, bhp, *dofVars); - - std::array surfaceRates; - computeSurfaceRates_(surfaceRates, resvRates, *dofVars); - - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - totalSurfaceRates[phaseIdx] += surfaceRates[phaseIdx]; - } - - resvRate += computeWeightedRate_(resvRates); - } - - BhpEval surfaceRate = computeWeightedRate_(totalSurfaceRates); - - // compute the residual of well equation. we currently use max(rateMax - rate, - // bhp - targetBhp) for producers and max(rateMax - rate, bhp - targetBhp) for - // injectors. (i.e., the target bottom hole pressure is an upper limit for - // injectors and a lower limit for producers.) Note that with this approach, one - // of the limits must always be reached to get the well equation to zero... - Valgrind::CheckDefined(maximumSurfaceRate_); - Valgrind::CheckDefined(maximumReservoirRate_); - Valgrind::CheckDefined(surfaceRate); - Valgrind::CheckDefined(resvRate); - - BhpEval result = 1e30; - - BhpEval maxSurfaceRate = maximumSurfaceRate_; - BhpEval maxResvRate = maximumReservoirRate_; - if (wellStatus() == Closed) { - // make the weight of the fluids on the surface equal and require that no - // fluids are produced on the surface... - maxSurfaceRate = 0.0; - surfaceRate = 0.0; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - surfaceRate += totalSurfaceRates[phaseIdx]; - } - - // don't care about the reservoir rate... - maxResvRate = 1e30; - } - - if (wellType_ == Injector) { - // for injectors the computed rates are positive and the target BHP is the - // maximum allowed pressure ... - result = BhpEvalToolbox::min(maxSurfaceRate - surfaceRate, result); - result = BhpEvalToolbox::min(maxResvRate - resvRate, result); - result = BhpEvalToolbox::min(1e-7*(targetBottomHolePressure_ - bhp), result); - } - else { - assert(wellType_ == Producer); - // ... for producers the rates are negative and the bottom hole pressure is - // is the minimum - result = BhpEvalToolbox::min(maxSurfaceRate + surfaceRate, result); - result = BhpEvalToolbox::min(maxResvRate + resvRate, result); - result = BhpEvalToolbox::min(1e-7*(bhp - targetBottomHolePressure_), result); - } - - const Scalar scalingFactor = 1e-3; - return scalingFactor*result; - } - - const Simulator& simulator_; - - std::string name_; - - std::vector > dofVarsStore_; - std::map dofVariables_; - - // the number of times beginIteration*() was called for the current time step - unsigned iterationIdx_; - - // the type of the well (injector, producer or undefined) - WellType wellType_; - - // Specifies whether the well is currently open, closed or shut. The difference - // between "closed" and "shut" is that for the former, the well is assumed to be - // closed above the reservoir so that cross-flow within the well is possible while - // the well is completely separated from the reservoir if it is shut. (i.e., no - // crossflow is possible in this case.) - WellStatus wellStatus_; - - // specifies the quantities which are controlled for (i.e., which - // should be assumed to be externally specified and which should - // be computed based on those) - ControlMode controlMode_; - - // the sum of the total volumes of all the degrees of freedoms that interact with the well - Scalar wellTotalVolume_; - - // the temperature assumed for the fluid (in the case of an injector well) - Scalar wellTemperature_; - - // The assumed bottom hole and tubing head pressures as specified by the user - Scalar bhpLimit_; - Scalar thpLimit_; - - // The bottom hole pressure to be targeted by the well model. This may be computed - // from the tubing head pressure (if the control mode is TubingHeadPressure), or it may be - // just the user-specified bottom hole pressure if the control mode is - // BottomHolePressure. - Scalar targetBottomHolePressure_; - - // The bottom hole pressure which is actually observed in the well - Scalar actualBottomHolePressure_; - - // The maximum weighted volumetric surface rates specified by the - // user. This is used to apply rate limits and it is to be read as - // the maximum absolute value of the rate, i.e., the well can - // produce or inject the given amount. - Scalar maximumSurfaceRate_; - - // The maximum weighted volumetric reservoir rates specified by - // the user. This is used to apply rate limits and it is to be - // read as the maximum absolute value of the rate, i.e., the well - // can produce or inject the given amount. - Scalar maximumReservoirRate_; - - // The volumetric surface rate which is actually observed in the well - Scalar actualWeightedSurfaceRate_; - std::array actualSurfaceRates_; - - // The volumetric reservoir rate which is actually observed in the well - Scalar actualWeightedResvRate_; - std::array actualResvRates_; - - // The relative weight of the volumetric rate of each fluid - std::array volumetricWeight_; - - // the reference depth for the bottom hole pressure. if not specified otherwise, this - // is the position of the _highest_ DOF in the well. - Scalar refDepth_; - - // The thermodynamic state of the fluid which gets injected - // - // The fact that this attribute is mutable is kind of an hack - // which can be avoided using a PressureOverlayFluidState, but - // then performance would be slightly worse... - mutable FluidState injectionFluidState_; - - unsigned injectedPhaseIdx_; -}; -} // namespace Opm - -#endif From 347ad2da6b99fa823151eef5b897ce2c991b9e63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 13 Oct 2021 11:10:48 +0200 Subject: [PATCH 5/5] Reinstate tests that were commented out due to EclWellManager issues. --- tests/test_ecl_output.cc | 68 +++++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/tests/test_ecl_output.cc b/tests/test_ecl_output.cc index d175b760fc..2409071a8e 100644 --- a/tests/test_ecl_output.cc +++ b/tests/test_ecl_output.cc @@ -41,6 +41,7 @@ #include #include #include +#include #if HAVE_DUNE_FEM #include @@ -70,10 +71,15 @@ namespace Opm::Properties { namespace TTag { struct TestEclOutputTypeTag { - using InheritsFrom = std::tuple; + using InheritsFrom = std::tuple; }; } +template +struct EnableTerminalOutput { + static constexpr bool value = true; +}; + template struct EnableGravity { static constexpr bool value = false; @@ -92,24 +98,24 @@ struct EclWellModel { } // namespace Opm::Properties namespace { -//std::unique_ptr readsum(const std::string& base) -//{ -// return std::make_unique(base); -//} -// -//double ecl_sum_get_field_var(const Opm::EclIO::ESmry* smry, -// const int timeIdx, -// const std::string& var) -//{ -// return smry->get(var)[timeIdx]; -//} -// -//double ecl_sum_get_general_var(const Opm::EclIO::ESmry* smry, -// const int timeIdx, -// const std::string& var) -//{ -// return smry->get(var)[timeIdx]; -//} +std::unique_ptr readsum(const std::string& base) +{ + return std::make_unique(base); +} + +double ecl_sum_get_field_var(const Opm::EclIO::ESmry* smry, + const int timeIdx, + const std::string& var) +{ + return smry->get(var)[timeIdx]; +} + +double ecl_sum_get_general_var(const Opm::EclIO::ESmry* smry, + const int timeIdx, + const std::string& var) +{ + return smry->get(var)[timeIdx]; +} template std::unique_ptr> @@ -132,14 +138,20 @@ initSimulator(const char *filename) struct EclOutputFixture { EclOutputFixture () { - int argc = boost::unit_test::framework::master_test_suite().argc; - char** argv = boost::unit_test::framework::master_test_suite().argv; + int argc = boost::unit_test::framework::master_test_suite().argc; + char** argv = boost::unit_test::framework::master_test_suite().argv; #if HAVE_DUNE_FEM - Dune::Fem::MPIManager::initialize(argc, argv); + Dune::Fem::MPIManager::initialize(argc, argv); #else - Dune::MPIHelper::instance(argc, argv); + Dune::MPIHelper::instance(argc, argv); #endif + Opm::EclGenericVanguard::setCommunication(std::make_unique()); using TypeTag = Opm::Properties::TTag::TestEclOutputTypeTag; + Opm::BlackoilModelParametersEbos::registerParameters(); + Opm::Parameters::registerParam("EnableTerminalOutput", + "EnableTerminalOutput", + Opm::getPropValue(), + "Dummy added for the well model to compile."); Opm::registerAllParameters_(); } }; @@ -150,15 +162,6 @@ BOOST_GLOBAL_FIXTURE(EclOutputFixture); BOOST_AUTO_TEST_CASE(Summary) { - /* - This test is commented out following commit: "Collect WellTestState on - IORank". As part of that commit the CollectTOIORank::collect() method will - access the method WellModel::wellTestState(). That method throws a not - implemented exception in the EclWellManager - and this test instantiates - and uses a EclWellManager deep down. - - ---------------------------------------------------------------------------- - using TypeTag = Opm::Properties::TTag::TestEclOutputTypeTag; const std::string filename = "SUMMARY_DECK_NON_CONSTANT_POROSITY.DATA"; const std::string casename = "SUMMARY_DECK_NON_CONSTANT_POROSITY"; @@ -231,7 +234,6 @@ BOOST_AUTO_TEST_CASE(Summary) // roip = sum_ (b * s * pv) // rs == 0; const double roip2 = ( (0.5 * 0.1 * 100 + 0.6 * 0.2 * 100) * (1 - 0.2) ); BOOST_CHECK_CLOSE(roip2, ecl_sum_get_general_var( resp, 1, "ROIP:2" ), 1e-1 ); - */ } BOOST_AUTO_TEST_CASE(readWriteWells)