Add Coats pressure-equation weights for blackoil and compositional mo…

[hnil]

this is based on a closed pullrequest made by copilot in opm-simulators main fork.

It is work in progess to add better cpr weights which also can work for compositional.

Agent-Logs-Url: https://github.com/OPM/opm-simulators/sessions/db94714d-c71d-48fb-87ba-6e9c13054f52

[Copilot]

Pull request overview

Adds two new CPR pressure-equation weight calculators — coatsblackoil and coatscompositional — based on Coats (1980), and exposes them through ISTLSolver's weightsType selector. The blackoil variant analytically inverts the 3-phase formation-volume/dissolution matrix (extended with Rsw/Rvw cross-terms), while the compositional variant builds saturation-weighted specific-volume weights for the ptflash model (plus a 1/ρ_w water weight, energy left at zero). The PR description notes this is work in progress.

Changes:

• New Amg::getCoatsWeightsBlackoil and Amg::getCoatsWeightsCompositional template functions in getQuasiImpesWeights.hpp.
• New coatsblackoil / coatscompositional branches in ISTLSolver::getWeightsCalculator, plus an updated error message listing the new options.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 2 comments.

File	Description
opm/simulators/linalg/ISTLSolver.hpp	Wires the two new weight calculators into the weightsType dispatch and extends the unknown-type error message.
opm/simulators/linalg/getQuasiImpesWeights.hpp	Adds the two new template functions implementing Coats-style pressure weights for blackoil and compositional models.

Comments suppressed due to low confidence (6)

opm/simulators/linalg/getQuasiImpesWeights.hpp:409

• The comments at lines 397-398 and 405-406 say the oil/gas phase is "undersaturated" in these branches, but when primaryVarsMeaningGas() == Rv the free oil phase is absent (oil saturation is zero, the variable encodes how much oil is vaporized in the gas phase). Likewise for == Rs, the free gas phase is absent. The actual code (setting rs=0 / rv=0) is correct, but the rationale should describe phase absence rather than undersaturation.

                        // When the gas primary variable represents Rv (oil vaporized in gas),
                        // the oil phase is undersaturated, so Rs should be treated as zero.
                        if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
                            rs = 0.0;
                        }
                    }
                    if (FluidSystem::enableVaporizedOil()) {
                        rv = Toolbox::template decay<double>(fs.Rv());
                        // When the gas primary variable represents Rs (gas dissolved in oil),
                        // the gas phase is undersaturated, so Rv should be treated as zero.
                        if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
                            rv = 0.0;
                        }

opm/simulators/linalg/getQuasiImpesWeights.hpp:454

• D = 1 - rs*rv - rvw*rsw is used as divisor without any safeguard. For degenerate / pathological PVT states (e.g. fully saturated systems near a singularity) this denominator can approach zero and produce extreme/NaN weights that are then propagated through the CPR preconditioner. Consider clamping D away from zero (similar to how trueimpesanalytic is guarded by phaseIsActive(oil)&&phaseIsActive(gas) before computing 1-rs*rv).

                // Denominator D = 1 - Rs*Rv - Rvw*Rsw
                const double D = 1.0 - rs * rv - rvw * rsw;

                // Compute weights by solving the linear system derived from the
                // requirement that sum_i w_i * M_i = total reservoir saturation.
                //
                // 3-phase (w, o, g) with full cross-terms:
                //   ww = [ Bw*(1 - Rs*Rv) - Rsw*(Bg - Rv*Bo) ] / D
                //   wo = [ Bo*(1 - Rvw*Rsw) - Rs*(Bg - Rvw*Bw) ] / D
                //   wg = ( Bg - Rvw*Bw - Rv*Bo ) / D
                //
                // Degenerate (fewer active phases) cases collapse correctly:
                //   water-only:    ww = Bw
                //   gas+oil:       wo = (Bo - Rs*Bg)/(1-Rs*Rv),  wg = (Bg - Rv*Bo)/(1-Rs*Rv)
                //   water+gas:     ww = (Bw - Rsw*Bg)/(1-Rsw*Rvw), wg = (Bg - Rvw*Bw)/(1-Rsw*Rvw)
                //   water+oil:     ww = Bw, wo = Bo  (no cross-dissolution)

                if (waterActive) {
                    const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(
                        FluidSystem::solventComponentIndex(FluidSystem::waterPhaseIdx));
                    bweights[activeCompIdx] = static_cast<LhsEval>(
                        (Bw * (1.0 - rs * rv) - rsw * (Bg - rv * Bo)) / D);
                }
                if (oilActive) {
                    const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(
                        FluidSystem::solventComponentIndex(FluidSystem::oilPhaseIdx));
                    bweights[activeCompIdx] = static_cast<LhsEval>(
                        (Bo * (1.0 - rvw * rsw) - rs * (Bg - rvw * Bw)) / D);
                }
                if (gasActive) {
                    const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(
                        FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx));
                    bweights[activeCompIdx] = static_cast<LhsEval>(
                        (Bg - rvw * Bw - rv * Bo) / D);

