diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4588c2dd9a8..f557480e6c0 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -574,9 +574,40 @@ list (APPEND TEST_DATA_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES + flow/flow_blackoil.hpp + flow/flow_blackoil_legacyassembly.hpp + flow/flow_brine.hpp + flow/flow_brine_precsalt_vapwat.hpp + flow/flow_brine_saltprecipitation.hpp + flow/flow_energy.hpp + flow/flow_extbo.hpp + flow/flow_foam.hpp + flow/flow_gasoil.hpp + flow/flow_gasoil_energy.hpp + flow/flow_gasoildiffuse.hpp + flow/flow_gaswater.hpp + flow/flow_gaswater_brine.hpp + flow/flow_gaswater_dissolution.hpp + flow/flow_gaswater_dissolution_diffuse.hpp + flow/flow_gaswater_energy.hpp + flow/flow_gaswater_saltprec_energy.hpp + flow/flow_gaswater_saltprec_vapwat.hpp + flow/flow_gaswater_solvent.hpp + flow/flow_micp.hpp + flow/flow_oilwater.hpp + flow/flow_oilwater_brine.hpp + flow/flow_oilwater_polymer.hpp + flow/flow_oilwater_polymer_injectivity.hpp + flow/flow_onephase.hpp + flow/flow_onephase_energy.hpp + flow/flow_polymer.hpp + flow/flow_solvent.hpp + flow/flow_solvent_foam.hpp + opm/models/blackoil/blackoilboundaryratevector.hh opm/models/blackoil/blackoilbrinemodules.hh opm/models/blackoil/blackoilbrineparams.hpp + opm/models/blackoil/blackoilconvectivemixingmodule.hh opm/models/blackoil/blackoildarcyfluxmodule.hh opm/models/blackoil/blackoildiffusionmodule.hh opm/models/blackoil/blackoildispersionmodule.hh @@ -905,6 +936,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/overlaptypes.hh opm/simulators/linalg/OwningBlockPreconditioner.hpp opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp + opm/simulators/linalg/MILU.hpp opm/simulators/linalg/parallelamgbackend.hh opm/simulators/linalg/parallelbasebackend.hh opm/simulators/linalg/parallelbicgstabbackend.hh @@ -912,10 +944,13 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/ParallelIstlInformation.hpp opm/simulators/linalg/ParallelOverlappingILU0.hpp opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp - opm/simulators/linalg/PressureSolverPolicy.hpp - opm/simulators/linalg/PressureTransferPolicy.hpp + opm/simulators/linalg/PreconditionerFactoryGPUIncludeWrapper.hpp opm/simulators/linalg/PreconditionerFactory.hpp + opm/simulators/linalg/PreconditionerFactory_impl.hpp opm/simulators/linalg/PreconditionerWithUpdate.hpp + opm/simulators/linalg/PressureBhpTransferPolicy.hpp + opm/simulators/linalg/PressureSolverPolicy.hpp + opm/simulators/linalg/PressureTransferPolicy.hpp opm/simulators/linalg/PropertyTree.hpp opm/simulators/linalg/residreductioncriterion.hh opm/simulators/linalg/SmallDenseMatrixUtils.hpp @@ -952,6 +987,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/gatherDeferredLogger.hpp opm/simulators/utils/moduleVersion.hpp opm/simulators/utils/phaseUsageFromDeck.hpp + opm/simulators/utils/ParallelCommunication.hpp + opm/simulators/utils/ParallelSerialization.hpp opm/simulators/utils/readDeck.hpp opm/simulators/utils/satfunc/RelpermDiagnostics.hpp opm/simulators/wells/ALQState.hpp @@ -962,6 +999,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/BlackoilWellModelGuideRates.hpp opm/simulators/wells/BlackoilWellModelRestart.hpp opm/simulators/wells/ConnFiltrateData.hpp + opm/simulators/wells/ConnFracStatistics.hpp opm/simulators/wells/FractionCalculator.hpp opm/simulators/wells/GasLiftCommon.hpp opm/simulators/wells/GasLiftGroupInfo.hpp @@ -992,6 +1030,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/RatioCalculator.hpp opm/simulators/wells/RegionAttributeHelpers.hpp opm/simulators/wells/RegionAverageCalculator.hpp + opm/simulators/wells/RunningStatistics.hpp opm/simulators/wells/SingleWellState.hpp opm/simulators/wells/StandardWell.hpp opm/simulators/wells/StandardWell_impl.hpp @@ -1015,14 +1054,16 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/WellGroupControls.hpp opm/simulators/wells/WellGroupHelpers.hpp opm/simulators/wells/WellHelpers.hpp - opm/simulators/wells/WellInterface.hpp + opm/simulators/wells/WellInterfaceFluidSystem.hpp opm/simulators/wells/WellInterfaceGeneric.hpp + opm/simulators/wells/WellInterface.hpp opm/simulators/wells/WellInterface_impl.hpp + opm/simulators/wells/WellInterfaceIndices.hpp opm/simulators/wells/WellProdIndexCalculator.hpp opm/simulators/wells/WellState.hpp opm/simulators/wells/WellTest.hpp opm/simulators/wells/WGState.hpp - ) +) if (USE_GPU_BRIDGE) list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/gpubridge/amgclSolverBackend.hpp diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 2de7bea2893..ea129abb717 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -56,8 +56,6 @@ #include #include -#include - namespace Opm { /*! * \ingroup BlackOilModel diff --git a/opm/models/discretization/common/fvbaseprimaryvariables.hh b/opm/models/discretization/common/fvbaseprimaryvariables.hh index dc900e20da4..33a647edda1 100644 --- a/opm/models/discretization/common/fvbaseprimaryvariables.hh +++ b/opm/models/discretization/common/fvbaseprimaryvariables.hh @@ -163,14 +163,14 @@ namespace Dune { }; - /** Specialization of FieldTraits for all PrimaryVariables derived from Opm::FvBasePrimaryVariables */ - template class EwomsPrimaryVariable> - struct FieldTraits< EwomsPrimaryVariable< TypeTag > > - : public FieldTraitsImpl< TypeTag, - std::is_base_of< Opm::FvBasePrimaryVariables< TypeTag >, - EwomsPrimaryVariable< TypeTag > > :: value > - { - }; + // /** Specialization of FieldTraits for all PrimaryVariables derived from Opm::FvBasePrimaryVariables */ + // template class EwomsPrimaryVariable> + // struct FieldTraits< EwomsPrimaryVariable< TypeTag > > + // : public FieldTraitsImpl< TypeTag, + // std::is_base_of< Opm::FvBasePrimaryVariables< TypeTag >, + // EwomsPrimaryVariable< TypeTag > > :: value > + // { + // }; } #endif diff --git a/opm/simulators/flow/ActionHandler.hpp b/opm/simulators/flow/ActionHandler.hpp index f8825c67368..1dd7a3bacd6 100644 --- a/opm/simulators/flow/ActionHandler.hpp +++ b/opm/simulators/flow/ActionHandler.hpp @@ -99,7 +99,6 @@ class ActionHandler void evalUDQAssignments(const unsigned episodeIdx, UDQState& udq_state); -private: /// Convey dynamic updates triggered by an action block back to the /// running simulator. /// @@ -125,6 +124,7 @@ class ActionHandler const TransFunc& updateTrans, bool& commit_wellstate); +private: /// Static properties such as permeability and transmissibility. EclipseState& ecl_state_; diff --git a/opm/simulators/flow/CpGridVanguard.hpp b/opm/simulators/flow/CpGridVanguard.hpp index 4ae4c747962..86fd0c648cb 100644 --- a/opm/simulators/flow/CpGridVanguard.hpp +++ b/opm/simulators/flow/CpGridVanguard.hpp @@ -227,6 +227,7 @@ class CpGridVanguard : public FlowBaseVanguard } this->doLoadBalance_(this->edgeWeightsMethod(), this->ownersFirst(), + this->addCorners(),this->numOverlap(), this->partitionMethod(), this->serialPartitioning(), this->enableDistributedWells(), this->allow_splitting_inactive_wells_, diff --git a/opm/simulators/flow/FlowBaseVanguard.hpp b/opm/simulators/flow/FlowBaseVanguard.hpp index fc0200d360f..11ddabc16a0 100644 --- a/opm/simulators/flow/FlowBaseVanguard.hpp +++ b/opm/simulators/flow/FlowBaseVanguard.hpp @@ -126,6 +126,8 @@ class FlowBaseVanguard : public BaseVanguard, ownersFirst_ = Parameters::Get(); #if HAVE_MPI + numOverlap_ = Parameters::Get(); + addCorners_ = Parameters::Get(); partitionMethod_ = Dune::PartitionMethod(Parameters::Get()); serialPartitioning_ = Parameters::Get(); imbalanceTol_ = Parameters::Get>(); diff --git a/opm/simulators/flow/FlowGenericVanguard.cpp b/opm/simulators/flow/FlowGenericVanguard.cpp index c9ced125299..da91dfbde4a 100644 --- a/opm/simulators/flow/FlowGenericVanguard.cpp +++ b/opm/simulators/flow/FlowGenericVanguard.cpp @@ -424,6 +424,10 @@ void FlowGenericVanguard::registerParameters_() Parameters::Register ("Order cells owned by rank before ghost/overlap cells."); #if HAVE_MPI + Parameters::Register + ("Add corners to partion."); + Parameters::Register + ("Numbers of layers overlap in parallel partition"); Parameters::Register ("Choose partitioning strategy: 0=simple, 1=Zoltan, 2=METIS."); Parameters::Register diff --git a/opm/simulators/flow/FlowGenericVanguard.hpp b/opm/simulators/flow/FlowGenericVanguard.hpp index 8e4935f80e4..92ec997a046 100644 --- a/opm/simulators/flow/FlowGenericVanguard.hpp +++ b/opm/simulators/flow/FlowGenericVanguard.hpp @@ -73,6 +73,8 @@ struct ActionParsingStrictness { static constexpr auto value = "normal"; }; // 0: simple, 1: Zoltan, 2: METIS, see GridEnums.hpp struct PartitionMethod { static constexpr int value = 1; }; +struct AddCorners { static constexpr bool value = false; }; +struct NumOverlap { static constexpr int value = 1; }; struct SchedRestart{ static constexpr bool value = false; }; struct SerialPartitioning{ static constexpr bool value = false; }; @@ -253,7 +255,14 @@ class FlowGenericVanguard { bool ownersFirst() const { return ownersFirst_; } + #if HAVE_MPI + bool addCorners() const + { return addCorners_; } + + int numOverlap() const + { return numOverlap_; } + /*! * \brief Parameter deciding which partition method to use */ @@ -349,6 +358,8 @@ class FlowGenericVanguard { bool ownersFirst_; #if HAVE_MPI + bool addCorners_; + int numOverlap_; Dune::PartitionMethod partitionMethod_; bool serialPartitioning_; double imbalanceTol_; diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index fc5353d84ad..42160aaf591 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -1676,9 +1676,8 @@ class FlowProblem : public GetPropType return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true); } - +protected: typename Vanguard::TransmissibilityType transmissibilities_; - std::shared_ptr materialLawManager_; std::shared_ptr thermalLawManager_; @@ -1694,6 +1693,8 @@ class FlowProblem : public GetPropType PffGridVector pffDofData_; TracerModel tracerModel_; + +private: template struct BCData { @@ -1726,8 +1727,9 @@ class FlowProblem : public GetPropType virtual void handleSolventBC(const BCProp::BCFace&, RateVector&) const = 0; virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0; - +protected: BCData bcindex_; +private: bool nonTrivialBoundaryConditions_ = false; bool explicitRockCompaction_ = false; bool first_step_ = true; diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp index 3540b1da090..8ef16d3c1e1 100644 --- a/opm/simulators/flow/FlowProblemBlackoil.hpp +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -80,8 +80,9 @@ template class FlowProblemBlackoil : public FlowProblem { // TODO: the naming of the Types might be able to be adjusted + public: using FlowProblemType = FlowProblem; - + private: using typename FlowProblemType::Scalar; using typename FlowProblemType::Simulator; using typename FlowProblemType::GridView; @@ -436,7 +437,10 @@ class FlowProblemBlackoil : public FlowProblem void endTimeStep() override { FlowProblemType::endTimeStep(); - + this->endStepApplyAction(); + } + void endStepApplyAction() + { // after the solution is updated, the values in output module needs also updated this->eclWriter()->mutableOutputModule().invalidateLocalData(); diff --git a/opm/simulators/flow/GenericCpGridVanguard.cpp b/opm/simulators/flow/GenericCpGridVanguard.cpp index 6b2145512db..ce795e4eb6b 100644 --- a/opm/simulators/flow/GenericCpGridVanguard.cpp +++ b/opm/simulators/flow/GenericCpGridVanguard.cpp @@ -144,6 +144,8 @@ template void GenericCpGridVanguard:: doLoadBalance_(const Dune::EdgeWeightMethod edgeWeightsMethod, const bool ownersFirst, + const bool addCorners, + const int numOverlap, const Dune::PartitionMethod partitionMethod, const bool serialPartitioning, const bool enableDistributedWells, @@ -205,7 +207,9 @@ doLoadBalance_(const Dune::EdgeWeightMethod edgeWeightsMethod, const auto& possibleFutureConnections = schedule.getPossibleFutureConnections(); // Distribute the grid and switch to the distributed view. if (mpiSize > 1) { - this->distributeGrid(edgeWeightsMethod, ownersFirst, partitionMethod, + this->distributeGrid(edgeWeightsMethod, ownersFirst, + addCorners, numOverlap, + partitionMethod, serialPartitioning, enableDistributedWells, imbalanceTol, loadBalancerSet != 0, faceTrans, wells, @@ -333,6 +337,8 @@ void GenericCpGridVanguard:: distributeGrid(const Dune::EdgeWeightMethod edgeWeightsMethod, const bool ownersFirst, + const bool addCorners, + const int numOverlap, const Dune::PartitionMethod partitionMethod, const bool serialPartitioning, const bool enableDistributedWells, @@ -347,7 +353,7 @@ distributeGrid(const Dune::EdgeWeightMethod edgeWeights if (auto* eclState = dynamic_cast(&eclState1); eclState != nullptr) { - this->distributeGrid(edgeWeightsMethod, ownersFirst, partitionMethod, + this->distributeGrid(edgeWeightsMethod, ownersFirst, addCorners, numOverlap, partitionMethod, serialPartitioning, enableDistributedWells, imbalanceTol, loadBalancerSet, faceTrans, wells, possibleFutureConnections, eclState, parallelWells); @@ -371,6 +377,8 @@ void GenericCpGridVanguard:: distributeGrid(const Dune::EdgeWeightMethod edgeWeightsMethod, const bool ownersFirst, + const bool addCorners, + const int numOverlap, const Dune::PartitionMethod partitionMethod, const bool serialPartitioning, const bool enableDistributedWells, @@ -389,8 +397,8 @@ distributeGrid(const Dune::EdgeWeightMethod edgeWeights *this->grid_, *eclState }; - const auto addCornerCells = false; - const auto overlapLayers = 1; + const auto addCornerCells = addCorners; + const auto overlapLayers = numOverlap; if (loadBalancerSet) { auto parts = isIORank diff --git a/opm/simulators/flow/GenericCpGridVanguard.hpp b/opm/simulators/flow/GenericCpGridVanguard.hpp index 4fd16aea9c6..9c6840f50c6 100644 --- a/opm/simulators/flow/GenericCpGridVanguard.hpp +++ b/opm/simulators/flow/GenericCpGridVanguard.hpp @@ -158,6 +158,8 @@ class GenericCpGridVanguard { #if HAVE_MPI void doLoadBalance_(const Dune::EdgeWeightMethod edgeWeightsMethod, const bool ownersFirst, + const bool addCorners, + const int numOverlap, const Dune::PartitionMethod partitionMethod, const bool serialPartitioning, const bool enableDistributedWells, @@ -176,6 +178,8 @@ class GenericCpGridVanguard { void distributeGrid(const Dune::EdgeWeightMethod edgeWeightsMethod, const bool ownersFirst, + const bool addCorners, + const int numOverlap, const Dune::PartitionMethod partitionMethod, const bool serialPartitioning, const bool enableDistributedWells, @@ -189,6 +193,8 @@ class GenericCpGridVanguard { void distributeGrid(const Dune::EdgeWeightMethod edgeWeightsMethod, const bool ownersFirst, + const bool addCorners, + const int numOverlap, const Dune::PartitionMethod partitionMethod, const bool serialPartitioning, const bool enableDistributedWells, diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 0a345e8ae52..70ea4998751 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -632,12 +632,24 @@ assignToSolution(data::Solution& sol) DataEntry{"STRAINXY", UnitSystem::measure::identity, strainXY_}, DataEntry{"STRAINXZ", UnitSystem::measure::identity, strainXZ_}, DataEntry{"STRAINYZ", UnitSystem::measure::identity, strainYZ_}, - DataEntry{"STRESSXX", UnitSystem::measure::length, stressXX_}, - DataEntry{"STRESSYY", UnitSystem::measure::length, stressYY_}, - DataEntry{"STRESSZZ", UnitSystem::measure::length, stressZZ_}, - DataEntry{"STRESSXY", UnitSystem::measure::length, stressXY_}, - DataEntry{"STRESSXZ", UnitSystem::measure::length, stressXZ_}, - DataEntry{"STRESSYZ", UnitSystem::measure::length, stressYZ_}, + DataEntry{"STRESSXX", UnitSystem::measure::pressure, stressXX_}, + DataEntry{"STRESSYY", UnitSystem::measure::pressure, stressYY_}, + DataEntry{"STRESSZZ", UnitSystem::measure::pressure, stressZZ_}, + DataEntry{"STRESSXY", UnitSystem::measure::pressure, stressXY_}, + DataEntry{"STRESSXZ", UnitSystem::measure::pressure, stressXZ_}, + DataEntry{"STRESSYZ", UnitSystem::measure::pressure, stressYZ_}, + DataEntry{"LINSTRXX", UnitSystem::measure::pressure, linstressXX_}, + DataEntry{"LINSTRYY", UnitSystem::measure::pressure, linstressYY_}, + DataEntry{"LINSTRZZ", UnitSystem::measure::pressure, linstressZZ_}, + DataEntry{"LINSTRXY", UnitSystem::measure::pressure, linstressXY_}, + DataEntry{"LINSTRXZ", UnitSystem::measure::pressure, linstressXZ_}, + DataEntry{"LINSTRYZ", UnitSystem::measure::pressure, linstressYZ_}, + DataEntry{"FRCSTRXX", UnitSystem::measure::pressure, fracstressXX_}, + DataEntry{"FRCSTRYY", UnitSystem::measure::pressure, fracstressYY_}, + DataEntry{"FRCSTRZZ", UnitSystem::measure::pressure, fracstressZZ_}, + DataEntry{"FRCSTRXY", UnitSystem::measure::pressure, fracstressXY_}, + DataEntry{"FRCSTRXZ", UnitSystem::measure::pressure, fracstressXZ_}, + DataEntry{"FRCSTRYZ", UnitSystem::measure::pressure, fracstressYZ_}, DataEntry{"TEMPPOTF", UnitSystem::measure::pressure, mechPotentialTempForce_}, DataEntry{"TMULT_RC", UnitSystem::measure::identity, rockCompTransMultiplier_}, DataEntry{"UREA", UnitSystem::measure::density, cUrea_}, @@ -1194,6 +1206,37 @@ doAllocBuffers(const unsigned bufferSize, rstKeywords["DELSTRXY"] = 0; this->delstressYZ_.resize(bufferSize,0.0); rstKeywords["DELSTRYZ"] = 0; + + + this->fracstressXX_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXX"] = 0; + this->fracstressYY_.resize(bufferSize,0.0); + rstKeywords["FRCSTRYY"] = 0; + this->fracstressZZ_.resize(bufferSize,0.0); + rstKeywords["FRCSTRZZ"] = 0; + this->fracstressXY_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXY"] = 0; + this->fracstressXZ_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXZ"] = 0; + this->fracstressXY_.resize(bufferSize,0.0); + rstKeywords["FRCSTRXY"] = 0; + this->fracstressYZ_.resize(bufferSize,0.0); + rstKeywords["FRCSTRYZ"] = 0; + + this->linstressXX_.resize(bufferSize,0.0); + rstKeywords["LINSTRXX"] = 0; + this->linstressYY_.resize(bufferSize,0.0); + rstKeywords["LINSTRYY"] = 0; + this->linstressZZ_.resize(bufferSize,0.0); + rstKeywords["LINSTRZZ"] = 0; + this->linstressXY_.resize(bufferSize,0.0); + rstKeywords["LINSTRXY"] = 0; + this->linstressXZ_.resize(bufferSize,0.0); + rstKeywords["LINSTRXZ"] = 0; + this->linstressXY_.resize(bufferSize,0.0); + rstKeywords["LINSTRXY"] = 0; + this->linstressYZ_.resize(bufferSize,0.0); + rstKeywords["LINSTRYZ"] = 0; } // If TEMP is set in RPTRST we output temperature even if THERMAL diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index b8f3afdcc8c..a6b3ca750fd 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -531,6 +531,18 @@ class GenericOutputBlackoilModule { ScalarBuffer delstressXY_; ScalarBuffer delstressXZ_; ScalarBuffer delstressYZ_; + ScalarBuffer linstressXX_; + ScalarBuffer linstressYY_; + ScalarBuffer linstressZZ_; + ScalarBuffer linstressXY_; + ScalarBuffer linstressXZ_; + ScalarBuffer linstressYZ_; + ScalarBuffer fracstressXX_; + ScalarBuffer fracstressYY_; + ScalarBuffer fracstressZZ_; + ScalarBuffer fracstressXY_; + ScalarBuffer fracstressXZ_; + ScalarBuffer fracstressYZ_; ScalarBuffer strainXX_; ScalarBuffer strainYY_; ScalarBuffer strainZZ_; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 9104e7debd2..52e7e6ef9e4 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -152,7 +152,7 @@ class OutputBlackOilModule : public GenericOutputBlackoilModuleforceDisableFipresvOutput_ = Parameters::Get(); - if (! Parameters::Get()) { + if (!Parameters::Get()){// && !Parameters::Get()) { const std::string msg = "The output code does not support --owner-cells-first=false."; if (collectToIORank.isIORank()) { OpmLog::error(msg); @@ -233,33 +233,55 @@ class OutputBlackOilModule : public GenericOutputBlackoilModulemechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx); this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx); - this->dispX_[globalDofIdx] = model.disp(globalDofIdx, 0); - this->dispY_[globalDofIdx] = model.disp(globalDofIdx, 1); - this->dispZ_[globalDofIdx] = model.disp(globalDofIdx, 2); - this->stressXX_[globalDofIdx] = model.stress(globalDofIdx, 0); - this->stressYY_[globalDofIdx] = model.stress(globalDofIdx, 1); - this->stressZZ_[globalDofIdx] = model.stress(globalDofIdx, 2); + auto disp = model.disp(globalDofIdx,/*include_fracture*/true); + this->dispX_[globalDofIdx] = disp[0]; + this->dispY_[globalDofIdx] = disp[1]; + this->dispZ_[globalDofIdx] = disp[2]; + //total stress is not stored but calulated result is voit notation + auto stress = model.stress(globalDofIdx,/*include_fracture*/true); + this->stressXX_[globalDofIdx] = stress[0]; + this->stressYY_[globalDofIdx] = stress[1]; + this->stressZZ_[globalDofIdx] = stress[2]; // voight notation - this->stressXY_[globalDofIdx] = model.stress(globalDofIdx, 5); - this->stressXZ_[globalDofIdx] = model.stress(globalDofIdx, 4); - this->stressYZ_[globalDofIdx] = model.stress(globalDofIdx, 3); - - this->strainXX_[globalDofIdx] = model.strain(globalDofIdx, 0); - this->strainYY_[globalDofIdx] = model.strain(globalDofIdx, 1); - this->strainZZ_[globalDofIdx] = model.strain(globalDofIdx, 2); - // voight notation - this->strainXY_[globalDofIdx] = model.strain(globalDofIdx, 5); - this->strainXZ_[globalDofIdx] = model.strain(globalDofIdx, 4); - this->strainYZ_[globalDofIdx] = model.strain(globalDofIdx, 3); - - - this->delstressXX_[globalDofIdx] = model.delstress(globalDofIdx, 0); - this->delstressYY_[globalDofIdx] = model.delstress(globalDofIdx, 1); - this->delstressZZ_[globalDofIdx] = model.delstress(globalDofIdx, 2); + this->stressXY_[globalDofIdx] = stress[5]; + this->stressXZ_[globalDofIdx] = stress[4]; + this->stressYZ_[globalDofIdx] = stress[3]; + + auto strain = model.strain(globalDofIdx,/*include_fracture*/true); + this->strainXX_[globalDofIdx] = strain[0]; + this->strainYY_[globalDofIdx] = strain[1]; + this->strainZZ_[globalDofIdx] = strain[2]; // voight notation - this->delstressXY_[globalDofIdx] = model.delstress(globalDofIdx, 5); - this->delstressXZ_[globalDofIdx] = model.delstress(globalDofIdx, 4); - this->delstressYZ_[globalDofIdx] = model.delstress(globalDofIdx, 3); + this->strainXY_[globalDofIdx] = strain[5]; + this->strainXZ_[globalDofIdx] = strain[4]; + this->strainYZ_[globalDofIdx] = strain[3]; + + auto delstress = model.delstress(globalDofIdx);//not including fracture + this->delstressXX_[globalDofIdx] = delstress[ 0]; + this->delstressYY_[globalDofIdx] = delstress[ 1]; + this->delstressZZ_[globalDofIdx] = delstress[ 2]; + // voight notation + this->delstressXY_[globalDofIdx] = delstress[ 5]; + this->delstressXZ_[globalDofIdx] = delstress[ 4]; + this->delstressYZ_[globalDofIdx] = delstress[ 3]; + + auto linstress = model.linstress(globalDofIdx);; + this->linstressXX_[globalDofIdx] = linstress[ 0]; + this->linstressYY_[globalDofIdx] = linstress[ 1]; + this->linstressZZ_[globalDofIdx] = linstress[ 2]; + // voight notation + this->linstressXY_[globalDofIdx] = linstress[ 5]; + this->linstressXZ_[globalDofIdx] = linstress[ 4]; + this->linstressYZ_[globalDofIdx] = linstress[ 3]; + + auto fracstress = model.fractureStress(globalDofIdx);//is the tresagii stress which make rock fracture + this->fracstressXX_[globalDofIdx] = fracstress[ 0]; + this->fracstressYY_[globalDofIdx] = fracstress[ 1]; + this->fracstressZZ_[globalDofIdx] = fracstress[ 2]; + // voight notation + this->fracstressXY_[globalDofIdx] = fracstress[ 5]; + this->fracstressXZ_[globalDofIdx] = fracstress[ 4]; + this->fracstressYZ_[globalDofIdx] = fracstress[ 3]; } } } @@ -822,6 +844,7 @@ class OutputBlackOilModule : public GenericOutputBlackoilModule>::value) { + const auto stress = problem.geoMechModel() + .stress(globalDofIdx, /*include_fracture*/ true); + + if (key.first == "BSTRSSXX") { val.second = stress[0]; } + else if (key.first == "BSTRSSYY") { val.second = stress[1]; } + else if (key.first == "BSTRSSZZ") { val.second = stress[2]; } + else if (key.first == "BSTRSSXY") { val.second = stress[5]; } + else if (key.first == "BSTRSSXZ") { val.second = stress[4]; } + else { val.second = stress[3]; } + } + else { + val.second = 0.0; + } + } else if (key.first == "BWKR" || key.first == "BKRW") val.second = getValue(intQuants.relativePermeability(waterPhaseIdx)); else if (key.first == "BGKR" || key.first == "BKRG") @@ -1253,6 +1298,17 @@ class OutputBlackOilModule : public GenericOutputBlackoilModule + using RemoveCVR = std::remove_cv_t>; + + template + struct HasGeoMech : public std::false_type {}; + + template + struct HasGeoMech< + Problem, std::void_t().geoMechModel())> + > : public std::true_type {}; + bool isDefunctParallelWell(std::string wname) const override { if (simulator_.gridView().comm().size() == 1) diff --git a/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp b/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp index 78968b717d5..21efa8ee8ea 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp @@ -41,6 +41,7 @@ template void ghost_last_bilu0_decomposition (M& A, std::size_t interiorSize) { // iterator types + assert(interiorSize <= A.N()); using rowiterator = typename M::RowIterator; using coliterator = typename M::ColIterator; using block = typename M::block_type; @@ -311,6 +312,7 @@ ParallelOverlappingILU0(const Matrix& A, { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. + assert(interiorSize <= A_->N()); update( ); } @@ -443,6 +445,7 @@ update() if (comm_) { interiorSize_ = detail::set_interiorSize(A_->N(), interiorSize_, *comm_); + assert(interiorSize_ <= A_->N()); } // create ILU-0 decomposition diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp index bc00271ceac..1c251368393 100644 --- a/opm/simulators/linalg/PreconditionerFactory_impl.hpp +++ b/opm/simulators/linalg/PreconditionerFactory_impl.hpp @@ -395,6 +395,7 @@ struct StandardPreconditioners { // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. if (ilulevel == 0) { const std::size_t num_interior = interiorIfGhostLast(comm); + assert(num_interior <= op.getmat().N()); return std::make_shared>( op.getmat(), comm, w, MILU_VARIANT::ILU, num_interior, redblack, reorder_spheres); } else { diff --git a/opm/simulators/linalg/PropertyTree.cpp b/opm/simulators/linalg/PropertyTree.cpp index 96d43e674a9..820f73887e6 100644 --- a/opm/simulators/linalg/PropertyTree.cpp +++ b/opm/simulators/linalg/PropertyTree.cpp @@ -72,6 +72,14 @@ void PropertyTree::write_json(std::ostream &os, bool pretty) const boost::property_tree::write_json(os, *tree_, pretty); } +// PropertyTree& +// PropertyTree::get_child(const std::string& key) const +// { +// auto& pt = tree_->get_child(key); + +// return PropertyTree(pt); +// } + PropertyTree PropertyTree::get_child(const std::string& key) const { @@ -113,6 +121,7 @@ template void PropertyTree::put(const std::string& key, const std:: template void PropertyTree::put(const std::string& key, const float& value); template void PropertyTree::put(const std::string& key, const double& value); template void PropertyTree::put(const std::string& key, const int& value); +template void PropertyTree::put(const std::string& key, const bool& value); } // namespace Opm diff --git a/opm/simulators/linalg/PropertyTree.hpp b/opm/simulators/linalg/PropertyTree.hpp index 23822531005..92e89bfc874 100644 --- a/opm/simulators/linalg/PropertyTree.hpp +++ b/opm/simulators/linalg/PropertyTree.hpp @@ -50,6 +50,8 @@ class PropertyTree { template T get(const std::string& key, const T& defValue) const; + //PropertyTree& get_child(const std::string& key) const; + PropertyTree get_child(const std::string& key) const; std::optional get_child_optional(const std::string& key) const; @@ -57,7 +59,10 @@ class PropertyTree { PropertyTree& operator=(const PropertyTree& tree); void write_json(std::ostream& os, bool pretty) const; - + boost::property_tree::ptree* getBoostParamPtr() const{ + //just to be able to use full boost interface + return tree_.get(); + } protected: PropertyTree(const boost::property_tree::ptree& tree); diff --git a/opm/simulators/linalg/WriteSystemMatrixHelper.hpp b/opm/simulators/linalg/WriteSystemMatrixHelper.hpp index 496c3790beb..0063c2202b0 100644 --- a/opm/simulators/linalg/WriteSystemMatrixHelper.hpp +++ b/opm/simulators/linalg/WriteSystemMatrixHelper.hpp @@ -30,10 +30,94 @@ namespace Opm { namespace Helper { + template + void writeVector(const SimulatorType& simulator, + const VectorType& rhs, + const std::string& mname, + [[maybe_unused]] const Communicator* comm) + { + std::string dir = simulator.problem().outputDir(); + if (dir == ".") { + dir = ""; + } else if (!dir.empty() && dir.back() != '/') { + dir += "/"; + } + namespace fs = ::std::filesystem; + fs::path output_dir(dir); + fs::path subdir("reports"); + output_dir = output_dir / subdir; + if (!(fs::exists(output_dir))) { + fs::create_directory(output_dir); + } + // Combine and return. + std::ostringstream oss; + oss << "prob_" << simulator.episodeIndex() << "_time_"; + oss << std::setprecision(15) << std::setw(12) << std::setfill('0') << simulator.time() << "_"; + int nit = simulator.model().newtonMethod().numIterations(); + oss << "_nit_" << nit << "_"; + std::string output_file(oss.str()); + fs::path full_path = output_dir / output_file; + std::string prefix = full_path.string(); + { + std::string filename = prefix + mname + "vector_istl"; +#if HAVE_MPI + if (comm != nullptr) { // comm is not set in serial runs + Dune::storeMatrixMarket(rhs, filename, *comm, true); + } else +#endif + { + Dune::storeMatrixMarket(rhs, filename + ".mm"); + } + } + } + + template + void writeMatrix(const SimulatorType& simulator, + const MatrixType& matrix, + const std::string& mname, + [[maybe_unused]] const Communicator* comm) + { + std::string dir = simulator.problem().outputDir(); + if (dir == ".") { + dir = ""; + } else if (!dir.empty() && dir.back() != '/') { + dir += "/"; + } + namespace fs = ::std::filesystem; + fs::path output_dir(dir); + fs::path subdir("reports"); + output_dir = output_dir / subdir; + if (!(fs::exists(output_dir))) { + fs::create_directory(output_dir); + } + // Combine and return. + std::ostringstream oss; + oss << "prob_" << simulator.episodeIndex() << "_time_"; + oss << std::setprecision(15) << std::setw(12) << std::setfill('0') << simulator.time() << "_"; + int nit = simulator.model().newtonMethod().numIterations(); + oss << "_nit_" << nit << "_"; + std::string output_file(oss.str()); + fs::path full_path = output_dir / output_file; + std::string prefix = full_path.string(); + { + std::string filename = prefix + mname + "matrix_istl"; +#if HAVE_MPI + if (comm != nullptr) { // comm is not set in serial runs + Dune::storeMatrixMarket(matrix, filename, *comm, true); + } else +#endif + { + Dune::storeMatrixMarket(matrix, filename + ".mm"); + } + } + } + + template void writeSystem(const SimulatorType& simulator, const MatrixType& matrix, const VectorType& rhs, + const std::string& mname, [[maybe_unused]] const Communicator* comm) { std::string dir = simulator.problem().outputDir(); @@ -59,7 +143,7 @@ namespace Helper fs::path full_path = output_dir / output_file; std::string prefix = full_path.string(); { - std::string filename = prefix + "matrix_istl"; + std::string filename = prefix + mname + "matrix_istl"; #if HAVE_MPI if (comm != nullptr) { // comm is not set in serial runs Dune::storeMatrixMarket(matrix, filename, *comm, true); @@ -70,7 +154,7 @@ namespace Helper } } { - std::string filename = prefix + "rhs_istl"; + std::string filename = prefix + mname + "rhs_istl"; #if HAVE_MPI if (comm != nullptr) { // comm is not set in serial runs Dune::storeMatrixMarket(rhs, filename, *comm, true); @@ -81,6 +165,23 @@ namespace Helper } } } + template + void writeMechSystem(const SimulatorType& simulator, + const MatrixType& matrix, + const VectorType& rhs, + [[maybe_unused]] const Communicator* comm) + { + writeSystem(simulator, matrix, rhs, "mech_", comm); + } + template + void writeSystem(const SimulatorType& simulator, + const MatrixType& matrix, + const VectorType& rhs, + [[maybe_unused]] const Communicator* comm) + { + writeSystem(simulator, matrix, rhs, "flow_", comm); + } + } // namespace Helper diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 7014c71e3ea..26429bb8e4c 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -729,34 +729,44 @@ namespace Opm { { this->closed_this_step_.clear(); - // time step is finished and we are not any more at the beginning of an report step + // Time step ends and we are no longer at the beginning of a report + // step. this->report_step_starts_ = false; - const int reportStepIdx = simulator_.episodeIndex(); + const int reportStepIdx = this->simulator_.episodeIndex(); - DeferredLogger local_deferredLogger; - for (const auto& well : well_container_) { - if (getPropValue() && well->isInjector()) { - well->updateWaterThroughput(dt, this->wellState()); + if (getPropValue()) { + for (auto& wellPtr : this->well_container_) { + if (wellPtr->isInjector()) { + wellPtr->updateWaterThroughput(dt, this->wellState()); + } } } - // update connection transmissibility factor and d factor (if applicable) in the wellstate - for (const auto& well : well_container_) { - well->updateConnectionTransmissibilityFactor(simulator_, this->wellState().well(well->indexOfWell())); - well->updateConnectionDFactor(simulator_, this->wellState().well(well->indexOfWell())); + + // Update connection transmissibility factor and d factor (if + // applicable) in the wellstate. + for (auto& wellPtr : this->well_container_) { + auto& ws = this->wellState().well(wellPtr->indexOfWell()); + + wellPtr->updateConnectionTransmissibilityFactor(this->simulator_, ws); + wellPtr->updateConnectionDFactor(this->simulator_, ws); } + DeferredLogger local_deferredLogger{}; if (Indices::waterEnabled) { this->updateFiltrationModelsPostStep(dt, FluidSystem::waterPhaseIdx, local_deferredLogger); } - // WINJMULT: At the end of the time step, update the inj_multiplier saved in WellState for later use + // WINJMULT: At the end of the time step, update the inj_multiplier + // saved in WellState for later use. this->updateInjMult(local_deferredLogger); - // report well switching - for (const auto& well : well_container_) { - well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger); + // Report well switching. + for (auto& wellPtr : this->well_container_) { + wellPtr->reportWellSwitching(this->wellState().well(wellPtr->indexOfWell()), + local_deferredLogger); } - // report group switching + + // Report group switching. if (this->terminal_output_) { for (const auto& [name, ctrls] : this->switched_prod_groups_) { @@ -796,40 +806,51 @@ namespace Opm { } } - // update the rate converter with current averages pressures etc in - rateConverter_->template defineState(simulator_); + // Update the rate converter with current averages pressures etc in. + this->rateConverter_->template defineState(simulator_); - // calculate the well potentials + // Calculate the well potentials. try { this->updateWellPotentials(reportStepIdx, - /*onlyAfterEvent*/false, + /* onlyAfterEvent = */ false, simulator_.vanguard().summaryConfig(), local_deferredLogger); - } catch ( std::runtime_error& e ) { + } + catch (const std::runtime_error& e) { const std::string msg = "A zero well potential is returned for output purposes. "; local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg); } updateWellTestState(simulationTime, this->wellTestState()); - // check group sales limits at the end of the timestep - const Group& fieldGroup = this->schedule_.getGroup("FIELD", reportStepIdx); - this->checkGEconLimits(fieldGroup, simulationTime, - simulator_.episodeIndex(), local_deferredLogger); - this->checkGconsaleLimits(fieldGroup, this->wellState(), - simulator_.episodeIndex(), local_deferredLogger); + // Check group sales limits at the end of the timestep. + { + const Group& fieldGroup = this->schedule_ + .getGroup("FIELD", reportStepIdx); + + this->checkGEconLimits(fieldGroup, + simulationTime, + reportStepIdx, + local_deferredLogger); + + this->checkGconsaleLimits(fieldGroup, + this->wellState(), + reportStepIdx, + local_deferredLogger); + } this->calculateProductivityIndexValues(local_deferredLogger); this->commitWGState(); - const Opm::Parallel::Communication& comm = grid().comm(); - DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm); - if (this->terminal_output_) { + if (auto global_deferredLogger = + gatherDeferredLogger(local_deferredLogger, grid().comm()); + this->terminal_output_) + { global_deferredLogger.logMessages(); } - //reporting output temperatures + // Reporting output temperatures. this->computeWellTemperature(); } diff --git a/opm/simulators/wells/ConnFracStatistics.hpp b/opm/simulators/wells/ConnFracStatistics.hpp new file mode 100644 index 00000000000..2521fee39b6 --- /dev/null +++ b/opm/simulators/wells/ConnFracStatistics.hpp @@ -0,0 +1,139 @@ +/* + Copyright 2024 Equinor ASA. + + 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 OPM_CONNFRACSTATISTICS_HPP +#define OPM_CONNFRACSTATISTICS_HPP + +#include + +#include +#include +#include + +namespace Opm { + +/// Collection of fracturing statistics measures at the connection level. +/// +/// \tparam Scalar Statistics element type. Typically a built-in arithmetic +/// type like \c float or \c double. +template +class ConnFracStatistics +{ +public: + /// Known quantities for which this collection provides statistics + /// measures. + enum class Quantity : std::size_t + { + /// Fracture pressure + Pressure, + + /// Fracture flow rate + FlowRate, + + /// Fracture width + Width, + + // ------------------------------------------------------------- + // Helper. Must be last enumerator. + NumQuantities, + }; + + /// Sample point representation. + /// + /// Client code must populate an object of this type in order to collect + /// statistics. + using SamplePoint = std::array + >(Quantity::NumQuantities)>; + + /// Convert between byte array and object representation. + /// + /// \tparam Serializer Byte array conversion protocol. + /// + /// \param[in,out] serializer Byte array conversion object. + template + void serializeOp(Serializer& serializer) + { + serializer(this->quantity_); + } + + /// Reset internal counters to prepare for calculating a new set of + /// sample statistics. + void reset() + { + for (auto& q : this->quantity_) { q.reset(); } + } + + /// Include new element into sample. + /// + /// Updates internal statistics counters. + /// + /// \param[in] samplePoint Collection of sample values for the + /// fracturing of the current well/reservoir connection. + void addSamplePoint(const SamplePoint& samplePoint) + { + for (auto qIdx = 0*samplePoint.size(); qIdx < samplePoint.size(); ++qIdx) { + this->quantity_[qIdx].addSamplePoint(samplePoint[qIdx]); + } + } + + /// Retrieve collection of sample statistics for a single quantity. + /// + /// \param[in] q Quantity for which to retrieve sample statistics. + /// + /// \return Sample statistics for quantity \p q. + const RunningStatistics& statistics(const Quantity q) const + { + return this->quantity_[ static_cast>(q) ]; + } + + /// Create a serialisation test object. + static ConnFracStatistics serializationTestObject() + { + auto stat = ConnFracStatistics{}; + + stat.quantity_ + .fill(RunningStatistics::serializationTestObject()); + + return stat; + } + + /// Equality predicate. + /// + /// \param[in] that Object against which \code *this \endcode will + /// be tested for equality. + /// + /// \return Whether or not \code *this \endcode is the same as \p that. + bool operator==(const ConnFracStatistics& that) const + { + return this->quantity_ == that.quantity_; + } + +private: + using StatArray = std::array< + RunningStatistics, + static_cast>(Quantity::NumQuantities) + >; + + /// Collection of connection level fracturing statistics. + StatArray quantity_{}; +}; + +} // namespace Opm + +#endif // OPM_CONNFRACSTATISTICS_HPP diff --git a/opm/simulators/wells/PerfData.cpp b/opm/simulators/wells/PerfData.cpp index 66176fb359a..d8d776c82b3 100644 --- a/opm/simulators/wells/PerfData.cpp +++ b/opm/simulators/wells/PerfData.cpp @@ -83,6 +83,7 @@ PerfData PerfData::serializationTestObject() result.skin_pressure = {27.0, 28.0}; result.water_velocity = {29.0, 30.0}; result.filtrate_data = ConnFiltrateData::serializationTestObject(); + result.connFracStatistics.assign(3, ConnFracStatistics::serializationTestObject()); return result; } @@ -124,6 +125,7 @@ bool PerfData::try_assign(const PerfData& other) this->skin_pressure = other.skin_pressure; this->water_velocity = other.water_velocity; this->filtrate_data = other.filtrate_data; + this->connFracStatistics = other.connFracStatistics; return true; } @@ -152,6 +154,7 @@ bool PerfData::operator==(const PerfData& rhs) const && (this->skin_pressure == rhs.skin_pressure) && (this->water_velocity == rhs.water_velocity) && (this->filtrate_data == rhs.filtrate_data) + && (this->connFracStatistics == rhs.connFracStatistics) ; } diff --git a/opm/simulators/wells/PerfData.hpp b/opm/simulators/wells/PerfData.hpp index b6b5cb1809b..1c00251077d 100644 --- a/opm/simulators/wells/PerfData.hpp +++ b/opm/simulators/wells/PerfData.hpp @@ -21,6 +21,7 @@ #define OPM_PERFDATA_HEADER_INCLUDED #include +#include #include #include @@ -71,6 +72,7 @@ class PerfData serializer(skin_pressure); serializer(water_velocity); serializer(filtrate_data); + serializer(connFracStatistics); } bool operator==(const PerfData&) const; @@ -105,6 +107,7 @@ class PerfData std::vector water_velocity{}; ConnFiltrateData filtrate_data{}; + std::vector> connFracStatistics{}; }; } // namespace Opm diff --git a/opm/simulators/wells/RunningStatistics.hpp b/opm/simulators/wells/RunningStatistics.hpp new file mode 100644 index 00000000000..4831c5f492c --- /dev/null +++ b/opm/simulators/wells/RunningStatistics.hpp @@ -0,0 +1,163 @@ +/* + Copyright 2024 Equinor ASA. + + 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 OPM_RUNNING_STATISTICS_HPP +#define OPM_RUNNING_STATISTICS_HPP + +#include +#include +#include + +namespace Opm { + +/// Facility for calculating simple sample statistics without having full +/// sample available. +/// +/// \tparam Scalar Sample element type. Typically a built-in arithmetic +/// type like \c float or \c double. +template +class RunningStatistics +{ +public: + /// Convert between byte array and object representation. + /// + /// \tparam Serializer Byte array conversion protocol. + /// + /// \param[in,out] serializer Byte array conversion object. + template + void serializeOp(Serializer& serializer) + { + serializer(this->sampleSize_); + serializer(this->min_); + serializer(this->max_); + serializer(this->mean_); + serializer(this->totalVariance_); + } + + /// Create a serialisation test object. + static RunningStatistics serializationTestObject() + { + auto stat = RunningStatistics{}; + + stat.sampleSize_ = 12; + stat.min_ = -static_cast(1); + stat.max_ = static_cast(2); + stat.mean_ = static_cast(0.03); + stat.totalVariance_ = static_cast(0.4); + + return stat; + } + + /// Equality predicate. + /// + /// \param[in] that Object against which \code *this \endcode will + /// be tested for equality. + /// + /// \return Whether or not \code *this \endcode is the same as \p that. + bool operator==(const RunningStatistics& that) const + { + return (this->sampleSize_ == that.sampleSize_) + && (this->min_ == that.min_) + && (this->max_ == that.max_) + && (this->mean_ == that.mean_) + && (this->totalVariance_ == that.totalVariance_) + ; + } + + /// Reset internal counters to prepare for calculating a new set of + /// sample statistics. + void reset() + { + this->sampleSize_ = 0; + this->min_ = std::numeric_limits::max(); + this->max_ = std::numeric_limits::lowest(); + this->mean_ = Scalar{}; + this->totalVariance_ = Scalar{}; + } + + /// Include new element into sample. + /// + /// Updates internal statistics counters. + /// + /// \param[in] x Sample point. + void addSamplePoint(const Scalar x) + { + if (x < this->min_) { this->min_ = x; } + if (x > this->max_) { this->max_ = x; } + + // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance + + ++this->sampleSize_; + + const auto d1 = x - this->mean(); + this->mean_ += d1 / this->sampleSize_; + + const auto d2 = x - this->mean(); + this->totalVariance_ += d1 * d2; + } + + /// Retrieve current sample size. + /// + /// Effectively returns the number of calls to addSamplePoint() since + /// object was constructed or since the previous call to reset(). + std::size_t sampleSize() const { return this->sampleSize_; } + + /// Retrieve smallest sample value seen so far. + Scalar min() const { return this->min_; } + + /// Retrieve largest sample value seen so far. + Scalar max() const { return this->max_; } + + /// Retrieve arithmetic average of all sample points seen so far. + Scalar mean() const { return this->mean_; } + + /// Retrieve unbiased standard deviation of all sample points seen so + /// far. + /// + /// Returns nullopt if number of sample points is less than two. + std::optional stdev() const + { + if (this->sampleSize_ < 2) { + return {}; + } + + using std::sqrt; + return sqrt(this->totalVariance_ / (this->sampleSize_ - 1)); + } + +private: + /// Current sample size. + std::size_t sampleSize_{}; + + /// Smallest sample value seen so far. + Scalar min_ { std::numeric_limits::max() }; + + /// Largest sample value seen so far. + Scalar max_ { std::numeric_limits::lowest() }; + + /// Arithmetic average of all sample points seen so far. + Scalar mean_{}; + + /// Variance measure. In particular, N-1 * Var{x_i}. + Scalar totalVariance_{}; +}; + +} // namespace Opm + +#endif // OPM_RUNNING_STATISTICS_HPP diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 26ffb2f36fd..8bfe2368b4f 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -20,6 +20,7 @@ */ #include + #include #include @@ -33,7 +34,9 @@ #include #include #include + #include + #include #include #include @@ -45,10 +48,13 @@ #include +#include #include #include #include #include +#include +#include namespace Opm { @@ -676,6 +682,35 @@ void WellInterfaceGeneric::resetWellOperability() this->operability_status_.resetOperability(); } +template +void WellInterfaceGeneric::addPerforations(const std::vector>& perfs) +{ + for (const auto& [cell, WI, depth] : perfs) { + auto it = std::find(well_cells_.begin(), well_cells_.end(), cell); + if (it != well_cells_.end()) { + // If perforation to cell already exists, just add contribution. + const auto ind = std::distance(well_cells_.begin(), it); + well_index_[ind] += WI; + } + else { + well_cells_.push_back(cell); + well_index_.push_back(WI); + perf_depth_.push_back(depth); + + // Not strictly needed. + const double nan = std::nan("1"); + perf_rep_radius_.push_back(nan); + perf_length_.push_back(nan); + bore_diameters_.push_back(nan); + + // For now use the saturation table for the first cell. + saturation_table_number_.push_back(saturation_table_number_[0]); + + number_of_perforations_ += 1; + } + } +} + template Scalar WellInterfaceGeneric::wmicrobes_() const { diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index f511eb09761..51a16c266e3 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -189,7 +189,6 @@ class WellInterfaceGeneric { void resetWellOperability(); - virtual std::vector getPrimaryVars() const { return {}; @@ -203,6 +202,8 @@ class WellInterfaceGeneric { virtual Scalar connectionDensity(const int globalConnIdx, const int openConnIdx) const = 0; + void addPerforations(const std::vector>& perfs); + protected: bool getAllowCrossFlow() const; diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 561546cc79c..02a8129b65c 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -23,22 +23,31 @@ #include #include + #include #include #include #include -#include +#include + +#include #include #include +#include + +#include + #include -#include #include #include +#include #include #include #include +#include +#include #include #include @@ -267,27 +276,32 @@ void WellState::init(const std::vector& cellPressures, const SummaryState& summary_state, const bool enableDistributedWells) { - // call init on base class - this->base_init(cellPressures, cellTemperatures, wells_ecl, parallel_well_info, + // Call init on base class. + this->base_init(cellPressures, cellTemperatures, + wells_ecl, parallel_well_info, well_perf_data, summary_state); + this->enableDistributedWells_ = enableDistributedWells; - this->global_well_info = std::make_optional>(schedule, - report_step, - wells_ecl); + + this->global_well_info.emplace(schedule, report_step, wells_ecl); + well_rates.clear(); + this->permanently_inactive_well_names_ = schedule.getInactiveWellNamesAtEnd(); - for (const auto& wname : schedule.wellNames(report_step)) - { + + for (const auto& wname : schedule.wellNames(report_step)) { well_rates.insert({wname, std::make_pair(false, std::vector(this->numPhases()))}); } - for (const auto& winfo: parallel_well_info) - { + + for (const auto& winfo : parallel_well_info) { well_rates[winfo.get().name()].first = winfo.get().isOwner(); } - const int nw = wells_ecl.size(); + if (wells_ecl.empty()) { + return; + } - if( nw == 0 ) return ; + const int nw = wells_ecl.size(); // Initialize perfphaserates_, which must be done here. const auto& pu = this->phaseUsage(); @@ -359,71 +373,75 @@ void WellState::init(const std::vector& cellPressures, // intialize wells that have been there before // order may change so the mapping is based on the well name - if (prevState && prevState->size() > 0) { + if ((prevState != nullptr) && (prevState->size() > 0)) { for (int w = 0; w < nw; ++w) { - const Well& well = wells_ecl[w]; - if (well.getStatus() == Well::Status::SHUT) { + if (wells_ecl[w].getStatus() == Well::Status::SHUT) { continue; } - auto& new_well = this->well(w); - const auto& old_index = prevState->index(well.name()); - if (old_index.has_value()) { - const auto& prev_well = prevState->well(old_index.value()); - new_well.init_timestep(prev_well); + const auto old_index = prevState->index(wells_ecl[w].name()); + if (! old_index.has_value()) { + continue; + } - if (prev_well.status == Well::Status::SHUT) { - // Well was shut in previous state, do not use its values. - continue; - } + const auto& prev_well = prevState->well(*old_index); - if (new_well.producer != prev_well.producer) - // Well changed to/from injector from/to producer, do not use its privious values. - continue; + auto& new_well = this->well(w); + new_well.init_timestep(prev_well); - // If new target is set using WCONPROD, WCONINJE etc. we use the new control - if (!new_well.events.hasEvent(WellState::event_mask)) { - new_well.injection_cmode = prev_well.injection_cmode; - new_well.production_cmode = prev_well.production_cmode; - } + if (prev_well.status == Well::Status::SHUT) { + // Well was shut in previous state, do not use its values. + continue; + } - new_well.surface_rates = prev_well.surface_rates; - new_well.reservoir_rates = prev_well.reservoir_rates; - new_well.well_potentials = prev_well.well_potentials; - - // perfPhaseRates - const int num_perf_old_well = prev_well.perf_data.size(); - const int num_perf_this_well = new_well.perf_data.size(); - const bool global_num_perf_same = (num_perf_this_well == num_perf_old_well); - - // copy perforation rates when the number of - // perforations is equal, otherwise initialize - // perfphaserates to well rates divided by the - // number of perforations. - // TODO: we might still need the values from the prev_well if the connection structure changes - if (global_num_perf_same) - { - auto& perf_data = new_well.perf_data; - const auto& prev_perf_data = prev_well.perf_data; - perf_data.try_assign( prev_perf_data ); - } else { - const int global_num_perf_this_well = well.getConnections().num_open(); - auto& perf_data = new_well.perf_data; - auto& target_rates = perf_data.phase_rates; - for (int perf_index = 0; perf_index < num_perf_this_well; perf_index++) { - for (int p = 0; p < np; ++p) { - target_rates[perf_index*np + p] = new_well.surface_rates[p] / Scalar(global_num_perf_this_well); - } + if (new_well.producer != prev_well.producer) { + // Well changed to/from injector from/to producer, do not + // use its previous values. + continue; + } + + // If new target is set using WCONPROD, WCONINJE etc. we use the new control + if (!new_well.events.hasEvent(WellState::event_mask)) { + new_well.injection_cmode = prev_well.injection_cmode; + new_well.production_cmode = prev_well.production_cmode; + } + + new_well.surface_rates = prev_well.surface_rates; + new_well.reservoir_rates = prev_well.reservoir_rates; + new_well.well_potentials = prev_well.well_potentials; + + // perfPhaseRates + // + // Copy perforation rates when the number of perforations is + // equal, otherwise initialize perfphaserates to well rates + // divided by the number of perforations. + // + // TODO: we might still need the values from the prev_well if + // the connection structure changes. + if (const auto num_perf_this_well = new_well.perf_data.size(); + num_perf_this_well == prev_well.perf_data.size()) + { + new_well.perf_data.try_assign(prev_well.perf_data); + } + else { + const auto global_num_perf_this_well = + static_cast(wells_ecl[w].getConnections().num_open()); + + auto target_rate = new_well.perf_data.phase_rates.begin(); + for (auto perf_index = 0*num_perf_this_well; perf_index < num_perf_this_well; ++perf_index) { + for (int p = 0; p < np; ++p, ++target_rate) { + *target_rate = new_well.surface_rates[p] / global_num_perf_this_well; } } + } - // Productivity index. - new_well.productivity_index = prev_well.productivity_index; + // Productivity index. + new_well.productivity_index = prev_well.productivity_index; - // if there is no valid VFP table associated, we set the THP value to be 0. - if (well.vfp_table_number() == 0) { - new_well.thp = 0.; - } + // If there is no valid VFP table associated, we set the THP + // value to zero. + if (wells_ecl[w].vfp_table_number() == 0) { + new_well.thp = Scalar{}; } } } @@ -578,8 +596,9 @@ WellState::report(const int* globalCellIdxMap, curr.inj = ws.injection_cmode; } - const auto& pwinfo = ws.parallel_info.get(); - if (pwinfo.communication().size() == 1) { + if (const auto& pwinfo = ws.parallel_info.get(); + pwinfo.communication().size() == 1) + { reportConnections(well.connections, pu, well_index, globalCellIdxMap); } else { @@ -601,84 +620,31 @@ WellState::report(const int* globalCellIdxMap, template void WellState::reportConnections(std::vector& connections, - const PhaseUsage &pu, - std::size_t well_index, + const PhaseUsage& pu, + const std::size_t well_index, const int* globalCellIdxMap) const { - using rt = data::Rates::opt; const auto& ws = this->well(well_index); + const auto& perf_data = ws.perf_data; - const int num_perf_well = perf_data.size(); + const auto num_perf_well = perf_data.size(); + connections.resize(num_perf_well); - const auto& perf_rates = perf_data.rates; - const auto& perf_pressure = perf_data.pressure; - const auto& perf_mixing_rates = perf_data.phase_mixing_rates; - for (int i = 0; i < num_perf_well; ++i) { - const auto active_index = perf_data.cell_index[i]; - auto& connection = connections[ i ]; - connection.index = globalCellIdxMap[active_index]; - connection.pressure = perf_pressure[i]; - connection.reservoir_rate = perf_rates[i]; - connection.trans_factor = perf_data.connection_transmissibility_factor[i]; - connection.d_factor = perf_data.connection_d_factor[i]; - connection.compact_mult = perf_data.connection_compaction_tmult[i]; - connection.rates.set(rt::dissolved_gas, perf_mixing_rates[i][ws.dissolved_gas]); - connection.rates.set(rt::vaporized_oil, perf_mixing_rates[i][ws.vaporized_oil]); - if (!ws.producer) { - const auto& filtrate_data = perf_data.filtrate_data; - auto& filtrate = connection.filtrate; - filtrate.rate = filtrate_data.rates[i]; - filtrate.total = filtrate_data.total[i]; - filtrate.skin_factor = filtrate_data.skin_factor[i]; - filtrate.thickness = filtrate_data.thickness[i]; - filtrate.poro = filtrate_data.poro[i]; - filtrate.perm = filtrate_data.perm[i]; - filtrate.radius = filtrate_data.radius[i]; - filtrate.area_of_flow = filtrate_data.area_of_flow[i]; - } - } - const int np = pu.num_phases; - std::vector< rt > phs( np ); - std::vector pi(np); - if (pu.phase_used[Water]) { - phs.at( pu.phase_pos[Water] ) = rt::wat; - pi .at( pu.phase_pos[Water] ) = rt::productivity_index_water; + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + connections[i].index = globalCellIdxMap[perf_data.cell_index[i]]; } - if (pu.phase_used[Oil]) { - phs.at( pu.phase_pos[Oil] ) = rt::oil; - pi .at( pu.phase_pos[Oil] ) = rt::productivity_index_oil; - } + this->reportConnectionFactors(well_index, connections); + this->reportConnectionPressuresAndRates(well_index, pu, connections); - if (pu.phase_used[Gas]) { - phs.at( pu.phase_pos[Gas] ) = rt::gas; - pi .at( pu.phase_pos[Gas] ) = rt::productivity_index_gas; + if (! ws.producer) { + this->reportConnectionFilterCake(well_index, connections); } - std::size_t local_conn_index = 0; - for (auto& comp : connections) { - const auto * rates = &perf_data.phase_rates[np * local_conn_index]; - const auto * connPI = &perf_data.prod_index[np * local_conn_index]; - - - for (int i = 0; i < np; ++i) { - comp.rates.set( phs[ i ], rates[i] ); - comp.rates.set( pi [ i ], connPI[i] ); - } - if (pu.has_polymer) { - const auto& perf_polymer_rate = perf_data.polymer_rates; - comp.rates.set( rt::polymer, perf_polymer_rate[local_conn_index]); - } - if (pu.has_brine) { - const auto& perf_brine_rate = perf_data.brine_rates; - comp.rates.set( rt::brine, perf_brine_rate[local_conn_index]); - } - if (pu.has_solvent) { - const auto& perf_solvent_rate = perf_data.solvent_rates; - comp.rates.set( rt::solvent, perf_solvent_rate[local_conn_index] ); - } - ++local_conn_index; + if (! perf_data.connFracStatistics.empty()) { + this->reportFractureStatistics(perf_data.connFracStatistics, + connections); } } @@ -1068,6 +1034,154 @@ WellState::reportSegmentResults(const int well_id, return seg_res; } +template +void WellState:: +reportConnectionFactors(const std::size_t well_index, + std::vector& connections) const +{ + const auto& perf_data = this->well(well_index).perf_data; + const auto num_perf_well = perf_data.size(); + + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + auto& connection = connections[i]; + + connection.trans_factor = perf_data.connection_transmissibility_factor[i]; + connection.d_factor = perf_data.connection_d_factor[i]; + connection.compact_mult = perf_data.connection_compaction_tmult[i]; + } +} + +template +void WellState:: +reportConnectionPressuresAndRates(const std::size_t well_index, + const PhaseUsage& pu, + std::vector& connections) const +{ + using rt = data::Rates::opt; + + const int np = pu.num_phases; + std::vector phs(np); + std::vector pi(np); + + if (pu.phase_used[Water]) { + phs.at(pu.phase_pos[Water]) = rt::wat; + pi .at(pu.phase_pos[Water]) = rt::productivity_index_water; + } + + if (pu.phase_used[Oil]) { + phs.at(pu.phase_pos[Oil]) = rt::oil; + pi .at(pu.phase_pos[Oil]) = rt::productivity_index_oil; + } + + if (pu.phase_used[Gas]) { + phs.at(pu.phase_pos[Gas]) = rt::gas; + pi .at(pu.phase_pos[Gas]) = rt::productivity_index_gas; + } + + const auto& ws = this->well(well_index); + const auto& perf_data = ws.perf_data; + const auto num_perf_well = perf_data.size(); + + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + auto& connection = connections[i]; + + { + const auto* rates = &perf_data.phase_rates[np * i]; + const auto* connPI = &perf_data.prod_index[np * i]; + + for (int p = 0; p < np; ++p) { + connection.rates.set(phs[p], rates [p]); + connection.rates.set(pi [p], connPI[p]); + } + } + + connection.pressure = perf_data.pressure[i]; + connection.reservoir_rate = perf_data.rates[i]; + + connection.rates.set(rt::dissolved_gas, perf_data.phase_mixing_rates[i][ws.dissolved_gas]); + connection.rates.set(rt::vaporized_oil, perf_data.phase_mixing_rates[i][ws.vaporized_oil]); + } + + if (pu.has_polymer) { + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + connections[i].rates.set(rt::polymer, perf_data.polymer_rates[i]); + } + } + + if (pu.has_brine) { + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + connections[i].rates.set(rt::brine, perf_data.brine_rates[i]); + } + } + + if (pu.has_solvent) { + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + connections[i].rates.set(rt::solvent, perf_data.solvent_rates[i]); + } + } +} + +template +void WellState:: +reportConnectionFilterCake(const std::size_t well_index, + std::vector& connections) const +{ + const auto& perf_data = this->well(well_index).perf_data; + const auto num_perf_well = perf_data.size(); + + const auto& filtrate_data = perf_data.filtrate_data; + + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + auto& filtrate = connections[i].filtrate; + + filtrate.rate = filtrate_data.rates[i]; + filtrate.total = filtrate_data.total[i]; + filtrate.skin_factor = filtrate_data.skin_factor[i]; + filtrate.thickness = filtrate_data.thickness[i]; + filtrate.poro = filtrate_data.poro[i]; + filtrate.perm = filtrate_data.perm[i]; + filtrate.radius = filtrate_data.radius[i]; + filtrate.area_of_flow = filtrate_data.area_of_flow[i]; + } +} + +template +void WellState:: +reportFractureStatistics(const std::vector>& stats, + std::vector& connections) const +{ + using Quantity = typename ConnFracStatistics::Quantity; + using StatResult = data::ConnectionFracturing; + + auto connIx = 0*connections.size(); + for (auto& connection : connections) { + for (const auto& [q, result] : { + std::pair { Quantity::Pressure, &StatResult::press }, + std::pair { Quantity::FlowRate, &StatResult::rate }, + std::pair { Quantity::Width , &StatResult::width }, + }) + { + const auto& stat = stats[connIx].statistics(q); + + if (stat.sampleSize() > 0) { + auto& x = connection.fract.*result; + + x.avg = stat.mean(); + x.min = stat.min(); + x.max = stat.max(); + + if (const auto stdev = stat.stdev(); stdev.has_value()) { + x.stdev = *stdev; + } + + connection.fract.numCells = stat.sampleSize(); + } + } + + ++connIx; + } +} + template bool WellState::wellIsOwned(std::size_t well_index, [[maybe_unused]] const std::string& wellName) const diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index fb5e4d06324..b7f1a00f46c 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -40,6 +40,9 @@ #include #include +#include +#include +#include #include #include #include @@ -52,6 +55,7 @@ namespace Opm template class ParallelWellInfo; template struct PerforationData; +template class ConnFracStatistics; class Schedule; enum class WellStatus; @@ -61,7 +65,10 @@ template class WellState { public: - static const uint64_t event_mask = ScheduleEvents::WELL_STATUS_CHANGE + ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::INJECTION_UPDATE; + static const std::uint64_t event_mask = ScheduleEvents::WELL_STATUS_CHANGE + | ScheduleEvents::PRODUCTION_UPDATE + | ScheduleEvents::INJECTION_UPDATE; + // TODO: same definition with WellInterface, eventually they should go to a common header file. static const int Water = BlackoilPhases::Aqua; static const int Oil = BlackoilPhases::Liquid; @@ -437,6 +444,18 @@ class WellState const int segment, std::vector& segment_rates); + void reportConnectionFactors(const std::size_t well_index, + std::vector& connections) const; + + void reportConnectionPressuresAndRates(const std::size_t well_index, + const PhaseUsage& pu, + std::vector& connections) const; + + void reportConnectionFilterCake(const std::size_t well_index, + std::vector& connections) const; + + void reportFractureStatistics(const std::vector>& stats, + std::vector& connections) const; }; } // namespace Opm