Numerical aquifers as auxiliary cells instead of grid cells

CMakeLists_files.cmake

@@ -1020,6 +1020,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/FlowGenericProblem_impl.hpp
   opm/simulators/flow/FlowGenericVanguard.hpp
   opm/simulators/flow/FlowMain.hpp
+  opm/simulators/flow/FlowAuxCellModule.hpp
   opm/simulators/flow/FlowProblem.hpp
   opm/simulators/flow/FlowProblemBlackoil.hpp
   opm/simulators/flow/FlowProblemBlackoilProperties.hpp
@@ -1095,6 +1096,8 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/aquifers/AquiferGridUtils.hpp
   opm/simulators/aquifers/AquiferInterface.hpp
   opm/simulators/aquifers/AquiferNumerical.hpp
+  opm/simulators/aquifers/AquiferNumericalAux.hpp
+  opm/simulators/aquifers/NumericalAquiferAuxCells.hpp
   opm/simulators/aquifers/BlackoilAquiferModel.hpp
   opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
   opm/simulators/aquifers/SupportsFaceTag.hpp

comparisonTests.cmake

@@ -31,3 +31,42 @@ add_test_compareSeparateECLFiles(
   MPI_PROCS
     1
 )
+
+
+###########################################################################
+# Numerical aquifers represented outside the grid
+#
+# Both representations solve the same discrete system -- same unknowns, pore volumes,
+# depths, regions and connection transmissibilities -- so these compare the two against
+# each other rather than against stored reference data, which is a sharper test than any
+# single-mode regression would be.  The driver pins the time steps and tightens the
+# convergence tolerances, so what is left to differ is arithmetic ordering.
+###########################################################################
+
+opm_set_test_driver(${PROJECT_SOURCE_DIR}/tests/run-numerical-aquifer-mode-comparison.sh "")
+
+# AQUNUM-02 asks for BPR at two of its own aquifer cells.  Those are grid cells in one
+# representation and not in the other, so the vectors exist in the grid-mode run alone and
+# the keyword sets cannot match; -x restricts the comparison to what both runs produce.
+# Serving block data at an aquifer cell from its auxiliary degree of freedom would remove
+# the exception.
+foreach(case AQUNUM-01 AQUNUM-02 AQUNUM-03 AQUNUM-04)
+  string(TOLOWER ${case} test)
+  set(aquifer_mode_extra_args "")
+  if(${case} STREQUAL AQUNUM-02)
+    set(aquifer_mode_extra_args -x)
+  endif()
+
+  opm_add_test(compareNumericalAquiferModes_flow+${test}
+    EXE_TARGET
+      flow
+    DRIVER_ARGS
+      -i ${OPM_TESTS_ROOT}/aquifers
+      -f ${case}
+      -r ${BASE_RESULT_PATH}/flow+aquifer_modes_${test}
+      -a ${abs_tol}
+      -t ${rel_tol}
+      -c $<TARGET_FILE:compareECL>
+      ${aquifer_mode_extra_args}
+  )
+endforeach()

opm/models/blackoil/blackoilintensivequantities.hh

@@ -135,6 +135,18 @@ class BlackOilIntensiveQuantities
     using BioeffectsIntQua = BlackOilBioeffectsIntensiveQuantities<TypeTag, enableBioeffects>;
 
 public:
+    /*!
+     * \brief Whether update() can be called without an ElementContext.
+     *
+     * The index-based overload of update() covers the plain black-oil equations and the
+     * energy module; the modules listed here still reach for the element context and have
+     * no index-based path.  Callers that have no element to offer -- an auxiliary degree
+     * of freedom has none -- have to know in advance whether they may ask.
+     */
+    static constexpr bool supportsElementContextFreeUpdate =
+        !enableSolvent && !enableExtbo && !enablePolymer && !enableFoam &&
+        !enableMICP && !enableBrine && !enableDiffusion && !enableDispersion;
+
     using FluidState = BlackOilFluidState<Evaluation,
                                           FluidSystem,
                                           energyModuleType != EnergyModules::NoTemperature,
@@ -764,20 +776,23 @@ public:
     {
         // This is the version of update() that does not use any ElementContext.
         // It is limited by some modules that are not yet adapted to that.
-        static_assert(!enableSolvent);
-        static_assert(!enableExtbo);
-        static_assert(!enablePolymer);
-        static_assert(!enableFoam);
-        static_assert(!enableMICP);
-        static_assert(!enableBrine);
-        static_assert(!enableDiffusion);
-        static_assert(!enableDispersion);
+        static_assert(supportsElementContextFreeUpdate);
 
         this->extrusionFactor_ = 1.0;// to avoid fixing parent update
         updateCommonPart<Args...>(problem, priVars, globalSpaceIdx, timeIdx);
         // Porosity requires separate calls so this can be instantiated with ReservoirProblem from the examples/ directory.
         updatePorosity(problem, priVars, globalSpaceIdx, timeIdx);
 
+        // The element-context update does this from updateCommonPart(); here it has to be
+        // called separately because that overload is shared with configurations built
+        // against a problem that has no energy at all.  Without it the fluid enthalpies,
+        // the rock internal energy and the thermal conductivity are never computed, so
+        // nothing in the cell depends on its temperature and its diagonal block comes out
+        // with an entirely zero temperature column.
+        if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
+            asImp_().updateEnergyQuantities_(problem, globalSpaceIdx, timeIdx);
+        }
+
         // TODO: Here we should do the parts for solvent etc. at the bottom of the other update() function.
     }
 