opm/simulators/linalg/getQuasiImpesWeights.hpp:464

• *std::ranges::max_element(bweights, |a|<|b|) returns the element with the largest |·| but the signed value. bweights = 0.0 initializes all entries, and when no phase contributes to a weight (e.g. compositional models where every den_i == 0 and rho_w <= 0), abs_max evaluates to 0.0 and the guard std::fabs(abs_max) > 0.0 prevents division-by-zero, but the resulting all-zero weight row will produce a degenerate pressure equation for that cell. Worth at least logging / asserting rather than silently returning a zero block, otherwise CPR will silently use a zero row and the solver behavior will be hard to debug.

                // Normalize so that the maximum absolute weight equals one.
                const double abs_max =
                    *std::ranges::max_element(bweights,
                                              [](double a, double b)
                                              { return std::fabs(a) < std::fabs(b); });
                if (std::fabs(abs_max) > 0.0) {
                    bweights /= std::fabs(abs_max);
                }

opm/simulators/linalg/getQuasiImpesWeights.hpp:578

• The water equation index in compositional ptflash is conti0EqIdx + numComponents, which collides with EnergyIndices::energyEqIdx (also PVOffset + numComponents) when energy is enabled. The existing ptflash residual at opm/models/ptflash/flashlocalresidual.hh:90 uses the same index for water storage, but neither this new code nor the residual code accounts for the ordering when both energy and water are enabled. If both are ever enabled simultaneously, this weight will be written into the wrong equation slot.

                if constexpr (waterEnabled) {
                    const int waterEqIdx = Indices::conti0EqIdx + numComponents;
                    const double rho_w =
                        Toolbox::template decay<LhsEval>(
                            fs.density(FluidSystem::waterPhaseIdx));
                    bweights[waterEqIdx] = (rho_w > 0.0)
                        ? static_cast<LhsEval>(1.0 / rho_w)
                        : static_cast<LhsEval>(0.0);
                }

opm/simulators/linalg/getQuasiImpesWeights.hpp:590

• The normalization block (max_element with |a|<|b| comparator + guarded divide by std::fabs(abs_max)) is now duplicated across getQuasiImpesWeights, getTrueImpesWeights, getCoatsWeightsBlackoil, and getCoatsWeightsCompositional. Consider extracting a small normalizeBlockWeights(bweights) helper to keep behavior consistent (and to centralize the zero-guard that only the two new functions currently apply).

                // Normalize so that the maximum absolute weight equals one.
                const double abs_max =
                    *std::ranges::max_element(bweights,
                                              [](double a, double b)
                                              { return std::fabs(a) < std::fabs(b); });
                if (std::fabs(abs_max) > 0.0) {
                    bweights /= std::fabs(abs_max);
                }

opm/simulators/linalg/getQuasiImpesWeights.hpp:344

• pressureVarIndex is silently ignored (/*pressureVarIndex*/) in both new functions, but ISTLSolver still passes it in. Either honor the index (so users with non-default pressIndex can actually target a different pressure row) or drop the parameter from the signature with a comment explaining why. As is, the function name suggests it computes the pressure-equation weights, yet pressIndex from the caller has no effect.

    void getCoatsWeightsBlackoil(int /*pressureVarIndex*/,
                                 Vector& weights,
                                 const ElementContext& elemCtx,
                                 const Model& model,
                                 const ElementChunksType& element_chunks,
                                 [[maybe_unused]] bool enable_thread_parallel)

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[hnil]

@copilot apply changes based on this feedback