Numerical aquifers as auxiliary cells instead of grid cells by hnil · Pull Request #7309 · OPM/opm-simulators

hnil · 2026-08-11T08:13:35Z

An alternative representation of numerical aquifers: instead of taking over a grid cell,
the aquifer gets a degree of freedom of its own, appended after the grid ones, with an
authored volume, depth and connection list. Opt in with --numerical-aquifer-mode=aux;
the default path is unchanged and byte-identical.

This is mainly useful where the grid has to stay a proper grid. AQUNUM currently
repurposes a real cell for what is pure flow bookkeeping, which is fine for flow but not
for anything that reads the grid as geometry — mechanics in particular, where such a cell
corrupts the assembly. Here the grid holds only physical rock. As a side effect AQUNUM
also works with EdgeConformal=true, which the NNC-based path cannot.

Getting there means making BaseAuxiliaryModule degrees of freedom that carry the model's
own conservation equations work at all: two-sided neighbour info in the TPFA linearizer,
per-DOF containers sized for them, Newton and convergence including them, and a CPR weight.
Most commits are that groundwork; the aquifer client is the last few. It also fixes a
latent out-of-bounds read of neighborInfo_ whenever an auxiliary module has DOFs.

A new CTest compares the two representations on the AQUNUM decks. Serial only for now —
auxiliary cells in parallel are refused at setup until they are part of the partitioning.

Needs OPM/opm-common#5286 and OPM/opm-grid#1059. Draft — CI not yet green, and I have not
run the full opm-tests suite.

🤖 Generated with Claude Code

Auxiliary modules may append degrees of freedom after the grid ones (BaseAuxiliaryModule::numDofs()), and the solution vector, the residual and the Jacobian are already sized for numTotalDof(). The per-DOF containers owned by the discretization were not: dofTotalVolume_, isLocalDof_, the storage cache and the intensive-quantity cache were all sized by numGridDof(), so an auxiliary DOF carrying the model's own equations had nowhere to store its volume or its intensive quantities. Size them by numTotalDof() instead. Auxiliary DOFs have no grid geometry to derive a volume from, so their entries are left for the auxiliary module to author; they are not shared through the grid's interior-border interface and are therefore local by construction. No existing configuration registers an auxiliary module with a non-zero DOF count, so numTotalDof() == numGridDof() throughout and this is a no-op for every current run. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

updateFlowsInfo() iterated over numTotalDof() while the containers it writes into -- flowsInfo_ and floresInfo_ -- have one row per grid element, and the block-flow lookup maps the DOF index back to a cartesian index. Both are grid-cell concepts. With auxiliary DOFs present the loop would run past the end of those tables and ask the vanguard for the cartesian index of a DOF that has none. Bound the loop by numGridDof(). Fluxes on auxiliary connections are reported by the client that owns them, not through the ECL flow tables. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The solution vector handed to solveJacobianSystem() was sized by numGridDof() while the Jacobian and residual span numTotalDof(), and dx_old_ -- which stabilizeNonlinearUpdate() compares against it -- had the same problem. With auxiliary DOFs present this is a size mismatch against the linear system. Size both by numTotalDof(). solUpd_ is deliberately left alone: it is only ever indexed through the element mapper, so it is a grid-only container. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

is_cell_perforated_ is indexed in computeTotalRatesForDof(), which the linearizer calls for every DOF of the model, but it was sized by the number of grid cells. With auxiliary DOFs present that is an out-of-bounds read on every source-term evaluation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

