Add PVDO table and PVTTableSet return type to PVTExperiments (#28)

* Initial plan

* Add PVDO table and PVTTableSet return type to PVTExperiments

Co-authored-by: moyner <454871+moyner@users.noreply.github.com>

---------

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Co-authored-by: moyner <454871+moyner@users.noreply.github.com>

Author: Copilot

src/PVTExperiments/PVTExperiments.jl

@@ -21,5 +21,5 @@ module PVTExperiments
     include("tables.jl")
     include("interface.jl")
 
-    export generate_pvt_tables
+    export generate_pvt_tables, PVTTableSet
 end

src/PVTExperiments/interface.jl

@@ -5,7 +5,7 @@ High-level interface that generates complete black oil PVT tables from a
 compositional fluid description.
 
 Goes from a fluid sample, reservoir temperature and surface conditions to
-PVTO/PVDG (or PVTG/PVDG) tables + surface densities.
+PVTO/PVDG/PVDO (or PVTG/PVDG/PVDO) tables + surface densities.
 
 # Arguments
 - `eos`: Equation of state (e.g., `GenericCubicEOS(mixture, PengRobinson())`)
@@ -25,13 +25,15 @@ PVTO/PVDG (or PVTG/PVDG) tables + surface densities.
 - `n_pvto`: Number of Rs levels for PVTO. Default: 15
 - `n_pvtg`: Number of pressure levels for PVTG. Default: 15
 - `n_pvdg`: Number of pressure points for PVDG. Default: 20
+- `n_pvdo`: Number of pressure points for PVDO. Default: 20
 - `n_undersaturated`: Undersaturated points per level. Default: 5
 
 # Returns
-A NamedTuple with fields:
+A `PVTTableSet` instance with fields:
 - `pvto`: PVTOTable or `nothing`
 - `pvtg`: PVTGTable or `nothing`
 - `pvdg`: PVDGTable
+- `pvdo`: PVDOTable
 - `surface_densities`: SurfaceDensities
 - `saturation_pressure`: Saturation pressure (Pa)
 - `is_bubblepoint`: Whether the fluid has a bubble point (oil) or dew point (gas)
@@ -46,6 +48,7 @@ function generate_pvt_tables(eos, z, T_res;
         n_pvto = 15,
         n_pvtg = 15,
         n_pvdg = 20,
+        n_pvdo = 20,
         n_undersaturated = 5
     )
     z = collect(Float64, z)
@@ -113,19 +116,38 @@ function generate_pvt_tables(eos, z, T_res;
         T_sc = T_sc
     )
 
+    # PVDO is always generated
+    # For oil systems, use the overall composition
+    # For gas systems, use the oil composition from the first condensation step
+    if is_bp
+        z_oil_dead = z
+    else
+        # Get oil composition near dew point
+        props_dp = flash_and_properties(eos, p_sat * 0.95, T_res, z)
+        z_oil_dead = props_dp.x
+    end
+
+    pvdo_result = pvdo_table(eos, z_oil_dead, T_res;
+        p_range = p_range,
+        n_points = n_pvdo,
+        p_sc = p_sc,
+        T_sc = T_sc
+    )
+
     # Surface densities
     sd = surface_densities(eos, z, T_res;
         p_sc = p_sc,
         T_sc = T_sc,
         separator_stages = separator_stages
     )
 
-    return (
-        pvto = pvto_result,
-        pvtg = pvtg_result,
-        pvdg = pvdg_result,
-        surface_densities = sd,
-        saturation_pressure = p_sat,
-        is_bubblepoint = is_bp
+    return PVTTableSet(
+        pvto_result,
+        pvtg_result,
+        pvdg_result,
+        pvdo_result,
+        sd,
+        p_sat,
+        is_bp
     )
 end

src/PVTExperiments/tables.jl

@@ -313,6 +313,47 @@ function pvtg_table(eos, z, T;
     return PVTGTable(p_values, Rv_sat_values, Rv_sub_table, Bg_table, mu_g_table)
 end
 
+"""
+    pvdo_table(eos, z, T; p_range, n_points, p_sc, T_sc)
+
+Generate a PVDO (dead oil) table.
+
+# Arguments
+- `eos`: Equation of state
+- `z`: Oil composition (mole fractions)
+- `T`: Temperature (K)
+- `p_range`: Pressure range. Default: (50e6, 1e5)
+- `n_points`: Number of pressure points. Default: 20
+- `p_sc`: Standard condition pressure (Pa). Default: 101325.0
+- `T_sc`: Standard condition temperature (K). Default: 288.706
+"""
+function pvdo_table(eos, z, T;
+        p_range = (50e6, 1e5),
+        n_points = 20,
+        p_sc = 101325.0,
+        T_sc = 288.706
+    )
+    z = collect(Float64, z)
+    z ./= sum(z)
+
+    pressures = collect(range(p_range[2], p_range[1], length = n_points))
+    sort!(pressures)
+
+    Bo_arr = zeros(n_points)
+    mu_o_arr = zeros(n_points)
+
+    props_sc = flash_and_properties(eos, p_sc, T_sc, z)
+    V_mol_sc = props_sc.V_mol_l
+
+    for (i, p) in enumerate(pressures)
+        props = flash_and_properties(eos, p, T, z)
+        Bo_arr[i] = props.V_mol_l / V_mol_sc
+        mu_o_arr[i] = props.μ_l
+    end
+
+    return PVDOTable(pressures, Bo_arr, mu_o_arr)
+end
+
 """
     surface_densities(eos, z, T; p_sc, T_sc, separator_stages)
 

src/PVTExperiments/types.jl

@@ -198,6 +198,22 @@ struct PVTGTable
     mu_g::Vector{Vector{Float64}}
 end
 
+"""
+    PVDOTable
+
+Black oil PVDO (dead oil) table.
+
+# Fields
+- `p`: Pressure values (Pa)
+- `Bo`: Oil formation volume factor (m³/m³)
+- `mu_o`: Oil viscosity (Pa·s)
+"""
+struct PVDOTable
+    p::Vector{Float64}
+    Bo::Vector{Float64}
+    mu_o::Vector{Float64}
+end
+
 """
     SurfaceDensities
 
@@ -211,3 +227,27 @@ struct SurfaceDensities
     oil::Float64
     gas::Float64
 end
+
+"""
+    PVTTableSet
+
+Collection of black oil PVT tables generated from a compositional fluid description.
+
+# Fields
+- `pvto`: Live oil table (PVTO) or `nothing` if not applicable
+- `pvtg`: Wet gas table (PVTG) or `nothing` if not applicable
+- `pvdg`: Dry gas table (PVDG)
+- `pvdo`: Dead oil table (PVDO)
+- `surface_densities`: Oil and gas densities at standard conditions
+- `saturation_pressure`: Bubble or dew point pressure (Pa)
+- `is_bubblepoint`: `true` if the fluid has a bubble point (oil), `false` if dew point (gas)
+"""
+struct PVTTableSet
+    pvto::Union{PVTOTable, Nothing}
+    pvtg::Union{PVTGTable, Nothing}
+    pvdg::PVDGTable
+    pvdo::PVDOTable
+    surface_densities::SurfaceDensities
+    saturation_pressure::Float64
+    is_bubblepoint::Bool
+end

src/PVTExperiments/utils.jl

@@ -328,3 +328,27 @@ function Base.show(io::IO, s::SurfaceDensities)
     @printf(io, "  Oil: %.2f kg/m³\n", s.oil)
     @printf(io, "  Gas: %.4f kg/m³\n", s.gas)
 end
+
+function Base.show(io::IO, t::PVDOTable)
+    println(io, "PVDO Table")
+    println(io, "=" ^ 60)
+    @printf(io, "%14s %12s %12s\n", "P (Pa)", "Bo", "μ_o (Pa·s)")
+    println(io, "-" ^ 60)
+    for i in eachindex(t.p)
+        @printf(io, "%14.4e %12.6f %12.4e\n", t.p[i], t.Bo[i], t.mu_o[i])
+    end
+end
+
+function Base.show(io::IO, s::PVTTableSet)
+    println(io, "PVTTableSet")
+    println(io, "=" ^ 60)
+    bp_str = s.is_bubblepoint ? "bubble point (oil)" : "dew point (gas)"
+    @printf(io, "Saturation pressure: %.4e Pa (%.2f bar) [%s]\n",
+        s.saturation_pressure, s.saturation_pressure / 1e5, bp_str)
+    println(io, "-" ^ 60)
+    println(io, "  pvto: ", isnothing(s.pvto) ? "nothing" : "PVTOTable ($(length(s.pvto.Rs)) Rs levels)")
+    println(io, "  pvtg: ", isnothing(s.pvtg) ? "nothing" : "PVTGTable ($(length(s.pvtg.p)) pressure levels)")
+    println(io, "  pvdg: PVDGTable ($(length(s.pvdg.p)) pressure points)")
+    println(io, "  pvdo: PVDOTable ($(length(s.pvdo.p)) pressure points)")
+    show(io, s.surface_densities)
+end

test/pvt_experiments.jl

@@ -147,6 +147,19 @@ import MultiComponentFlash.PVTExperiments as PVTExp
         @test contains(output, "PVDG")
     end
 
+    @testset "PVDO Table" begin
+        table = PVTExp.pvdo_table(eos, z_oil, T_res; n_points = 10)
+        @test table isa PVTExp.PVDOTable
+        @test length(table.p) == 10
+        @test all(table.Bo .> 0)
+        @test all(table.mu_o .> 0)
+        # Test printing
+        io = IOBuffer()
+        show(io, table)
+        output = String(take!(io))
+        @test contains(output, "PVDO")
+    end
+
     @testset "PVTG Table" begin
         table = PVTExp.pvtg_table(eos, z_gas, T_res; n_rv = 5, n_undersaturated = 2)
         @test table isa PVTExp.PVTGTable
@@ -178,18 +191,27 @@ import MultiComponentFlash.PVTExperiments as PVTExp
 
     @testset "High-Level Interface - Oil" begin
         tables = generate_pvt_tables(eos, z_oil, T_res;
-            n_pvto = 5, n_pvdg = 10, n_undersaturated = 2)
+            n_pvto = 5, n_pvdg = 10, n_pvdo = 10, n_undersaturated = 2)
+        @test tables isa PVTExp.PVTTableSet
         @test tables.pvto !== nothing
         @test tables.pvdg !== nothing
+        @test tables.pvdo !== nothing
         @test tables.surface_densities isa PVTExp.SurfaceDensities
         @test tables.saturation_pressure > 0
         @test tables.is_bubblepoint == true
+        # Test printing
+        io = IOBuffer()
+        show(io, tables)
+        output = String(take!(io))
+        @test contains(output, "PVTTableSet")
     end
 
     @testset "High-Level Interface - Gas" begin
         tables = generate_pvt_tables(eos, z_gas, T_res;
-            n_pvtg = 5, n_pvdg = 10)
+            n_pvtg = 5, n_pvdg = 10, n_pvdo = 10)
+        @test tables isa PVTExp.PVTTableSet
         @test tables.pvdg !== nothing
+        @test tables.pvdo !== nothing
         @test tables.surface_densities isa PVTExp.SurfaceDensities
         @test tables.saturation_pressure > 0
     end
@@ -201,9 +223,11 @@ import MultiComponentFlash.PVTExperiments as PVTExp
         ]
         tables = generate_pvt_tables(eos, z_oil, T_res;
             separator_stages = stages,
-            n_pvto = 5, n_pvdg = 10, n_undersaturated = 2)
+            n_pvto = 5, n_pvdg = 10, n_pvdo = 10, n_undersaturated = 2)
+        @test tables isa PVTExp.PVTTableSet
         @test tables.pvto !== nothing
         @test tables.pvdg !== nothing
+        @test tables.pvdo !== nothing
         @test tables.surface_densities isa PVTExp.SurfaceDensities
     end
 end