opm/models/blackoil/blackoilmodel.hh

@@ -499,8 +499,9 @@ public:
     Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
     {
         // do not care about the auxiliary equations as they are supposed to scale
-        // themselves
-        if (globalDofIdx >= this->numGridDof()) {
+        // themselves -- but an auxiliary cell holds this model's own primary variables
+        // and needs this model's scaling
+        if (!this->dofCarriesModelEquations(globalDofIdx)) {
             return 1.0;
         }
 
@@ -563,8 +564,9 @@ public:
     Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
     {
         // do not care about the auxiliary equations as they are supposed to scale
-        // themselves
-        if (globalDofIdx >= this->numGridDof()) {
+        // themselves -- but an auxiliary cell carries this model's own equations and
+        // needs this model's scaling
+        if (!this->dofCarriesModelEquations(globalDofIdx)) {
             return 1.0;
         }
 

opm/models/discretization/common/baseauxiliarymodule.hh

@@ -63,6 +63,27 @@ protected:
     using NeighborSet = std::set<unsigned>;
 
 public:
+    /*!
+     * \brief A flux connection between an auxiliary degree of freedom and another
+     *        degree of freedom of the model.
+     *
+     * This is for auxiliary modules whose degrees of freedom carry the *model's own*
+     * conservation equations -- an auxiliary "cell" with an authored volume and an
+     * authored connection list, as opposed to a genuinely different unknown such as a
+     * well's bottom hole pressure. Such a module only has to declare which pairs of
+     * degrees of freedom exchange fluxes; the discretization then assembles them with
+     * the same local residual it uses for the grid, reading the transmissibility (and
+     * the thermal/diffusive counterparts) from the problem exactly as it does for a
+     * geometric face. Both endpoints are stored as plain degree-of-freedom indices, so
+     * a connection may join two auxiliary degrees of freedom or an auxiliary one and a
+     * grid cell.
+     */
+    struct AuxiliaryConnection
+    {
+        unsigned dof1{};
+        unsigned dof2{};
+    };
+
     virtual ~BaseAuxiliaryModule() = default;
 
     /*!
@@ -82,9 +103,37 @@ public:
      * \brief Return the offset in the global system of equations for the first degree of
      *        freedom of this auxiliary module.
      */
-    int dofOffset()
+    int dofOffset() const
     { return dofOffset_; }
 
+    /*!
+     * \brief Whether this module's degrees of freedom carry the model's own
+     *        conservation equations.
+     *
+     * False by default, which describes a module whose unknowns are of a different
+     * kind -- a well's bottom hole pressure, a mortar multiplier -- and which
+     * assembles and scales its own equations in linearize().  Such degrees of freedom
+     * are deliberately kept out of the model's error norm and primary-variable
+     * switching.
+     *
+     * A module which returns true is declaring the opposite: its degrees of freedom are
+     * cells as far as the model is concerned, with the model's own unknowns, and they
+     * take part in the Newton update, the convergence measures and the variable
+     * switching exactly like a grid cell.
+     */
+    virtual bool carriesModelEquations() const
+    { return false; }
+
+    /*!
+     * \brief The volume associated with one of this module's degrees of freedom.
+     *
+     * Zero unless the module's degrees of freedom are cells.  A grid degree of freedom
+     * takes this from the geometry of its entity; an auxiliary cell has no entity, so it
+     * has to state the volume itself.
+     */
+    virtual Scalar dofVolume(unsigned /*localDofIdx*/) const
+    { return 0.0; }
+
     /*!
      * \brief Given a degree of freedom relative to the current auxiliary equation,
      *        return the corresponding index in the global system of equations.
@@ -101,6 +150,20 @@ public:
      */
     virtual void addNeighbors(std::vector<NeighborSet>& neighbors) const = 0;
 
+    /*!
+     * \brief Append this module's flux connections, if any.
+     *
+     * Modules whose degrees of freedom do not carry the model's conservation equations
+     * -- the well models, for instance, which assemble their own equations in
+     * linearize() -- leave this empty, which is the default.
+     *
+     * The discretization inserts each reported connection into the sparsity pattern in
+     * both directions and assembles it from both endpoints, so a connection must be
+     * reported exactly once, not once per endpoint.
+     */
+    virtual void addConnections(std::vector<AuxiliaryConnection>&) const
+    {}
+
     /*!
      * \brief Set the initial condition of the auxiliary module in the solution vector.
      */