diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 89b30df5615..541e01b3806 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/comparisonTests.cmake b/comparisonTests.cmake index 1e95e42d272..66cf864622d 100644 --- a/comparisonTests.cmake +++ b/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 $ + ${aquifer_mode_extra_args} + ) +endforeach() diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index a766922b72e..d03f224ce90 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -135,6 +135,18 @@ class BlackOilIntensiveQuantities using BioeffectsIntQua = BlackOilBioeffectsIntensiveQuantities; 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 = BlackOilFluidStateextrusionFactor_ = 1.0;// to avoid fixing parent update updateCommonPart(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. } diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh index e0a08bf4efc..cacbae41b47 100644 --- a/opm/models/blackoil/blackoilmodel.hh +++ b/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; } diff --git a/opm/models/discretization/common/baseauxiliarymodule.hh b/opm/models/discretization/common/baseauxiliarymodule.hh index 98a85971ca2..b38504318f6 100644 --- a/opm/models/discretization/common/baseauxiliarymodule.hh +++ b/opm/models/discretization/common/baseauxiliarymodule.hh @@ -63,6 +63,27 @@ protected: using NeighborSet = std::set; 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& 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&) const + {} + /*! * \brief Set the initial condition of the auxiliary module in the solution vector. */ diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh index eae22639c53..8be1d2122ef 100644 --- a/opm/models/discretization/common/fvbasediscretization.hh +++ b/opm/models/discretization/common/fvbasediscretization.hh @@ -454,8 +454,36 @@ public: */ void finishInit() { + // Give the problem the chance to register auxiliary modules which introduce + // degrees of freedom, before anything below is sized from numTotalDof(). + simulator_.problem().registerAuxiliaryCellModules(); + + // An auxiliary module whose degrees of freedom carry the model's own equations + // needs a linearizer that assembles over degrees of freedom rather than over grid + // elements: an auxiliary DOF has no element, so an element-driven linearizer walks + // straight past it and leaves its row empty. That failure is entirely silent -- + // an all-zero row is a singular matrix at best and a spurious 0 = 0 equation at + // worst -- so refuse the combination here. + if constexpr (!Linearizer::assemblesAuxiliaryDofEquations) { + for (const auto* auxMod : auxEqModules_) { + if ((auxMod->numDofs() > 0) && auxMod->carriesModelEquations()) { + throw std::logic_error("An auxiliary module introduces degrees of freedom " + "carrying the model's conservation equations, but " + "the linearizer in use assembles element by element " + "and cannot reach them. Use the TPFA linearizer."); + } + } + } + // initialize the volume of the finite volumes to zero - const std::size_t numDof = asImp_().numGridDof(); + // + // Auxiliary modules may introduce degrees of freedom which are appended after + // the grid ones, so the per-DOF containers are sized for the total number of + // DOFs. The entries beyond the grid DOFs are authored by the auxiliary modules + // themselves: they have no grid geometry to derive a volume from. If no + // auxiliary DOFs exist (numTotalDof() == numGridDof()), this is unchanged. + const std::size_t numGridDof = asImp_().numGridDof(); + const std::size_t numDof = numTotalDof(); dofTotalVolume_.resize(numDof); std::ranges::fill(dofTotalVolume_, 0.0); @@ -489,10 +517,24 @@ public: // local process grid partition: those which do not have a non-zero volume // before taking the peer processes into account... isLocalDof_.resize(numDof); - for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) { + for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) { isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0); } + // Auxiliary DOFs have no grid entity to take a volume from, so the modules state + // it themselves. They are also not shared with peer processes via the grid's + // interior-border interface, hence local by construction. + for (const auto* auxMod : auxEqModules_) { + for (unsigned localIdx = 0; localIdx < auxMod->numDofs(); ++localIdx) { + const auto globalIdx = static_cast(auxMod->localToGlobalDof(localIdx)); + dofTotalVolume_[globalIdx] = auxMod->dofVolume(localIdx); + } + } + + for (std::size_t dofIdx = numGridDof; dofIdx < numDof; ++dofIdx) { + isLocalDof_[dofIdx] = true; + } + // add the volumes of the DOFs on the process boundaries const auto sumHandle = GridCommHandleFactory::template sumHandle(dofTotalVolume_, @@ -512,6 +554,10 @@ public: resizeAndResetIntensiveQuantitiesCache_(); newtonMethod_.finishInit(); + + // from here on the per-DOF containers are sized, so a module which introduces + // degrees of freedom can no longer be registered (see addAuxiliaryModule()) + finishInitCalled_ = true; } /*! @@ -556,6 +602,18 @@ public: } } + // Let auxiliary modules which introduce degrees of freedom set their initial + // condition. Their applyInitial() already ran when they were registered, but + // that was before the solution vector was zeroed above, so whatever they wrote + // then has just been erased. This has to happen before the history copy below + // and before the checkDefined() sweep at the end, both of which span the whole + // solution vector. + for (unsigned auxModIdx = 0; auxModIdx < numAuxiliaryModules(); ++auxModIdx) { + if (auxiliaryModule(auxModIdx)->numDofs() > 0) { + auxiliaryModule(auxModIdx)->applyInitial(); + } + } + // synchronize the ghost DOFs (if necessary) asImp_().syncOverlap(); @@ -1832,6 +1890,19 @@ public: */ void addAuxiliaryModule(BaseAuxiliaryModule* auxMod) { + // A module which introduces degrees of freedom changes numTotalDof(), and + // finishInit() sizes every per-DOF container from it. Registering such a module + // afterwards would leave all of them short -- silently, and in a way that only + // shows up much later as garbage in the intensive quantities or as an + // out-of-bounds write. Modules which declare no degrees of freedom (the well + // models, which assemble their own equations) are unaffected and may still be + // registered at any point, as they are today. + if (auxMod->numDofs() > 0 && finishInitCalled_) { + throw std::logic_error("An auxiliary module which introduces degrees of freedom must " + "be registered before the model's finishInit(), so that the " + "per-DOF containers are sized for its degrees of freedom"); + } + auxMod->setDofOffset(numTotalDof()); auxEqModules_.push_back(auxMod); @@ -1846,7 +1917,13 @@ public: solution(timeIdx).resize(numDof); } - auxMod->applyInitial(); + // A module which introduces degrees of freedom is initialised later, from + // applyInitialSolution(): at registration time the solution vector has not been + // written yet, and anything set here would be erased when it is zeroed there. + // Modules without degrees of freedom keep being initialised on the spot. + if (auxMod->numDofs() == 0) { + auxMod->applyInitial(); + } } /*! @@ -1867,6 +1944,32 @@ public: std::size_t numAuxiliaryModules() const { return auxEqModules_.size(); } + /*! + * \brief Whether the given degree of freedom carries this model's own conservation + * equations. + * + * True for every grid degree of freedom. For an auxiliary one it is up to the + * module that owns it: a well's unknowns are of a different kind and scale + * themselves, while an auxiliary *cell* carries the model's unknowns and has to + * take part in the error norm, the primary-variable switching and the convergence + * measures exactly like a grid cell. + */ + bool dofCarriesModelEquations(unsigned globalDofIdx) const + { + if (globalDofIdx < asImp_().numGridDof()) { + return true; + } + + for (const auto* auxMod : auxEqModules_) { + const auto begin = static_cast(auxMod->dofOffset()); + if (globalDofIdx >= begin && globalDofIdx < begin + auxMod->numDofs()) { + return auxMod->carriesModelEquations(); + } + } + + return false; + } + /*! * \brief Returns a given module for auxiliary equations */ @@ -1915,9 +2018,13 @@ public: protected: void resizeAndResetIntensiveQuantitiesCache_() { + // Auxiliary DOFs which carry the model's own equations need a storage cache and + // intensive quantities just like grid DOFs, so both caches are sized for the + // total number of DOFs. Without auxiliary DOFs this is unchanged. + // allocate the storage cache if (enableStorageCache()) { - const std::size_t numDof = asImp_().numGridDof(); + const std::size_t numDof = numTotalDof(); for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) { storageCache_[timeIdx].resize(numDof); storageCacheUpToDate_[timeIdx].resize(numDof, /*value=*/0); @@ -1926,7 +2033,7 @@ protected: // allocate the intensive quantities cache if (storeIntensiveQuantities()) { - const std::size_t numDof = asImp_().numGridDof(); + const std::size_t numDof = numTotalDof(); cachedIntensiveQuantityHistorySize_ = simulator_.problem().intensiveQuantityHistorySize(); const unsigned intensiveHistorySize = cachedIntensiveQuantityHistorySize_; @@ -1994,6 +2101,9 @@ protected: // a vector with all auxiliary equations to be considered std::vector*> auxEqModules_; + // guards the registration order of auxiliary modules with degrees of freedom + bool finishInitCalled_{false}; + NewtonMethod newtonMethod_; Timer prePostProcessTimer_; diff --git a/opm/models/discretization/common/fvbaselinearizer.hh b/opm/models/discretization/common/fvbaselinearizer.hh index c8aadcd1cfe..cb95027ed96 100644 --- a/opm/models/discretization/common/fvbaselinearizer.hh +++ b/opm/models/discretization/common/fvbaselinearizer.hh @@ -118,6 +118,15 @@ class FvBaseLinearizer //! \endcond public: + /*! + * \brief Whether the linearizer assembles the model's equations on auxiliary DOFs. + * + * It does not: the assembly is driven by an element loop, and an auxiliary DOF has no + * element. Auxiliary modules here are expected to assemble their rows themselves, + * from linearizeAuxiliaryEquations(), the way the well model does. + */ + static constexpr bool assemblesAuxiliaryDofEquations = false; + FvBaseLinearizer() = default; // copying the linearizer is not a good idea diff --git a/opm/models/discretization/common/fvbaseproblem.hh b/opm/models/discretization/common/fvbaseproblem.hh index aa907a66907..2d9b97c046c 100644 --- a/opm/models/discretization/common/fvbaseproblem.hh +++ b/opm/models/discretization/common/fvbaseproblem.hh @@ -286,6 +286,25 @@ public: void finishInit() {} + /*! + * \brief Called by the model before it sizes anything, so that the problem can + * register auxiliary modules which introduce degrees of freedom. + * + * This exists because of the order in which the simulator brings things up: the + * model and the problem are constructed, then the model's finishInit() runs, then + * the problem's. Every per-DOF container is sized during the model's finishInit() + * from numTotalDof(), so a module that adds degrees of freedom has to be registered + * before that -- which is earlier than the problem's own finishInit(), and hence + * needs a seam of its own. The grid and the deck are already available at this + * point, which is what such a module needs in order to know how many degrees of + * freedom it has. + * + * Modules that declare no degrees of freedom (the well models) do not need this and + * may keep registering themselves whenever they like. + */ + void registerAuxiliaryCellModules() + {} + /*! * \brief Allows to improve the performance by prefetching all data which is * associated with a given element. diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 4d77b1fb367..ba5fef5fbbc 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -63,6 +63,7 @@ #include #include #include +#include #include #include #include @@ -178,6 +179,15 @@ class TpfaLinearizer //! \endcond public: + /*! + * \brief Whether the linearizer assembles the model's equations on auxiliary DOFs. + * + * It does: the cell loop runs over every degree of freedom, and an auxiliary DOF's + * connections are in neighborInfo_ alongside the geometric ones, so the storage and + * flux terms are computed for it exactly as for a grid cell. + */ + static constexpr bool assemblesAuxiliaryDofEquations = true; + TpfaLinearizer() { simulatorPtr_ = nullptr; @@ -515,6 +525,87 @@ private: unsigned numCells = model.numTotalDof(); neighborInfo_.reserve(numCells, 6 * numCells); // Expect ~6 neighbors per cell std::vector loc_nbinfo; + + // Collect the flux connections contributed by auxiliary modules whose degrees of + // freedom carry the model's own conservation equations. + // + // A connection is reported once by its module, but it has to appear in the + // neighbour list of *both* of its endpoints: linearize_cell() only ever writes + // the row of the degree of freedom it is iterating -- residual[globI], the + // diagonal block (globI,globI) and the off-diagonal block (globJ,globI) -- so + // the equal and opposite flux into the partner, and the partner's own diagonal + // contribution, are produced only when the loop reaches the partner's row. A + // connection entered on one side alone would let the auxiliary cell exchange + // mass with a reservoir cell that never sees it in return: the Newton iteration + // would still converge, on a system that does not conserve mass. + // The GPU assembly path copies neighborInfo_, the domain and the residual to the + // device and assembles there, while an auxiliary module's linearize() runs on + // the host after the copy back. Auxiliary degrees of freedom would therefore be + // assembled inconsistently rather than wrongly-but-visibly, so refuse the + // combination outright until the device path handles them. + if constexpr (runAssemblyOnGpu) { + if (model.numAuxiliaryDof() > 0) { + throw std::logic_error("Auxiliary degrees of freedom are not supported by the GPU " + "assembly path; run the assembly on the CPU instead"); + } + } + + std::vector> auxNeighbors(numCells); + { + std::vector::AuxiliaryConnection> auxConns; + const std::size_t numAuxModules = model.numAuxiliaryModules(); + for (unsigned auxModIdx = 0; auxModIdx < numAuxModules; ++auxModIdx) { + model.auxiliaryModule(auxModIdx)->addConnections(auxConns); + } + for (const auto& conn : auxConns) { + if (conn.dof1 >= numCells || conn.dof2 >= numCells || conn.dof1 == conn.dof2) { + throw std::logic_error("Auxiliary module reported an invalid connection (" + + std::to_string(conn.dof1) + ", " + + std::to_string(conn.dof2) + ") for a model with " + + std::to_string(numCells) + " degrees of freedom"); + } + auxNeighbors[conn.dof1].push_back(conn.dof2); + auxNeighbors[conn.dof2].push_back(conn.dof1); + sparsityPattern[conn.dof1].insert(conn.dof2); + sparsityPattern[conn.dof2].insert(conn.dof1); + } + // every auxiliary degree of freedom needs a diagonal block, including one + // that currently has no connections at all + for (unsigned dofIdx = model.numGridDof(); dofIdx < numCells; ++dofIdx) { + sparsityPattern[dofIdx].insert(dofIdx); + } + } + + // Build the neighbour info for a connection that has no geometric face. The + // transmissibility carries the whole of the geometry, so the face area is unity + // and there is no face direction. + const auto makeAuxNeighborInfo = [this, gravity](unsigned myIdx, unsigned neighborIdx) { + ResidualNBInfo nbinfo{problem_().transmissibility(myIdx, neighborIdx), + 1.0, + problem_().thresholdPressure(myIdx, neighborIdx), + problem_().thresholdPressure(neighborIdx, myIdx), + (problem_().dofCenterDepth(myIdx) - + problem_().dofCenterDepth(neighborIdx)) * gravity, + FaceDir::DirEnum::Unknown, + problem_().model().dofTotalVolume(myIdx), + problem_().model().dofTotalVolume(neighborIdx), + {}, + {}, + {}, + {}}; + if constexpr (enableFullyImplicitThermal) { + nbinfo.inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx); + nbinfo.outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx); + } + if constexpr (enableDiffusion) { + nbinfo.diffusivity = problem_().diffusivity(myIdx, neighborIdx); + } + if constexpr (enableDispersion) { + nbinfo.dispersivity = problem_().dispersivity(myIdx, neighborIdx); + } + return NeighborInfoCPU{neighborIdx, nbinfo, nullptr}; + }; + for (const auto& elem : elements(gridView_())) { stencil.update(elem); @@ -565,6 +656,10 @@ private: loc_nbinfo[dofIdx - 1] = NeighborInfoCPU{neighborIdx, nbinfo, nullptr}; } } + // this grid cell's side of any auxiliary connection attached to it + for (const unsigned auxNeighbor : auxNeighbors[myIdx]) { + loc_nbinfo.push_back(makeAuxNeighborInfo(myIdx, auxNeighbor)); + } neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end()); if (problem_().nonTrivialBoundaryConditions()) { for (unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) { @@ -601,6 +696,29 @@ private: model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern); } + // Append one row per auxiliary degree of freedom. The rows of neighborInfo_ are + // addressed by degree-of-freedom index, and the grid loop above has produced + // exactly the first numGridDof() of them, so appending the auxiliary rows here + // puts each one at its own index. + for (unsigned auxDofIdx = model.numGridDof(); auxDofIdx < numCells; ++auxDofIdx) { + loc_nbinfo.clear(); + for (const unsigned neighborIdx : auxNeighbors[auxDofIdx]) { + loc_nbinfo.push_back(makeAuxNeighborInfo(auxDofIdx, neighborIdx)); + } + neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end()); + } + + // neighborInfo_ is indexed up to numTotalDof() by linearize_cell(), + // updateStoredTransmissibilities() and the loop below, so a short table is an + // out-of-bounds read rather than a missing contribution. SparseTable only + // guards this with an assert, so check it unconditionally. + if (static_cast(neighborInfo_.size()) != numCells) { + throw std::logic_error("The neighbor info table has " + + std::to_string(neighborInfo_.size()) + + " rows but the model has " + std::to_string(numCells) + + " degrees of freedom"); + } + // allocate raw matrix jacobian_ = std::make_unique(simulator_()); diagMatAddress_.resize(numCells); @@ -794,7 +912,11 @@ public: if (!enableFlows && !enableFlores && blockFlows.empty()) { return; } - const unsigned int numCells = model_().numTotalDof(); + // Flow reporting is a grid-cell concept: the rows of flowsInfo_/floresInfo_ are + // built per grid element, and the block-flow lookup below maps the DOF back to a + // cartesian index. Auxiliary DOFs have neither, so they are not visited here; + // fluxes on auxiliary connections are reported through their own client. + const unsigned int numCells = model_().numGridDof(); #ifdef _OPENMP #pragma omp parallel for #endif diff --git a/opm/models/nonlinear/newtonmethod.hh b/opm/models/nonlinear/newtonmethod.hh index f81ce17c5f6..2c4772eeeda 100644 --- a/opm/models/nonlinear/newtonmethod.hh +++ b/opm/models/nonlinear/newtonmethod.hh @@ -576,8 +576,13 @@ protected: // the solution's residual error_ = 0; for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) { - // do not consider auxiliary DOFs for the error - if (dofIdx >= model().numGridDof() || model().dofTotalVolume(dofIdx) <= 0.0) { + // Do not consider auxiliary DOFs for the error -- unless they carry the + // model's own equations, in which case they are cells as far as the model + // is concerned and their residual has to be measured like any other. + // A dormant auxiliary cell has no volume and is skipped by the second test. + if (!model().dofCarriesModelEquations(dofIdx) || + model().dofTotalVolume(dofIdx) <= 0.0) + { continue; } @@ -705,8 +710,20 @@ protected: // update the DOFs of the auxiliary equations std::size_t numDof = model().numTotalDof(); for (std::size_t dofIdx = numGridDof; dofIdx < numDof; ++dofIdx) { - nextSolution[dofIdx] = currentSolution[dofIdx]; - nextSolution[dofIdx] -= solutionUpdate[dofIdx]; + if (model().dofCarriesModelEquations(dofIdx)) { + // an auxiliary cell holds the model's own primary variables, so it needs + // the model's update -- without it the variable switching never runs and + // the cell cannot cross a phase-appearance boundary + asImp_().updatePrimaryVariables_(dofIdx, + nextSolution[dofIdx], + currentSolution[dofIdx], + solutionUpdate[dofIdx], + currentResidual[dofIdx]); + } + else { + nextSolution[dofIdx] = currentSolution[dofIdx]; + nextSolution[dofIdx] -= solutionUpdate[dofIdx]; + } } } diff --git a/opm/simulators/aquifers/AquiferNumericalAux.hpp b/opm/simulators/aquifers/AquiferNumericalAux.hpp new file mode 100644 index 00000000000..96f339498a4 --- /dev/null +++ b/opm/simulators/aquifers/AquiferNumericalAux.hpp @@ -0,0 +1,248 @@ +/* + Copyright (C) 2026 SINTEF Digital + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_AQUIFERNUMERICALAUX_HEADER_INCLUDED +#define OPM_AQUIFERNUMERICALAUX_HEADER_INCLUDED + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +namespace Opm { + +/*! + * \brief Reporting for a numerical aquifer represented as auxiliary cells. + * + * The counterpart of AquiferNumerical, which finds its cells in the grid and reads their + * fluxes off the element context. Neither is possible here -- the aquifer is not in the + * grid -- so both quantities are formed from the degrees of freedom directly: + * + * - the aquifer pressure is the water-volume-weighted water pressure of its cells, the + * same average the grid-cell representation forms over the same cells; + * - the influx is the water flux across the connections that reach into the reservoir, + * computed with the very local residual that assembled them, so the number reported is + * the number the equations used rather than a second opinion about it. + * + * The intra-aquifer chain is deliberately excluded: what the aquifer reports is what it + * gives the reservoir, not what moves inside itself. + */ +template +class AquiferNumericalAux : public AquiferInterface +{ + using Scalar = GetPropType; + using Simulator = GetPropType; + using FluidSystem = GetPropType; + using Indices = GetPropType; + using LocalResidual = GetPropType; + using Linearizer = GetPropType; + using AuxCells = NumericalAquiferAuxCells; + + using typename AquiferInterface::RateVector; + + static constexpr bool conserveSurfaceVolume = + getPropValue(); + +public: + AquiferNumericalAux(const int aquiferId, + const AuxCells& cells, + const Simulator& simulator) + : AquiferInterface(aquiferId, simulator) + , cells_(cells) + , dofs_(cells.aquiferDofs(aquiferId)) + , connections_(cells.reservoirConnections(aquiferId)) + , init_pressure_(cells.initialPressure(aquiferId)) + {} + + void initFromRestart(const data::Aquifers&) override + { + // Restart together with auxiliary cells is refused at setup; there is nothing + // sensible to restore here and pretending otherwise would hide that. + } + + void initialSolutionApplied() override + { + this->pressure_ = this->aquiferPressure(); + this->flux_rate_ = 0.0; + this->cumulative_flux_ = 0.0; + } + + void beginTimeStep() override {} + + void endTimeStep() override + { + this->pressure_ = this->aquiferPressure(); + this->flux_rate_ = this->aquiferFluxRate(); + this->cumulative_flux_ += this->flux_rate_ * this->simulator_.timeStepSize(); + } + + void addToSource(RateVector&, const unsigned, const unsigned) override + { + // The aquifer is a set of degrees of freedom with their own equations, not a + // source term on a reservoir cell. + } + + data::AquiferData aquiferData() const override + { + data::AquiferData data; + data.aquiferID = this->aquiferID(); + data.pressure = this->pressure_; + data.fluxRate = this->flux_rate_; + data.volume = this->cumulative_flux_; + + auto* aquNum = data.typeData.template create(); + aquNum->initPressure = this->init_pressure_; + + return data; + } + + void computeFaceAreaFraction(const std::vector&) override {} + + Scalar totalFaceArea() const override + { return 1.0; } + + //! Pore volume of the aquifer, for the field totals it would otherwise drop out of. + Scalar poreVolume() const + { + const auto& model = this->simulator_.model(); + + Scalar pv = 0.0; + for (const auto dof : this->dofs_) { + const auto& iq = model.intensiveQuantities(dof, /*timeIdx=*/0); + pv += model.dofTotalVolume(dof) * getValue(iq.porosity()); + } + + return pv; + } + + Scalar cumulativeFlux() const + { return this->cumulative_flux_; } + +private: + Scalar aquiferPressure() const + { + const auto& model = this->simulator_.model(); + const auto waterPos = this->phaseIdx_(); + + Scalar sumPressureWaterVolume = 0.0; + Scalar sumWaterVolume = 0.0; + + for (const auto dof : this->dofs_) { + const auto& iq = model.intensiveQuantities(dof, /*timeIdx=*/0); + const auto& fs = iq.fluidState(); + + const Scalar waterVolume = model.dofTotalVolume(dof) + * getValue(iq.porosity()) + * getValue(fs.saturation(waterPos)); + + sumPressureWaterVolume += waterVolume * getValue(fs.pressure(waterPos)); + sumWaterVolume += waterVolume; + } + + return (sumWaterVolume > 0.0) + ? sumPressureWaterVolume / sumWaterVolume + : Scalar{0}; + } + + /*! + * \brief Surface water rate from the aquifer into the reservoir. + * + * Formed with the same LocalResidual::computeFlux() the linearizer used, given the + * same neighbour information it cached, so the reported influx is the flux the + * equations were assembled with -- upwinding, gravity, threshold pressures and all -- + * rather than a reimplementation that would drift from it. + */ + Scalar aquiferFluxRate() const + { + if constexpr (! Linearizer::assemblesAuxiliaryDofEquations) { + // Cannot happen: auxiliary DOFs carrying equations are refused on such a + // linearizer when they are registered. + return Scalar{0}; + } + else { + const auto& model = this->simulator_.model(); + const auto& problem = this->simulator_.problem(); + const auto& neighborInfo = model.linearizer().getNeighborInfo(); + const auto waterPos = this->phaseIdx_(); + + Scalar rate = 0.0; + for (const auto& [aquiferDof, reservoirDof] : this->connections_) { + for (const auto& nbInfo : neighborInfo[aquiferDof]) { + if (nbInfo.neighbor != reservoirDof) { + continue; + } + + const auto& iqIn = model.intensiveQuantities(aquiferDof, 0); + const auto& iqEx = model.intensiveQuantities(reservoirDof, 0); + + RateVector flux(0.0); + RateVector darcy(0.0); + LocalResidual::computeFlux(flux, darcy, + aquiferDof, reservoirDof, + iqIn, iqEx, + nbInfo.res_nbinfo, + problem.moduleParams()); + + const auto& fsys = iqIn.fluidState().fluidSystem(); + const auto waterEqIdx = Indices::conti0EqIdx + + fsys.canonicalToActiveCompIdx(fsys.solventComponentIndex(waterPos)); + + // computeFlux() reports a flux per unit area, which the linearizer + // scales by the face area; an authored connection carries the whole + // geometry in its transmissibility, so that area is unity. + Scalar connRate = getValue(flux[waterEqIdx]) * nbInfo.res_nbinfo.faceArea; + + // ... and it is written in mass unless the model conserves surface + // volume, while the aquifer influx is reported as a surface rate. + if constexpr (! conserveSurfaceVolume) { + connRate /= fsys.referenceDensity(waterPos, + problem.pvtRegionIndex(aquiferDof)); + } + + rate += connRate; + break; + } + } + + return rate; + } + } + + const AuxCells& cells_; + std::vector dofs_{}; + std::vector> connections_{}; + std::vector init_pressure_{}; + + Scalar pressure_{0.0}; + Scalar flux_rate_{0.0}; + Scalar cumulative_flux_{0.0}; +}; + +} // namespace Opm + +#endif // OPM_AQUIFERNUMERICALAUX_HEADER_INCLUDED diff --git a/opm/simulators/aquifers/BlackoilAquiferModel.hpp b/opm/simulators/aquifers/BlackoilAquiferModel.hpp index 660cebf4ddd..7ab3329ae3b 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel.hpp @@ -99,6 +99,7 @@ class BlackoilAquiferModel void createDynamicAquifers(const int episode_index); void initializeStaticAquifers(); + void createAuxiliaryCellAquifers(); void initializeRestartDynamicAquifers(); bool needRestartDynamicAquifers() const; diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index abed4597805..9e09403a2c4 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -27,6 +27,7 @@ #endif #include +#include #include @@ -55,6 +56,7 @@ template void BlackoilAquiferModel::initialSolutionApplied() { + this->createAuxiliaryCellAquifers(); this->computeConnectionAreaFraction(); for (auto& aquifer : this->aquifers) { @@ -212,6 +214,7 @@ serializeOp(Serializer& serializer) auto* ct = dynamic_cast*>(aiPtr.get()); auto* fetp = dynamic_cast*>(aiPtr.get()); auto* num = dynamic_cast*>(aiPtr.get()); + auto* numAux = dynamic_cast*>(aiPtr.get()); auto* flux = dynamic_cast*>(aiPtr.get()); if (ct) { serializer(*ct); @@ -219,6 +222,9 @@ serializeOp(Serializer& serializer) serializer(*fetp); } else if (num) { serializer(*num); + } else if (numAux) { + // Nothing to serialize: restart is refused while auxiliary cells are live, + // and the reported values are recomputed from the solution at every step. } else if (flux) { serializer(*flux); } else { @@ -236,6 +242,42 @@ void BlackoilAquiferModel::initializeRestartDynamicAquifers() this->createDynamicAquifers(rstStep); } +/*! + * \brief Attach reporting to the numerical aquifers that live outside the grid. + * + * Deferred until the initial solution is in place: the aquifers proper are created from + * the problem's constructor, at which point neither the problem nor the auxiliary modules + * it owns exist yet, and the reported quantities are read off the degrees of freedom. + */ +template +void BlackoilAquiferModel::createAuxiliaryCellAquifers() +{ + const auto& aquifer = this->simulator_.vanguard().eclState().aquifer(); + if (! aquifer.hasNumericalAquifer()) { + return; + } + + if (this->simulator_.vanguard().eclState().numericalAquiferMode() != + NumericalAquiferMode::AuxiliaryCells) + { + return; + } + + for (const auto& module : this->simulator_.problem().auxCellModules()) { + const auto* auxAquifers = + dynamic_cast*>(module.get()); + if (auxAquifers == nullptr) { + continue; + } + + for (const auto id : auxAquifers->aquiferIds()) { + this->aquifers.push_back + (std::make_unique> + (id, *auxAquifers, this->simulator_)); + } + } +} + template void BlackoilAquiferModel::initializeStaticAquifers() { @@ -272,7 +314,15 @@ void BlackoilAquiferModel::initializeStaticAquifers() } } - if (aquifer.hasNumericalAquifer()) { + // AquiferNumerical is a reporting shim over aquifer cells that live in the grid: it + // locates them through the grid and post-computes their pressure and influx from the + // neighbouring cells' fluxes. None of that applies when the aquifer is represented + // outside the grid -- the cells it would look for are not there. + const bool numAquifersInGrid = + this->simulator_.vanguard().eclState().numericalAquiferMode() == + NumericalAquiferMode::GridCells; + + if (aquifer.hasNumericalAquifer() && numAquifersInGrid) { for (const auto& aquNum : aquifer.numericalAquifers().aquifers()) { auto aquNumPtr = std::make_unique> (aquNum.second, this->simulator_); @@ -280,6 +330,10 @@ void BlackoilAquiferModel::initializeStaticAquifers() this->aquifers.push_back(std::move(aquNumPtr)); } } + // The auxiliary-cell counterpart cannot be built here: this runs from the problem's + // own constructor, so neither the problem nor the auxiliary modules exist yet. See + // createAuxiliaryCellAquifers(). + } template diff --git a/opm/simulators/aquifers/NumericalAquiferAuxCells.hpp b/opm/simulators/aquifers/NumericalAquiferAuxCells.hpp new file mode 100644 index 00000000000..693fe2f6a91 --- /dev/null +++ b/opm/simulators/aquifers/NumericalAquiferAuxCells.hpp @@ -0,0 +1,457 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::NumericalAquiferAuxCells + */ +#ifndef OPM_NUMERICAL_AQUIFER_AUX_CELLS_HPP +#define OPM_NUMERICAL_AQUIFER_AUX_CELLS_HPP + +#include + +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + +/*! + * \brief Numerical aquifers (AQUNUM/AQUCON) represented as auxiliary cells. + * + * The alternative to letting each aquifer cell take over a grid cell. The aquifer + * satisfies the same flow equations either way, and with the same coefficients: the pore + * volume, depth and regions come from the AQUNUM record exactly as they would when + * overriding a grid cell's field properties, and the connections are the very NNCs the + * grid-cell representation would have generated -- reused here rather than reimplemented, + * so that the two representations are the same discrete system and can be compared + * directly. + * + * What differs is only where the unknown lives. The grid keeps the shape the deck gives + * it, which means no cell is resurrected to host an aquifer, no field property is + * overridden, and no non-neighbour connection is created. + */ +template +class NumericalAquiferAuxCells : public FlowAuxCellModule +{ + using ParentType = FlowAuxCellModule; + using Scalar = GetPropType; + using Simulator = GetPropType; + using SparseMatrixAdapter = GetPropType; + using GlobalEqVector = GetPropType; + using FluidSystem = GetPropType; + using GridView = GetPropType; + + enum { dimWorld = GridView::dimensionworld }; + +public: + using Connection = typename ParentType::Connection; + + explicit NumericalAquiferAuxCells(Simulator& simulator) + : simulator_(simulator) + { + const auto& eclState = simulator_.vanguard().eclState(); + const auto& aquifers = eclState.aquifer().numericalAquifers(); + + // Enumerate aquifer by aquifer and, within an aquifer, in the order its cells + // were declared. Not allAquiferCells(), which is an unordered map: the degree of + // freedom numbering would then depend on the hash order, and the chain structure + // -- which cell of an aquifer carries the reservoir connections -- would be lost. + for (const auto& [id, aquifer] : aquifers.aquifers()) { + const auto first = this->cells_.size(); + + for (std::size_t i = 0; i < aquifer.numCells(); ++i) { + const auto* cell = aquifer.getCellPrt(i); + this->cells_.push_back(cell); + this->cartesianToLocal_.emplace(cell->global_index, this->cells_.size() - 1); + } + + this->aquiferRange_.emplace(id, std::make_pair(first, this->cells_.size())); + } + + // The cells these records name were deactivated when the grid was built, which + // is what makes room for the aquifer to live outside it. + this->checkAquiferCellsAreNotInterior(); + } + + unsigned numDofs() const override + { return static_cast(this->cells_.size()); } + + Scalar poreVolume(unsigned localIdx) const override + { return this->cells_.at(localIdx)->poreVolume(); } + + Scalar bulkVolume(unsigned localIdx) const override + { return this->cells_.at(localIdx)->cellVolume(); } + + Scalar depth(unsigned localIdx) const override + { return this->cells_.at(localIdx)->depth; } + + unsigned pvtRegionIndex(unsigned localIdx) const override + { return static_cast(this->cells_.at(localIdx)->pvttable) - 1; } + + unsigned satRegionIndex(unsigned localIdx) const override + { return static_cast(this->cells_.at(localIdx)->sattable) - 1; } + + int hostCartesianIndex(unsigned localIdx) const override + { return static_cast(this->cells_.at(localIdx)->global_index); } + + /*! + * \brief The reservoir cell this aquifer cell hangs off. + * + * Used to start the aquifer from a sensible fluid state. All AQUCON connections of + * an aquifer attach to its first cell, so the whole chain is initialised from that + * cell's first reservoir neighbour. + */ + unsigned initialisationPartner(unsigned localIdx) const override + { return this->initialisationPartner_.at(localIdx); } + + void connections(std::vector& conns) const override + { conns.insert(conns.end(), this->connections_.begin(), this->connections_.end()); } + + /*! + * \brief Build the connection list. + * + * Deferred out of the constructor because it needs the grid's cartesian-to-compressed + * mapping, which is only meaningful once the grid has been distributed. + */ + void buildConnections() + { + const auto& vanguard = simulator_.vanguard(); + const auto& eclState = vanguard.eclState(); + const auto& aquifers = eclState.aquifer().numericalAquifers(); + + this->connections_.clear(); + this->initialisationPartner_.assign(this->cells_.size(), 0); + this->hasReservoirConnection_.assign(this->cells_.size(), false); + + // The grid-cell representation puts these connections into the input NNC list, + // where the region-based multipliers are applied to them alongside the deck's own + // NNCs (Transmissibility::applyMultRegTToInputNncTrans_). A MULTREGT record can + // reduce an aquifer's connection to the reservoir by orders of magnitude, or shut + // it off entirely, so the same multiplier has to be applied here -- otherwise the + // aquifer is connected on one path and not on the other. + const auto& transMult = eclState.getTransMult(); + const auto multiplier = [&transMult](const std::size_t cell1, const std::size_t cell2) { + return static_cast(transMult.getRegionMultiplierNNC(cell1, cell2)); + }; + + // Aquifer cell to aquifer cell: the chain within one aquifer. Both endpoints + // are auxiliary cells. + for (const auto& nnc : aquifers.aquiferCellNNCs()) { + const auto dof1 = this->auxDofOf(nnc.cell1); + const auto dof2 = this->auxDofOf(nnc.cell2); + const auto trans = static_cast(nnc.trans) * multiplier(nnc.cell1, nnc.cell2); + this->connections_.push_back({dof1, dof2, trans, 0.0, 0.0}); + } + + // Aquifer cell to reservoir cell. The second endpoint is a real grid cell, so it + // has to be translated into the local compressed numbering -- and skipped when + // this rank does not own it. + // NNCdata normalises the order of its two cells, so which endpoint is the + // aquifer cell is not fixed; identify it by membership rather than by position. + std::size_t skipped = 0; + for (const auto& nnc : aquifers.aquiferConnectionNNCs(eclState.getInputGrid(), + eclState.fieldProps())) + { + const bool firstIsAquifer = this->cartesianToLocal_.count(nnc.cell1) > 0; + const auto aquiferCartesian = firstIsAquifer ? nnc.cell1 : nnc.cell2; + const auto reservoirCartesian = firstIsAquifer ? nnc.cell2 : nnc.cell1; + + const auto dof1 = this->auxDofOf(aquiferCartesian); + const int reservoirCell = vanguard.compressedIndexForInterior(reservoirCartesian); + if (reservoirCell < 0) { + ++skipped; + continue; + } + + const auto dof2 = static_cast(reservoirCell); + const auto trans = static_cast(nnc.trans) + * multiplier(aquiferCartesian, reservoirCartesian); + this->connections_.push_back({dof1, dof2, trans, 0.0, 0.0}); + + const auto localIdx = this->localOf(aquiferCartesian); + if (!this->hasReservoirConnection_[localIdx]) { + this->initialisationPartner_[localIdx] = dof2; + this->hasReservoirConnection_[localIdx] = true; + } + } + + // Only the cell that carries the AQUCON connections has a reservoir neighbour of + // its own; the rest of the chain hangs off it. Give them all that cell's + // neighbour to start from -- per aquifer, so that two aquifers cannot borrow each + // other's. Note this is a fallback for initialisation only: it says where a cell + // takes its initial state from, not what it is connected to. + for (const auto& [id, range] : this->aquiferRange_) { + unsigned connected = 0; + bool found = false; + for (auto i = range.first; i < range.second; ++i) { + if (this->hasReservoirConnection_.at(i)) { + connected = this->initialisationPartner_[i]; + found = true; + break; + } + } + + if (!found) { + OPM_THROW(std::runtime_error, + fmt::format("Numerical aquifer {} has no connection to any " + "reservoir cell owned by this process", id)); + } + + for (auto i = range.first; i < range.second; ++i) { + if (!this->hasReservoirConnection_.at(i)) { + this->initialisationPartner_[i] = connected; + } + } + } + + if (skipped > 0) { + OpmLog::debug(fmt::format("Numerical aquifer: {} connection(s) target cells " + "not owned by this rank", skipped)); + } + } + + /*! + * \brief Set the initial state of the aquifer cells. + * + * They cannot be equilibrated the ordinary way -- that needs the cell's geometry -- + * so each one starts from the state of the reservoir cell it is connected to, with + * the phase pressures carried to its own depth and the cell filled with water. That + * is the same thing the grid-cell representation arranges for its aquifer cells, and + * an explicit initial pressure on the AQUNUM record overrides it either way. + */ + void applyInitial() override + { + auto& solution = simulator_.model().solution(/*timeIdx=*/0); + const auto& problem = simulator_.problem(); + const auto gravity = problem.gravity()[dimWorld - 1]; + + for (unsigned localIdx = 0; localIdx < this->numDofs(); ++localIdx) { + const auto globalIdx = static_cast(this->localToGlobalDof(localIdx)); + const auto partner = this->initialisationPartner_.at(localIdx); + + auto fs = problem.initialFluidState(partner); + + const auto waterPos = FluidSystem::waterPhaseIdx; + const auto rho = getValue(fs.density(waterPos)); + const auto dz = this->depth(localIdx) - problem.dofCenterDepth(partner); + + const auto& cell = *this->cells_.at(localIdx); + for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) { + if (!FluidSystem::phaseIsActive(phase)) { + continue; + } + + fs.setSaturation(phase, (phase == waterPos) ? 1.0 : 0.0); + fs.setPressure(phase, cell.init_pressure.has_value() + ? cell.init_pressure.value() + : getValue(fs.pressure(phase)) + rho * gravity * dz); + } + + // Order matters, as it does on the grid path: assignNaive() decides what the + // primary variables mean and needs the PVT region to do it. + solution[globalIdx].setPvtRegionIndex(this->pvtRegionIndex(localIdx)); + solution[globalIdx].assignNaive(fs); + } + } + + void linearize(SparseMatrixAdapter&, GlobalEqVector&) override + { + // Nothing to add: the aquifer cells are assembled by the model's own local + // residual, like grid cells, and they are always active so there are no dormant + // rows to condition. + } + + //! The aquifer identifiers this module represents, in declaration order. + std::vector aquiferIds() const + { + std::vector ids; + ids.reserve(this->aquiferRange_.size()); + for (const auto& [id, range] : this->aquiferRange_) { + static_cast(range); + ids.push_back(static_cast(id)); + } + + return ids; + } + + //! The global degrees of freedom carrying one aquifer, in declaration order. + std::vector aquiferDofs(const int aquiferId) const + { + std::vector dofs; + + auto pos = this->aquiferRange_.find(static_cast(aquiferId)); + if (pos == this->aquiferRange_.end()) { + return dofs; + } + + dofs.reserve(pos->second.second - pos->second.first); + for (auto localIdx = pos->second.first; localIdx < pos->second.second; ++localIdx) { + dofs.push_back(static_cast(this->localToGlobalDof(localIdx))); + } + + return dofs; + } + + /*! + * \brief The connections of one aquifer that reach into the grid. + * + * The intra-aquifer chain is left out: what an aquifer reports as its influx is what + * crosses into the reservoir, not what moves inside itself. Each entry is + * (aquifer degree of freedom, reservoir degree of freedom). + */ + std::vector> reservoirConnections(const int aquiferId) const + { + std::vector> conns; + + const auto dofs = this->aquiferDofs(aquiferId); + const auto isMine = [&dofs](const unsigned dof) { + return std::find(dofs.begin(), dofs.end(), dof) != dofs.end(); + }; + const auto isGridDof = [this](const unsigned dof) { + return dof < this->simulator_.model().numGridDof(); + }; + + for (const auto& conn : this->connections_) { + if (isMine(conn.dof1) && isGridDof(conn.dof2)) { + conns.emplace_back(conn.dof1, conn.dof2); + } + else if (isMine(conn.dof2) && isGridDof(conn.dof1)) { + conns.emplace_back(conn.dof2, conn.dof1); + } + } + + return conns; + } + + //! Initial pressure of one aquifer's cells, for the restart record. + std::vector initialPressure(const int aquiferId) const + { + std::vector pressures; + + auto pos = this->aquiferRange_.find(static_cast(aquiferId)); + if (pos == this->aquiferRange_.end()) { + return pressures; + } + + pressures.reserve(pos->second.second - pos->second.first); + for (auto localIdx = pos->second.first; localIdx < pos->second.second; ++localIdx) { + const auto& cell = *this->cells_.at(localIdx); + pressures.push_back(cell.init_pressure.has_value() + ? static_cast(cell.init_pressure.value()) + : Scalar{0}); + } + + return pressures; + } + +private: + unsigned auxDofOf(std::size_t cartesianIndex) const + { return static_cast(this->localToGlobalDof(this->localOf(cartesianIndex))); } + + unsigned localOf(std::size_t cartesianIndex) const + { + auto pos = this->cartesianToLocal_.find(cartesianIndex); + if (pos == this->cartesianToLocal_.end()) { + OPM_THROW(std::logic_error, + fmt::format("Numerical aquifer connection names cell {} which is " + "not an aquifer cell", cartesianIndex)); + } + + return static_cast(pos->second); + } + + /*! + * \brief Refuse an aquifer cell buried inside the model. + * + * Deactivating the cell an AQUNUM record names is only sound where that cell is at + * the edge of the model. With live rock both above and below it, removing it opens + * a hole in the middle of a column: the neighbours lose a connection they would + * otherwise have, and with PINCH active the grid processing may bridge across the + * gap and connect them to each other instead. Neither is what the deck describes, + * and neither matches what the grid-cell representation does, so the two modes would + * quietly stop being comparable. + * + * Placing an aquifer cell there is questionable modelling to begin with -- a + * numerical aquifer is meant to hang off the model, not to sit inside it -- so this + * is refused rather than approximated. The proper answer is for a numerical aquifer + * to be defined independently of the grid and then connected to it, which is the + * direction the auxiliary-cell representation is going; the restriction can be + * lifted once the aquifer no longer needs a cell to name at all. + */ + void checkAquiferCellsAreNotInterior() const + { + const auto& grid = simulator_.vanguard().eclState().getInputGrid(); + const auto nz = grid.getNZ(); + + for (const auto* cell : this->cells_) { + if ((cell->K == 0) || (cell->K + 1 >= nz)) { + continue; // at the top or the bottom of the model + } + + const auto above = grid.getGlobalIndex(cell->I, cell->J, cell->K - 1); + const auto below = grid.getGlobalIndex(cell->I, cell->J, cell->K + 1); + + if (grid.cellActive(above) && grid.cellActive(below)) { + OPM_THROW(std::runtime_error, + fmt::format("AQUNUM record for aquifer {} names cell " + "({},{},{}), which has active cells both above and " + "below it. Representing numerical aquifers outside " + "the grid removes the cell they name, which would " + "open a hole inside the model; place the aquifer " + "cell at the edge of the model, or run with the " + "grid-cell representation.", + cell->aquifer_id, + cell->I + 1, cell->J + 1, cell->K + 1)); + } + } + } + + Simulator& simulator_; + + //! Aquifer cells, in auxiliary-DOF order. + std::vector cells_{}; + + //! Cartesian index named by an AQUNUM record -> auxiliary-DOF-local index. + std::unordered_map cartesianToLocal_{}; + + std::vector connections_{}; + std::vector initialisationPartner_{}; + std::vector hasReservoirConnection_{}; + + //! aquifer id -> [first, last) range of its cells in the local numbering + std::map> aquiferRange_{}; +}; + +} // namespace Opm + +#endif // OPM_NUMERICAL_AQUIFER_AUX_CELLS_HPP diff --git a/opm/simulators/flow/EclWriter.hpp b/opm/simulators/flow/EclWriter.hpp index 434a87e764a..2fa02f1a85f 100644 --- a/opm/simulators/flow/EclWriter.hpp +++ b/opm/simulators/flow/EclWriter.hpp @@ -845,6 +845,20 @@ class EclWriter : public EclGenericWriter this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume); } + + // Degrees of freedom introduced by an auxiliary module hold real fluid in a + // real pore volume, so they belong in the field and region totals just as the + // grid cells do -- and they have to be, for those totals to stay comparable + // with a run that represents the same thing inside the grid. They are + // reached by index here exactly as the grid cells are; what they are not part + // of is the per-cell output, which stops at the grid. + const auto& model = simulator_.model(); + for (unsigned dofIdx = model.numGridDof(); dofIdx < model.numTotalDof(); ++dofIdx) { + const auto& intQuants = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/0); + + this->outputModule_->updateFluidInPlace(dofIdx, intQuants, + model.dofTotalVolume(dofIdx)); + } } this->outputModule_->validateLocalData(); diff --git a/opm/simulators/flow/FIBlackoilModel.hpp b/opm/simulators/flow/FIBlackoilModel.hpp index 5a56ed016f8..9a1fd0bf636 100644 --- a/opm/simulators/flow/FIBlackoilModel.hpp +++ b/opm/simulators/flow/FIBlackoilModel.hpp @@ -87,21 +87,22 @@ class FIBlackOilModel : public BlackOilModel if constexpr (gridIsUnchanging) { if constexpr (avoidElementContext) { updateCachedIntQuants(timeIdx); - return; } - OPM_BEGIN_PARALLEL_TRY_CATCH(); + else { + OPM_BEGIN_PARALLEL_TRY_CATCH(); #ifdef _OPENMP #pragma omp parallel for #endif - for (const auto& chunk : element_chunks_) { - ElementContext elemCtx(this->simulator_); - for (const auto& elem : chunk) { - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(timeIdx); + for (const auto& chunk : element_chunks_) { + ElementContext elemCtx(this->simulator_); + for (const auto& elem : chunk) { + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(timeIdx); + } } + OPM_END_PARALLEL_TRY_CATCH("invalidateAndUpdateIntensiveQuantities: state error", + this->simulator_.vanguard().grid().comm()); } - OPM_END_PARALLEL_TRY_CATCH("invalidateAndUpdateIntensiveQuantities: state error", - this->simulator_.vanguard().grid().comm()); } else { // Grid is possibly refined or otherwise changed between calls. ElementContext elemCtx(this->simulator_); @@ -110,6 +111,50 @@ class FIBlackOilModel : public BlackOilModel elemCtx.updatePrimaryIntensiveQuantities(timeIdx); } } + + updateAuxiliaryIntQuants(timeIdx); + } + + /*! + * \brief Update the intensive quantities of the degrees of freedom introduced by + * auxiliary modules. + * + * These are not reachable through the grid, so none of the loops above visit them. + * The update itself needs no element context: it is driven entirely by the DOF + * index and the problem's index-based accessors, which is what allows an auxiliary + * cell to carry the model's own equations. + * + * The index-based update does not cover every module -- solvent, extbo, polymer, + * foam, MICP, brine, diffusion and dispersion still need an element -- so it is + * instantiated only where the intensive quantities say it is available. Note that + * this is a weaker condition than AvoidElementContext, which is about how the *grid* + * cells are updated: a configuration may well drive the grid through element contexts + * and still be able to update an auxiliary DOF without one. A configuration that has + * auxiliary DOFs and genuinely cannot update them says so rather than leaving them + * uninitialised. + */ + void updateAuxiliaryIntQuants(const unsigned timeIdx) const + { + if (this->numTotalDof() == this->numGridDof()) { + return; + } + + if constexpr (!IntensiveQuantities::supportsElementContextFreeUpdate) { + throw std::logic_error("Auxiliary degrees of freedom need intensive quantities " + "updated without an element context, which this model " + "configuration does not support"); + } + else { + if (!this->storeIntensiveQuantities()) { + return; + } + + const unsigned numGridDof = this->numGridDof(); + const unsigned numTotalDof = this->numTotalDof(); + for (unsigned globalIdx = numGridDof; globalIdx < numTotalDof; ++globalIdx) { + this->updateSingleCachedIntQuantUnchecked(globalIdx, timeIdx); + } + } } void invalidateAndUpdateIntensiveQuantitiesOverlap(unsigned timeIdx) const diff --git a/opm/simulators/flow/FIPContainer.cpp b/opm/simulators/flow/FIPContainer.cpp index e2284e69e0c..fdce2e24817 100644 --- a/opm/simulators/flow/FIPContainer.cpp +++ b/opm/simulators/flow/FIPContainer.cpp @@ -39,6 +39,7 @@ template bool FIPContainer:: allocate(const std::size_t bufferSize, + const std::size_t gridSize, const SummaryConfig& summaryConfig, const bool forceAlloc, std::map& rstKeywords) @@ -64,6 +65,7 @@ allocate(const std::size_t bufferSize, bool computeFip = false; bufferSize_ = bufferSize; + gridSize_ = gridSize; for (const auto& phase : Inplace::phases()) { if (forceAlloc || summaryConfig.require3DField(Inplace::EclString(phase))) { this->add(phase); @@ -506,9 +508,21 @@ outputRestart(data::Solution& sol) }); } + // The restart file is written per grid cell. Where auxiliary degrees of freedom + // extend these buffers past the grid they contribute to the region sums, which is + // what they are for, but they have no cell to be written against. + const auto perCell = [this](const Inplace::Phase phase) { + auto& v = this->fip_[phase]; + if (v.size() <= this->gridSize_) { + return std::move(v); + } + + return std::vector(v.begin(), v.begin() + this->gridSize_); + }; + for (const auto& [mnemonic, unit, phase] : fipArrays) { if (! this->fip_[phase].empty()) { - sol.insert(mnemonic, unit, std::move(this->fip_[phase]), + sol.insert(mnemonic, unit, perCell(phase), data::TargetType::RESTART_SOLUTION); } } @@ -517,7 +531,7 @@ outputRestart(data::Solution& sol) if (! this->fip_[phase].empty()) { sol.insert(Inplace::EclString(phase), UnitSystem::measure::volume, - std::move(this->fip_[phase]), + perCell(phase), data::TargetType::SUMMARY); } } diff --git a/opm/simulators/flow/FIPContainer.hpp b/opm/simulators/flow/FIPContainer.hpp index b1b1251b6d0..c45f7449e80 100644 --- a/opm/simulators/flow/FIPContainer.hpp +++ b/opm/simulators/flow/FIPContainer.hpp @@ -54,7 +54,12 @@ class FIPContainer { static constexpr auto oilPhaseIdx = FluidSystem::oilPhaseIdx; static constexpr auto waterPhaseIdx = FluidSystem::waterPhaseIdx; + //! \brief Allocate the per-DOF buffers. + //! \param bufferSize Number of degrees of freedom, auxiliary ones included. + //! \param gridSize Number of grid cells; the restart arrays stop there, since a + //! degree of freedom outside the grid has no cell to be written against. bool allocate(const std::size_t bufferSize, + const std::size_t gridSize, const SummaryConfig& summaryConfig, const bool forceAlloc, std::map& rstKeywords); @@ -161,6 +166,8 @@ class FIPContainer { return this->noPrefix || this->surface || this->reservoir; } } outputRestart_{}; + + std::size_t gridSize_{}; }; } // namespace Opm diff --git a/opm/simulators/flow/FlowAuxCellModule.hpp b/opm/simulators/flow/FlowAuxCellModule.hpp new file mode 100644 index 00000000000..e32df7b93c9 --- /dev/null +++ b/opm/simulators/flow/FlowAuxCellModule.hpp @@ -0,0 +1,200 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::FlowAuxCellModule + */ +#ifndef OPM_FLOW_AUX_CELL_MODULE_HPP +#define OPM_FLOW_AUX_CELL_MODULE_HPP + +#include + +#include +#include + +namespace Opm { + +/*! + * \brief Base class for auxiliary modules whose degrees of freedom are *cells*. + * + * An auxiliary cell satisfies the same conservation equations as a grid cell, and is + * assembled by the same local residual. What it does not have is geometry: its pore + * volume, depth and regions are authored rather than derived from a grid entity, and so + * is its connection list. A numerical aquifer is the simplest example -- it is pure + * bookkeeping, a volume and a list of connections -- and fracture flow cells are the + * dynamic one, where the same quantities are recomputed as the aperture changes. + * + * This is deliberately not the shape of the well models, which are also auxiliary + * modules. Their unknowns are of a different kind and they assemble and scale their own + * equations in linearize(); they leave carriesModelEquations() false. A module derived + * from this class declares the opposite, and in exchange has to supply, per degree of + * freedom, everything the model would otherwise read off the grid. + * + * Indices passed to the accessors below are *module-local*, in [0, numDofs()); + * localToGlobalDof() converts to the model's numbering. The connections reported by + * connections() use global numbering, because they may name grid cells. + */ +template +class FlowAuxCellModule : public BaseAuxiliaryModule +{ + using Scalar = GetPropType; + using ParentType = BaseAuxiliaryModule; + +protected: + using NeighborSet = typename ParentType::NeighborSet; + +public: + /*! + * \brief A flux connection authored by this module. + * + * Either endpoint may be a grid cell or an auxiliary cell; both are global degree of + * freedom indices. The transmissibility is the whole of the geometry -- there is no + * face area to multiply by. + */ + struct Connection + { + unsigned dof1{}; + unsigned dof2{}; + Scalar trans{}; + + //! Thermal half transmissibilities, dof1->dof2 and dof2->dof1. Only consulted + //! when the model solves an energy equation. + Scalar thermalHalfTrans12{}; + Scalar thermalHalfTrans21{}; + }; + + //! The degrees of freedom of an auxiliary cell module are cells. + bool carriesModelEquations() const override + { return true; } + + //! A cell's volume is its bulk volume. + Scalar dofVolume(unsigned localIdx) const override + { return this->bulkVolume(localIdx); } + + /*! + * \brief The connections of this module, each reported exactly once. + * + * The discretization enters every connection in the neighbour list of *both* of its + * endpoints and assembles it from both sides, so reporting one per pair is correct; + * reporting it twice would double the flux. + */ + virtual void connections(std::vector& conns) const = 0; + + /*! + * \brief Pore volume of an auxiliary cell. + */ + virtual Scalar poreVolume(unsigned localIdx) const = 0; + + /*! + * \brief Bulk volume of an auxiliary cell. + * + * The model divides the pore volume by this to obtain a porosity, so the two have to + * be authored consistently. A module with no meaningful bulk volume should report + * the pore volume and accept a porosity of one. + */ + virtual Scalar bulkVolume(unsigned localIdx) const = 0; + + /*! + * \brief Datum depth of an auxiliary cell, used for the gravity head. + */ + virtual Scalar depth(unsigned localIdx) const = 0; + + //! Zero-based PVT region of an auxiliary cell. + virtual unsigned pvtRegionIndex(unsigned localIdx) const = 0; + + //! Zero-based saturation-function region of an auxiliary cell. + virtual unsigned satRegionIndex(unsigned localIdx) const = 0; + + /*! + * \brief The cell of the input grid whose field properties describe this one. + * + * A numerical aquifer names a cell, and takes that cell's region numbers wherever the + * AQUNUM record does not override them -- so the reporting regions it belongs to are + * that cell's, whether or not the cell is part of the simulation grid. Auxiliary + * cells with no such cell to name return -1 and fall back to their initialisation + * partner's regions. + */ + virtual int hostCartesianIndex(unsigned /*localIdx*/) const + { return -1; } + + /*! + * \brief A grid cell whose initial state this auxiliary cell is derived from. + * + * Auxiliary cells cannot be equilibrated by the ordinary machinery, which needs the + * cell's geometry. Reporting a partner lets the problem start from that cell's + * initial fluid state; a module which sets its own initial state in applyInitial() + * may ignore this. + */ + virtual unsigned initialisationPartner(unsigned localIdx) const = 0; + + /*! + * \brief Whether this auxiliary cell currently takes part in the flow problem. + * + * A module may preallocate degrees of freedom it does not use yet -- fracture cells + * that have not opened. Such a cell has no volume and no connections, so its row + * would be empty; the module is responsible for conditioning it in linearize(). + */ + virtual bool isActive(unsigned /*localIdx*/) const + { return true; } + + /*! + * \brief Report the sparsity contribution of this module's connections. + * + * Both directions, plus the diagonal, so that a cell with no connections still has a + * block to be conditioned in. + */ + void addNeighbors(std::vector& neighbors) const override + { + std::vector conns; + this->connections(conns); + + for (const auto& conn : conns) { + neighbors[conn.dof1].insert(conn.dof2); + neighbors[conn.dof2].insert(conn.dof1); + } + + for (unsigned localIdx = 0; localIdx < this->numDofs(); ++localIdx) { + const auto globalIdx = static_cast(this->localToGlobalDof(localIdx)); + neighbors[globalIdx].insert(globalIdx); + } + } + + /*! + * \brief Hand the same connections to the discretization, which builds the + * neighbour info from them. + */ + void addConnections(std::vector& conns) const override + { + std::vector own; + this->connections(own); + + conns.reserve(conns.size() + own.size()); + for (const auto& conn : own) { + conns.push_back({conn.dof1, conn.dof2}); + } + } +}; + +} // namespace Opm + +#endif // OPM_FLOW_AUX_CELL_MODULE_HPP diff --git a/opm/simulators/flow/FlowGenericProblem_impl.hpp b/opm/simulators/flow/FlowGenericProblem_impl.hpp index e15d5e351c6..45dae74a9ef 100644 --- a/opm/simulators/flow/FlowGenericProblem_impl.hpp +++ b/opm/simulators/flow/FlowGenericProblem_impl.hpp @@ -211,7 +211,10 @@ readRockParameters_(const std::vector& cellCenterDepths, for (std::size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) { unsigned tableIdx = 0; if (!rockTableIdx_.empty()) { - tableIdx = rockTableIdx_[elemIdx]; + // Auxiliary DOFs have no entry here; they take the first ROCK region. + if (elemIdx < rockTableIdx_.size()) { + tableIdx = rockTableIdx_[elemIdx]; + } } overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true); @@ -327,7 +330,10 @@ rockCompressibility(unsigned globalSpaceIdx) const unsigned tableIdx = 0; if (!this->rockTableIdx_.empty()) { - tableIdx = this->rockTableIdx_[globalSpaceIdx]; + // Auxiliary DOFs have no entry here; they take the first ROCK region. + if (globalSpaceIdx < this->rockTableIdx_.size()) { + tableIdx = this->rockTableIdx_[globalSpaceIdx]; + } } return this->rockParams_[tableIdx].compressibility; } diff --git a/opm/simulators/flow/FlowGenericVanguard.cpp b/opm/simulators/flow/FlowGenericVanguard.cpp index 0a4b308f846..3ad176b3131 100644 --- a/opm/simulators/flow/FlowGenericVanguard.cpp +++ b/opm/simulators/flow/FlowGenericVanguard.cpp @@ -455,6 +455,12 @@ void FlowGenericVanguard::registerParameters_() "enable restart of OPM simulators from these files"); Parameters::Register ("List of Eclipse keywords which should be ignored. As a ':' separated string."); + Parameters::Register + ("How numerical aquifers are represented. Available options are " + "grid (each aquifer cell takes over the grid cell its AQUNUM record names) and " + "aux (the aquifers are kept out of the grid and represented as degrees of " + "freedom of their own, which leaves the grid as the deck describes it and " + "generates no non-neighbour connections)"); Parameters::Register ("Set strictness of parsing process. Available options are " "normal (stop for critical errors), " diff --git a/opm/simulators/flow/FlowGenericVanguard.hpp b/opm/simulators/flow/FlowGenericVanguard.hpp index 3c13587e6ce..74264160233 100644 --- a/opm/simulators/flow/FlowGenericVanguard.hpp +++ b/opm/simulators/flow/FlowGenericVanguard.hpp @@ -69,6 +69,10 @@ struct NumJacobiBlocks { static constexpr int value = 0; }; #endif // HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA struct OwnerCellsFirst { static constexpr bool value = true; }; + +/// How numerical aquifers (AQUNUM/AQUCON) are represented: "grid" lets each aquifer cell +/// take over the grid cell it names, "aux" keeps them out of the grid entirely. +struct NumericalAquiferMode { static constexpr auto value = "grid"; }; struct ParsingStrictness { static constexpr auto value = "normal"; }; struct ActionParsingStrictness { static constexpr auto value = "normal"; }; diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index 87a1ea6fc7d..42aeeffe37c 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -58,6 +58,7 @@ #include #include // TODO: maybe we can name it FlowProblemProperties.hpp +#include #include #include #include @@ -730,6 +731,12 @@ class FlowProblem : public GetPropType */ Scalar dofCenterDepth(unsigned globalSpaceIdx) const { + // An auxiliary cell has no grid entity to take a depth from; it states its own, + // and that is what the gravity head between it and its neighbours is built from. + if (globalSpaceIdx >= this->model().numGridDof()) { + return this->auxCellDepth_(globalSpaceIdx); + } + return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx); } @@ -798,7 +805,10 @@ class FlowProblem : public GetPropType unsigned tableIdx = 0; if (!this->rockTableIdx_.empty()) { - tableIdx = this->rockTableIdx_[globalSpaceIdx]; + // Auxiliary DOFs have no entry here; they take the first ROCK region. + if (globalSpaceIdx < this->rockTableIdx_.size()) { + tableIdx = this->rockTableIdx_[globalSpaceIdx]; + } } return this->rockParams_[tableIdx].referencePressure; } @@ -817,12 +827,12 @@ class FlowProblem : public GetPropType const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const { - return materialLawManager_->materialLawParams(globalDofIdx); + return materialLawManager_->materialLawParams(auxCellSaturationProxy_(globalDofIdx)); } const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const { - return materialLawManager_->materialLawParams(globalDofIdx, facedir); + return materialLawManager_->materialLawParams(auxCellSaturationProxy_(globalDofIdx), facedir); } /*! @@ -975,13 +985,13 @@ class FlowProblem : public GetPropType solidEnergyLawParams(unsigned globalSpaceIdx, unsigned /*timeIdx*/) const { - return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx); + return this->thermalLawManager_->solidEnergyLawParams(auxCellSaturationProxy_(globalSpaceIdx)); } const ThermalConductionLawParams & thermalConductionLawParams(unsigned globalSpaceIdx, unsigned /*timeIdx*/)const { - return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx); + return this->thermalLawManager_->thermalConductionLawParams(auxCellSaturationProxy_(globalSpaceIdx)); } /*! @@ -1151,7 +1161,10 @@ class FlowProblem : public GetPropType unsigned tableIdx = 0; if (!this->rockTableIdx_.empty()) - tableIdx = this->rockTableIdx_[elementIdx]; + // Auxiliary DOFs have no entry here; they take the first ROCK region. + if (elementIdx < this->rockTableIdx_.size()) { + tableIdx = this->rockTableIdx_[elementIdx]; + } const auto& fs = intQuants.fluidState(); LhsEval effectivePressure = decay(fs.pressure(refPressurePhaseIdx_())); @@ -1291,6 +1304,15 @@ class FlowProblem : public GetPropType return drift_; } + /*! + * \brief The auxiliary cell modules this problem owns. + * + * Exposed so that the parts of the simulator which have to know what lives outside + * the grid -- reporting, chiefly -- can find them without a second registry. + */ + const std::vector>>& auxCellModules() const + { return auxCellModules_; } + private: Implementation& asImp_() { return *static_cast(this); } @@ -1457,6 +1479,7 @@ class FlowProblem : public GetPropType this->updatePlmixnum_(); OPM_END_PARALLEL_TRY_CATCH("Invalid region numbers: ", vanguard.gridView().comm()); + this->authorAuxCellRegions_(); //////////////////////////////// // porosity updateReferencePorosity_(); @@ -1473,6 +1496,12 @@ class FlowProblem : public GetPropType // fluid-matrix interactions (saturation functions; relperm/capillary pressure) materialLawManager_ = std::make_shared(); materialLawManager_->initFromState(eclState); + // NOTE: sized by the grid degrees of freedom on purpose. Extending this over the + // auxiliary cells does not work: initParamsForElements walks the elements into the + // endpoint-scaling and hysteresis machinery, which is keyed on grid entities + // throughout, so raising the count merely pushes that machinery off the end of its + // own arrays. Auxiliary cells are given saturation-function parameters by + // redirecting materialLawParams() instead. materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(), this-> template fieldPropIntTypeOnLeafAssigner_(), this-> lookupIdxOnLevelZeroAssigner_()); @@ -1504,7 +1533,10 @@ class FlowProblem : public GetPropType std::size_t numDof = this->model().numGridDof(); - this->referencePorosity_[/*timeIdx=*/0].resize(numDof); + // Auxiliary cells have their own authored pore and bulk volume; sizing for them + // here keeps referencePorosity() answerable for every degree of freedom. + this->referencePorosity_[/*timeIdx=*/0].resize(this->model().numTotalDof()); + this->authorAuxCellPorosity_(); const auto& fp = eclState.fieldProps(); const std::vector porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV"); @@ -1533,7 +1565,11 @@ class FlowProblem : public GetPropType const auto& eclState = vanguard.eclState(); std::size_t numDof = this->model().numGridDof(); - this->rockFraction_[/*timeIdx=*/0].resize(numDof); + this->rockFraction_[/*timeIdx=*/0].resize(this->model().numTotalDof(), 0.0); + + // An auxiliary cell is a bookkeeping volume, not rock: it stores no heat of its + // own. Leaving its rock fraction at zero is what expresses that, since the + // energy storage term is rockFraction times the rock's internal energy. // For the energy equation, we need the volume of the rock. // The volume of the rock is computed by rockFraction * geometric volume of the element. // The reference porosity is defined as porosity * ntg * pore-volume-multiplier. @@ -1617,6 +1653,181 @@ class FlowProblem : public GetPropType } protected: + /*! + * \brief Take ownership of an auxiliary cell module and register it with the model. + * + * Registration has to happen before the model sizes its per-DOF containers, which is + * why this is only ever called from registerAuxiliaryCellModules(). The module's + * connections are built afterwards, because they are expressed in global degree of + * freedom indices and the offset is only assigned on registration. + */ + template + Module& registerAuxCellModule_(std::unique_ptr module) + { + // Which rank owns a degree of freedom that has no geometry is a choice, and one + // that has to be made where the partition is made: every rank holding a cell the + // auxiliary cell connects to needs it at least as a copy, with its intensive + // quantities refreshed, and the owner has to be the only one contributing its + // accumulation. None of that is in place yet -- the auxiliary degrees of freedom + // are not in the communication index set -- so refuse the combination rather than + // let each rank solve its own version of the aquifer. + if (this->simulator().gridView().comm().size() > 1) { + OPM_THROW(std::runtime_error, + "Degrees of freedom outside the grid (numerical aquifers " + "represented as auxiliary cells) are not supported in parallel " + "yet. Run on one process, or use --numerical-aquifer-mode=grid."); + } + + auto& ref = *module; + this->simulator().model().addAuxiliaryModule(&ref); + this->auxCellModules_.push_back(std::move(module)); + ref.buildConnections(); + return ref; + } + + /*! + * \brief Map an auxiliary cell onto a grid cell for the per-element rock tables. + * + * The saturation-function machinery is built per grid element and stays that way: + * endpoint scaling and hysteresis are keyed on grid entities throughout, so the + * parameter tables cannot simply be extended over degrees of freedom that have no + * entity. Instead an auxiliary cell borrows the parameters of a grid cell -- the one + * it is connected to and initialised from. + * + * That is exact when the two share a saturation region, which is the ordinary case + * (and the only one an aquifer is likely to present, being water-filled). Where they + * differ, the borrowed curves are the wrong ones; giving auxiliary cells parameters of + * their own region is the follow-up, and needs the saturation-function tables to be + * addressable by region rather than by element. + */ + unsigned auxCellSaturationProxy_(unsigned globalDofIdx) const + { + if (globalDofIdx < this->model().numGridDof()) { + return globalDofIdx; + } + + for (const auto& module : this->auxCellModules_) { + const auto begin = static_cast(module->dofOffset()); + if ((globalDofIdx >= begin) && (globalDofIdx < begin + module->numDofs())) { + return module->initialisationPartner(globalDofIdx - begin); + } + } + + return 0; + } + + //! Depth of an auxiliary cell, by global degree of freedom index. + Scalar auxCellDepth_(unsigned globalSpaceIdx) const + { + for (const auto& module : this->auxCellModules_) { + const auto begin = static_cast(module->dofOffset()); + if ((globalSpaceIdx >= begin) && (globalSpaceIdx < begin + module->numDofs())) { + return module->depth(globalSpaceIdx - begin); + } + } + + return 0.0; + } + + /*! + * \brief Give the auxiliary cells their region numbers. + * + * The region arrays are built from the field properties over the grid, which leaves + * them one entry short of every degree of freedom the moment an auxiliary module adds + * any. They are read by degree-of-freedom index -- pvtRegionIndex() is asked for one + * on the way to the PVT tables -- so a short array is an out-of-bounds read that ends + * up indexing the tables with whatever was next in memory. An empty array means "one + * region", which needs no extending. + */ + void authorAuxCellRegions_() + { + const auto numTotalDof = this->model().numTotalDof(); + if (numTotalDof == this->model().numGridDof()) { + return; + } + + const auto extend = [numTotalDof](auto& numbers, auto&& regionOf) { + if (numbers.empty()) { + return; + } + + numbers.resize(numTotalDof, 0); + regionOf(numbers); + }; + + extend(this->pvtnum_, [this](auto& numbers) { + for (const auto& module : this->auxCellModules_) { + for (unsigned localIdx = 0; localIdx < module->numDofs(); ++localIdx) { + numbers[module->localToGlobalDof(localIdx)] = module->pvtRegionIndex(localIdx); + } + } + }); + + extend(this->satnum_, [this](auto& numbers) { + for (const auto& module : this->auxCellModules_) { + for (unsigned localIdx = 0; localIdx < module->numDofs(); ++localIdx) { + numbers[module->localToGlobalDof(localIdx)] = module->satRegionIndex(localIdx); + } + } + }); + + // The solvent and polymer models have no auxiliary-cell story yet; region zero is + // the only defensible placeholder, and it is what an absent array would give. + extend(this->miscnum_, [](auto&) {}); + extend(this->plmixnum_, [](auto&) {}); + } + + /*! + * \brief Fill in the reference porosity of the auxiliary cells. + * + * The model reads a porosity and multiplies it by the degree of freedom's total + * volume to recover a pore volume, so the two have to be authored consistently with + * the volume the module reports. + */ + void authorAuxCellPorosity_() + { + for (const auto& module : this->auxCellModules_) { + for (unsigned localIdx = 0; localIdx < module->numDofs(); ++localIdx) { + const auto globalIdx = static_cast(module->localToGlobalDof(localIdx)); + const auto bulkVolume = module->bulkVolume(localIdx); + + this->referencePorosity_[/*timeIdx=*/0][globalIdx] = (bulkVolume > 0.0) + ? module->poreVolume(localIdx) / bulkVolume + : 0.0; + } + } + } + + /*! + * \brief Publish the auxiliary modules' connection transmissibilities. + * + * The transmissibility store is keyed on the degree of freedom pair and does not care + * whether a connection came from a face, so an authored connection simply goes in + * alongside the geometric ones and problem.transmissibility() answers for it. Called + * once the grid's own transmissibilities are final; a module whose values change with + * the solution refreshes them the same way. + */ + void applyAuxCellTransmissibilities_() + { + using ConnectionVector = std::vector::Connection>; + + for (const auto& module : this->auxCellModules_) { + ConnectionVector conns; + module->connections(conns); + + for (const auto& conn : conns) { + this->transmissibilities_.setTransmissibility(conn.dof1, conn.dof2, conn.trans); + + if constexpr (enableFullyImplicitThermal) { + this->transmissibilities_ + .setThermalHalfTrans(conn.dof1, conn.dof2, conn.thermalHalfTrans12); + this->transmissibilities_ + .setThermalHalfTrans(conn.dof2, conn.dof1, conn.thermalHalfTrans21); + } + } + } + } + struct PffDofData_ { ConditionalStorage thermalHalfTransIn; @@ -1800,7 +2011,10 @@ class FlowProblem : public GetPropType unsigned tableIdx = 0; if (!this->rockTableIdx_.empty()) - tableIdx = this->rockTableIdx_[elementIdx]; + // Auxiliary DOFs have no entry here; they take the first ROCK region. + if (elementIdx < this->rockTableIdx_.size()) { + tableIdx = this->rockTableIdx_[elementIdx]; + } const auto& fs = intQuants.fluidState(); LhsEval effectivePressure = obtain(fs.pressure(refPressurePhaseIdx_())); @@ -1831,6 +2045,11 @@ class FlowProblem : public GetPropType typename Vanguard::TransmissibilityType transmissibilities_; + //! Auxiliary modules whose degrees of freedom are cells (numerical aquifers, and + //! later fracture flow cells). Owned here because their authored data feeds the + //! problem's own per-DOF tables. + std::vector>> auxCellModules_; + std::shared_ptr materialLawManager_; std::shared_ptr thermalLawManager_; diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp index bb5a4348a2a..e85383434e6 100644 --- a/opm/simulators/flow/FlowProblemBlackoil.hpp +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -49,6 +49,7 @@ #include #include +#include #include #include #include @@ -321,6 +322,29 @@ class FlowProblemBlackoil : public FlowProblem /*! * \copydoc FvBaseProblem::finishInit */ + /*! + * \brief Create the auxiliary cell modules this deck asks for. + * + * Called by the model before it sizes anything, which is the only point at which a + * module that introduces degrees of freedom may still be registered. + */ + void registerAuxiliaryCellModules() + { + const auto& eclState = this->simulator().vanguard().eclState(); + + if ((eclState.numericalAquiferMode() == NumericalAquiferMode::AuxiliaryCells) && + eclState.aquifer().hasNumericalAquifer()) + { + const auto& aquifers = this->registerAuxCellModule_( + std::make_unique>(this->simulator())); + + if (this->simulator().gridView().comm().rank() == 0) { + OpmLog::info(fmt::format("Numerical aquifers represented as {} auxiliary " + "cells outside the grid", aquifers.numDofs())); + } + } + } + void finishInit() { // TODO: there should be room to remove duplication for this @@ -404,7 +428,9 @@ class FlowProblemBlackoil : public FlowProblem if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0); + // By total degree of freedom count: maxOilSaturation() is asked by + // degree-of-freedom index, and an auxiliary cell asks the same way. + this->maxOilSaturation_.resize(this->model().numTotalDof(), 0.0); } this->readRockParameters_(simulator.vanguard().cellCenterDepths(), @@ -427,7 +453,21 @@ class FlowProblemBlackoil : public FlowProblem finishTransmissibilities(); + // The auxiliary cells' connections are authored, not geometric, so they are + // published once the grid's own transmissibilities are final. + this->applyAuxCellTransmissibilities_(); + const auto& initconfig = eclState.getInitConfig(); + + // Auxiliary cells do not appear in the restart file: it is written per grid + // element, so their state would come back undefined rather than merely stale. + if (initconfig.restartRequested() && !this->auxCellModules_.empty()) { + OPM_THROW(std::runtime_error, + "Restart is not supported together with auxiliary cells " + "(numerical aquifers represented outside the grid): their state " + "is not written to the restart file."); + } + this->tracerModel_.init(initconfig.restartRequested()); if (initconfig.restartRequested()) { this->readEclRestartSolution_(); @@ -453,7 +493,8 @@ class FlowProblemBlackoil : public FlowProblem this->computeAndSetEqWeights_(); if (this->enableDriftCompensation_ || this->enableDriftCompensationTemp_) { - this->drift_.resize(this->model().numGridDof()); + // Sized like the residual it compensates, which spans the auxiliary DOFs. + this->drift_.resize(this->model().numTotalDof()); this->drift_ = 0.0; } @@ -487,7 +528,10 @@ class FlowProblemBlackoil : public FlowProblem // TODO: move to the end for later refactoring of the function finishInit() // // deal with DRSDT - this->mixControls_.init(this->model().numGridDof(), + // Sized over every degree of freedom, auxiliary ones included: the intensive + // quantities ask for the dissolution limits by degree-of-freedom index, and an + // auxiliary cell reaches this the same way a grid cell does. + this->mixControls_.init(this->model().numTotalDof(), this->episodeIndex(), eclState.runspec().tabdims().getNumPVTTables()); @@ -1175,7 +1219,7 @@ class FlowProblemBlackoil : public FlowProblem if (fip_init) { this->updateReferencePorosity_(); - this->mixControls_.init(this->model().numGridDof(), + this->mixControls_.init(this->model().numTotalDof(), this->episodeIndex(), eclState.runspec().tabdims().getNumPVTTables()); } diff --git a/opm/simulators/flow/GenericOutputModule.cpp b/opm/simulators/flow/GenericOutputModule.cpp index b249e1074a6..6d9232fb4a5 100644 --- a/opm/simulators/flow/GenericOutputModule.cpp +++ b/opm/simulators/flow/GenericOutputModule.cpp @@ -641,8 +641,14 @@ doAllocBuffers(const unsigned bufferSize, const bool isRestart, const EclHysteresisConfig* hysteresisConfig, const unsigned numOutputNnc, - std::map rstKeywords) + std::map rstKeywords, + const unsigned auxDofCount) { + // The fluid-in-place buffers span the auxiliary degrees of freedom as well, so that + // what lives outside the grid still counts towards the field and region totals; every + // other buffer here is written per grid cell and stays that size. + const unsigned fipBufferSize = bufferSize + auxDofCount; + if (rstKeywords.empty()) { rstKeywords = schedule_.rst_keywords(reportStepNum); } @@ -668,7 +674,8 @@ doAllocBuffers(const unsigned bufferSize, rstKeywords["PRES"] = 0; // Fluid in place - this->computeFip_ = this->fipC_.allocate(bufferSize, + this->computeFip_ = this->fipC_.allocate(fipBufferSize, + bufferSize, summaryConfig_, !substep, rstKeywords); @@ -687,15 +694,15 @@ doAllocBuffers(const unsigned bufferSize, if (needPoreVolume) { this->fipC_.add(Inplace::Phase::PoreVolume); this->fipC_.add(Inplace::Phase::DynamicPoreVolume); - this->hydrocarbonPoreVolume_.resize(bufferSize, 0.0); + this->hydrocarbonPoreVolume_.resize(fipBufferSize, 0.0); } else { this->hydrocarbonPoreVolume_.clear(); } if (needAvgPress) { - this->pressureTimesPoreVolume_.resize(bufferSize, 0.0); - this->pressureTimesHydrocarbonVolume_.resize(bufferSize, 0.0); + this->pressureTimesPoreVolume_.resize(fipBufferSize, 0.0); + this->pressureTimesHydrocarbonVolume_.resize(fipBufferSize, 0.0); } else { this->pressureTimesPoreVolume_.clear(); @@ -1048,6 +1055,36 @@ update(Inplace& inplace, inplace.add(phase, sum); } +template +void GenericOutputModule:: +extendRegionsForAuxiliaryDofs(const std::vector& hostCartesianIndex) +{ + if (hostCartesianIndex.empty()) { + return; + } + + const auto& fp = this->eclState_.fieldProps(); + + for (auto& [name, region] : this->regions_) { + if (region.empty()) { + continue; + } + + // The global array still describes the cell an AQUNUM record names even when that + // cell has been left out of the simulation grid to make room for the aquifer. + const auto& global = fp.get_global_int(name); + + const auto firstAux = region.size(); + region.resize(firstAux + hostCartesianIndex.size(), 0); + for (std::size_t i = 0; i < hostCartesianIndex.size(); ++i) { + const auto host = hostCartesianIndex[i]; + region[firstAux + i] = (host < 0) + ? 0 // no cell to take a region from: no region + : global[static_cast(host)]; + } + } +} + template void GenericOutputModule:: makeRegionSum(Inplace& inplace, diff --git a/opm/simulators/flow/GenericOutputModule.hpp b/opm/simulators/flow/GenericOutputModule.hpp index 8815096b36b..4f13369d749 100644 --- a/opm/simulators/flow/GenericOutputModule.hpp +++ b/opm/simulators/flow/GenericOutputModule.hpp @@ -380,7 +380,18 @@ class GenericOutputModule { const bool isRestart, const EclHysteresisConfig* hysteresisConfig, unsigned numOutputNnc = 0, - std::map rstKeywords = {}); + std::map rstKeywords = {}, + unsigned auxDofCount = 0); + + /*! + * \brief Extend the reporting regions over the auxiliary degrees of freedom. + * + * The region arrays come from the field properties over the grid. An auxiliary cell + * belongs to a region all the same -- a numerical aquifer takes the regions of the + * cell its AQUNUM record names -- and it has to, or its contribution to the field and + * region totals has nowhere to go. + */ + void extendRegionsForAuxiliaryDofs(const std::vector& hostCartesianIndex); /// Allocate the buffers a derived module owns. Called while the restart /// keywords are being handled, so that a keyword consumed here is marked diff --git a/opm/simulators/flow/Main.cpp b/opm/simulators/flow/Main.cpp index 95b80928ee3..b6ca3d26e14 100644 --- a/opm/simulators/flow/Main.cpp +++ b/opm/simulators/flow/Main.cpp @@ -330,6 +330,7 @@ void Main::readDeck(const std::string& deckFilename, const std::size_t numThreads, const int output_param, const bool slaveMode, + const std::string& numericalAquiferMode, const std::string& parameters, std::string_view moduleVersion, std::string_view compileTimestamp) @@ -350,6 +351,19 @@ void Main::readDeck(const std::string& deckFilename, if (output_param >= 0) outputInterval = output_param; + // How the numerical aquifers of the deck are to be represented. This has to be + // decided before the EclipseState is built, because taking over a grid cell reshapes + // the grid and the field properties as the state is constructed. + auto aquiferMode = NumericalAquiferMode::GridCells; + if (numericalAquiferMode == "aux") { + aquiferMode = NumericalAquiferMode::AuxiliaryCells; + } + else if (numericalAquiferMode != "grid") { + OPM_THROW(std::runtime_error, + "Unknown numerical aquifer mode '" + numericalAquiferMode + + "'; valid choices are 'grid' and 'aux'"); + } + Opm::readDeck(FlowGenericVanguard::comm(), deckFilename, eclipseState_, @@ -366,7 +380,8 @@ void Main::readDeck(const std::string& deckFilename, outputCout_, keepKeywords, outputInterval, - slaveMode); + slaveMode, + aquiferMode); verifyValidCellGeometry(FlowGenericVanguard::comm(), *this->eclipseState_); diff --git a/opm/simulators/flow/Main.hpp b/opm/simulators/flow/Main.hpp index 483bd7ff9b8..a3dbd08e2c5 100644 --- a/opm/simulators/flow/Main.hpp +++ b/opm/simulators/flow/Main.hpp @@ -327,6 +327,7 @@ class Main getNumThreads(), Parameters::Get(), Parameters::Get(), + Parameters::Get(), cmdline_params, Opm::moduleVersion(), Opm::compileTimestamp()); @@ -555,6 +556,7 @@ class Main const std::size_t numThreads, const int output_param, const bool slaveMode, + const std::string& numericalAquiferMode, const std::string& parameters, std::string_view moduleVersion, std::string_view compileTimestamp); diff --git a/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp b/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp index e39e61b4df0..086bac62d46 100644 --- a/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp +++ b/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp @@ -252,7 +252,10 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer, perfTimer.start(); report.total_newton_iterations = 1; - const unsigned nc = this->simulator_.model().numGridDof(); + // The Jacobian and the residual are sized for the total number of DOFs, i.e. + // they include the rows contributed by auxiliary modules. The solution vector + // handed to the linear solver has to match. + const unsigned nc = this->simulator_.model().numTotalDof(); BVector x(nc); linear_solve_setup_time_ = 0.0; @@ -616,6 +619,28 @@ localConvergenceData(std::vector& R_sum, OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::localConvergenceData() failed: ", this->grid_.comm()); + // Auxiliary cells carry the same equations but are not reachable through the grid, + // so the loop above never sees them. Their residual has to enter the convergence + // measures like any other cell's, otherwise the Newton iteration would be declared + // converged while their mass balance is still violated. + const unsigned numGridDof = model.numGridDof(); + const unsigned numTotalDof = model.numTotalDof(); + for (unsigned cell_idx = numGridDof; cell_idx < numTotalDof; ++cell_idx) { + if (!model.dofCarriesModelEquations(cell_idx)) { + continue; + } + + const auto& intQuants = model.intensiveQuantities(cell_idx, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) * + model.dofTotalVolume(cell_idx); + pvSumLocal += pvValue; + + this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue, + B_avg, R_sum, maxCoeff, maxCoeffCell); + } + // compute local average in terms of global number of elements const int bSize = B_avg.size(); for (int i = 0; i < bSize; ++i) { diff --git a/opm/simulators/flow/NonlinearSystem_impl.hpp b/opm/simulators/flow/NonlinearSystem_impl.hpp index 2c202c99503..9f2b7fa8faa 100644 --- a/opm/simulators/flow/NonlinearSystem_impl.hpp +++ b/opm/simulators/flow/NonlinearSystem_impl.hpp @@ -144,7 +144,9 @@ NonlinearSystem(Simulator& simulator, , param_(param) , well_model_(wellModel) , current_relaxation_(1.0) - , dx_old_(simulator_.model().numGridDof()) + // sized like the linear-system solution vector it is compared against in + // stabilizeNonlinearUpdate(), which spans the auxiliary DOFs as well + , dx_old_(simulator_.model().numTotalDof()) {} template diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 93b7e251643..1c9593a5d25 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -254,6 +254,8 @@ class OutputBlackOilModule : public GenericOutputModulesimulator_.problem(); + const auto& model = this->simulator_.model(); + const auto auxDofCount = model.numTotalDof() - model.numGridDof(); this->doAllocBuffers(bufferSize, reportStepNum, @@ -261,7 +263,33 @@ class OutputBlackOilModule : public GenericOutputModulehysteresisConfig(), - problem.eclWriter().getOutputNnc().front().size()); + problem.eclWriter().getOutputNnc().front().size(), + /*rstKeywords=*/{}, + auxDofCount); + + if (auxDofCount > 0) { + this->extendRegionsForAuxiliaryDofs(this->auxCellHostCartesianIndices_()); + } + } + + /*! + * \brief The input-grid cell each auxiliary degree of freedom takes its regions from. + * + * In auxiliary-DOF order, so that it lines up with the tail of the per-DOF buffers. + */ + std::vector auxCellHostCartesianIndices_() const + { + const auto& model = this->simulator_.model(); + std::vector hosts(model.numTotalDof() - model.numGridDof(), -1); + + for (const auto& module : this->simulator_.problem().auxCellModules()) { + for (unsigned localIdx = 0; localIdx < module->numDofs(); ++localIdx) { + const auto globalIdx = static_cast(module->localToGlobalDof(localIdx)); + hosts[globalIdx - model.numGridDof()] = module->hostCartesianIndex(localIdx); + } + } + + return hosts; } //! \brief Setup list of active element-level data extractors diff --git a/opm/simulators/flow/Transmissibility.hpp b/opm/simulators/flow/Transmissibility.hpp index 9bc3d6c62bc..81a08f7854a 100644 --- a/opm/simulators/flow/Transmissibility.hpp +++ b/opm/simulators/flow/Transmissibility.hpp @@ -79,6 +79,25 @@ class Transmissibility { */ Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const; + /*! + * \brief Set the transmissibility of a connection which is not a geometric face. + * + * This is for connections authored by something other than the grid -- an auxiliary + * cell's connection list, for instance -- where the transmissibility is computed by + * the author rather than from face geometry. It may be called again to update the + * value of an existing connection, which is how a connection whose transmissibility + * depends on the solution (a fracture aperture, say) is refreshed between + * iterations without touching the sparsity pattern. + */ + void setTransmissibility(unsigned elemIdx1, unsigned elemIdx2, Scalar value); + + /*! + * \brief Set the thermal half transmissibility of a non-geometric connection. + * + * Directional: call it once for each ordering of the two degrees of freedom. + */ + void setThermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx, Scalar value); + /*! * \brief Return the transmissibility for a given boundary segment. */ diff --git a/opm/simulators/flow/Transmissibility_impl.hpp b/opm/simulators/flow/Transmissibility_impl.hpp index c268d2eed91..b488bd77339 100644 --- a/opm/simulators/flow/Transmissibility_impl.hpp +++ b/opm/simulators/flow/Transmissibility_impl.hpp @@ -124,6 +124,20 @@ transmissibility(unsigned elemIdx1, unsigned elemIdx2) const return trans_.at(details::isId(elemIdx1, elemIdx2)); } +template +void Transmissibility:: +setTransmissibility(unsigned elemIdx1, unsigned elemIdx2, Scalar value) +{ + trans_[details::isId(elemIdx1, elemIdx2)] = value; +} + +template +void Transmissibility:: +setThermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx, Scalar value) +{ + thermalHalfTrans_[details::directionalIsId(insideElemIdx, outsideElemIdx)] = value; +} + template Scalar Transmissibility:: transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const diff --git a/opm/simulators/linalg/ISTLSolver.cpp b/opm/simulators/linalg/ISTLSolver.cpp index c7d7031a25b..ee8daa04760 100644 --- a/opm/simulators/linalg/ISTLSolver.cpp +++ b/opm/simulators/linalg/ISTLSolver.cpp @@ -127,8 +127,20 @@ void FlexibleSolverInfo::create(const Matrix& matrix, this->solver_ = std::move(sol); } else { using ParOperatorType = WellModelGhostLastMatrixAdapter; + // The rows this rank owns are the interior grid cells and, behind the ghost + // rows, the degrees of freedom that have no grid cell at all. Passed as + // bands rather than as a single count, because the second group sits after + // the ghosts and would otherwise be projected out of the Krylov operator -- + // silently, as a 0 = 0 equation. + auto ownedRowBands = std::vector> + {{0, interiorCellNum_}}; + if (numAuxiliaryDof_ > 0) { + ownedRowBands.emplace_back(matrix.N() - numAuxiliaryDof_, matrix.N()); + } + auto pop = std::make_unique(matrix, *wellOperator_, - interiorCellNum_); + std::move(ownedRowBands), + matrix.N()); using FlexibleSolverType = Dune::FlexibleSolver; auto sol = std::make_unique(*pop, *comm, prm, weightsCalculator, diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 109d52b90cb..85ab3ea0305 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -118,6 +118,9 @@ struct FlexibleSolverInfo std::unique_ptr> wellOperator_; AbstractPreconditionerType* pre_ = nullptr; std::size_t interiorCellNum_ = 0; + //! Degrees of freedom appended after the grid rows; owned by this rank by + //! construction, but not part of the interior prefix. + std::size_t numAuxiliaryDof_ = 0; }; @@ -511,6 +514,11 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp); } } + // Only known once the model has been built, which is after this solver is + // constructed -- hence here rather than in initialize(). + flexibleSolver_[activeSolverNum_].numAuxiliaryDof_ = + simulator_.model().numTotalDof() - simulator_.model().numGridDof(); + std::function weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex); OPM_TIMEBLOCK(flexibleSolverCreate); flexibleSolver_[activeSolverNum_].create(getMatrix(), diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp index f11c080b72f..ca204247323 100644 --- a/opm/simulators/linalg/WellOperators.hpp +++ b/opm/simulators/linalg/WellOperators.hpp @@ -33,6 +33,8 @@ #include #include +#include +#include namespace Opm { @@ -315,22 +317,53 @@ class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator& wellOper, const std::size_t interiorSize ) - : A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize) + : WellModelGhostLastMatrixAdapter(A, wellOper, {{0, interiorSize}}, A.N()) + {} + + /*! + * \brief Constructor taking the owned rows as a list of bands. + * + * The rows this rank owns need not be a prefix of the matrix. Degrees of freedom + * that have no grid cell -- a numerical aquifer or a fracture represented outside + * the grid -- are appended after the grid rows, which puts them behind the ghost + * rows; and a grid that gains or loses cells during the run cannot keep its owned + * rows at the front without renumbering. + * + * Describing the owned rows explicitly is what keeps this operator honest in either + * case. With a single band it is exactly the prefix it always was. + * + * \param ownedRowBands Half-open [begin, end) ranges of rows this rank owns, in + * increasing order and not overlapping. + * \param numRows Total number of rows; everything outside the bands is projected out. + */ + WellModelGhostLastMatrixAdapter (const M& A, + const LinearOperatorExtra& wellOper, + std::vector> ownedRowBands, + const std::size_t numRows) + : A_( A ) + , wellOper_( wellOper ) + , ownedRowBands_(std::move(ownedRowBands)) + , numRows_(numRows) + , interiorSize_(ownedRowBands_.empty() ? 0 : ownedRowBands_.front().second) {} void apply(const X& x, Y& y) const override { OPM_TIMEBLOCK(apply); - for (auto row = A_.begin(); row.index() < interiorSize_; ++row) - { - y[row.index()]=0; - auto endc = (*row).end(); - for (auto col = (*row).begin(); col != endc; ++col) - (*col).umv(x[col.index()], y[row.index()]); + for (const auto& [first, last] : ownedRowBands_) { + for (auto row = A_.begin() + first; row.index() < last; ++row) + { + y[row.index()]=0; + auto endc = (*row).end(); + for (auto col = (*row).begin(); col != endc; ++col) + (*col).umv(x[col.index()], y[row.index()]); + } } // add well model modification to y @@ -343,11 +376,13 @@ class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator& wellOper_; + std::vector> ownedRowBands_; + std::size_t numRows_; + //! First band's end, kept for the derived operators which still reason about a prefix. std::size_t interiorSize_; }; diff --git a/opm/simulators/linalg/getQuasiImpesWeights.hpp b/opm/simulators/linalg/getQuasiImpesWeights.hpp index 5c288396701..6714f8bea32 100644 --- a/opm/simulators/linalg/getQuasiImpesWeights.hpp +++ b/opm/simulators/linalg/getQuasiImpesWeights.hpp @@ -61,6 +61,78 @@ namespace Details namespace Amg { + /*! + * \brief The quasi-IMPES weight of one matrix row. + * + * Built from the row's diagonal block alone, so it asks nothing of the grid. + */ + template + VectorBlockType quasiImpesWeightForRow(const Matrix& A, + const int rowIdx, + const int pressureVarIndex, + const bool transpose) + { + using MatrixBlockType = typename Matrix::block_type; + + VectorBlockType rhs(0.0); + rhs[pressureVarIndex] = 1.0; + + MatrixBlockType diag_block(0.0); + const auto row_it = A.begin() + rowIdx; + const auto endj = (*row_it).end(); + for (auto j = (*row_it).begin(); j != endj; ++j) { + if (row_it.index() == j.index()) { + diag_block = (*j); + break; + } + } + + VectorBlockType bweights; + if (transpose) { + diag_block.solve(bweights, rhs); + } else { + MatrixBlockType diag_block_transpose = Details::transposeDenseMatrix(diag_block); + diag_block_transpose.solve(bweights, rhs); + } + + const double abs_max = + *std::ranges::max_element(bweights, + [](double a, double b) + { return std::fabs(a) < std::fabs(b); }); + bweights /= std::fabs(abs_max); + + return bweights; + } + + /*! + * \brief Give the auxiliary degrees of freedom a CPR weight. + * + * The true-IMPES weights are built by walking the grid, so they leave the auxiliary + * degrees of freedom -- which have no element -- untouched. That is not a small + * inaccuracy: the weight multiplies the whole row on its way into the coarse pressure + * system, so a zero weight deletes the row and makes that system singular. + * + * They get the quasi-IMPES weight instead, which needs only the assembled diagonal + * block. Mixing the two is a compromise on scaling, not on correctness; an auxiliary + * cell's true-IMPES weight needs the storage term evaluated without an element + * context, which is the same thing the TPFA linearizer does and is the natural + * follow-up. + */ + template + void getAuxiliaryDofWeights(const Matrix& matrix, + const int firstAuxiliaryRow, + const int pressureVarIndex, + const bool transpose, + Vector& weights) + { + using VectorBlockType = typename Vector::block_type; + + for (int rowIdx = firstAuxiliaryRow; rowIdx < static_cast(matrix.N()); ++rowIdx) { + weights[rowIdx] = quasiImpesWeightForRow(matrix, rowIdx, + pressureVarIndex, transpose); + } + } + template void getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, @@ -254,10 +326,14 @@ namespace Amg } } OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm()); + + getAuxiliaryDofWeights(model.linearizer().jacobian().istlMatrix(), + static_cast(model.numGridDof()), + pressureVarIndex, /*transpose=*/false, weights); } template - void getTrueImpesWeightsAnalytic(int /*pressureVarIndex*/, + void getTrueImpesWeightsAnalytic(int pressureVarIndex, Vector& weights, const ElementContext& elemCtx, const Model& model, @@ -343,6 +419,10 @@ namespace Amg } } OPM_END_PARALLEL_TRY_CATCH("getTrueImpesAnalyticWeights() failed: ", elemCtx.simulator().vanguard().grid().comm()); + + getAuxiliaryDofWeights(model.linearizer().jacobian().istlMatrix(), + static_cast(model.numGridDof()), + pressureVarIndex, /*transpose=*/false, weights); } } // namespace Amg diff --git a/opm/simulators/utils/ParallelEclipseState.cpp b/opm/simulators/utils/ParallelEclipseState.cpp index 95134b7f733..2b64c9e4e08 100644 --- a/opm/simulators/utils/ParallelEclipseState.cpp +++ b/opm/simulators/utils/ParallelEclipseState.cpp @@ -252,14 +252,15 @@ ParallelEclipseState::ParallelEclipseState(Parallel::Communication comm) } -ParallelEclipseState::ParallelEclipseState(const Deck& deck) - : EclipseState(deck) +ParallelEclipseState::ParallelEclipseState(const Deck& deck, NumericalAquiferMode aquiferMode) + : EclipseState(deck, aquiferMode) , m_fieldProps(field_props) { } -ParallelEclipseState::ParallelEclipseState(const Deck& deck, Parallel::Communication comm) - : EclipseState(deck) +ParallelEclipseState::ParallelEclipseState(const Deck& deck, Parallel::Communication comm, + NumericalAquiferMode aquiferMode) + : EclipseState(deck, aquiferMode) , m_fieldProps(field_props, comm) , m_comm(comm) { diff --git a/opm/simulators/utils/ParallelEclipseState.hpp b/opm/simulators/utils/ParallelEclipseState.hpp index f9243ddb885..e8d0817390c 100644 --- a/opm/simulators/utils/ParallelEclipseState.hpp +++ b/opm/simulators/utils/ParallelEclipseState.hpp @@ -160,15 +160,19 @@ class ParallelEclipseState : public EclipseState { //! \brief Construct from a deck instance. //! \param deck The deck to construct from + //! \param aquiferMode How numerical aquifers are represented //! \details Only called on root process - ParallelEclipseState(const Deck& deck); + ParallelEclipseState(const Deck& deck, + NumericalAquiferMode aquiferMode = NumericalAquiferMode::GridCells); //! EXPERIMENTAL FUNCTION TO ADD COMM AS INPUT. //! \brief Construct from a deck instance. //! \param deck The deck to construct from //! \param comm Parallel communicator + //! \param aquiferMode How numerical aquifers are represented //! \details Only called on root process - ParallelEclipseState(const Deck& deck, Parallel::Communication comm); + ParallelEclipseState(const Deck& deck, Parallel::Communication comm, + NumericalAquiferMode aquiferMode = NumericalAquiferMode::GridCells); //! \brief Switch to global field properties. //! \details Called on root process to use the global field properties diff --git a/opm/simulators/utils/readDeck.cpp b/opm/simulators/utils/readDeck.cpp index a8c22b96cf1..37732b3f253 100644 --- a/opm/simulators/utils/readDeck.cpp +++ b/opm/simulators/utils/readDeck.cpp @@ -249,12 +249,13 @@ namespace { std::shared_ptr createEclipseState([[maybe_unused]] Opm::Parallel::Communication comm, - const Opm::Deck& deck) + const Opm::Deck& deck, + const Opm::NumericalAquiferMode aquiferMode) { #if HAVE_MPI - return std::make_shared(deck, comm); + return std::make_shared(deck, comm, aquiferMode); #else - return std::make_shared(deck); + return std::make_shared(deck, aquiferMode); #endif } @@ -360,7 +361,8 @@ namespace { const bool keepKeywords, const std::optional& outputInterval, Opm::ErrorGuard& errorGuard, - const bool slaveMode) + const bool slaveMode, + const Opm::NumericalAquiferMode aquiferMode) { OPM_TIMEBLOCK(readDeck); @@ -380,7 +382,7 @@ namespace { if (eclipseState == nullptr) { OPM_TIMEBLOCK(createEclState); - eclipseState = createEclipseState(comm, deck); + eclipseState = createEclipseState(comm, deck, aquiferMode); } if (eclipseState->getInitConfig().restartRequested()) { @@ -783,7 +785,8 @@ void Opm::readDeck(Opm::Parallel::Communication comm, const bool checkDeck, const bool keepKeywords, const std::optional& outputInterval, - const bool slaveMode) + const bool slaveMode, + const NumericalAquiferMode aquiferMode) { auto errorGuard = std::make_unique(); int parseSuccess = 1; // > 0 is success @@ -812,7 +815,8 @@ void Opm::readDeck(Opm::Parallel::Communication comm, eclipseState, schedule, udqState, actionState, wtestState, summaryConfig, std::move(python), initFromRestart, checkDeck, treatCriticalAsNonCritical, lowActionParsingStrictness, - keepKeywords, outputInterval, *errorGuard, slaveMode); + keepKeywords, outputInterval, *errorGuard, slaveMode, + aquiferMode); // Update schedule so that re-parsing after actions use same strictness assert(schedule); diff --git a/opm/simulators/utils/readDeck.hpp b/opm/simulators/utils/readDeck.hpp index dc7692ad471..c46b3e35481 100644 --- a/opm/simulators/utils/readDeck.hpp +++ b/opm/simulators/utils/readDeck.hpp @@ -22,6 +22,8 @@ #ifndef OPM_READDECK_HEADER_INCLUDED #define OPM_READDECK_HEADER_INCLUDED +#include + #include #include @@ -100,7 +102,8 @@ void readDeck(Parallel::Communication comm, bool checkDeck, bool keepKeywords, const std::optional& outputInterval, - bool slaveMode); + bool slaveMode, + NumericalAquiferMode aquiferMode = NumericalAquiferMode::GridCells); void verifyValidCellGeometry(Parallel::Communication comm, const EclipseState& eclipseState); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index d4051b58dbb..75535a04242 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -172,7 +172,10 @@ namespace Opm { // add the eWoms auxiliary module for the wells to the list simulator_.model().addAuxiliaryModule(this); - is_cell_perforated_.resize(local_num_cells_, false); + // Indexed by DOF index in computeTotalRatesForDof(), which the linearizer calls + // for every DOF of the model -- including those contributed by auxiliary + // modules -- so it has to span the whole DOF range, not just the grid cells. + is_cell_perforated_.resize(simulator_.model().numTotalDof(), false); } @@ -1503,7 +1506,10 @@ namespace Opm { const bool use_well_weights) const { int nw = this->numLocalWellsEnd(); - int rdofs = local_num_cells_; + // The well rows sit behind every reservoir row of the pressure system, auxiliary + // ones included -- which is how the wells themselves index them, from the size of + // the weight vector (StandardWellEquations::extractCPRPressureMatrix). + int rdofs = simulator_.model().numTotalDof(); for ( int i = 0; i < nw; i++ ) { int wdof = rdofs + i; jacobian[wdof][wdof] = 1.0;// better scaling ? @@ -1553,7 +1559,12 @@ namespace Opm { addWellPressureEquationsStruct(PressureMatrix& jacobian) const { int nw = this->numLocalWellsEnd(); - int rdofs = local_num_cells_; + // Same numbering as addWellPressureEquations(): behind every reservoir row, + // auxiliary ones included. Placing the well rows at the grid row count instead + // would have them land on top of the auxiliary rows, which is not merely a wrong + // value -- the extra columns break the assumption that the coarse pressure matrix + // has the fine matrix's sparsity pattern row for row. + int rdofs = simulator_.model().numTotalDof(); const auto wellconnections = this->getMaxWellConnections(); for (int i = 0; i < nw; ++i) { int wdof = rdofs + i; diff --git a/tests/run-numerical-aquifer-mode-comparison.sh b/tests/run-numerical-aquifer-mode-comparison.sh new file mode 100755 index 00000000000..8deb614b81f --- /dev/null +++ b/tests/run-numerical-aquifer-mode-comparison.sh @@ -0,0 +1,78 @@ +#!/bin/bash + +# Runs one deck twice -- once with each representation of its numerical aquifers -- and +# compares the two against each other. +# +# The aquifer is the same discrete system either way: the same unknowns, pore volumes, +# depths, regions and connection transmissibilities. Only where the unknown lives +# changes, so the two runs should agree to the level of arithmetic-ordering roundoff, and +# anything looser is a real difference that wants explaining rather than a wider band. +# +# That comparison is only meaningful with the two runs pinned to the same time steps and +# converged well past the reporting precision. Left to itself the adaptive stepper takes +# different substeps in the two runs -- the DOF ordering alone is enough to change where +# each Newton iteration stops inside the tolerance band -- and the difference grows to a +# few parts in a thousand for reasons that have nothing to do with the aquifer. + +if test $# -eq 0 +then + echo -e "Usage:\t$0 -- [additional simulator options]" + echo -e "\tMandatory options:" + echo -e "\t\t -i Path to read deck from" + echo -e "\t\t -f Deck file name" + echo -e "\t\t -r Path to store results in" + echo -e "\t\t -a Absolute tolerance in comparison" + echo -e "\t\t -t Relative tolerance in comparison" + echo -e "\t\t -c Path to comparison tool" + echo -e "\t\t -e Simulator binary to use" + echo -e "\tOptional options:" + echo -e "\t\t -x Compare only the summary vectors both runs produce. Needed" + echo -e "\t\t for a deck that asks for block data at an aquifer cell: that" + echo -e "\t\t cell is a grid cell in one representation and not in the" + echo -e "\t\t other, so the vector exists in one run only." + exit 1 +fi + +IGNORE_ONE_SIDED="" +OPTIND=1 +while getopts "i:f:r:a:t:c:e:x" OPT +do + case "${OPT}" in + i) INPUT_DATA_PATH=${OPTARG} ;; + f) FILENAME=${OPTARG} ;; + r) RESULT_PATH=${OPTARG} ;; + a) ABS_TOL=${OPTARG} ;; + t) REL_TOL=${OPTARG} ;; + c) COMPARE_ECL_COMMAND=${OPTARG} ;; + e) EXE_NAME=${OPTARG} ;; + x) IGNORE_ONE_SIDED="-y" ;; + esac +done +shift $(($OPTIND-1)) +TEST_ARGS="$@" + +PINNED_ARGS="--enable-adaptive-time-stepping=false \ + --tolerance-cnv=1e-8 \ + --tolerance-mb=1e-12 \ + --newton-min-iterations=2" + +mkdir -p ${RESULT_PATH} +for MODE in grid aux +do + rm -rf ${RESULT_PATH}/${MODE} + mkdir -p ${RESULT_PATH}/${MODE} + "${EXE_NAME}" ${INPUT_DATA_PATH}/${FILENAME} ${TEST_ARGS} ${PINNED_ARGS} \ + --numerical-aquifer-mode=${MODE} \ + --output-dir=${RESULT_PATH}/${MODE} + test $? -eq 0 || exit 1 +done + +echo "=== Comparing the summary of the two numerical-aquifer representations ===" +if test -n "${IGNORE_ONE_SIDED}" +then + echo " (comparing only the vectors both runs produce -- see -x)" +fi +${COMPARE_ECL_COMMAND} -t SMRY ${IGNORE_ONE_SIDED} -a ${RESULT_PATH}/grid/${FILENAME} \ + ${RESULT_PATH}/aux/${FILENAME} \ + ${ABS_TOL} ${REL_TOL} +exit $?