compositional equilibration with ZMFVD

CMakeLists_files.cmake

@@ -1083,7 +1083,9 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/equil/EquilibrationHelpers.hpp
   opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp
   opm/simulators/flow/equil/InitStateEquil.hpp
+  opm/simulators/flow/equil/InitStateEquilComp.hpp
   opm/simulators/flow/equil/InitStateEquil_impl.hpp
+  opm/simulators/flow/equil/PressureFunction.hpp
   opm/simulators/flow/rescoup/ReservoirCouplingEnabled.hpp
   opm/simulators/wells/SegmentState.hpp
   opm/simulators/wells/WellContainer.hpp

opm/simulators/flow/FlowProblemComp.hpp

@@ -1,7 +1,7 @@
 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
 // vi: set et ts=4 sw=4 sts=4:
 /*
-  Copyright 2024 SINTEF Digital
+  Copyright 2024, 2026 SINTEF Digital
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -34,6 +34,7 @@
 #include <opm/simulators/flow/FlowProblem.hpp>
 #include <opm/simulators/flow/FlowThresholdPressure.hpp>
 #include <opm/simulators/flow/OutputCompositionalModule.hpp>
+#include <opm/simulators/flow/equil/InitStateEquilComp.hpp>
 
 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
 
@@ -141,17 +142,26 @@ class FlowProblemComp : public FlowProblem<TypeTag>
                 [&vg = this->simulator().vanguard()](const unsigned int it) { return vg.gridIdxToEquilGridIdx(it); });
             updated = true;
         };
-        // TODO: we might need to do the same with FlowProblemBlackoil for parallel
 
         finishTransmissibilities();
 
         if (enableEclOutput_) {
-            eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
+            // The output of TRANX, TRANY, TRANZ and NNC is on the whole grid: the
+            // I/O rank needs the global transmissibilities when running in parallel.
+            if (simulator.vanguard().grid().comm().size() > 1) {
+                if (simulator.vanguard().grid().comm().rank() == 0) {
+                    eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
+                }
+            }
+            else {
+                eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
+            }
             std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
                 return simulator.vanguard().gridEquilIdxToGridIdx(i);
             };
             eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
         }
+        simulator.vanguard().releaseGlobalTransmissibilities();
 
         const auto& eclState = simulator.vanguard().eclState();
         const auto& schedule = simulator.vanguard().schedule();
@@ -450,7 +460,29 @@ class FlowProblemComp : public FlowProblem<TypeTag>
 
     void readEquilInitialCondition_() override
     {
-        throw std::logic_error("Equilibration is not supported by compositional modeling yet");
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
+        const auto& eclState = vanguard.eclState();
+
+        // Zero-based equilibration region of every cell (EQLNUM, or region 0).
+        std::vector<int> eqlnum(this->model().numGridDof(), 0);
+        if (eclState.fieldProps().has_int("EQLNUM")) {
+            const auto& e = eclState.fieldProps().get_int("EQLNUM");
+            std::ranges::transform(e, eqlnum.begin(), [](const int r) { return r - 1; });
+        }
+
+        EQUIL::Comp::InitialStateComputer<FluidSystem> initialState(
+            eclState,
+            getEosType(),
+            vanguard.cellCenterDepths(),
+            eqlnum,
+            vanguard.gridView().comm(),
+            this->gravity()[dimWorld - 1],
+            this->numPressurePointsEquil());
+
+        initialFluidStates_ = std::move(initialState.fluidStates());
+        // The primary variables are formed from the total composition; see initial().
+        zmf_initialization_ = true;
     }
 
     void readEclRestartSolution_()

opm/simulators/flow/NonlinearSystemCompositional.hpp

@@ -83,6 +83,8 @@ class NonlinearSystemCompositional : public NonlinearSystem<TypeTag>
 
     void solveJacobianSystem(BVector& x);
 
+    void updateSolution(const BVector& dx);
+
     bool hasNlddSolver() const
     { return false; }
 

opm/simulators/flow/NonlinearSystemCompositional_impl.hpp

@@ -274,6 +274,32 @@ relativeChange() const
     return resultDenom > 0.0 ? resultDelta / resultDenom : 0.0;
 }
 
+template <class TypeTag>
+void
+NonlinearSystemCompositional<TypeTag>::
+updateSolution(const BVector& dx)
+{
+    OPM_TIMEBLOCK(updateSolution);
+
+    auto& model = this->simulator_.model();
+    auto& solution = model.solution(/*timeIdx=*/0);
+
+    model.newtonMethod().applyUpdate(/*nextSolution=*/solution,
+                                     /*curSolution=*/solution,
+                                     /*update=*/dx,
+                                     /*resid=*/dx);
+
+    // The linear solver leaves the rows of ghost cells untouched: fetch their
+    // updated primary variables from the owning processes before the intensive
+    // quantities are recomputed.
+    model.syncOverlap();
+
+    {
+        OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
+        model.invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
+    }
+}
+
 template <class TypeTag>
 void
 NonlinearSystemCompositional<TypeTag>::
