diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4588c2dd9a8..c91dbdfa775 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 @@ -1015,14 +1052,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/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;