neighborInfo_ had one row per grid element while linearize_cell(), updateStoredTransmissibilities() and the block-address loop in createMatrix_() all index it up to numTotalDof(). With an auxiliary module that declares degrees of freedom this reads past the end of the table; SparseTable::operator[] only guards it with an assert, so a release build reads a garbage row rather than failing. Let an auxiliary module report its flux connections through the new BaseAuxiliaryModule::addConnections(), and build the missing rows from them. Such a connection has no geometric face, so the transmissibility carries the geometry, the face area is unity and there is no face direction; everything else -- threshold pressure, gravity head, thermal half-transmissibility, diffusivity, dispersivity -- is read from the problem exactly as it is for a grid face. This is what lets an auxiliary cell be assembled by the model's own local residual instead of a hand-written copy of it, and it makes the values follow updateStoredTransmissibilities() for free. Each connection is reported once but entered on both endpoints. linearize_cell() writes only the row it is iterating -- residual[globI] and the blocks (globI,globI) and (globJ,globI) -- so the counter-flux and the partner's diagonal contribution come from the partner's own row. Entering a connection on one side alone would let an auxiliary cell push mass into a grid cell that never sees it in return, and the run would converge on a system that does not conserve mass. The default addConnections() is empty, so modules that assemble their own equations (the well models) are unaffected, and with no auxiliary DOFs the table is built exactly as before. A hard check on the row count replaces the assert, since the failure is silent otherwise. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Two ordering problems that only bite once an auxiliary module declares degrees of freedom. finishInit() sizes every per-DOF container from numTotalDof(), so a module registered afterwards leaves all of them short. That is exactly what happens today: BlackoilWellModel registers itself from initialSolutionApplied(), long after finishInit(). It is harmless there because the well model declares no degrees of freedom, but it would be a silent, hard-to-trace corruption for one that does. Reject that registration order instead of letting it through. addAuxiliaryModule() also calls applyInitial() at registration time, while applyInitialSolution() starts by zeroing the whole solution vector -- so anything an auxiliary module wrote is erased before the run begins. Give the modules a second applyInitial() after the wipe, before the history copy and the checkDefined() sweep, both of which span the whole solution vector. Only modules with degrees of freedom are called, so this is a no-op for the existing ones. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Every path in invalidateAndUpdateIntensiveQuantities() walks the grid, so degrees of freedom introduced by auxiliary modules never have their intensive quantities refreshed. Add a loop over them. No element context is needed: the blackoil intensive quantities update from the problem's index-based accessors alone, which is what lets an auxiliary cell carry the model's own equations rather than a hand-written copy of them. The early return in the avoid-element-context branch becomes an else so there is a single exit to hook onto; the behaviour of that branch is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The device path copies the neighbor info, the domain and the residual to the GPU and assembles there, while an auxiliary module's linearize() runs on the host after the copy back. Auxiliary degrees of freedom would be assembled inconsistently rather than visibly wrongly, so fail at setup instead of producing plausible numbers. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Auxiliary degrees of freedom are treated one of two ways depending on what kind of unknown they hold, and until now the distinction was made by index -- everything past numGridDof() was assumed to scale itself and was kept out of the error norm and the primary-variable switching. That is right for a well's bottom hole pressure, and wrong for an auxiliary cell holding the model's own primary variables. Add BaseAuxiliaryModule::carriesModelEquations(), false by default, and a model-level dofCarriesModelEquations(dofIdx) that answers it for any degree of freedom. The consumers follow in the next commits. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Use dofCarriesModelEquations() in place of the index test, at the three places that treated every auxiliary degree of freedom as somebody else's problem: the Newton error norm, the primary-variable and equation weights of the blackoil model, and the update loop. The update is the substantive one. Auxiliary degrees of freedom got a plain u -= du, bypassing updatePrimaryVariables_() and therefore the variable-meaning switching -- so an auxiliary cell could never cross a phase-appearance boundary. Dispatch on the same predicate: cells get the model's update, everything else keeps the plain one. A dormant auxiliary cell has zero volume and is still skipped by the existing volume test in the error norm. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

trans_ and thermalHalfTrans_ are hash maps keyed on the packed DOF pair, so a connection that does not come from face geometry fits the existing storage exactly; only a way to write it was missing. Calling the setter again updates an existing connection, which is how a transmissibility that depends on the solution -- a fracture aperture -- can be refreshed between iterations without disturbing the sparsity pattern. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

localConvergenceData() walks the grid, so auxiliary cells never entered the CNV/MB measures. The iteration could then be declared converged with their mass balance still violated. Add them, with their pore volume, on the same footing as grid cells. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The simulator constructs the model and the problem, then runs the model's finishInit(), then the problem's. Every per-DOF container is sized during the model's finishInit() from numTotalDof(), so a module that introduces degrees of freedom has to be registered before that -- earlier than the problem's own finishInit(), which is where a problem would naturally do such things. Add Problem::registerAuxiliaryCellModules(), empty by default, called at the top of the model's finishInit(). The grid and the deck are available by then, which is what such a module needs to know its DOF count. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

FvBaseLinearizer assembles by looping over grid elements. An auxiliary degree of freedom has no element, so the loop walks past it and its row is left empty -- which for a module that expects the model's own conservation equations to be assembled there is a silent failure: the row is all zeros, and the system is either singular or, once anything conditions the diagonal, a spurious 0 = 0. TpfaLinearizer loops over degrees of freedom instead and keeps the auxiliary connections in neighborInfo_ next to the geometric ones, so it does assemble them. State the difference as a trait on the two linearizers and check it where the auxiliary modules have just been registered. Only modules that both introduce DOFs and declare that those DOFs carry model equations are affected; the well model and the other existing auxiliary modules are not. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

