diff --git a/CMakeLists.txt b/CMakeLists.txt index e803b3cac..03efb9ecf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -345,7 +345,7 @@ set_property(TARGET moduleVersion PROPERTY POSITION_INDEPENDENT_CODE ON) add_dependencies(moduleVersion opmsimulators) set(COMMON_MODELS brine energy extbo foam gasoil gaswater oilwater oilwater_polymer polymer solvent) -set(FLOW_MODELS blackoil oilwater_brine oilwater_polymer_injectivity) +set(FLOW_MODELS blackoil oilwater_brine oilwater_polymer_injectivity micp) set(FLOW_TGTS) foreach(OBJ ${COMMON_MODELS} ${FLOW_MODELS}) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index b0da7681e..c047b6517 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -44,6 +44,8 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/FlexibleSolver2.cpp opm/simulators/linalg/FlexibleSolver3.cpp opm/simulators/linalg/FlexibleSolver4.cpp + opm/simulators/linalg/FlexibleSolver5.cpp + opm/simulators/linalg/FlexibleSolver6.cpp opm/simulators/linalg/PropertyTree.cpp opm/simulators/linalg/setupPropertyTree.cpp opm/simulators/utils/PartiallySupportedFlowKeywords.cpp diff --git a/ebos/ebos_gasoil.cc b/ebos/ebos_gasoil.cc index 41fd44bda..ab02322fe 100644 --- a/ebos/ebos_gasoil.cc +++ b/ebos/ebos_gasoil.cc @@ -57,7 +57,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::waterCompIdx> type; + /*disabledCompIdx=*/FluidSystem::waterCompIdx, + getPropValue()> type; }; } // namespace Opm::Properties diff --git a/ebos/ebos_gaswater.cc b/ebos/ebos_gaswater.cc index c54c9fc65..8d9958e65 100644 --- a/ebos/ebos_gaswater.cc +++ b/ebos/ebos_gaswater.cc @@ -60,7 +60,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::oilCompIdx> type; + /*disabledCompIdx=*/FluidSystem::oilCompIdx, + getPropValue()> type; }; } // namespace Opm::Properties diff --git a/ebos/ebos_oilwater.cc b/ebos/ebos_oilwater.cc index fdb636f61..807c485f9 100644 --- a/ebos/ebos_oilwater.cc +++ b/ebos/ebos_oilwater.cc @@ -56,7 +56,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::gasCompIdx> type; + /*disabledCompIdx=*/FluidSystem::gasCompIdx, + getPropValue()> type; }; } // namespace Opm::Properties diff --git a/ebos/ebos_oilwater_polymer.cc b/ebos/ebos_oilwater_polymer.cc index 655eea51b..3a5deb387 100644 --- a/ebos/ebos_oilwater_polymer.cc +++ b/ebos/ebos_oilwater_polymer.cc @@ -61,7 +61,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::gasCompIdx> type; + /*disabledCompIdx=*/FluidSystem::gasCompIdx, + getPropValue()> type; }; } // namespace Opm::Properties diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc index fd0d3e073..1aa44d18f 100644 --- a/ebos/eclgenericoutputblackoilmodule.cc +++ b/ebos/eclgenericoutputblackoilmodule.cc @@ -2,20 +2,16 @@ // 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. @@ -79,7 +75,8 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState, bool enablePolymer, bool enableFoam, bool enableBrine, - bool enableExtbo) + bool enableExtbo, + bool enableMICP) : eclState_(eclState) , schedule_(schedule) , summaryConfig_(summaryConfig) @@ -91,6 +88,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState, , enableFoam_(enableFoam) , enableBrine_(enableBrine) , enableExtbo_(enableExtbo) + , enableMICP_(enableMICP) { const auto& fp = eclState_.fieldProps(); @@ -592,15 +590,19 @@ assignToSolution(data::Solution& sol) {"1OVERBG", UnitSystem::measure::gas_inverse_formation_volume_factor, data::TargetType::RESTART_AUXILIARY, invB_[gasPhaseIdx]}, {"1OVERBO", UnitSystem::measure::oil_inverse_formation_volume_factor, data::TargetType::RESTART_AUXILIARY, invB_[oilPhaseIdx]}, {"1OVERBW", UnitSystem::measure::water_inverse_formation_volume_factor, data::TargetType::RESTART_AUXILIARY, invB_[waterPhaseIdx]}, + {"BIOFILM", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, cBiofilm_}, + {"CALCITE", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, cCalcite_}, {"FOAM", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, cFoam_}, {"GASKR", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, relativePermeability_[gasPhaseIdx]}, {"GAS_DEN", UnitSystem::measure::density, data::TargetType::RESTART_AUXILIARY, density_[gasPhaseIdx]}, {"GAS_VISC", UnitSystem::measure::viscosity, data::TargetType::RESTART_AUXILIARY, viscosity_[gasPhaseIdx]}, {"KRNSW_GO", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, krnSwMdcGo_}, {"KRNSW_OW", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, krnSwMdcOw_}, + {"MICROBES", UnitSystem::measure::density, data::TargetType::RESTART_SOLUTION, cMicrobes_}, {"OILKR", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, relativePermeability_[oilPhaseIdx]}, {"OIL_DEN", UnitSystem::measure::density, data::TargetType::RESTART_AUXILIARY, density_[oilPhaseIdx]}, {"OIL_VISC", UnitSystem::measure::viscosity, data::TargetType::RESTART_AUXILIARY, viscosity_[oilPhaseIdx]}, + {"OXYGEN", UnitSystem::measure::density, data::TargetType::RESTART_SOLUTION, cOxygen_}, {"PBUB", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, bubblePointPressure_}, {"PCSWM_GO", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, pcSwMdcGo_}, {"PCSWM_OW", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, pcSwMdcOw_}, @@ -625,6 +627,7 @@ assignToSolution(data::Solution& sol) {"STD_GAS", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, mFracGas_}, {"STD_OIL", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, mFracOil_}, {"SWMAX", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, swMax_}, + {"UREA", UnitSystem::measure::density, data::TargetType::RESTART_SOLUTION, cUrea_}, {"TMULT_RC", UnitSystem::measure::identity, data::TargetType::RESTART_SOLUTION, rockCompTransMultiplier_}, {"WATKR", UnitSystem::measure::identity, data::TargetType::RESTART_AUXILIARY, relativePermeability_[waterPhaseIdx]}, {"WAT_DEN", UnitSystem::measure::density, data::TargetType::RESTART_AUXILIARY, density_[waterPhaseIdx]}, @@ -734,6 +737,16 @@ setRestart(const data::Solution& sol, krnSwMdcGo_[elemIdx] = sol.data("KRNSW_GO")[globalDofIndex]; if (!ppcw_.empty() && sol.has("PPCW")) ppcw_[elemIdx] = sol.data("PPCW")[globalDofIndex]; + if (!cMicrobes_.empty() && sol.has("MICROBES")) + cMicrobes_[elemIdx] = sol.data("MICROBES")[globalDofIndex]; + if (!cOxygen_.empty() && sol.has("OXYGEN")) + cOxygen_[elemIdx] = sol.data("OXYGEN")[globalDofIndex]; + if (!cUrea_.empty() && sol.has("UREA")) + cUrea_[elemIdx] = sol.data("UREA")[globalDofIndex]; + if (!cBiofilm_.empty() && sol.has("BIOFILM")) + cBiofilm_[elemIdx] = sol.data("BIOFILM")[globalDofIndex]; + if (!cCalcite_.empty() && sol.has("CALCITE")) + cCalcite_[elemIdx] = sol.data("CALCITE")[globalDofIndex]; } template @@ -913,6 +926,13 @@ doAllocBuffers(unsigned bufferSize, mFracGas_.resize(bufferSize, 0.0); mFracCo2_.resize(bufferSize, 0.0); } + if (enableMICP_){ + cMicrobes_.resize(bufferSize, 0.0); + cOxygen_.resize(bufferSize, 0.0); + cUrea_.resize(bufferSize, 0.0); + cBiofilm_.resize(bufferSize, 0.0); + cCalcite_.resize(bufferSize, 0.0); + } if (vapparsActive) soMax_.resize(bufferSize, 0.0); @@ -1169,14 +1189,14 @@ outputResvFluidInPlace_(std::unordered_map cipr, const i } ss << ":---------:---------------:---------------:---------------:---------------:---------------:\n" << ": REGION : TOTAL PORE : PORE VOLUME : PORE VOLUME : PORE VOLUME : PORE VOLUME :\n" - << ": : VOLUME : CONTAINING : CONTAINING : CONTAINING : CONTAINING :\n" - << ": : : OIL : WATER : GAS : HYDRO-CARBON :\n" + << ": : VOLUME : CONTAINING : CONTAINING : CONTAINING : CONTAINING :\n" + << ": : : OIL : WATER : GAS : HYDRO-CARBON :\n" << ":---------:---------------:---------------:---------------:---------------:---------------\n"; - } - else { + } + else { ss << std::right << std::fixed << std::setprecision(0) << ":" << std::setw (9) << reg << ":" << std::setw(15) << cipr[Inplace::Phase::DynamicPoreVolume] << ":" << std::setw(15) << cipr[Inplace::Phase::OilResVolume] << ":" << std::setw(15) << cipr[Inplace::Phase::WaterResVolume] << ":" << std::setw(15) << cipr[Inplace::Phase::GasResVolume] << ":" << std::setw(15) << cipr[Inplace::Phase::OilResVolume] + cipr[Inplace::Phase::GasResVolume] << ":\n" << ":---------:---------------:---------------:---------------:---------------:---------------:\n"; - } + } OpmLog::note(ss.str()); } diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh index bc375e0aa..caab55e2f 100644 --- a/ebos/eclgenericoutputblackoilmodule.hh +++ b/ebos/eclgenericoutputblackoilmodule.hh @@ -2,20 +2,16 @@ // 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. @@ -132,6 +128,46 @@ public: return 0; } + Scalar getMicrobialConcentration(unsigned elemIdx) const + { + if (cMicrobes_.size() > elemIdx) + return cMicrobes_[elemIdx]; + + return 0; + } + + Scalar getOxygenConcentration(unsigned elemIdx) const + { + if (cOxygen_.size() > elemIdx) + return cOxygen_[elemIdx]; + + return 0; + } + + Scalar getUreaConcentration(unsigned elemIdx) const + { + if (cUrea_.size() > elemIdx) + return cUrea_[elemIdx]; + + return 0; + } + + Scalar getBiofilmConcentration(unsigned elemIdx) const + { + if (cBiofilm_.size() > elemIdx) + return cBiofilm_[elemIdx]; + + return 0; + } + + Scalar getCalciteConcentration(unsigned elemIdx) const + { + if (cCalcite_.size() > elemIdx) + return cCalcite_[elemIdx]; + + return 0; + } + const std::map& getWBPData() const { return this->wbpData_; @@ -168,7 +204,8 @@ protected: bool enablePolymer, bool enableFoam, bool enableBrine, - bool enableExtbo); + bool enableExtbo, + bool enableMICP); struct WellProdDataType { @@ -324,6 +361,7 @@ protected: bool enableFoam_; bool enableBrine_; bool enableExtbo_; + bool enableMICP_; bool forceDisableFipOutput_; bool forceDisableFipresvOutput_; @@ -376,6 +414,11 @@ protected: ScalarBuffer minimumOilPressure_; ScalarBuffer saturatedOilFormationVolumeFactor_; ScalarBuffer rockCompTransMultiplier_; + ScalarBuffer cMicrobes_; + ScalarBuffer cOxygen_; + ScalarBuffer cUrea_; + ScalarBuffer cBiofilm_; + ScalarBuffer cCalcite_; std::array saturation_; std::array invB_; diff --git a/ebos/eclgenericproblem.cc b/ebos/eclgenericproblem.cc index f7986acc0..f7ab2bb86 100644 --- a/ebos/eclgenericproblem.cc +++ b/ebos/eclgenericproblem.cc @@ -407,7 +407,8 @@ checkDeckCompatibility_(const Deck& deck, int numPhases, bool indicesGasEnabled, bool indicesOilEnabled, - bool indicesWaterEnabled) const + bool indicesWaterEnabled, + bool enableMICP) const { if (enableApiTracking) throw std::logic_error("API tracking is not yet implemented but requested at compile time."); @@ -424,6 +425,11 @@ checkDeckCompatibility_(const Deck& deck, else if (!enablePolymer && deck.hasKeyword("POLYMER")) throw std::runtime_error("The deck enables the polymer option, but the simulator is compiled without it."); + if (enableMICP && !deck.hasKeyword("MICP")) + throw std::runtime_error("The simulator requires the MICP option to be enabled, but the deck does not."); + else if (!enableMICP && deck.hasKeyword("MICP")) + throw std::runtime_error("The deck enables the MICP option, but the simulator is compiled without it."); + if (enableExtbo && !deck.hasKeyword("PVTSOL")) throw std::runtime_error("The simulator requires the extendedBO option to be enabled, but the deck does not."); else if (!enableExtbo && deck.hasKeyword("PVTSOL")) @@ -476,7 +482,8 @@ void EclGenericProblem:: readBlackoilExtentionsInitialConditions_(size_t numDof, bool enableSolvent, bool enablePolymer, - bool enablePolymerMolarWeight) + bool enablePolymerMolarWeight, + bool enableMICP) { if (enableSolvent) { if (eclState_.fieldProps().has_double("SSOL")) @@ -498,6 +505,29 @@ readBlackoilExtentionsInitialConditions_(size_t numDof, else polymerMoleWeight_.resize(numDof, 0.0); } + + if (enableMICP) { + if (eclState_.fieldProps().has_double("SMICR")) + microbialConcentration_ = eclState_.fieldProps().get_double("SMICR"); + else + microbialConcentration_.resize(numDof, 0.0); + if (eclState_.fieldProps().has_double("SOXYG")) + oxygenConcentration_ = eclState_.fieldProps().get_double("SOXYG"); + else + oxygenConcentration_.resize(numDof, 0.0); + if (eclState_.fieldProps().has_double("SUREA")) + ureaConcentration_ = eclState_.fieldProps().get_double("SUREA"); + else + ureaConcentration_.resize(numDof, 0.0); + if (eclState_.fieldProps().has_double("SBIOF")) + biofilmConcentration_ = eclState_.fieldProps().get_double("SBIOF"); + else + biofilmConcentration_.resize(numDof, 0.0); + if (eclState_.fieldProps().has_double("SCALC")) + calciteConcentration_ = eclState_.fieldProps().get_double("SCALC"); + else + calciteConcentration_.resize(numDof, 0.0); +} } @@ -561,6 +591,56 @@ polymerMolecularWeight(const unsigned elemIdx) const return polymerMoleWeight_[elemIdx]; } +template +Scalar EclGenericProblem:: +microbialConcentration(unsigned elemIdx) const +{ + if (microbialConcentration_.empty()) + return 0; + + return microbialConcentration_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +oxygenConcentration(unsigned elemIdx) const +{ + if (oxygenConcentration_.empty()) + return 0; + + return oxygenConcentration_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +ureaConcentration(unsigned elemIdx) const +{ + if (ureaConcentration_.empty()) + return 0; + + return ureaConcentration_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +biofilmConcentration(unsigned elemIdx) const +{ + if (biofilmConcentration_.empty()) + return 0; + + return biofilmConcentration_[elemIdx]; +} + +template +Scalar EclGenericProblem:: +calciteConcentration(unsigned elemIdx) const +{ + if (calciteConcentration_.empty()) + return 0; + + return calciteConcentration_[elemIdx]; +} + template unsigned EclGenericProblem:: pvtRegionIndex(unsigned elemIdx) const diff --git a/ebos/eclgenericproblem.hh b/ebos/eclgenericproblem.hh index 8ca65ec9d..123aa7400 100644 --- a/ebos/eclgenericproblem.hh +++ b/ebos/eclgenericproblem.hh @@ -145,6 +145,31 @@ public: // TODO: remove this function if not called Scalar polymerMolecularWeight(const unsigned elemIdx) const; + /*! + * \brief Returns the initial microbial concentration for a given a cell index + */ + Scalar microbialConcentration(unsigned elemIdx) const; + + /*! + * \brief Returns the initial oxygen concentration for a given a cell index + */ + Scalar oxygenConcentration(unsigned elemIdx) const; + + /*! + * \brief Returns the initial urea concentration for a given a cell index + */ + Scalar ureaConcentration(unsigned elemIdx) const; + + /*! + * \brief Returns the initial biofilm concentration for a given a cell index + */ + Scalar biofilmConcentration(unsigned elemIdx) const; + + /*! + * \brief Returns the initial calcite concentration for a given a cell index + */ + Scalar calciteConcentration(unsigned elemIdx) const; + /*! * \brief Returns the index the relevant PVT region given a cell index */ @@ -231,7 +256,8 @@ protected: int numPhases, bool indicesGasEnabled, bool indicesOilEnabled, - bool indicesWaterEnabled) const; + bool indicesWaterEnabled, + bool enableMICP) const; void readRockParameters_(const std::vector& cellCenterDepths); @@ -240,7 +266,8 @@ protected: void readBlackoilExtentionsInitialConditions_(size_t numDof, bool enableSolvent, bool enablePolymer, - bool enablePolymerMolarWeight); + bool enablePolymerMolarWeight, + bool enableMICP); void updatePvtnum_(); void updateSatnum_(); @@ -274,6 +301,11 @@ protected: std::vector polymerConcentration_; std::vector polymerMoleWeight_; // polymer molecular weight std::vector solventSaturation_; + std::vector microbialConcentration_; + std::vector oxygenConcentration_; + std::vector ureaConcentration_; + std::vector biofilmConcentration_; + std::vector calciteConcentration_; std::vector lastRv_; std::vector maxDRv_; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 4c5910967..cfedd1bc1 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -135,7 +135,8 @@ public: getPropValue(), getPropValue(), getPropValue(), - getPropValue()) + getPropValue(), + getPropValue()) , simulator_(simulator) { const SummaryConfig summaryConfig = simulator_.vanguard().summaryConfig(); @@ -359,6 +360,26 @@ public: this->mFracCo2_[globalDofIdx] = stdVolCo2*rhoCO2/stdMassTotal; } + if (!this->cMicrobes_.empty()) { + this->cMicrobes_[globalDofIdx] = intQuants.microbialConcentration().value(); + } + + if (!this->cOxygen_.empty()) { + this->cOxygen_[globalDofIdx] = intQuants.oxygenConcentration().value(); + } + + if (!this->cUrea_.empty()) { + this->cUrea_[globalDofIdx] = 10 * intQuants.ureaConcentration().value(); //Reescaling back the urea concentration (see WellInterface_impl.hpp) + } + + if (!this->cBiofilm_.empty()) { + this->cBiofilm_[globalDofIdx] = intQuants.biofilmConcentration().value(); + } + + if (!this->cCalcite_.empty()) { + this->cCalcite_[globalDofIdx] = intQuants.calciteConcentration().value(); + } + if (!this->bubblePointPressure_.empty()) { try { this->bubblePointPressure_[globalDofIdx] = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex())); @@ -666,7 +687,7 @@ private: this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) * pv; this->pressureTimesHydrocarbonVolume_[globalDofIdx] = this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon; } else if (FluidSystem::phaseIsActive(waterPhaseIdx)) { - this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)) * pv; + this->pressureTimesPoreVolume_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx)) * pv; } } @@ -680,7 +701,7 @@ private: continue; const double b = getValue(fs.invB(phaseIdx)); - const double s = getValue(fs.saturation(phaseIdx)); + const double s = getValue(fs.saturation(phaseIdx)); fipr[phaseIdx] = s * pv; fip[phaseIdx] = b * fipr[phaseIdx]; } @@ -698,8 +719,8 @@ private: this->fip_[Inplace::Phase::GasResVolume][globalDofIdx] = fipr[gasPhaseIdx]; if (FluidSystem::phaseIsActive(waterPhaseIdx) && !this->fip_[Inplace::Phase::WaterResVolume].empty()) this->fip_[Inplace::Phase::WaterResVolume][globalDofIdx] = fipr[waterPhaseIdx]; - - if (FluidSystem::phaseIsActive(waterPhaseIdx) && !this->fip_[Inplace::Phase::SALT].empty()) + + if (FluidSystem::phaseIsActive(waterPhaseIdx) && !this->fip_[Inplace::Phase::SALT].empty()) this->fip_[Inplace::Phase::SALT][globalDofIdx] = fipr[waterPhaseIdx] * fs.saltConcentration().value(); // Store the pure oil and gas Fip diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 2a3792314..cb885ff22 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -539,6 +539,10 @@ template struct EnableExtbo { static constexpr bool value = false; }; +template +struct EnableMICP { + static constexpr bool value = false; +}; // disable thermal flux boundaries by default template @@ -627,6 +631,7 @@ class EclProblem : public GetPropType enum { enableDiffusion = getPropValue() }; enum { enableThermalFluxBoundaries = getPropValue() }; enum { enableApiTracking = getPropValue() }; + enum { enableMICP = getPropValue() }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; @@ -658,6 +663,7 @@ class EclProblem : public GetPropType using FoamModule = BlackOilFoamModule; using BrineModule = BlackOilBrineModule; using ExtboModule = BlackOilExtboModule; + using MICPModule= BlackOilMICPModule; using InitialFluidState = typename EclEquilInitializer::ScalarFluidState; @@ -789,6 +795,7 @@ public: FoamModule::initFromState(vanguard.eclState()); BrineModule::initFromState(vanguard.eclState()); ExtboModule::initFromState(vanguard.eclState()); + MICPModule::initFromState(vanguard.eclState()); // create the ECL writer eclWriter_.reset(new EclWriterType(simulator)); @@ -908,7 +915,8 @@ public: Indices::numPhases, Indices::gasEnabled, Indices::oilEnabled, - Indices::waterEnabled); + Indices::waterEnabled, + enableMICP); } catch(const std::exception& e) { @@ -1161,6 +1169,16 @@ public: schedule, simulator.vanguard().actionState(), simulator.vanguard().summaryState()); + + // deal with "clogging" for the MICP model + if constexpr (enableMICP){ + auto& model = this->model(); + const auto& residual = this->model().linearizer().residual(); + for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) { + auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx]; + MICPModule::checkCloggingMICP(model, phi, globalDofIdx); + } + } } /*! @@ -1805,6 +1823,14 @@ public: if constexpr (enableBrine) values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration(); + if constexpr (enableMICP){ + values[Indices::microbialConcentrationIdx]= this->microbialConcentration_[globalDofIdx]; + values[Indices::oxygenConcentrationIdx]= this->oxygenConcentration_[globalDofIdx]; + values[Indices::ureaConcentrationIdx]= this->ureaConcentration_[globalDofIdx]; + values[Indices::calciteConcentrationIdx]= this->calciteConcentration_[globalDofIdx]; + values[Indices::biofilmConcentrationIdx]= this->biofilmConcentration_[globalDofIdx]; + } + values.checkDefined(); } @@ -2297,11 +2323,12 @@ private: else readExplicitInitialCondition_(); - if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight) + if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP) this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(), enableSolvent, enablePolymer, - enablePolymerMolarWeight); + enablePolymerMolarWeight, + enableMICP); //initialize min/max values size_t numElems = this->model().numGridDof(); @@ -2372,6 +2399,14 @@ private: this->polymerMoleWeight_.resize(numElems, 0.0); } + if constexpr (enableMICP){ + this->microbialConcentration_.resize(numElems, 0.0); + this->oxygenConcentration_.resize(numElems, 0.0); + this->ureaConcentration_.resize(numElems, 0.0); + this->biofilmConcentration_.resize(numElems, 0.0); + this->calciteConcentration_.resize(numElems, 0.0); + } + for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { auto& elemFluidState = initialFluidStates_[elemIdx]; elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx)); @@ -2405,6 +2440,13 @@ private: if constexpr (enablePolymer) this->polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx); + if constexpr (enableMICP){ + this->microbialConcentration_[elemIdx] = eclWriter_->eclOutputModule().getMicrobialConcentration(elemIdx); + this->oxygenConcentration_[elemIdx] = eclWriter_->eclOutputModule().getOxygenConcentration(elemIdx); + this->ureaConcentration_[elemIdx] = eclWriter_->eclOutputModule().getUreaConcentration(elemIdx); + this->biofilmConcentration_[elemIdx] = eclWriter_->eclOutputModule().getBiofilmConcentration(elemIdx); + this->calciteConcentration_[elemIdx] = eclWriter_->eclOutputModule().getCalciteConcentration(elemIdx); + } // if we need to restart for polymer molecular weight simulation, we need to add related here } diff --git a/flow/flow_ebos_gasoil.cpp b/flow/flow_ebos_gasoil.cpp index 8d1ce8773..c9a5b4364 100644 --- a/flow/flow_ebos_gasoil.cpp +++ b/flow/flow_ebos_gasoil.cpp @@ -51,14 +51,15 @@ private: using FluidSystem = GetPropType; public: - typedef BlackOilTwoPhaseIndices(), - getPropValue(), - getPropValue(), - getPropValue(), - getPropValue(), - getPropValue(), - /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::waterCompIdx> type; + typedef BlackOilTwoPhaseIndices(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*disabledCompIdx=*/FluidSystem::waterCompIdx, + getPropValue()> type; }; }} diff --git a/flow/flow_ebos_gaswater.cpp b/flow/flow_ebos_gaswater.cpp index 237cb8a36..d7c314e62 100644 --- a/flow/flow_ebos_gaswater.cpp +++ b/flow/flow_ebos_gaswater.cpp @@ -61,7 +61,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::oilCompIdx> type; + /*disabledCompIdx=*/FluidSystem::oilCompIdx, + getPropValue()> type; }; }} diff --git a/flow/flow_ebos_micp.cpp b/flow/flow_ebos_micp.cpp new file mode 100644 index 000000000..797f90bd3 --- /dev/null +++ b/flow/flow_ebos_micp.cpp @@ -0,0 +1,103 @@ +/* + 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 3 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 . +*/ +#include "config.h" + +#include + +#include +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlowMICPProblem { + using InheritsFrom = std::tuple; +}; +} +template +struct EnableMICP { + static constexpr bool value = true; +}; +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef BlackOilOnePhaseIndices(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*enabledCompIdx=*/FluidSystem::waterCompIdx, + 5> type; //Five MICP components +}; +}} + +namespace Opm { +void flowEbosMICPSetDeck(double setupTime, std::shared_ptr deck, + std::shared_ptr eclState, + std::shared_ptr schedule, + std::shared_ptr summaryConfig) +{ + using TypeTag = Properties::TTag::EclFlowMICPProblem; + using Vanguard = GetPropType; + + Vanguard::setExternalSetupTime(setupTime); + Vanguard::setExternalDeck(std::move(deck)); + Vanguard::setExternalEclState(std::move(eclState)); + Vanguard::setExternalSchedule(std::move(schedule)); + Vanguard::setExternalSummaryConfig(std::move(summaryConfig)); +} + +// ----------------- Main program ----------------- +int flowEbosMICPMain(int argc, char** argv, bool outputCout, bool outputFiles) +{ + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + resetLocale(); + +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + + FlowMainEbos + mainfunc {argc, argv, outputCout, outputFiles}; + return mainfunc.execute(); +} + +} diff --git a/flow/flow_ebos_micp.hpp b/flow/flow_ebos_micp.hpp new file mode 100644 index 000000000..2921b7c19 --- /dev/null +++ b/flow/flow_ebos_micp.hpp @@ -0,0 +1,37 @@ +/* + 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 3 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 . +*/ +#ifndef FLOW_EBOS_MICP_HPP +#define FLOW_EBOS_MICP_HPP + +#include + +namespace Opm { + +class Deck; +class EclipseState; +class Schedule; +class SummaryConfig; + +void flowEbosMICPSetDeck(double setupTime, std::shared_ptr deck, + std::shared_ptr eclState, + std::shared_ptr schedule, + std::shared_ptr summaryConfig); +int flowEbosMICPMain(int argc, char** argv, bool outputCout, bool outputFiles); + +} + +#endif // FLOW_EBOS_MICP_HPP diff --git a/flow/flow_ebos_oilwater.cpp b/flow/flow_ebos_oilwater.cpp index 0c73f9443..1fe28d262 100644 --- a/flow/flow_ebos_oilwater.cpp +++ b/flow/flow_ebos_oilwater.cpp @@ -58,7 +58,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::gasCompIdx> type; + /*disabledCompIdx=*/FluidSystem::gasCompIdx, + getPropValue()> type; }; }} diff --git a/flow/flow_ebos_oilwater_brine.cpp b/flow/flow_ebos_oilwater_brine.cpp index ec56dff26..91da841bb 100644 --- a/flow/flow_ebos_oilwater_brine.cpp +++ b/flow/flow_ebos_oilwater_brine.cpp @@ -61,7 +61,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::gasCompIdx> type; + /*disabledCompIdx=*/FluidSystem::gasCompIdx, + getPropValue()> type; }; }} diff --git a/flow/flow_ebos_oilwater_polymer.cpp b/flow/flow_ebos_oilwater_polymer.cpp index 5bedc4ef3..586a680fd 100644 --- a/flow/flow_ebos_oilwater_polymer.cpp +++ b/flow/flow_ebos_oilwater_polymer.cpp @@ -61,7 +61,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::gasCompIdx> type; + /*disabledCompIdx=*/FluidSystem::gasCompIdx, + getPropValue()> type; }; }} diff --git a/flow/flow_ebos_oilwater_polymer_injectivity.cpp b/flow/flow_ebos_oilwater_polymer_injectivity.cpp index e0639eaa6..7484ea6ad 100644 --- a/flow/flow_ebos_oilwater_polymer_injectivity.cpp +++ b/flow/flow_ebos_oilwater_polymer_injectivity.cpp @@ -67,7 +67,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*disabledCompIdx=*/FluidSystem::gasCompIdx> type; + /*disabledCompIdx=*/FluidSystem::gasCompIdx, + getPropValue()> type; }; }} diff --git a/flow/flow_ebos_onephase.cpp b/flow/flow_ebos_onephase.cpp index 4176b92f3..1da67791a 100644 --- a/flow/flow_ebos_onephase.cpp +++ b/flow/flow_ebos_onephase.cpp @@ -42,6 +42,7 @@ private: using BaseTypeTag = TTag::EclFlowProblem; using FluidSystem = GetPropType; +public: public: using type = BlackOilOnePhaseIndices(), getPropValue(), @@ -50,7 +51,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*enabledCompIdx=*/FluidSystem::waterCompIdx>; + /*enabledCompIdx=*/FluidSystem::waterCompIdx, + getPropValue()>; }; } // namespace Opm::Properties diff --git a/flow/flow_ebos_onephase_energy.cpp b/flow/flow_ebos_onephase_energy.cpp index e674810c9..505aa9f27 100644 --- a/flow/flow_ebos_onephase_energy.cpp +++ b/flow/flow_ebos_onephase_energy.cpp @@ -53,7 +53,8 @@ public: getPropValue(), getPropValue(), /*PVOffset=*/0, - /*enebledCompIdx=*/FluidSystem::waterCompIdx> + /*enebledCompIdx=*/FluidSystem::waterCompIdx, + getPropValue()> type; }; diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index c1697f773..b315db07d 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -126,6 +126,10 @@ template struct EnableBrine { static constexpr bool value = false; }; +template +struct EnableMICP { + static constexpr bool value = false; +}; template struct EclWellModel { @@ -172,6 +176,11 @@ namespace Opm { static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx; static const int contiFoamEqIdx = Indices::contiFoamEqIdx; static const int contiBrineEqIdx = Indices::contiBrineEqIdx; + static const int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx; + static const int contiOxygenEqIdx = Indices::contiOxygenEqIdx; + static const int contiUreaEqIdx = Indices::contiUreaEqIdx; + static const int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx; + static const int contiCalciteEqIdx = Indices::contiCalciteEqIdx; static const int solventSaturationIdx = Indices::solventSaturationIdx; static const int zFractionIdx = Indices::zFractionIdx; static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx; @@ -179,6 +188,11 @@ namespace Opm { static const int temperatureIdx = Indices::temperatureIdx; static const int foamConcentrationIdx = Indices::foamConcentrationIdx; static const int saltConcentrationIdx = Indices::saltConcentrationIdx; + static const int microbialConcentrationIdx = Indices::microbialConcentrationIdx; + static const int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx; + static const int ureaConcentrationIdx = Indices::ureaConcentrationIdx; + static const int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx; + static const int calciteConcentrationIdx = Indices::calciteConcentrationIdx; typedef Dune::FieldVector VectorBlockType; typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; @@ -452,13 +466,13 @@ namespace Opm { Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationNew = 1.0; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::numActivePhases() > 1) { saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx]; } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { assert(Indices::compositionSwitchIdx >= 0 ); saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx]; @@ -489,7 +503,7 @@ namespace Opm { } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { assert(Indices::compositionSwitchIdx >= 0 ); @@ -726,7 +740,29 @@ namespace Opm { maxCoeff[ contiEnergyEqIdx ] = std::max( maxCoeff[ contiEnergyEqIdx ], std::abs( R2 ) / pvValue ); } + if constexpr (has_micp_) { + B_avg[ contiMicrobialEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R1 = ebosResid[cell_idx][contiMicrobialEqIdx]; + R_sum[ contiMicrobialEqIdx ] += R1; + maxCoeff[ contiMicrobialEqIdx ] = std::max( maxCoeff[ contiMicrobialEqIdx ], std::abs( R1 ) / pvValue ); + B_avg[ contiOxygenEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R2 = ebosResid[cell_idx][contiOxygenEqIdx]; + R_sum[ contiOxygenEqIdx ] += R2; + maxCoeff[ contiOxygenEqIdx ] = std::max( maxCoeff[ contiOxygenEqIdx ], std::abs( R2 ) / pvValue ); + B_avg[ contiUreaEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R3 = ebosResid[cell_idx][contiUreaEqIdx]; + R_sum[ contiUreaEqIdx ] += R3; + maxCoeff[ contiUreaEqIdx ] = std::max( maxCoeff[ contiUreaEqIdx ], std::abs( R3 ) / pvValue ); + B_avg[ contiBiofilmEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R4 = ebosResid[cell_idx][contiBiofilmEqIdx]; + R_sum[ contiBiofilmEqIdx ] += R4; + maxCoeff[ contiBiofilmEqIdx ] = std::max( maxCoeff[ contiBiofilmEqIdx ], std::abs( R4 ) / pvValue ); + B_avg[ contiCalciteEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); + const auto R5 = ebosResid[cell_idx][contiCalciteEqIdx]; + R_sum[ contiCalciteEqIdx ] += R5; + maxCoeff[ contiCalciteEqIdx ] = std::max( maxCoeff[ contiCalciteEqIdx ], std::abs( R5 ) / pvValue ); } + } OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm()); @@ -802,7 +838,7 @@ namespace Opm { // max_strict_iter_ is 8. Hence only iteration chooses // whether to use relaxed or not. // To activate only fraction use fraction below 1 and iter 0. - const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.max_strict_iter_; + const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.max_strict_iter_; const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_; // Finish computation @@ -849,6 +885,13 @@ namespace Opm { if constexpr (has_brine_) { compNames[saltConcentrationIdx] = "Brine"; } + if constexpr (has_micp_) { + compNames[microbialConcentrationIdx] = "Microbes"; + compNames[oxygenConcentrationIdx] = "Oxygen"; + compNames[ureaConcentrationIdx] = "Urea"; + compNames[biofilmConcentrationIdx] = "Biofilm"; + compNames[calciteConcentrationIdx] = "Calcite"; + } } // Create convergence report. @@ -994,6 +1037,7 @@ namespace Opm { static constexpr bool has_energy_ = getPropValue(); static constexpr bool has_foam_ = getPropValue(); static constexpr bool has_brine_ = getPropValue(); + static constexpr bool has_micp_ = getPropValue(); ModelParameters param_; SimulatorReportSingle failureReport_; diff --git a/opm/simulators/flow/Main.hpp b/opm/simulators/flow/Main.hpp index fa6232746..5f1ade84b 100644 --- a/opm/simulators/flow/Main.hpp +++ b/opm/simulators/flow/Main.hpp @@ -37,6 +37,7 @@ # include # include # include +# include # endif #include @@ -288,6 +289,19 @@ namespace Opm if (false) {} #ifndef FLOW_BLACKOIL_ONLY + // Single-phase case + else if( eclipseState_->runspec().micp() ) { + // micp + if ( !phases.active( Phase::WATER) || phases.size() > 2) { + if (outputCout_) + std::cerr << "No valid configuration is found for MICP simulation, the only valid option is " + << "water + MICP" << std::endl; + return EXIT_FAILURE; + } + flowEbosMICPSetDeck( + setupTime_, deck_, eclipseState_, schedule_, summaryConfig_); + return flowEbosMICPMain(argc_, argv_, outputCout_, outputFiles_); + } // Twophase cases else if (phases.size() == 2) { return this->runTwoPhase(phases); diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp index 64d492685..9de0a127f 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp @@ -83,6 +83,7 @@ public: typedef AdaptiveTimeSteppingEbos TimeStepper; typedef BlackOilPolymerModule PolymerModule; + typedef BlackOilMICPModule MICPModule; typedef BlackoilModelEbos Model; typedef NonlinearSolverEbos Solver; diff --git a/opm/simulators/linalg/FlexibleSolver5.cpp b/opm/simulators/linalg/FlexibleSolver5.cpp new file mode 100644 index 000000000..c8b80a95a --- /dev/null +++ b/opm/simulators/linalg/FlexibleSolver5.cpp @@ -0,0 +1,24 @@ +/* + Copyright 2019, 2020 SINTEF Digital, Mathematics and Cybernetics. + + 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 3 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 . +*/ + +#include "config.h" + +#include + +INSTANTIATE_FLEXIBLESOLVER(5); diff --git a/opm/simulators/linalg/FlexibleSolver6.cpp b/opm/simulators/linalg/FlexibleSolver6.cpp new file mode 100644 index 000000000..cc2242baa --- /dev/null +++ b/opm/simulators/linalg/FlexibleSolver6.cpp @@ -0,0 +1,24 @@ +/* + Copyright 2019, 2020 SINTEF Digital, Mathematics and Cybernetics. + + 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 3 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 . +*/ + +#include "config.h" + +#include + +INSTANTIATE_FLEXIBLESOLVER(6); diff --git a/opm/simulators/linalg/bda/BILU0.cpp b/opm/simulators/linalg/bda/BILU0.cpp index 7b06f9c8d..ed1b0466e 100644 --- a/opm/simulators/linalg/bda/BILU0.cpp +++ b/opm/simulators/linalg/bda/BILU0.cpp @@ -687,9 +687,9 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS } // end namespace bda - - diff --git a/opm/simulators/linalg/bda/BdaBridge.cpp b/opm/simulators/linalg/bda/BdaBridge.cpp index d88cae9b0..e889fb950 100644 --- a/opm/simulators/linalg/bda/BdaBridge.cpp +++ b/opm/simulators/linalg/bda/BdaBridge.cpp @@ -292,9 +292,9 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS } // namespace Opm - - diff --git a/opm/simulators/linalg/bda/BlockedMatrix.cpp b/opm/simulators/linalg/bda/BlockedMatrix.cpp index 811d1b484..d9119cfd3 100644 --- a/opm/simulators/linalg/bda/BlockedMatrix.cpp +++ b/opm/simulators/linalg/bda/BlockedMatrix.cpp @@ -478,6 +478,8 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS diff --git a/opm/simulators/linalg/bda/FPGABILU0.cpp b/opm/simulators/linalg/bda/FPGABILU0.cpp index 6246140e4..cb7755018 100644 --- a/opm/simulators/linalg/bda/FPGABILU0.cpp +++ b/opm/simulators/linalg/bda/FPGABILU0.cpp @@ -407,6 +407,8 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS diff --git a/opm/simulators/linalg/bda/FPGASolverBackend.cpp b/opm/simulators/linalg/bda/FPGASolverBackend.cpp index ccdfb9092..de76e2101 100644 --- a/opm/simulators/linalg/bda/FPGASolverBackend.cpp +++ b/opm/simulators/linalg/bda/FPGASolverBackend.cpp @@ -726,6 +726,8 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS diff --git a/opm/simulators/linalg/bda/Reorder.cpp b/opm/simulators/linalg/bda/Reorder.cpp index 1834540c8..78c478290 100644 --- a/opm/simulators/linalg/bda/Reorder.cpp +++ b/opm/simulators/linalg/bda/Reorder.cpp @@ -376,6 +376,8 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS diff --git a/opm/simulators/linalg/bda/amgclSolverBackend.cpp b/opm/simulators/linalg/bda/amgclSolverBackend.cpp index 8914dcb25..f2ac39900 100644 --- a/opm/simulators/linalg/bda/amgclSolverBackend.cpp +++ b/opm/simulators/linalg/bda/amgclSolverBackend.cpp @@ -376,8 +376,9 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS } // namespace bda - diff --git a/opm/simulators/linalg/bda/amgclSolverBackend.cu b/opm/simulators/linalg/bda/amgclSolverBackend.cu index a19901c01..12c3414b1 100644 --- a/opm/simulators/linalg/bda/amgclSolverBackend.cu +++ b/opm/simulators/linalg/bda/amgclSolverBackend.cu @@ -82,6 +82,8 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS diff --git a/opm/simulators/linalg/bda/cusparseSolverBackend.cu b/opm/simulators/linalg/bda/cusparseSolverBackend.cu index 3a29cd573..67da1203e 100644 --- a/opm/simulators/linalg/bda/cusparseSolverBackend.cu +++ b/opm/simulators/linalg/bda/cusparseSolverBackend.cu @@ -506,6 +506,8 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS diff --git a/opm/simulators/linalg/bda/openclSolverBackend.cpp b/opm/simulators/linalg/bda/openclSolverBackend.cpp index edbafe56b..e2228baf6 100644 --- a/opm/simulators/linalg/bda/openclSolverBackend.cpp +++ b/opm/simulators/linalg/bda/openclSolverBackend.cpp @@ -173,7 +173,7 @@ openclSolverBackend::openclSolverBackend(int verbosity_, int maxit_, context = std::make_shared(devices[0]); queue.reset(new cl::CommandQueue(*context, devices[0], 0, &err)); - + } catch (const cl::Error& error) { std::ostringstream oss; oss << "OpenCL Error: " << error.what() << "(" << error.err() << ")\n"; @@ -815,8 +815,9 @@ INSTANTIATE_BDA_FUNCTIONS(1); INSTANTIATE_BDA_FUNCTIONS(2); INSTANTIATE_BDA_FUNCTIONS(3); INSTANTIATE_BDA_FUNCTIONS(4); +INSTANTIATE_BDA_FUNCTIONS(5); +INSTANTIATE_BDA_FUNCTIONS(6); #undef INSTANTIATE_BDA_FUNCTIONS } // namespace bda - diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index c1dc044f1..727cf3774 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -119,6 +119,7 @@ namespace Opm { static constexpr bool has_solvent_ = getPropValue(); static constexpr bool has_polymer_ = getPropValue(); static constexpr bool has_energy_ = getPropValue(); + static constexpr bool has_micp_ = getPropValue(); // TODO: where we should put these types, WellInterface or Well Model? // or there is some other strategy, like TypeTag @@ -128,6 +129,7 @@ namespace Opm { typedef Dune::FieldMatrix MatrixBlockType; typedef BlackOilPolymerModule PolymerModule; + typedef BlackOilMICPModule MICPModule; // For the conversion between the surface volume rate and resrevoir voidage rate using RateConverterType = RateConverter:: diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 4ff322d8b..afa00fb42 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1852,26 +1852,27 @@ addWellContribution(WellContributions& wellContribs) const template class MultisegmentWellEval,__VA_ARGS__,double>; // One phase -INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) // Two phase -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) // Blackoil -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,1u,0u>) } // namespace Opm diff --git a/opm/simulators/wells/PerfData.cpp b/opm/simulators/wells/PerfData.cpp index 45c9b0b03..aee3a0dd7 100644 --- a/opm/simulators/wells/PerfData.cpp +++ b/opm/simulators/wells/PerfData.cpp @@ -37,6 +37,7 @@ PerfData::PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases) , connection_transmissibility_factor(num_perf) , satnum_id(num_perf) , ecl_index(num_perf) + , micp_rates(num_perf) { if (injector) { this->water_throughput.resize(num_perf); @@ -70,8 +71,8 @@ bool PerfData::try_assign(const PerfData& other) { this->skin_pressure = other.skin_pressure; this->water_velocity = other.water_velocity; this->prod_index = other.prod_index; + this->micp_rates = other.micp_rates; return true; } } - diff --git a/opm/simulators/wells/PerfData.hpp b/opm/simulators/wells/PerfData.hpp index 70497b3d9..bd8920dd9 100644 --- a/opm/simulators/wells/PerfData.hpp +++ b/opm/simulators/wells/PerfData.hpp @@ -45,6 +45,7 @@ public: std::vector polymer_rates; std::vector brine_rates; std::vector prod_index; + std::vector micp_rates; std::vector cell_index; std::vector connection_transmissibility_factor; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index fbf414c91..35935ebaa 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include @@ -95,6 +96,7 @@ namespace Opm using Base::has_foam; using Base::has_brine; using Base::has_energy; + using Base::has_micp; using PolymerModule = BlackOilPolymerModule; using FoamModule = BlackOilFoamModule; diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 7f79043a6..f91993d2f 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -88,13 +88,13 @@ relaxationFactorFractionsProducer(const std::vector& primary_variables, double relaxation_factor = 1.0; if (FluidSystem::numActivePhases() > 1) { - if constexpr (has_wfrac_variable) { + if constexpr (has_wfrac_variable) { const double relaxation_factor_w = StandardWellGeneric:: relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); } - if constexpr (has_gfrac_variable) { + if constexpr (has_gfrac_variable) { const double relaxation_factor_g = StandardWellGeneric:: relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); @@ -329,11 +329,11 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log // this will happen. } else if (baseif_.isProducer()) { // producers // TODO: the following are not addressed for the solvent case yet - if constexpr (has_wfrac_variable) { + if constexpr (has_wfrac_variable) { primary_variables_[WFrac] = 1.0 / np; } - - if constexpr (has_gfrac_variable) { + + if constexpr (has_gfrac_variable) { primary_variables_[GFrac] = 1.0 / np; } } else { @@ -534,12 +534,12 @@ processFractions() const F[pu.phase_pos[Oil]] = 0.0; } } - - if constexpr (has_wfrac_variable) { + + if constexpr (has_wfrac_variable) { primary_variables_[WFrac] = F[pu.phase_pos[Water]]; } - - if constexpr (has_gfrac_variable) { + + if constexpr (has_gfrac_variable) { primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; } if constexpr (Indices::enableSolvent) { @@ -732,7 +732,7 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells, : 1.0; // update the second and third well variable (The flux fractions) - + if constexpr (has_wfrac_variable) { const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit); @@ -740,7 +740,7 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells, primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; } - if constexpr (has_gfrac_variable) { + if constexpr (has_gfrac_variable) { const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit); primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; @@ -1110,26 +1110,27 @@ addWellContribution(WellContributions& wellContribs) const template class StandardWellEval,__VA_ARGS__,double>; // One phase -INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) // Two phase -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) // Blackoil -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,1u,0u>) } diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 02a85d8c9..cde3c5e8a 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -694,6 +694,31 @@ namespace Opm connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm); } + if constexpr (has_micp) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + EvalWell cq_s_microbe = cq_s[waterCompIdx]; + if (this->isInjector()) { + cq_s_microbe *= this->wmicrobes(); + } else { + cq_s_microbe *= this->extendEval(intQuants.microbialConcentration()); + } + connectionRates[perf][Indices::contiMicrobialEqIdx] = Base::restrictEval(cq_s_microbe); + EvalWell cq_s_oxygen = cq_s[waterCompIdx]; + if (this->isInjector()) { + cq_s_oxygen *= this->woxygen(); + } else { + cq_s_oxygen *= this->extendEval(intQuants.oxygenConcentration()); + } + connectionRates[perf][Indices::contiOxygenEqIdx] = Base::restrictEval(cq_s_oxygen); + EvalWell cq_s_urea = cq_s[waterCompIdx]; + if (this->isInjector()) { + cq_s_urea *= this->wurea(); + } else { + cq_s_urea *= this->extendEval(intQuants.ureaConcentration()); + } + connectionRates[perf][Indices::contiUreaEqIdx] = Base::restrictEval(cq_s_urea); + } + // Store the perforation pressure for later usage. perf_data.pressure[perf] = ws.bhp + this->perf_pressure_diffs_[perf]; } @@ -1323,7 +1348,7 @@ namespace Opm { // the following implementation assume that the polymer is always after the w-o-g phases // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells - assert((int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction); + assert((int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp); std::vector res; ConvergenceReport report = this->StdWellEval::getWellConvergence(well_state, diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 28a7042b6..f1fa4110c 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -112,8 +112,9 @@ public: static constexpr bool has_polymermw = getPropValue(); static constexpr bool has_foam = getPropValue(); static constexpr bool has_brine = getPropValue(); + static constexpr bool has_micp = getPropValue(); - // For the conversion between the surface volume rate and reservoir voidage rate + // For the conversion between the surface volume rate and reservoir voidage rate using FluidState = BlackOilFluidState, \ double>; // One phase -INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) // Two phase -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u,0u>) // Blackoil -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,2u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,1u>) -INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,1u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,2u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,1u,0u>) +INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,1u,0u>) } // namespace Opm diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 1e41bb778..1e26d2263 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -150,6 +150,74 @@ namespace Opm return 0.0; } + template + double + WellInterface:: + wmicrobes() const + { + if constexpr (has_micp) { + auto injectorType = this->well_ecl_.injectorType(); + + if (injectorType == InjectorType::WATER) { + WellMICPProperties microbes = this->well_ecl_.getMICPProperties(); + const double microbial_injection_concentration = microbes.m_microbialConcentration; + return microbial_injection_concentration; + } else { + // Not a water injection well => no microbes. + return 0.0; + } + } + + return 0.0; + } + + template + double + WellInterface:: + woxygen() const + { + if constexpr (has_micp) { + auto injectorType = this->well_ecl_.injectorType(); + + if (injectorType == InjectorType::WATER) { + WellMICPProperties oxygen = this->well_ecl_.getMICPProperties(); + const double oxygen_injection_concentration = oxygen.m_oxygenConcentration; + return oxygen_injection_concentration; + } else { + // Not a water injection well => no oxygen. + return 0.0; + } + } + + return 0.0; + } + + // The urea injection concentration is scaled down by a factor of 10, since its value + // can be much bigger than 1 (not doing this slows the simulations). The + // corresponding values are scaled accordingly in blackoilmicpmodules.hh when computing + // the reactions and also when writing the output files (vtk and eclipse format, i.e., + // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively). + + template + double + WellInterface:: + wurea() const + { + if constexpr (has_micp) { + auto injectorType = this->well_ecl_.injectorType(); + + if (injectorType == InjectorType::WATER) { + WellMICPProperties urea = this->well_ecl_.getMICPProperties(); + const double urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10 + return urea_injection_concentration; + } else { + // Not a water injection well => no urea. + return 0.0; + } + } + + return 0.0; + } template bool