diff --git a/opm/models/blackoil/blackoilnewtonmethod.hh b/opm/models/blackoil/blackoilnewtonmethod.hh index 70cf68146..381e66c91 100644 --- a/opm/models/blackoil/blackoilnewtonmethod.hh +++ b/opm/models/blackoil/blackoilnewtonmethod.hh @@ -222,8 +222,42 @@ public: succeeded = 1; } catch (...) { - std::cout << "Newton update threw an exception on rank " - << comm.rank() << "\n"; + succeeded = 0; + } + succeeded = comm.min(succeeded); + + if (!succeeded) + throw NumericalProblem("A process did not succeed in adapting the primary variables"); + + numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_); + } + + template + void update_(SolutionVector& nextSolution, + const SolutionVector& currentSolution, + const GlobalEqVector& solutionUpdate, + const GlobalEqVector& currentResidual, + const DofIndices& dofIndices) + { + const auto& comm = this->simulator_.gridView().comm(); + + int succeeded; + try { + auto zero = solutionUpdate[0]; + zero = 0.0; + for (auto dofIdx : dofIndices) { + if (solutionUpdate[dofIdx] == zero) { + continue; + } + updatePrimaryVariables_(dofIdx, + nextSolution[dofIdx], + currentSolution[dofIdx], + solutionUpdate[dofIdx], + currentResidual[dofIdx]); + } + succeeded = 1; + } + catch (...) { succeeded = 0; } succeeded = comm.min(succeeded); @@ -326,15 +360,14 @@ protected: delta = currentValue[Indices::compositionSwitchIdx]; } } - else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) + else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) { // solvent saturation updates are also subject to the Appleyard chop delta *= satAlpha; + } else if (enableExtbo && pvIdx == Indices::zFractionIdx) { // z fraction updates are also subject to the Appleyard chop - if (delta > currentValue[Indices::zFractionIdx]) - delta = currentValue[Indices::zFractionIdx]; - if (delta < currentValue[Indices::zFractionIdx]-1.0) - delta = currentValue[Indices::zFractionIdx]-1.0; + const auto& curr = currentValue[Indices::zFractionIdx]; // or currentValue[pvIdx] given the block condition + delta = std::clamp(delta, curr - 1.0, curr); } else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) { const double sign = delta >= 0. ? 1. : -1.; diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh index 65f4d7566..aa482ce3a 100644 --- a/opm/models/discretization/common/fvbasediscretization.hh +++ b/opm/models/discretization/common/fvbasediscretization.hh @@ -72,6 +72,7 @@ #endif #include +#include #include #include #include @@ -809,6 +810,36 @@ public: } } + template + void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const + { + // loop over all elements... + ThreadedEntityIterator threadedElemIt(gridView); +#ifdef _OPENMP +#pragma omp parallel +#endif + { + + ElementContext elemCtx(simulator_); + auto elemIt = threadedElemIt.beginParallel(); + for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) { + if (elemIt->partitionType() != Dune::InteriorEntity) { + continue; + } + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + // Mark cache for this element as invalid. + const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx); + for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) { + const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx); + setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false); + } + // Update for this element. + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + } + } + } + /*! * \brief Move the intensive quantities for a given time index to the back. * diff --git a/opm/models/discretization/common/fvbaselinearizer.hh b/opm/models/discretization/common/fvbaselinearizer.hh index 1313df1b0..989809d77 100644 --- a/opm/models/discretization/common/fvbaselinearizer.hh +++ b/opm/models/discretization/common/fvbaselinearizer.hh @@ -152,6 +152,7 @@ public: delete *it; } elementCtx_.resize(0); + fullDomain_ = std::make_unique(simulator.gridView()); } /*! @@ -192,6 +193,12 @@ public: * represented by the model object. */ void linearizeDomain() + { + linearizeDomain(*fullDomain_); + } + + template + void linearizeDomain(const SubDomainType& domain) { OPM_TIMEBLOCK(linearizeDomain); // we defer the initialization of the Jacobian matrix until here because the @@ -200,9 +207,17 @@ public: if (!jacobian_) initFirstIteration_(); + // Called here because it is no longer called from linearize_(). + if (static_cast(domain.view.size(0)) == model_().numTotalDof()) { + // We are on the full domain. + resetSystem_(); + } else { + resetSystem_(domain); + } + int succeeded; try { - linearize_(); + linearize_(domain); succeeded = 1; } catch (const std::exception& e) @@ -219,7 +234,7 @@ public: << "\n" << std::flush; succeeded = 0; } - succeeded = gridView_().comm().min(succeeded); + succeeded = simulator_().gridView().comm().min(succeeded); if (!succeeded) throw NumericalProblem("A process did not succeed in linearizing the system"); @@ -315,6 +330,41 @@ public: const auto& getFloresInfo() const {return floresInfo_;} + template + void resetSystem_(const SubDomainType& domain) + { + if (!jacobian_) { + initFirstIteration_(); + } + + // loop over selected elements + using GridViewType = decltype(domain.view); + ThreadedEntityIterator threadedElemIt(domain.view); +#ifdef _OPENMP +#pragma omp parallel +#endif + { + unsigned threadId = ThreadManager::threadId(); + auto elemIt = threadedElemIt.beginParallel(); + MatrixBlock zeroBlock; + zeroBlock = 0.0; + for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) { + const Element& elem = *elemIt; + ElementContext& elemCtx = *elementCtx_[threadId]; + elemCtx.updatePrimaryStencil(elem); + // Set to zero the relevant residual and jacobian parts. + for (unsigned primaryDofIdx = 0; + primaryDofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); + ++primaryDofIdx) + { + unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, /*timeIdx=*/0); + residual_[globI] = 0.0; + jacobian_->clearRow(globI, 0.0); + } + } + } + } + private: Simulator& simulator_() { return *simulatorPtr_; } @@ -363,8 +413,8 @@ private: // for the main model, find out the global indices of the neighboring degrees of // freedom of each primary degree of freedom - using NeighborSet = std::set< unsigned >; - std::vector sparsityPattern(model.numTotalDof()); + sparsityPattern_.clear(); + sparsityPattern_.resize(model.numTotalDof()); for (const auto& elem : elements(gridView_())) { stencil.update(elem); @@ -374,7 +424,7 @@ private: for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) { unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); - sparsityPattern[myIdx].insert(neighborIdx); + sparsityPattern_[myIdx].insert(neighborIdx); } } } @@ -383,13 +433,13 @@ private: // equations size_t numAuxMod = model.numAuxiliaryModules(); for (unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx) - model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern); + model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern_); // allocate raw matrix jacobian_.reset(new SparseMatrixAdapter(simulator_())); // create matrix structure based on sparsity pattern - jacobian_->reserve(sparsityPattern); + jacobian_->reserve(sparsityPattern_); } // reset the global linear system of equations. @@ -446,11 +496,15 @@ private: } } - // linearize the whole system - void linearize_() + // linearize the whole or part of the system + template + void linearize_(const SubDomainType& domain) { OPM_TIMEBLOCK(linearize_); - resetSystem_(); + + // We do not call resetSystem_() here, since that will set + // the full system to zero, not just our part. + // Instead, that must be called before starting the linearization. // before the first iteration of each time step, we need to update the // constraints. (i.e., we assume that constraints can be time dependent, but they @@ -470,13 +524,14 @@ private: std::exception_ptr exceptionPtr = nullptr; // relinearize the elements... - ThreadedEntityIterator threadedElemIt(gridView_()); + using GridViewType = decltype(domain.view); + ThreadedEntityIterator threadedElemIt(domain.view); #ifdef _OPENMP #pragma omp parallel #endif { - ElementIterator elemIt = threadedElemIt.beginParallel(); - ElementIterator nextElemIt = elemIt; + auto elemIt = threadedElemIt.beginParallel(); + auto nextElemIt = elemIt; try { for (; !threadedElemIt.isFinished(elemIt); elemIt = nextElemIt) { // give the model and the problem a chance to prefetch the data required @@ -492,7 +547,7 @@ private: } } - const Element& elem = *elemIt; + const auto& elem = *elemIt; if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity) continue; @@ -525,8 +580,10 @@ private: applyConstraintsToLinearization_(); } + // linearize an element in the interior of the process' grid partition - void linearizeElement_(const Element& elem) + template + void linearizeElement_(const ElementType& elem) { unsigned threadId = ThreadManager::threadId(); @@ -628,6 +685,19 @@ private: LinearizationType linearizationType_; std::mutex globalMatrixMutex_; + + std::vector> sparsityPattern_; + + struct FullDomain + { + explicit FullDomain(const GridView& v) : view (v) {} + GridView view; + std::vector interior; // Should remain empty. + }; + // Simple domain object used for full-domain linearization, it allows + // us to have the same interface for sub-domain and full-domain work. + // Pointer since it must defer construction, due to GridView member. + std::unique_ptr fullDomain_; }; } // namespace Opm diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 644227e28..438c87672 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -50,6 +50,7 @@ #include #include // current_exception, rethrow_exception #include +#include namespace Opm::Properties { template @@ -187,6 +188,22 @@ public: * represented by the model object. */ void linearizeDomain() + { + linearizeDomain(fullDomain_); + } + + /*! + * \brief Linearize the part of the non-linear system of equations that is associated + * with a part of the spatial domain. + * + * That means that the Jacobian of the residual is assembled and the residual + * is evaluated for the current solution, on the domain passed in as argument. + * + * The current state of affairs (esp. the previous and the current solutions) is + * represented by the model object. + */ + template + void linearizeDomain(const SubDomainType& domain) { OPM_TIMEBLOCK(linearizeDomain); // we defer the initialization of the Jacobian matrix until here because the @@ -195,9 +212,17 @@ public: if (!jacobian_) initFirstIteration_(); + // Called here because it is no longer called from linearize_(). + if (domain.cells.size() == model_().numTotalDof()) { + // We are on the full domain. + resetSystem_(); + } else { + resetSystem_(domain); + } + int succeeded; try { - linearize_(); + linearize_(domain); succeeded = 1; } catch (const std::exception& e) @@ -214,7 +239,7 @@ public: << "\n" << std::flush; succeeded = 0; } - succeeded = gridView_().comm().min(succeeded); + succeeded = simulator_().gridView().comm().min(succeeded); if (!succeeded) throw NumericalProblem("A process did not succeed in linearizing the system"); @@ -314,6 +339,18 @@ public: const std::map constraintsMap() const { return {}; } + template + void resetSystem_(const SubDomainType& domain) + { + if (!jacobian_) { + initFirstIteration_(); + } + for (int globI : domain.cells) { + residual_[globI] = 0.0; + jacobian_->clearRow(globI, 0.0); + } + } + private: Simulator& simulator_() { return *simulatorPtr_; } @@ -434,6 +471,10 @@ private: nbInfo.matBlockAddress = jacobian_->blockAddress(nbInfo.neighbor, globI); } } + + // Create dummy full domain. + fullDomain_.cells.resize(numCells); + std::iota(fullDomain_.cells.begin(), fullDomain_.cells.end(), 0); } // reset the global linear system of equations. @@ -530,7 +571,8 @@ public: } private: - void linearize_() + template + void linearize_(const SubDomainType& domain) { // This check should be removed once this is addressed by // for example storing the previous timesteps' values for @@ -542,16 +584,23 @@ private: } OPM_TIMEBLOCK(linearize); - resetSystem_(); - unsigned numCells = model_().numTotalDof(); + + // We do not call resetSystem_() here, since that will set + // the full system to zero, not just our part. + // Instead, that must be called before starting the linearization. + const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows(); const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores(); + const unsigned int numCells = domain.cells.size(); + const bool on_full_domain = (numCells == model_().numTotalDof()); + #ifdef _OPENMP #pragma omp parallel for #endif - for (unsigned globI = 0; globI < numCells; globI++) { + for (unsigned ii = 0; ii < numCells; ++ii) { OPM_TIMEBLOCK_LOCAL(linearizationForEachCell); - const auto& nbInfos = neighborInfo_[globI]; // this is a set but should maybe be changed + const unsigned globI = domain.cells[ii]; + const auto& nbInfos = neighborInfo_[globI]; VectorBlock res(0.0); MatrixBlock bMat(0.0); ADVectorBlock adres(0.0); @@ -617,7 +666,15 @@ private: if (problem_().recycleFirstIterationStorage()) { // Assumes nothing have changed in the system which // affects masses calculated from primary variables. - model_().updateCachedStorage(globI, /*timeIdx=*/1, res); + if (on_full_domain) { + // This is to avoid resetting the start-of-step storage + // to incorrect numbers when we do local solves, where the iteration + // number will start from 0, but the starting state may not be identical + // to the start-of-step state. + // Note that a full assembly must be done before local solves + // otherwise this will be left un-updated. + model_().updateCachedStorage(globI, /*timeIdx=*/1, res); + } } else { Dune::FieldVector tmp; IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1); @@ -751,6 +808,12 @@ private: }; std::vector boundaryInfo_; bool separateSparseSourceTerms_ = false; + struct FullDomain + { + std::vector cells; + std::vector interior; + }; + FullDomain fullDomain_; }; } // namespace Opm