…y exist The element-context-free IntensiveQuantities::update() overload is implemented for the plain black-oil equations alone; solvent, extbo, polymer, foam, MICP, brine, diffusion and dispersion all static_assert against it. Calling it unconditionally from the auxiliary tail loop instantiated it for every type tag, including flow_brine and flow_extbo, which then failed to compile. Instantiate it under AvoidElementContext, which is the property that already decides whether that overload is usable, and leave a configuration that does have auxiliary DOFs but cannot update them with an error rather than with silently uninitialised intensive quantities. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The CPR pressure system appends one row per well after the reservoir rows. The well model put the first of them at the grid cell count, while the wells themselves derive the same index from the size of the weight vector -- which is the number of rows of the fine matrix. The two agree only as long as there are no auxiliary degrees of freedom. With auxiliary DOFs the well rows land on top of them: the well's coupling columns are inserted into an auxiliary row, so the coarse pressure matrix no longer has the fine matrix's sparsity pattern row for row, and the well rows proper are left empty. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The true-IMPES weights are computed by walking the grid, so they never reach a degree of freedom introduced by an auxiliary module. The weight is not a refinement: it multiplies the row on its way into the coarse pressure system, so an unset one deletes the row and leaves that system singular. Fill them in from the assembled diagonal block -- the quasi-IMPES weight, which needs no element context. The per-row part of the quasi-IMPES calculation is factored out rather than duplicated. Nothing changes where there are no auxiliary DOFs: the added loop is then empty. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

…an element The index-based IntensiveQuantities::update() covers the plain black-oil equations and the energy module; solvent, extbo, polymer, foam, MICP, brine, diffusion and dispersion still reach for an element context. That was expressed only as a wall of static_asserts inside the overload, so a caller with no element to offer had no way to ask in advance whether it may call it. Name the condition and use it in both places. The auxiliary intensive-quantity loop was previously gated on AvoidElementContext, which is a stricter and different question -- that property is about how the *grid* cells are updated, and a configuration may well drive the grid through element contexts and still be perfectly able to update an auxiliary degree of freedom without one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

…pdate The element-context update() calls updateEnergyQuantities_() for a fully implicit thermal model; the index-based overload did not. Nothing exercised the combination before -- AvoidElementContext is never set on a fully-implicit- thermal target -- but an auxiliary degree of freedom has no element and so is updated through the index-based path whatever the grid cells do. The symptom is quiet and total: the fluid enthalpies, the rock internal energy and the thermal conductivity are left uncomputed, so nothing in the cell depends on its temperature, its diagonal block has an all-zero temperature column, and the ILU decomposition reports a singular matrix. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

