Merge pull request #682 from atgeirr/aspin

Linearization support for nonlinear domain decomposition methods

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -222,8 +222,42 @@ public:
             succeeded = 1;
         }
         catch (...) {
-            std::cout << "Newton update threw an exception on rank "
+            succeeded = 0;
-                      << comm.rank() << "\n";
+        }
+        succeeded = comm.min(succeeded);
+
+        if (!succeeded)
+            throw NumericalProblem("A process did not succeed in adapting the primary variables");
+
+        numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
+    }
+
+    template <class DofIndices>
+    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])
+                const auto& curr = currentValue[Indices::zFractionIdx]; // or currentValue[pvIdx] given the block condition
-                        delta = currentValue[Indices::zFractionIdx];
+                delta = std::clamp(delta, curr - 1.0, curr);
-                if (delta < currentValue[Indices::zFractionIdx]-1.0)
-                        delta = currentValue[Indices::zFractionIdx]-1.0;
             }
             else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
                 const double sign = delta >= 0. ? 1. : -1.;

opm/models/discretization/common/fvbasediscretization.hh

@@ -72,6 +72,7 @@
 #endif
 
 #include <algorithm>
+#include <cstddef>
 #include <limits>
 #include <list>
 #include <stdexcept>
@@ -809,6 +810,36 @@ public:
         }
     }
 
+    template <class GridViewType>
+    void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
+    {
+        // loop over all elements...
+        ThreadedEntityIterator<GridViewType, /*codim=*/0> 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.
      *

opm/models/discretization/common/fvbaselinearizer.hh

@@ -152,6 +152,7 @@ public:
             delete *it;
         }
         elementCtx_.resize(0);
+        fullDomain_ = std::make_unique<FullDomain>(simulator.gridView());
     }
 
     /*!
@@ -192,6 +193,12 @@ public:
      * represented by the model object.
      */
     void linearizeDomain()
+    {
+        linearizeDomain(*fullDomain_);
+    }
+
+    template <class SubDomainType>
+    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<std::size_t>(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 <class SubDomainType>
+    void resetSystem_(const SubDomainType& domain)
+    {
+        if (!jacobian_) {
+            initFirstIteration_();
+        }
+
+        // loop over selected elements
+        using GridViewType = decltype(domain.view);
+        ThreadedEntityIterator<GridViewType, /*codim=*/0> 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 >;
+        sparsityPattern_.clear();
-        std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
+        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
+    // linearize the whole or part of the system
-    void linearize_()
+    template <class SubDomainType>
+    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<GridView, /*codim=*/0> threadedElemIt(gridView_());
+        using GridViewType = decltype(domain.view);
+        ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(domain.view);
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
         {
-            ElementIterator elemIt = threadedElemIt.beginParallel();
+            auto elemIt = threadedElemIt.beginParallel();
-            ElementIterator nextElemIt = elemIt;
+            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 <class ElementType>
+    void linearizeElement_(const ElementType& elem)
     {
         unsigned threadId = ThreadManager::threadId();
 
@@ -628,6 +685,19 @@ private:
     LinearizationType linearizationType_;
 
     std::mutex globalMatrixMutex_;
+
+    std::vector<std::set<unsigned int>> sparsityPattern_;
+
+    struct FullDomain
+    {
+        explicit FullDomain(const GridView& v) : view (v) {}
+        GridView view;
+        std::vector<bool> 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> fullDomain_;
 };
 
 } // namespace Opm

opm/models/discretization/common/tpfalinearizer.hh

@@ -50,6 +50,7 @@
 #include <set>
 #include <exception>   // current_exception, rethrow_exception
 #include <mutex>
+#include <numeric>
 
 namespace Opm::Properties {
     template<class TypeTag, class MyTypeTag>
@@ -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 <class SubDomainType>
+    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<unsigned, Constraints> constraintsMap() const
     { return {}; }
 
+    template <class SubDomainType>
+    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 <class SubDomainType>
+    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.
+                        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<Scalar, numEq> tmp;
                         IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1);
@@ -751,6 +808,12 @@ private:
     };
     std::vector<BoundaryInfo> boundaryInfo_;
     bool separateSparseSourceTerms_ = false;
+    struct FullDomain
+    {
+        std::vector<int> cells;
+        std::vector<bool> interior;
+    };
+    FullDomain fullDomain_;
 };
 
 } // namespace Opm