Attempts at improving the zero og case for BO

[moyner]

No description provided.

[b11111010000]

Thanks for your work, Olav! This does seems to work for some well configurations but not all of them. I'll share an reproducible example case and generated report below.

---

MSW BlackOil injection deadlock — PR #226 verification

Follow-up to jutuldarcy_msw_blackoil_bug_report.md. Reports what we
observed when testing sintefmath/JutulDarcy.jl PR #226 (bo-wat-tests,
HEAD 2757dded) as a candidate fix for the MultiSegmentWell BlackOil
injection deadlock.

TL;DR — a minimal reproducer using only the registered JutulDarcy
API still exhibits the full deadlock (injector auto-disabled, zero flux)
on both stock v0.3.5 and the PR #226 branch. PR #226 does not resolve
this case.

Shareable reproducer

msw_blackoil_repro.jl (single file, ~90 lines) reproduces the deadlock
using only Jutul, JutulDarcy, MultiComponentFlash, HYPRE, LinearAlgebra.

Setup summary:

• 20 × 1 × 20 unstructured Cartesian mesh, top at 3500 m depth
• 2-phase BlackOil (LiquidPhase + VaporPhase) with dissolved gas, PVT
  tables generated from a CO₂/H₂O Peng–Robinson EOS
• One vertical injector penetrating all nz cells
• Pressure equilibrium: 360 bar at top, pure brine (GOC above domain, Rs=0)
• Injector BHP target 500 bar, pure CO₂ composition [0.0, 1.0]
• Simulate 2 × 2-day steps

Toggle simple_well at the top of the script for the working comparison.

Observed output

Version	simple_well	Injector BHP (MPa)	gas_mass_rate (kg/s)
v0.3.5 (stock)	true	50.0, 50.0	0.034, 0.030
v0.3.5 (stock)	false	36.04, 36.04	0.0, 0.0
PR #226 (bo-wat-tests)	true	50.0, 50.0	0.109, 0.0006
PR #226 (bo-wat-tests)	false	36.04, 36.04	0.0, 0.0

MSW injector pinned at hydrostatic (36.04 MPa ≈ 360 bar) on both
branches — DisabledControl engaged by the facility on the first
iteration because the perforation flux evaluates to zero.

How to run

mkdir repro && cd repro
julia --project=. -e '
using Pkg
Pkg.add(PackageSpec(url="https://github.com/sintefmath/JutulDarcy.jl", rev="bo-wat-tests"))
Pkg.add(["Jutul", "MultiComponentFlash", "HYPRE"])
'
# (copy msw_blackoil_repro.jl here)
julia --project=. msw_blackoil_repro.jl

For stock v0.3.5, replace the first Pkg.add with
Pkg.add(PackageSpec(name="JutulDarcy", version="0.3.5")).

What PR #226 changes

From the diff (src/blackoil/variables/varswitch.jl, types.jl,
utils.jl):

1. blackoil_unknown_init — water-filled cells (sw ≈ 1) are now
   classified as OilAndGas phase instead of OilOnly.
2. eps_s default lowered from 1e-10 to MINIMUM_COMPOSITIONAL_SATURATION.
3. ensure_endpoints! call in PhaseRelativePermeability commented out.
4. optimization_level plumbed through setup_reservoir_model.

None of these touch src/blackoil/wells.jl, where the actual MSW
perforation flux formula lives (q = λ_t * dp_v * sG * bG at lines
50–54). Since the formula still multiplies by the well-segment gas
saturation, and the segment is initialised from a brine-filled reservoir
(Sg = 0), the flux is still identically zero on the first Newton step
in cases where varswitch's revised initialisation doesn't bootstrap sG
above zero fast enough.

Suggested next step

Follow the pattern used by SimpleWell and the immiscible wells: accumulate
injection flux from all phases using total mobility, then distribute by
well mass fractions instead of segment saturations. Snippet from the
original bug report:

# Instead of (MSW, current):
q = λ_t * dp_v * sG * bG

# Consider (analogous to SimpleWell):
if dp_v < 0.0
    q_inj += λ_t * dp_v * ρ_w[v, wc]
end
# then distribute q_inj using state_well.MassFractions

This decouples injection composition from the well's internal saturation
state, which is the root cause of the deadlock.

[b11111010000]

----- Julia Script:

#
# Minimal reproducer: MultiSegmentWell BlackOil injection deadlock
#
# Related to sintefmath/JutulDarcy.jl PR #226. Self-contained — depends only on
# registered packages.
#
# Behaviour we observe on the `bo-wat-tests` branch of PR #226:
#   a MultiSegmentWell gas injector stays at reservoir pressure and produces
#   zero surface gas rate. A SimpleWell with otherwise identical setup injects
#   normally. The difference comes from how MSW BlackOil computes the
#   perforation flux: it multiplies by the well-segment gas saturation
#   (`sG = s_w[v, wc]`), which is zero for a well initialised from a
#   brine-filled reservoir.
#
# Flip `simple_well` at the top to compare.
#
# Run: `julia --project=. msw_blackoil_repro.jl`
#