So that the parallel state can be built in either mode, like the serial one. Defaulted, so existing callers are unaffected. Note that the mode is a member of EclipseState and is serialized with it, so the ranks which receive the broadcast state agree with rank 0 about how the aquifers are represented rather than having to be told separately. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Two classes, no consumer yet; the problem-side wiring follows. FlowAuxCellModule is the base for auxiliary modules whose degrees of freedom are cells: they satisfy the model's own conservation equations and are assembled by its local residual, but their pore volume, depth, regions and connection list are authored rather than read off a grid entity. That is deliberately not the shape of the well models, which are also auxiliary modules but whose unknowns are of a different kind and which assemble their own equations. NumericalAquiferAuxCells is the first client. Its data comes straight from the AQUNUM records, and -- this is the point -- its connections are the very NNCs the grid-cell representation generates, reused rather than reimplemented: the intra-aquifer chain from aquiferCellNNCs() and the AQUCON connections from aquiferConnectionNNCs(). The two representations are therefore the same discrete system with the same coefficients, and can be compared directly rather than merely eyeballed. Connections whose reservoir cell this rank does not own are skipped; the build is deferred out of the constructor because it needs the cartesian-to-compressed mapping. An AQUNUM record on an active grid cell is rejected: the grid-cell representation quietly repurposes such a cell, overwriting its volume, depth and regions and zeroing its permeability, and there is no faithful reading of that when the aquifer lives outside the grid. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Create the modules from registerAuxiliaryCellModules(), which the model calls before it sizes anything, and give them what a cell needs but an auxiliary degree of freedom cannot read off the grid: - total volume, stated by the module, filled where the grid cells take theirs from their entity's geometry; - reference porosity, as the authored pore volume over that volume; - rock fraction left at zero -- an auxiliary cell is a bookkeeping volume, not rock, and the energy storage term is rockFraction times the rock's internal energy, so zero is what says "stores no heat"; - depth, answered by dofCenterDepth() so that the gravity head between an auxiliary cell and its neighbours is built from it; - connection transmissibilities, published into the transmissibility store once the grid's own are final. It is keyed on the degree of freedom pair and does not care whether a connection came from a face. Restart together with auxiliary cells is refused: the restart file is written per grid element, so their state would come back undefined rather than merely stale. No deck can select the auxiliary representation yet, so this changes nothing for any run; the flag comes with the last piece. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Threads the choice from the command line down to where EclipseState is built, which is the only point at which it can take effect: the representation is decided as the state is constructed. Held back until now on purpose. With no client able to represent an aquifer outside the grid, selecting that mode would have dropped the aquifer silently and produced a plausible wrong answer rather than an error. Also stops creating AquiferNumerical in the auxiliary-cell mode. It 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 -- and none of that applies when the cells are not in the grid. Its AAQ* reporting has to be rebuilt from the auxiliary block; that is still to do. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Enough for the auxiliary aquifer cells to be assembled and solved. They still do not converge; see below. - Numbering. Enumerate aquifer by aquifer and, within an aquifer, in declaration order, rather than through allAquiferCells(), which is an unordered map. The hash order was making the degree of freedom numbering arbitrary and losing the chain structure -- which cell carries the reservoir connections. It showed up as the connections landing on the wrong cell of the chain. - Initial state. Auxiliary cells cannot be equilibrated the ordinary way, which needs the cell's geometry, so each starts from the state of the reservoir cell it is connected to, carried to its own depth and filled with water; an explicit AQUNUM initial pressure overrides that. Only the cell carrying the AQUCON connections has a reservoir neighbour, so the rest of the chain starts from the same place -- per aquifer, so two aquifers cannot borrow each other's. - Per-element rock tables. The saturation-function and thermal parameter tables are built per grid element and keyed on grid entities throughout, so they cannot be extended over degrees of freedom that have no entity. An auxiliary cell borrows the parameters of the grid cell it hangs off instead; exact while the two share a region, which is the ordinary case for a water-filled aquifer. - The ROCK table lookup is guarded: auxiliary DOFs have no entry and take the first region, where before it read off the end of the array. - applyInitial() is no longer called at registration time for modules with degrees of freedom. It runs from applyInitialSolution() instead, once the solution vector and the initial fluid states exist; the early call happened before either and threw. Known remaining problem: the linear solver reports a singular diagonal block for an auxiliary cell. The likely cause is the initial state -- copying the partner cell's fluid state and overriding the saturations can leave the primary variable meaning inconsistent with the composition, so an equation ends up with an identically zero row. Reproduce with flow_energy_geomech data/SIMPLE_MECH_NX_11_NY_11_NZ_25_AQUNUM \ --numerical-aquifer-mode=aux Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The mixing-rate controls, the maximum oil saturation and the drift-compensation vector are all addressed by degree-of-freedom index -- the intensive quantities reach maxGasDissolutionFactor() and maxOilSaturation() the same way for a grid cell and for an auxiliary one, and the source-term loop covers every row of the residual. Sizing them for the grid alone leaves an auxiliary DOF reading past the end, which is an out-of-bounds access rather than a missing contribution. A no-op where there are no auxiliary DOFs: numTotalDof() is then numGridDof(). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

The region arrays come from the field properties over the grid, so they stop one entry short of the auxiliary degrees of freedom. They are read by degree-of-freedom index -- the energy quantities ask pvtRegionIndex() on the way to the PVT tables -- which makes a short array an out-of-bounds read that then indexes the tables with whatever happened to follow in memory. Extend them and fill in what the auxiliary module declares. An empty array already means "one region for everything", so it is left alone. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

AquiferNumerical locates its cells through the grid and reads their fluxes off an element context, so it cannot say anything about an aquifer represented as auxiliary cells -- and with it absent, ANQP, ANQR and ANQT simply read zero. Form the same two quantities from the degrees of freedom instead. The pressure is the water-volume-weighted water pressure over the aquifer's cells, which is the average the grid-cell representation takes over the same cells. The influx is the water flux across the connections that reach into the reservoir, computed by handing the linearizer's own cached neighbour information back to LocalResidual::computeFlux() -- so what is reported is the flux the equations were assembled with, upwinding, gravity and threshold pressures included, rather than a second implementation that would drift away from it. The intra-aquifer chain is left out on purpose: an aquifer's influx is what it gives the reservoir, not what moves inside itself. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