@@ -306,7 +332,12 @@ reservoirResidualMetrics() const
 
     std::vector<Scalar> residualMetrics(numEq, 0.0);
 
-    for (unsigned dofIdx = 0; dofIdx < residual.size(); ++dofIdx) {
+    // Only the interior cells: the residual of a ghost cell misses the flux
+    // contributions of neighbors outside the overlap layer and does not
+    // converge, and its converged value lives on the owning process.
+    const auto& elemMapper = model.elementMapper();
+    for (const auto& elem : elements(this->simulator_.gridView(), Dune::Partitions::interior)) {
+        const unsigned dofIdx = elemMapper.index(elem);
         if (dofIdx >= model.numGridDof() || model.dofTotalVolume(dofIdx) <= 0.0) {
             continue;
         }

opm/simulators/flow/equil/InitStateEquil.hpp

@@ -34,6 +34,7 @@
 #include <opm/material/common/Tabulated1DFunction.hpp>
 #include <opm/material/fluidstates/SimpleModularFluidState.hpp>
 
+#include <opm/simulators/flow/equil/PressureFunction.hpp>
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
 #include <array>
@@ -75,25 +76,6 @@ template<class Scalar> class EquilReg;
 namespace Miscibility { template<class Scalar> class RsFunction; }
 
 namespace Details {
-template <class Scalar, class RHS>
-class RK4IVP
-{
-public:
-    RK4IVP(const RHS& f,
-           const std::array<Scalar,2>& span,
-           const Scalar y0,
-           const int N);
-
-    Scalar operator()(const Scalar x) const;
-
-private:
-    int N_;
-    std::array<Scalar,2> span_;
-    std::vector<Scalar>  y_;
-    std::vector<Scalar>  f_;
-
-    Scalar stepsize() const;
-};
 
 namespace PhasePressODE {
 template <class FluidSystem>
@@ -262,38 +244,7 @@ class PressureTable
 
 private:
     template <class ODE>
-    class PressureFunction
-    {
-    public:
-        struct InitCond {
-            Scalar depth;
-            Scalar pressure;
-        };
-
-        explicit PressureFunction(const ODE&      ode,
-                                  const InitCond& ic,
-                                  const int       nsample,
-                                  const VSpan&    span);
-
-        PressureFunction(const PressureFunction& rhs);
-
-        PressureFunction(PressureFunction&& rhs) = default;
-
-        PressureFunction& operator=(const PressureFunction& rhs);
-
-        PressureFunction& operator=(PressureFunction&& rhs);
-
-        Scalar value(const Scalar depth) const;
-
-    private:
-        enum Direction : std::size_t { Up, Down, NumDir };
-
-        using Distribution = Details::RK4IVP<Scalar,ODE>;
-        using DistrPtr = std::unique_ptr<Distribution>;
-
-        InitCond initial_;
-        std::array<DistrPtr, Direction::NumDir> value_;
-    };
+    using PressureFunction = Details::PressureFunction<Scalar, ODE>;
 
     using OilPressODE = PhasePressODE::Oil<
         FluidSystem, typename Region::CalcDissolution
@@ -773,7 +724,7 @@ class InitialStateComputer
                                PhaseSat&               psat);
 
      template<class CellRange, class PressTable, class PhaseSat>
-     void equilibrateTiltedFaultBlock(const CellRange& cells, 
+     void equilibrateTiltedFaultBlock(const CellRange& cells,
                             const EquilReg<Scalar>& eqreg,
                             const GridView& gridView, const int numLevels,
                             const PressTable& ptable, PhaseSat& psat);