using Jutul, JutulDarcy, MultiComponentFlash, LinearAlgebra

simple_well = false   # set to true for the working SimpleWell comparison

# ---------------- Mesh & rock ----------------
# Deep reservoir with non-trivial hydrostatic column (matches the setup in
# the original bug report: ~360 bar hydrostatic, injection at ~500 bar).
nx, ny, nz = 20, 1, 20
phys = (500.0, 10.0, 200.0)     # 3500 m deep top
z_top = 3500.0
mesh = UnstructuredMesh(CartesianMesh((nx, ny, nz), phys), z_is_depth = true)
for (i, pt) in enumerate(mesh.node_points)
    mesh.node_points[i] = pt + [0.0, 0.0, z_top]
end
nc   = number_of_cells(mesh)

Darcy, bar, day, kg, meter = si_units(:darcy, :bar, :day, :kilogram, :meter)

domain = reservoir_domain(mesh;
    permeability = fill(0.1 * Darcy, 3, nc),
    porosity = fill(0.2, nc))

# ---------------- PVT: oil (brine) + gas (CO2) with dissolved gas --------
# Generate a blackoil PVTO/PVTG pair from a CO2/H2O EOS sample so the
# resulting system is a StandardBlackOilSystem with dissolved-gas switching.
mix  = MultiComponentMixture(["CO2", "H2O"])
eos  = GenericCubicEOS(mix, PengRobinson())
T_res = convert_to_si(50.0, :Celsius)
z_oil = normalize([0.01, 0.99], 1)

tables = generate_pvt_tables(eos, z_oil, T_res;
    n_pvto = 10, n_pvdg = 10, n_undersaturated = 10)

# ---------------- Single vertical MSW injector ----------------
x_inj = (nx ÷ 2) * phys[1] / nx
y_inj = 0.5 * phys[2] / ny
injector_trajectory = [
    x_inj  y_inj  z_top;
    x_inj  y_inj  z_top + phys[3]
]
Inj = setup_well_from_trajectory(domain, injector_trajectory;
    name = :Inj, simple_well = simple_well)

# ---------------- Reservoir model: 2-phase (oil+gas), no water -------------
model = JutulDarcy.setup_reservoir_model_from_blackoil_tables(
    domain, tables; wells = [Inj], disgas = true, water = false)

import JutulDarcy: brooks_corey_relperm
sg   = collect(range(0.0, 1.0, length = 20))
krg  = brooks_corey_relperm.(sg,        n = 2.0, residual = 0.0, residual_total = 0.0)
krog = brooks_corey_relperm.(1 .- sg,  n = 2.0, residual = 0.0, residual_total = 0.0)
set_relative_permeability!(model; sgof = hcat(sg, krg, krog))

# ---------------- Equilibrium: brine-filled, Sg=0 everywhere --------------
# Put GOC above the domain so every cell is liquid; set Rs=0 (undersaturated
# oil with zero dissolved gas) so the initial state has Sg identically zero.
# ~360 bar at the reservoir top (matches original bug report)
eql = EquilibriumRegion(model, 360.0*bar, z_top;
    goc = z_top - 10.0, rs = 0.0)
state0 = setup_reservoir_state(model, eql)

# ---------------- Schedule: 2 report steps, BHP injector ~2× hydrostatic ---
dt = fill(2.0 * day, 2)
forces = map(dt) do _
    ctrl = Dict(
        :Inj => setup_injector_control(
            500.0*bar, :bhp, [0.0, 1.0];    # pure gas mixture, well above hydrostatic
            density = 700.0),
    )
    setup_reservoir_forces(model; control = ctrl)
end

# ---------------- Simulate and print well diagnostics ----------------
sim = simulate_reservoir(JutulCase(model, dt, forces; state0 = state0);
    info_level = -1)

ws   = sim.wells
bhp  = get(ws[:Inj], :bhp, nothing)
grat = get(ws[:Inj], :grat, nothing)
gmr  = get(ws[:Inj], :Vapor_mass_rate, nothing)

println("simple_well = ", simple_well)
println("  bhp           = ", bhp)
println("  grat          = ", grat)
println("  gas_mass_rate = ", gmr)
if gmr !== nothing
    println("  any |gas_mass_rate| > 1e-6 : ", any(abs.(gmr) .> 1e-6))
end