An aquifer's connections enter the grid-cell representation as input NNCs, where Transmissibility::applyMultRegTToInputNncTrans_() scales them by the MULTREGT multiplier before they are used. The auxiliary-cell representation took the raw value, so a deck that damps or shuts an aquifer connection with MULTREGT had the aquifer connected on one path and effectively isolated on the other. Visible on opm-tests/aquifers/AQUNUM-03, whose MULTREGT of 0.001 across the region boundary the AQUCON connection crosses reduces the influx by three orders of magnitude: ANQR reads 0.049 in grid mode against 44.8 without the multiplier. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

… other The two representations are meant to be the same discrete system -- same unknowns, pore volumes, depths, regions and connection transmissibilities, with only the place the unknown lives differing -- so running one deck both ways and comparing the results is a sharper test than a regression against stored data for either mode on its own. It needs no reference data, and it is what caught the missing MULTREGT multiplier on the auxiliary cells' connections. The driver pins the time steps and tightens the convergence tolerances first. Without that the adaptive stepper takes different substeps in the two runs -- the degree-of-freedom ordering alone changes where each Newton iteration stops inside the tolerance band -- and the two drift apart by a few parts in a thousand for reasons that have nothing to do with the aquifer. Pinned, they agree to the precision the summary file stores. AQUNUM-02 asks for BPR at two of its own aquifer cells, which are grid cells in one representation and not in the other; that case compares the vectors both runs produce, and says so. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

A numerical aquifer inside the grid contributes its pore volume and its fluid to FPR, the FIP report and every region total. Moving it out of the grid took it out of those numbers as well, so the same deck reported a different field pressure depending only on how the aquifer was represented -- on one deck 215.4 against 212.2 bar, with the well and mechanics answers identical. Let the fluid-in-place buffers span the auxiliary degrees of freedom and fill them through the same index-based path the grid cells use, so what is accumulated is what a grid cell would have accumulated rather than a separate estimate. The reporting regions are extended alongside, from the field properties of the cell each auxiliary cell names -- which is where an aquifer cell's regions come from in the grid representation too. Everything written per grid cell stays grid-sized; the restart arrays stop at the grid, since a degree of freedom outside it has no cell to be written against. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

WellModelGhostLastMatrixAdapter took the owned rows to be the first N of the matrix: apply() stops there and ghostLastProject() zeroes everything after it. That holds while every row belongs to a grid cell, but a degree of freedom which has none -- a numerical aquifer or a fracture represented outside the grid -- is appended after the grid rows, and so lands behind the ghosts. The operator would then skip its row and the projection would zero it: the auxiliary equations disappear from the Krylov operator and are solved as 0 = 0, with nothing thrown and no residual looking wrong. Take the owned rows as a list of half-open bands instead. One band is the prefix it always was, so nothing changes where there are no auxiliary degrees of freedom; the projection zeroes the gaps between bands rather than everything past a single count. The same assumption is what dynamic grid refinement runs into -- adding or removing cells reshuffles the interior/ghost split, and "the owned rows are the first N" has to be re-established on every topology change -- so describing them explicitly is shared groundwork rather than a detour. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

…ioned Which rank owns a degree of freedom that has no geometry is a choice, and it 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 live intensive quantities, and exactly one rank may contribute its accumulation -- which means entering it in the communication index set as owner there and copy elsewhere. None of that exists yet. What happens instead is that every rank builds the same auxiliary cells and then finds it has connections to only some of them, which surfaces as a confusing complaint about an aquifer with no reservoir connection on this process. Say what is actually wrong, and where the way out is. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

hnil force-pushed the pr/numerical-aquifer-auxcells branch from 8742bc2 to cd67499 Compare August 11, 2026 12:05

hnil added the manual:enhancement This is an enhancement/improvent that needs to be documented in the manual label Aug 12, 2026

hnil and others added 28 commits August 14, 2026 16:29

hnil and others added 4 commits August 14, 2026 16:29

hnil force-pushed the pr/numerical-aquifer-auxcells branch from cd67499 to cdf3acd Compare August 14, 2026 14:31

Navigation Menu

Uh oh!

Numerical aquifers as auxiliary cells instead of grid cells - #7309

Numerical aquifers as auxiliary cells instead of grid cells#7309
hnil wants to merge 32 commits into
OPM:masterfrom
hnil:pr/numerical-aquifer-auxcells

hnil commented Aug 11, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

hnil commented Aug 11, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant