From f2d879b8049e4e1425363cb2208805ff0e239fb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 12:24:55 +0100 Subject: [PATCH 01/23] Some keywords --- src/InputParser/deckinput/keywords/props.jl | 2 +- src/InputParser/deckinput/keywords/runspec.jl | 38 +++++++++++++++++++ src/InputParser/deckinput/utils.jl | 4 ++ 3 files changed, 43 insertions(+), 1 deletion(-) diff --git a/src/InputParser/deckinput/keywords/props.jl b/src/InputParser/deckinput/keywords/props.jl index 91e2ad16..2801820e 100644 --- a/src/InputParser/deckinput/keywords/props.jl +++ b/src/InputParser/deckinput/keywords/props.jl @@ -179,7 +179,7 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTG}) end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTW}) - tdims = [NaN, NaN, NaN, NaN, NaN] + tdims = [NaN, NaN, NaN, NaN, 0.0] utypes = (:pressure, :liquid_formation_volume_factor, :compressibility, :viscosity, :compressibility) nreg = number_of_tables(outer_data, :pvt) out = [] diff --git a/src/InputParser/deckinput/keywords/runspec.jl b/src/InputParser/deckinput/keywords/runspec.jl index 2c2a76c8..0397d925 100644 --- a/src/InputParser/deckinput/keywords/runspec.jl +++ b/src/InputParser/deckinput/keywords/runspec.jl @@ -52,6 +52,31 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DISGAS}) data["DISGAS"] = true end +function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:STONE1}, Val{:STONE2}, Val{:BAKER1}, Val{:BAKER2}}) + k = unpack_val(v) + parser_message(cfg, outer_data, "$k", PARSER_JUTULDARCY_MISSING_SUPPORT) + data["$k"] = true +end + +function parse_keyword!(data, outer_data, units, cfg, f, v::Val{:MISCSTR}) + parser_message(cfg, outer_data, "MISCSTR", PARSER_JUTULDARCY_MISSING_SUPPORT) + rec = read_record(f) + tdims = [NaN, NaN, NaN]; + l = parse_defaulted_line(rec, tdims) + @assert !isnan(l[1]) + for i in 2:3 + if isnan(l[i]) + l[i] = l[1] + end + end + cm = si_unit(:centi)*si_unit(:meter) + d = si_unit(:dyne) + for i in eachindex(l) + l[i] *= d/cm + end + data["MISCSTR"] = l +end + function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AIM}) data["AIM"] = true end @@ -152,6 +177,12 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELLDIMS}) data["WELLDIMS"] = parse_defaulted_line(rec, tdims) end +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WSEGDIMS}) + rec = read_record(f) + tdims = [0, 1, 1, 0] + data["WSEGDIMS"] = parse_defaulted_line(rec, tdims) +end + function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:VFPPDIMS}) rec = read_record(f) tdims = [0, 0, 0, 0, 0, 0] @@ -169,3 +200,10 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUDIMS}) tdims = [1, 1, 1, 36, 1, 1, 0, 0] data["AQUDIMS"] = parse_defaulted_line(rec, tdims) end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MISCIBLE}) + rec = read_record(f) + tdims = [1, 20, "NONE"] + parser_message(cfg, outer_data, "MISCIBLE", PARSER_JUTULDARCY_MISSING_SUPPORT) + data["MISCIBLE"] = parse_defaulted_line(rec, tdims) +end diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index 71e9fad6..d1d4a352 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -399,6 +399,8 @@ function parse_keyword!(data, outer_data, units, cfg, f, v::Val{T}) where T skip_kw!(:MULTSAVE, 1) skip_kw!(:VECTABLE, 1) skip_kw!(:MULTSAVE, 1) + skip_kw!(:MEMORY, 1) + skip_kw!(:OPTIONS3, 1) skip_kw!(:MULTOUT, 0) skip_kw!(:NOSIM, 0) @@ -407,6 +409,7 @@ function parse_keyword!(data, outer_data, units, cfg, f, v::Val{T}) where T skip_kw!(:CO2STORE, 0, PARSER_JUTULDARCY_MISSING_SUPPORT) skip_kw!(:CO2STOR, 0, PARSER_JUTULDARCY_MISSING_SUPPORT) skip_kw!(:DIFFUSE, 0, PARSER_JUTULDARCY_MISSING_SUPPORT) + skip_kw!(:UNIFSAVE, 0) skip_kw!(:SATOPTS, 1, PARSER_MISSING_SUPPORT) skip_kw!(:EQLOPTS, 1, PARSER_MISSING_SUPPORT) @@ -421,6 +424,7 @@ function parse_keyword!(data, outer_data, units, cfg, f, v::Val{T}) where T skip_kw!(:WPAVE, 1, PARSER_MISSING_SUPPORT) skip_kw!(:VAPPARS, 1, PARSER_MISSING_SUPPORT) skip_kw!(:NETBALAN, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:JFUNC, 1, PARSER_MISSING_SUPPORT) skip_kw!(:TRACER, Inf, PARSER_MISSING_SUPPORT) skip_kw!(:THPRES, Inf, PARSER_MISSING_SUPPORT) From d14e1e1cac8981c12b8579a144ac6fc67356f41d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 12:43:53 +0100 Subject: [PATCH 02/23] Another kw --- .../deckinput/keywords/solution.jl | 21 ++++++++++++++++--- src/InputParser/deckinput/units.jl | 2 +- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/InputParser/deckinput/keywords/solution.jl b/src/InputParser/deckinput/keywords/solution.jl index d342351e..57a7438f 100644 --- a/src/InputParser/deckinput/keywords/solution.jl +++ b/src/InputParser/deckinput/keywords/solution.jl @@ -118,7 +118,22 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RSVD}) data["RSVD"] = out end -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RTEMPVD}) +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ZMFVD}) + n = number_of_tables(outer_data, :equil) + out = [] + ncomp = compositional_number_of_components(outer_data) + for i = 1:n + zmfvd = parse_deck_vector(f) + zmfvd = collect(reshape(zmfvd, ncomp+1, :)') + for j in 1:size(zmfvd, 1) + zmfvd[j, 1] = swap_unit_system(zmfvd[j, 1], units, :length) + end + push!(out, zmfvd) + end + data["ZMFVD"] = out +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Union{Val{:TEMPVD}, Val{:RTEMPVD}}) n = number_of_tables(outer_data, :equil) out = [] for i = 1:n @@ -143,8 +158,8 @@ end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EQUIL}) n = number_of_tables(outer_data, :equil) - def = [0.0, NaN, 0.0, 0.0, 0.0, 0.0, 0, 0, 0] - eunits = (:length, :pressure, :length, :pressure, :length, :pressure, :id, :id, :id) + def = [0.0, NaN, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 1, 0] + eunits = (:length, :pressure, :length, :pressure, :length, :pressure, :id, :id, :id, :id, :id) out = [] for i = 1:n rec = read_record(f) diff --git a/src/InputParser/deckinput/units.jl b/src/InputParser/deckinput/units.jl index f1bcfde1..198ce269 100644 --- a/src/InputParser/deckinput/units.jl +++ b/src/InputParser/deckinput/units.jl @@ -191,7 +191,7 @@ end function swap_unit_system_axes!(x::AbstractVector, systems, eachunit) @assert eltype(eachunit)<:Symbol - @assert length(x) == length(eachunit) + @assert length(x) == length(eachunit) "Recieved vector of length $(length(x)) but units were $(length(eachunit)) long." for i in eachindex(x) x[i] = swap_unit_system(x[i], systems, eachunit[i]) end From edd0b9f8a4cb2ee186c63a3b6e0899749be6cdc1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 13:31:45 +0100 Subject: [PATCH 03/23] FIELDSEP --- .../deckinput/keywords/solution.jl | 26 +++++++++++++++++++ src/InputParser/deckinput/units.jl | 4 +++ 2 files changed, 30 insertions(+) diff --git a/src/InputParser/deckinput/keywords/solution.jl b/src/InputParser/deckinput/keywords/solution.jl index 57a7438f..b2057e14 100644 --- a/src/InputParser/deckinput/keywords/solution.jl +++ b/src/InputParser/deckinput/keywords/solution.jl @@ -169,3 +169,29 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EQUIL}) end data["EQUIL"] = out end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FIELDSEP}) + u_type = deck_unit_system_label(units.from) + if u_type == :field + temp = 60.0 + else + temp = 15.56 + end + # TODO: Can be improved. + p = si_unit(:atm) + + n = number_of_tables(outer_data, :equil) + def = [1, temp, p, 0, 0, 0, 0, 1, NaN, NaN] + eunits = (:id, :relative_temperature, :pressure, :id, :id, :id, :id, :id, :relative_temperature, :pressure) + out = [] + while true + rec = read_record(f) + if length(rec) == 0 + break + end + result = parse_defaulted_line(rec, def) + swap_unit_system_axes!(result, units, eunits) + push!(out, result) + end + data["FIELDSEP"] = out +end diff --git a/src/InputParser/deckinput/units.jl b/src/InputParser/deckinput/units.jl index 198ce269..564141d8 100644 --- a/src/InputParser/deckinput/units.jl +++ b/src/InputParser/deckinput/units.jl @@ -174,6 +174,10 @@ function DeckUnitSystem(sys::Symbol, T = Float64) ) end +function deck_unit_system_label(::DeckUnitSystem{S, T}) where {S, T} + return S +end + function swap_unit_system_axes!(x::AbstractMatrix, systems, eachunit; dim = 2) @assert eltype(eachunit)<:Symbol @assert size(x, dim) == length(eachunit) From 76ad344b392c36d48bcaf24c29d1d19bf807c32b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 14:13:28 +0100 Subject: [PATCH 04/23] Update schedule.jl --- .../deckinput/keywords/schedule.jl | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/src/InputParser/deckinput/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl index d630135b..998130c6 100644 --- a/src/InputParser/deckinput/keywords/schedule.jl +++ b/src/InputParser/deckinput/keywords/schedule.jl @@ -350,3 +350,61 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUANCON}) parser_message(cfg, outer_data, "AQUANCON", PARSER_MISSING_SUPPORT) skip_record(f) end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELSEGS}) + d = "Default" + # TODO: Last two entries for heat capacity / thermal conductivity are not + # properly handled w.r.t. units. + defaults = [d, NaN, 0.0, 1e-5, d, "HFA", "HO", 0.0, 0.0, 0.0, 0.0, 0.0] + utypes = [:id, :length, :length, :volume, :id, :id, :id, :length, :length, :area, :id, :id] + wells = get_wells(outer_data) + rec = read_record(f) + wheader = parse_defaulted_line(rec, defaults) + swap_unit_system_axes!(wheader, units, utypes) + + wname = wheader[1] + wrec = wheader[2:end] + + # Now follows the segments + defaults = [-1, -1, -1, -1, 0.0] + utypes = [:id, :id, :id, :id, :length] + + defaults = [] + utypes = Symbol[] + + function add_entry!(v, u = :id) + push!(defaults, v) + push!(utypes, u) + end + + add_entry!(-1) + add_entry!(-1) + add_entry!(-1) + add_entry!(-1) + add_entry!(0.0, :length) + add_entry!(0.0, :length) + add_entry!(0.0, :length) + add_entry!(0.0, :area) + add_entry!(0.0, :volume) + add_entry!(0.0, :length) + add_entry!(0.0, :length) + add_entry!(0.0, :length) + add_entry!(0.0, :area) + add_entry!(0.0, :id) + add_entry!(0.0, :id) + + segments = [] + while true + rec = read_record(f) + if length(rec) == 0 + break + end + l = parse_defaulted_line(rec, defaults) + swap_unit_system_axes!(l, units, utypes) + push!(segments, l) + end + if !haskey(data, "WELSEGS") + data["WELSEGS"] = Dict{String, Any}() + end + data["WELSEGS"][wname] = (header = wheader, segments = segments) +end From 5b5b95fcda10a057ed09e23289aea2583bb2b78a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 14:30:36 +0100 Subject: [PATCH 05/23] A bunch of smaller parser fixes --- src/InputParser/deckinput/keywords/props.jl | 2 +- .../deckinput/keywords/schedule.jl | 49 +++++++++++++++++-- src/input_simulation/data_input.jl | 10 ++-- 3 files changed, 50 insertions(+), 11 deletions(-) diff --git a/src/InputParser/deckinput/keywords/props.jl b/src/InputParser/deckinput/keywords/props.jl index 2801820e..906a21d1 100644 --- a/src/InputParser/deckinput/keywords/props.jl +++ b/src/InputParser/deckinput/keywords/props.jl @@ -19,7 +19,7 @@ end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EOS}) rec = read_record(f) - data["EOS"] = only(rec) + data["EOS"] = only(parse_defaulted_line(rec, ["PR"])) end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NCOMPS}) diff --git a/src/InputParser/deckinput/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl index 998130c6..1706eaac 100644 --- a/src/InputParser/deckinput/keywords/schedule.jl +++ b/src/InputParser/deckinput/keywords/schedule.jl @@ -48,7 +48,7 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COMPORD}) out = outer_data["SCHEDULE"]["COMPORD"] for co in compord well, val = co - out[name] = val + out[well] = val end end @@ -218,7 +218,7 @@ end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WEFAC}) parser_message(cfg, outer_data, "WEFAC", PARSER_JUTULDARCY_PARTIAL_SUPPORT) - defaults = ["Default", 1.0] + defaults = ["Default", 1.0, "YES"] wells = get_wells(outer_data) parsed = parse_defaulted_group_well(f, defaults, wells, 1) push_and_create!(data, "WEFAC", parsed) @@ -363,12 +363,10 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELSEGS}) swap_unit_system_axes!(wheader, units, utypes) wname = wheader[1] + @assert wname != d wrec = wheader[2:end] # Now follows the segments - defaults = [-1, -1, -1, -1, 0.0] - utypes = [:id, :id, :id, :id, :length] - defaults = [] utypes = Symbol[] @@ -408,3 +406,44 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELSEGS}) end data["WELSEGS"][wname] = (header = wheader, segments = segments) end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COMPSEGS}) + rec = read_record(f) + wname = only(parse_defaulted_line(rec, ["Default"])) + @assert wname != "Default" + + defaults = [] + utypes = Symbol[] + + function add_entry!(v, u = :id) + push!(defaults, v) + push!(utypes, u) + end + + add_entry!(-1) + add_entry!(-1) + add_entry!(-1) + add_entry!(-1) + add_entry!(0.0, :length) + add_entry!(0.0, :length) + add_entry!("Default") + add_entry!(-1) + add_entry!(0.0, :length) + add_entry!(0.0, :length) + add_entry!(0) + + segments = [] + while true + rec = read_record(f) + if length(rec) == 0 + break + end + l = parse_defaulted_line(rec, defaults) + swap_unit_system_axes!(l, units, utypes) + push!(segments, l) + end + if !haskey(data, "COMPSEGS") + data["COMPSEGS"] = Dict{String, Any}() + end + data["COMPSEGS"][wname] = segments +end diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl index 25d8590c..00d3432c 100644 --- a/src/input_simulation/data_input.jl +++ b/src/input_simulation/data_input.jl @@ -712,15 +712,15 @@ function parse_physics_types(datafile) mp = MolecularProperty.(mw, p_c, T_c, V_c, acf) mixture = MultiComponentMixture(mp, A_ij = A_ij, names = cnames) if haskey(props, "EOS") - eos_str = lowercase(props["EOS"]) - if eos_str == "pr" + eos_str = uppercase(props["EOS"]) + if eos_str == "PR" eos_type = PengRobinson() - elseif eos_str == "srk" + elseif eos_str == "SRK" eos_type = SoaveRedlichKwong() - elseif eos_str == "rk" + elseif eos_str == "RK" eos_type = RedlichKwong() else - @assert eos_str == "zj" + @assert eos_str == "ZJ" "Unexpected EOS $eos_str: Should be one of PR, SRK, RK, ZJ" eos_type = ZudkevitchJoffe() end else From 900743008691cbfa14fb4b3990b50db025ee0458 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 14:42:53 +0100 Subject: [PATCH 06/23] Handle compositional water injectors --- src/input_simulation/data_input.jl | 53 +++++++++++++++++------------- 1 file changed, 30 insertions(+), 23 deletions(-) diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl index 00d3432c..b65fec1d 100644 --- a/src/input_simulation/data_input.jl +++ b/src/input_simulation/data_input.jl @@ -1123,35 +1123,42 @@ function select_injector_mixture_spec(sys::CompositionalSystem, name, streams, t props = eos.mixture.properties rho_s = JutulDarcy.reference_densities(sys) phases = JutulDarcy.get_phases(sys) - mix = Float64[] - rho = 0.0 - # Well stream will be molar fracitons. offset = Int(has_other_phase(sys)) ncomp = number_of_components(sys) + # Well stream will be molar fractions. mix = zeros(Float64, ncomp) - stream_id = streams.wells[name] - stream = streams.streams[stream_id] + if uppercase(type) == "WATER" + has_other_phase(sys) || throw(ArgumentError("Cannot have WATER injector without water phase.")) + mix[1] = 1.0 + rho = rho_s[1] + else + if !haskey(streams.wells, name) + throw(ArgumentError("Well $name does not have a stream declared with well type $type.")) + end + stream_id = streams.wells[name] + stream = streams.streams[stream_id] - ϵ = MultiComponentFlash.MINIMUM_COMPOSITION - z = map(z_i -> max(z_i, ϵ), stream.mole_fractions) - z /= sum(z) - cond = stream.cond + ϵ = MultiComponentFlash.MINIMUM_COMPOSITION + z = map(z_i -> max(z_i, ϵ), stream.mole_fractions) + z /= sum(z) + cond = stream.cond + + z_mass = map( + (z_i, prop) -> max(z_i*prop.mw, ϵ), + z, props + ) + z_mass /= sum(z_mass) + for i in 1:ncomp + mix[i+offset] = z_mass[i] + end + @assert sum(mix) ≈ 1.0 "Sum of mixture was $(sum(mix)) != 1 for mole mixture $(z) as mass $z_mass" - z_mass = map( - (z_i, prop) -> max(z_i*prop.mw, ϵ), - z, props - ) - z_mass /= sum(z_mass) - for i in 1:ncomp - mix[i+offset] = z_mass[i] + flash_cond = (p = cond.p, T = cond.T, z = z) + flash = MultiComponentFlash.flashed_mixture_2ph(eos, flash_cond) + rho_l, rho_v = MultiComponentFlash.mass_densities(eos, cond.p, cond.T, flash) + S_l, S_v = MultiComponentFlash.phase_saturations(eos, cond.p, cond.T, flash) + rho = S_l*rho_l + S_v*rho_v end - @assert sum(mix) ≈ 1.0 "Sum of mixture was $(sum(mix)) != 1 for mole mixture $(z) as mass $z_mass" - - flash_cond = (p = cond.p, T = cond.T, z = z) - flash = MultiComponentFlash.flashed_mixture_2ph(eos, flash_cond) - rho_l, rho_v = MultiComponentFlash.mass_densities(eos, cond.p, cond.T, flash) - S_l, S_v = MultiComponentFlash.phase_saturations(eos, cond.p, cond.T, flash) - rho = S_l*rho_l + S_v*rho_v return (rho, mix) end From f2ea1b0abb4a108812ec2491ccb511ad27dd6be1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 19:09:58 +0100 Subject: [PATCH 07/23] Improve behavior of special/edit kw --- src/InputParser/deckinput/keywords/special.jl | 43 ++++++++++++++----- 1 file changed, 33 insertions(+), 10 deletions(-) diff --git a/src/InputParser/deckinput/keywords/special.jl b/src/InputParser/deckinput/keywords/special.jl index 369140a4..c7b37c53 100644 --- a/src/InputParser/deckinput/keywords/special.jl +++ b/src/InputParser/deckinput/keywords/special.jl @@ -119,10 +119,13 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:OPERATE}) op_prm2 = parsed[11] IJK = get_box_indices(outer_data, il, iu, jl, ju, kl, ku) - if !haskey(data, target) + operation_target = get_operation_section(outer_data, target) + operation_source = get_operation_section(outer_data, source) + if ismissing(operation_target) data[target] = zeros(dims) + operation_target = data[target] end - apply_operate!(data[target], data[source], IJK, op, op_prm1, op_prm2) + apply_operate!(operation_target, operation_source, IJK, op, op_prm1, op_prm2) # On to the next one. rec = read_record(f) end @@ -186,11 +189,21 @@ function apply_multiply!(vals, src, IX, dims) vals[I, J, K] .*= src end +function get_operation_section(outer_data, kw) + for (k, data) in pairs(outer_data) + if data isa AbstractDict && haskey(data, kw) + return data[kw] + end + end + return missing +end + function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTIPLY}) # TODO: Merge shared code with COPY rec = read_record(f) - l, u = outer_data["GRID"]["CURRENT_BOX"] - dims = outer_data["GRID"]["cartDims"] + grid = outer_data["GRID"] + l, u = grid["CURRENT_BOX"] + dims = grid["cartDims"] il = l[1] iu = u[1] @@ -213,8 +226,16 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTIPLY}) ju = parsed[6] kl = parsed[7] ku = parsed[8] - @assert haskey(data, dst) "Unable to apply MULTIPLY to non-declared field $dst" - apply_multiply!(data[dst], factor, (il, iu), (jl, ju), (kl, ku), dims) + target = get_operation_section(outer_data, dst) + if ismissing(target) + if dst == "PORV" && haskey(grid, "PORO") + # TODO: Bit of a hack + target = grid["PORO"] + else + throw(ArgumentError("Unable to apply MULTIPLY to non-declared field $dst")) + end + end + apply_multiply!(target, factor, (il, iu), (jl, ju), (kl, ku), dims) rec = read_record(f) end end @@ -254,7 +275,11 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EQUALS}) dst = parsed[1] constval = parsed[2] @assert dst != "Default" - if haskey(data, dst) + target = get_operation_section(outer_data, dst) + if ismissing(target) + # TODO: Different keywords go in different spots... + data[dst] = fill(constval, dims...) + else # Box can be kept. il = parsed[3] iu = parsed[4] @@ -262,9 +287,7 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EQUALS}) ju = parsed[6] kl = parsed[7] ku = parsed[8] - apply_equals!(data[dst], constval, (il, iu), (jl, ju), (kl, ku), dims) - else - data[dst] = fill(constval, dims...) + apply_equals!(target, constval, (il, iu), (jl, ju), (kl, ku), dims) end rec = read_record(f) end From 1852c8812b2b65b87fb17bd0912e40e2335b5356 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 19:20:25 +0100 Subject: [PATCH 08/23] Update utils.jl --- src/InputParser/deckinput/utils.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index d1d4a352..2e4de022 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -469,7 +469,10 @@ end function next_keyword!(f) m = nothing while isnothing(m) && !eof(f) - line = readline(f) + line = strip(readline(f)) + if line == "/" + continue + end m = keyword_start(line) end return m From eca9a05cc4827670dfd6cadf3f9031f31ac69ee8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 24 Jan 2024 19:28:21 +0100 Subject: [PATCH 09/23] RTEMP varies with EOSNUM --- src/InputParser/deckinput/keywords/solution.jl | 13 +++++++++---- src/InputParser/deckinput/utils.jl | 9 +++++++-- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/InputParser/deckinput/keywords/solution.jl b/src/InputParser/deckinput/keywords/solution.jl index b2057e14..09a0295c 100644 --- a/src/InputParser/deckinput/keywords/solution.jl +++ b/src/InputParser/deckinput/keywords/solution.jl @@ -91,10 +91,15 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PRESSURE}) end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RTEMP}) - rec = read_record(f) - result = parse_defaulted_line(rec, [NaN]) - swap_unit_system!(result, units, :relative_temperature) - data["RTEMP"] = result + nreg = number_of_tables(outer_data, :eos) + out = Float64[] + for i in 1:nreg + rec = read_record(f) + result = parse_defaulted_line(rec, [NaN]) + swap_unit_system!(result, units, :relative_temperature) + push!(out, only(result)) + end + data["RTEMP"] = out end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RS}) diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index 2e4de022..99efa25d 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -495,9 +495,14 @@ function number_of_tables(outer_data, t::Symbol) else return 1 end - else - error(":$t is not known") + elseif t == :eos + if haskey(rs, "TABDIMS") + return rs["TABDIMS"][9] + else + return 1 + end end + error(":$t is not known") end function compositional_number_of_components(outer_data) From 66a9e6eaf6bb77911fac7a41ecb25d19f4f7028b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 10:00:53 +0100 Subject: [PATCH 10/23] A couple of small parser fixes --- src/InputParser/deckinput/keywords/runspec.jl | 18 ++++++++++++++++++ src/InputParser/deckinput/parser.jl | 8 ++++++-- src/InputParser/deckinput/utils.jl | 4 ++++ test/parser.jl | 1 + 4 files changed, 29 insertions(+), 2 deletions(-) diff --git a/src/InputParser/deckinput/keywords/runspec.jl b/src/InputParser/deckinput/keywords/runspec.jl index 0397d925..e9648928 100644 --- a/src/InputParser/deckinput/keywords/runspec.jl +++ b/src/InputParser/deckinput/keywords/runspec.jl @@ -134,6 +134,24 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TABDIMS}) data["TABDIMS"] = parse_defaulted_line(rec, tdims) end +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:UDADIMS}) + rec = read_record(f) + tdims = [0, 0, 100] + data["UDADIMS"] = parse_defaulted_line(rec, tdims) +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:UDQDIMS}) + rec = read_record(f) + tdims = [16, 16, 0, 0, 0, 0, 0, 0, 0, 0, "N"] + data["UDQDIMS"] = parse_defaulted_line(rec, tdims) +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EXTRAPMS}) + rec = read_record(f) + tdims = [0] + data["EXTRAPMS"] = only(parse_defaulted_line(rec, tdims)) +end + function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FAULTDIM}) rec = read_record(f) tdims = [0]; diff --git a/src/InputParser/deckinput/parser.jl b/src/InputParser/deckinput/parser.jl index 3f14f2e8..7bb5ddf8 100644 --- a/src/InputParser/deckinput/parser.jl +++ b/src/InputParser/deckinput/parser.jl @@ -107,6 +107,7 @@ function parse_data_file!(outer_data, filename, data = outer_data; units::Union{Symbol, Nothing} = :si, warn_parsing = true, warn_feature = true, + basedir = missing, silent = false, is_outer = true, input_units::Union{Symbol, Nothing} = nothing, @@ -131,7 +132,9 @@ function parse_data_file!(outer_data, filename, data = outer_data; unit_systems = get_unit_system_pair(input_units, units) end filename = realpath(filename) - basedir = dirname(filename) + if ismissing(basedir) + basedir = dirname(filename) + end f = open(filename, "r") cfg = ParserVerbosityConfig( @@ -176,7 +179,7 @@ function parse_data_file!(outer_data, filename, data = outer_data; readline(f) end include_path = clean_include_path(basedir, next) - parser_message(cfg, outer_data, "$m", "Including file: $include_path") + parser_message(cfg, outer_data, "$m", "Including file: $include_path. Basedir: $basedir with INCLUDE = $next)") parse_data_file!( outer_data, include_path, data, verbose = verbose, @@ -185,6 +188,7 @@ function parse_data_file!(outer_data, filename, data = outer_data; silent = silent, sections = sections, skip = skip, + basedir = basedir, skip_mode = skip_mode, is_outer = false, units = units, diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index 99efa25d..6dbe9484 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -553,6 +553,10 @@ function clean_include_path(basedir, include_file_name) if startswith(include_file_name, "./") include_file_name = include_file_name[3:end] end + pos = findlast(" /", include_file_name) + if !isnothing(pos) + include_file_name = include_file_name[1:pos[1]-1] + end include_path = joinpath(basedir, include_file_name) return include_path end diff --git a/test/parser.jl b/test/parser.jl index 00e701f3..c831ee1e 100644 --- a/test/parser.jl +++ b/test/parser.jl @@ -6,6 +6,7 @@ import JutulDarcy.InputParser: clean_include_path, parse_defaulted_line @test clean_include_path("/some/path", " 'MYFILE'") == joinpath("/some/path", "MYFILE") @test clean_include_path("/some/path", " ./MYFILE") == joinpath("/some/path", "MYFILE") @test clean_include_path("/some/path", " './MYFILE'") == joinpath("/some/path", "MYFILE") + @test clean_include_path("/some/path", " 'INCLUDE/file.txt' / (Some comment)") == "/some/path/INCLUDE/file.txt" @test parse_defaulted_line("3.0 2* 7", [1.0, 2, 3, 4]) == [3.0, 2, 3, 7] @test parse_defaulted_line("2.0", [1.0, 2, 3, 4]) == [2.0, 2, 3, 4] From b1cb9fa82f916baa78d99f0e835f44a962e12b17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 10:07:27 +0100 Subject: [PATCH 11/23] Another parser fix --- src/InputParser/deckinput/utils.jl | 5 +++++ test/parser.jl | 1 + 2 files changed, 6 insertions(+) diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index 6dbe9484..bae128b4 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -544,6 +544,11 @@ function table_region(outer_data, t::Symbol; active = nothing) end function clean_include_path(basedir, include_file_name) + m = match(r"'[^']*'", include_file_name) + if !isnothing(m) + # Strip away anything that isn't ' + include_file_name = m.match[2:end-1] + end include_file_name = strip(include_file_name) if startswith(include_file_name, "./") include_file_name = include_file_name[3:end] diff --git a/test/parser.jl b/test/parser.jl index c831ee1e..9aed7c8a 100644 --- a/test/parser.jl +++ b/test/parser.jl @@ -7,6 +7,7 @@ import JutulDarcy.InputParser: clean_include_path, parse_defaulted_line @test clean_include_path("/some/path", " ./MYFILE") == joinpath("/some/path", "MYFILE") @test clean_include_path("/some/path", " './MYFILE'") == joinpath("/some/path", "MYFILE") @test clean_include_path("/some/path", " 'INCLUDE/file.txt' / (Some comment)") == "/some/path/INCLUDE/file.txt" + @test clean_include_path("/some/path", " 'INCLUDE/file.txt'/(Some comment)") == "/some/path/INCLUDE/file.txt" @test parse_defaulted_line("3.0 2* 7", [1.0, 2, 3, 4]) == [3.0, 2, 3, 7] @test parse_defaulted_line("2.0", [1.0, 2, 3, 4]) == [2.0, 2, 3, 4] From a48e7780a191906c311890d50842793958333b97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 10:23:37 +0100 Subject: [PATCH 12/23] Update grid.jl --- src/InputParser/deckinput/keywords/grid.jl | 28 ++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/src/InputParser/deckinput/keywords/grid.jl b/src/InputParser/deckinput/keywords/grid.jl index 4ea3ed01..b5975feb 100644 --- a/src/InputParser/deckinput/keywords/grid.jl +++ b/src/InputParser/deckinput/keywords/grid.jl @@ -59,9 +59,33 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GRIDUNIT}) partial_parse!(data, outer_data, units, cfg, f, :GRIDUNIT) end +function check_unit(unit_str, units, kw) + ref = uppercase("$(deck_unit_system_label(units.from))") + u = uppercase(unit_str) + if u != ref + @warn "Unit mismatch in $kw: Was $u but RUNSPEC declared $ref" + end +end + function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FILEUNIT}) - # TODO: This needs to be handled - partial_parse!(data, outer_data, units, cfg, f, :FILEUNIT) + rec = strip(only(read_record(f))) + check_unit(rec, units, "FILEUNIT") + data["FILEUNIT"] = rec +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GRIDUNIT}) + rec = read_record(f) + tdims = ["Default", "MAP"] + v = parse_defaulted_line(rec, tdims) + check_unit(v[1], units, "GRIDUNIT") + data["GRIDUNIT"] = v +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MAPUNITS}) + rec = read_record(f) + tdims = ["Default"] + v = parse_defaulted_line(rec, tdims) + data["MAPUNITS"] = only(v) end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GDORIENT}) From 50331ec027388f0da64bde2fce17a0825fb1c477 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 10:29:51 +0100 Subject: [PATCH 13/23] MAXVALUE/MINVALUE --- src/InputParser/deckinput/keywords/special.jl | 49 +++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/src/InputParser/deckinput/keywords/special.jl b/src/InputParser/deckinput/keywords/special.jl index c7b37c53..be8388ff 100644 --- a/src/InputParser/deckinput/keywords/special.jl +++ b/src/InputParser/deckinput/keywords/special.jl @@ -302,3 +302,52 @@ function apply_equals!(target, constval, I, J, K, dims) end end end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MAXVALUE}) + edit_apply_clamping!(f, outer_data, min) +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MINVALUE}) + edit_apply_clamping!(f, outer_data, max) +end + +function edit_apply_clamping!(f, outer_data, FUNCTION) + rec = read_record(f) + l, u = outer_data["GRID"]["CURRENT_BOX"] + dims = outer_data["GRID"]["cartDims"] + + il = l[1] + iu = u[1] + jl = l[2] + ju = u[2] + kl = l[3] + ku = u[3] + + while length(rec) > 0 + d = "Default" + parsed = parse_defaulted_line(rec, [d, 1.0, il, iu, jl, ju, kl, ku]) + dst = parsed[1] + lim = parsed[2] + @assert dst != "Default" + target = get_operation_section(outer_data, dst) + # Box can be kept. + il = parsed[3] + iu = parsed[4] + jl = parsed[5] + ju = parsed[6] + kl = parsed[7] + ku = parsed[8] + apply_equals!(target, lim, (il, iu), (jl, ju), (kl, ku), dims) + rec = read_record(f) + end +end + +function apply_minmax!(target, lim, I, J, K, dims, F) + for i in I[1]:I[2] + for j in J[1]:J[2] + for k in K[1]:K[2] + target[i, j, k] = F(target[i, j, k], lim) + end + end + end +end From aa79faee15596891f0b03314ebeda1983e1ee26f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 10:50:41 +0100 Subject: [PATCH 14/23] Some more keyword fixes --- src/InputParser/deckinput/keywords/grid.jl | 3 ++- src/InputParser/deckinput/keywords/props.jl | 18 ++++++++++++++++++ src/InputParser/deckinput/keywords/schedule.jl | 5 +++++ src/InputParser/deckinput/keywords/solution.jl | 11 +++++++++++ src/InputParser/deckinput/utils.jl | 8 +++++++- 5 files changed, 43 insertions(+), 2 deletions(-) diff --git a/src/InputParser/deckinput/keywords/grid.jl b/src/InputParser/deckinput/keywords/grid.jl index b5975feb..5613beac 100644 --- a/src/InputParser/deckinput/keywords/grid.jl +++ b/src/InputParser/deckinput/keywords/grid.jl @@ -63,7 +63,8 @@ function check_unit(unit_str, units, kw) ref = uppercase("$(deck_unit_system_label(units.from))") u = uppercase(unit_str) if u != ref - @warn "Unit mismatch in $kw: Was $u but RUNSPEC declared $ref" + # Commented out due to missing logic (e.g. METRIC should equal METRES) + # @warn "Unit mismatch in $kw: Was $u but RUNSPEC declared $ref" end end diff --git a/src/InputParser/deckinput/keywords/props.jl b/src/InputParser/deckinput/keywords/props.jl index 906a21d1..9dc0b3f3 100644 --- a/src/InputParser/deckinput/keywords/props.jl +++ b/src/InputParser/deckinput/keywords/props.jl @@ -302,3 +302,21 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ROCKOPTS}) parser_message(cfg, outer_data, "ROCKOPTS", PARSER_JUTULDARCY_MISSING_SUPPORT) data["ROCKOPTS"] = parsed end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUTAB}) + num = outer_data["RUNSPEC"]["AQUDIMS"][3]-1 + for i in 1:num + skip_record(f) + end + parser_message(cfg, outer_data, "AQUTAB", PARSER_MISSING_SUPPORT) +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUCT}) + skip_record(f) + parser_message(cfg, outer_data, "AQUCT", PARSER_MISSING_SUPPORT) +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUANCON}) + skip_record(f) + parser_message(cfg, outer_data, "AQUANCON", PARSER_MISSING_SUPPORT) +end diff --git a/src/InputParser/deckinput/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl index 1706eaac..97375053 100644 --- a/src/InputParser/deckinput/keywords/schedule.jl +++ b/src/InputParser/deckinput/keywords/schedule.jl @@ -447,3 +447,8 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COMPSEGS}) end data["COMPSEGS"][wname] = segments end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:UDQ}) + skip_record(f) + parser_message(cfg, outer_data, "UDQ", PARSER_MISSING_SUPPORT) +end diff --git a/src/InputParser/deckinput/keywords/solution.jl b/src/InputParser/deckinput/keywords/solution.jl index 09a0295c..546cdb20 100644 --- a/src/InputParser/deckinput/keywords/solution.jl +++ b/src/InputParser/deckinput/keywords/solution.jl @@ -123,6 +123,17 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RSVD}) data["RSVD"] = out end +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PBVD}) + n = number_of_tables(outer_data, :equil) + out = [] + for i = 1:n + rs = parse_deck_matrix(f) + swap_unit_system_axes!(rs, units, (:length, :pressure)) + push!(out, rs) + end + data["PBVD"] = out +end + function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ZMFVD}) n = number_of_tables(outer_data, :equil) out = [] diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index bae128b4..9ee7ac9a 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -173,7 +173,13 @@ function parse_deck_matrix_line!(data::Vector{T}, seg, n) where T @assert m == n "Expected $n was $m" end for d in seg - push!(data, Parsers.parse(T, d)) + if d == raw"1*" + # Defaulted... + @assert T == Float64 + push!(data, NaN) + else + push!(data, Parsers.parse(T, d)) + end end return n end From a6b5ec5063bd6f32ae26ed17f12e34f3bc8d0782 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 10:57:25 +0100 Subject: [PATCH 15/23] Handle defaulted saturation tables --- src/input_simulation/mrst_input.jl | 1 + src/io.jl | 10 ++++++++++ src/types.jl | 1 + 3 files changed, 12 insertions(+) diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl index 320f1390..c9daa740 100644 --- a/src/input_simulation/mrst_input.jl +++ b/src/input_simulation/mrst_input.jl @@ -435,6 +435,7 @@ function deck_pc(props; oil, water, gas, satnum = nothing) pc = -pc[ix] s = s[ix] end + s, pc = saturation_table_handle_defaults(s, pc) interp_ow = get_1d_interpolator(s, pc) push!(PC, interp_ow) end diff --git a/src/io.jl b/src/io.jl index 30b1c7f2..5d375dad 100644 --- a/src/io.jl +++ b/src/io.jl @@ -10,6 +10,16 @@ function table_to_relperm(swof; swcon = 0.0, first_label = :w, second_label = :o return (krw, krow) end +function saturation_table_handle_defaults(s, f) + if any(isnan, f) + # NaN values are removed due to INPUT file shenanigans + ix = findall(!isnan, f) + s = s[ix] + f = f[ix] + end + return (s, f) +end + function add_missing_endpoints(s, kr) if s[1] > 0.0 s = vcat(0.0, s) diff --git a/src/types.jl b/src/types.jl index 6cdbfd48..4803b6df 100644 --- a/src/types.jl +++ b/src/types.jl @@ -238,6 +238,7 @@ specified. """ function PhaseRelativePermeability(s, k; label = :w, connate = s[1], epsilon = 1e-16) msg(i) = "k = $(k[i]) at entry $i corresponding to saturation $(s[i])" + s, k = saturation_table_handle_defaults(s, k) for i in eachindex(s) if i == 1 if s[1] == 0.0 From 5c4738045ba9ec9eb6c0cdc9321051f862a29295 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 14:29:12 +0100 Subject: [PATCH 16/23] Improve Rs/Rv init --- src/init/init.jl | 40 ++++++++++++++++++++++++++-------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/src/init/init.jl b/src/init/init.jl index 59b76968..a945c4d9 100644 --- a/src/init/init.jl +++ b/src/init/init.jl @@ -30,18 +30,23 @@ function equilibriate_state(model, contacts, datum_depth = missing, datum_pressu sat = init[:Saturations] pressure = init[:Pressure] nph, nc = size(sat) - if ismissing(rs) - Rs = zeros(nc) - else - Rs = rs.(pts) + + if has_disgas(sys) + if ismissing(rs) + Rs = zeros(nc) + else + Rs = rs.(pts) + end + init[:Rs] = Rs end - if ismissing(rv) - Rv = zeros(nc) - else - Rv = rv.(pts) + if has_vapoil(sys) + if ismissing(rv) + Rv = zeros(nc) + else + Rv = rv.(pts) + end + init[:Rv] = Rv end - init[:Rs] = Rs - init[:Rv] = Rv end return init end @@ -255,10 +260,17 @@ function parse_state0_equil(model, datafile) end if disgas - @assert haskey(sol, "RSVD") - rsvd = sol["RSVD"][ereg] - z = rsvd[:, 1] - Rs = rsvd[:, 2] + if haskey(sol, "RSVD") + rsvd = sol["RSVD"][ereg] + z = rsvd[:, 1] + Rs = rsvd[:, 2] + else + @assert haskey(sol, "PBVD") + pbvd = sol["PBVD"][ereg] + z = pbvd[:, 1] + pb = vec(pbvd[:, 2]) + Rs = sys.rs_max.(pb) + end rs = Jutul.LinearInterpolant(z, Rs) else rs = missing From 670a99d897908b34c7321e7a7aa7460ed4e1d4b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 25 Jan 2024 15:17:12 +0100 Subject: [PATCH 17/23] Improve a warning message --- src/utils.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 01a3ee99..2751f3df 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1121,7 +1121,9 @@ function reservoir_transmissibility(d::DataDomain; version = :xyz) T_hf[i] = abs(T_hf_i) end if neg_count > 0 - @warn "Replaced $neg_count negative half-transmissibilities (out of $(length(T_hf))) with their absolute value." + tran_tot = length(T_hf) + perc = round(100*neg_count/tran_tot, digits = 2) + @warn "Replaced $neg_count negative half-transmissibilities (out of $tran_tot, $perc%) with their absolute value." end if haskey(d, :net_to_gross) # Net to gross applies to vertical trans only From 2f23f85433c1cfe6e6e38c46718f22a83349e309 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 26 Jan 2024 12:51:35 +0100 Subject: [PATCH 18/23] Add WIP support for KValue flash (depends o MultiComponentFlash.jl branch) --- src/multicomponent/variables/flash.jl | 61 ++++++++++++++++++++++++--- src/types.jl | 4 +- 2 files changed, 58 insertions(+), 7 deletions(-) diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index 11405d89..74b26612 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -1,5 +1,5 @@ mutable struct InPlaceFlashBuffer - z + z::Vector{Float64} forces function InPlaceFlashBuffer(n) z = zeros(n) @@ -75,16 +75,17 @@ function initialize_variable_ad!(state, model, pvar::FlashResults, symb, npartia end @jutul_secondary function update_flash!(flash_results, fr::FlashResults, model, Pressure, Temperature, OverallMoleFractions, ix) - flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, nothing, ix) + eos = model.system.equation_of_state + flash_entity_loop!(flash_results, fr, model, eos, Pressure, Temperature, OverallMoleFractions, nothing, ix) end @jutul_secondary function update_flash!(flash_results, fr::FlashResults, model::LVCompositionalModel3Phase, Pressure, Temperature, OverallMoleFractions, ImmiscibleSaturation, ix) - flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, ImmiscibleSaturation, ix) + eos = model.system.equation_of_state + flash_entity_loop!(flash_results, fr, model, eos, Pressure, Temperature, OverallMoleFractions, ImmiscibleSaturation, ix) end -function flash_entity_loop!(flash_results, fr, model, Pressure, Temperature, OverallMoleFractions, sw, ix) +function flash_entity_loop!(flash_results, fr, model, eos, Pressure, Temperature, OverallMoleFractions, sw, ix) storage, m, buffers = fr.storage, fr.method, fr.update_buffer - eos = model.system.equation_of_state S, buf = thread_buffers(storage, buffers) update_flash_buffer!(buf, eos, Pressure, Temperature, OverallMoleFractions) @@ -232,3 +233,53 @@ two_phase_pre!(S, P, T, Z, x, y, V, eos, c) = V return (convert(AD, Z_L), convert(AD, Z_V), convert(AD, V), phase_state) end + +function flash_entity_loop!(flash_results, fr, model, eos::KValuesEOS, Pressure, Temperature, OverallMoleFractions, sw, ix) + storage, m, buffers = fr.storage, fr.method, fr.update_buffer + S, buf = thread_buffers(storage, buffers) + z_buf = buf.z + for i in ix + P = Pressure[i] + T = Temperature[i] + Z = @view OverallMoleFractions[:, i] + + fr = flash_results[i] + flash_results[i] = k_value_flash!(fr, eos, P, T, Z, z_buf) + end +end + +function k_value_flash!(result::FR, eos, P, T, Z, z) where FR + Num_t = Base.promote_type(typeof(P), typeof(T), eltype(Z)) + @. z = max(value(Z), 1e-8) + # Conditions + c = (p = value(P), T = value(T), z = z) + c_ad = (p = P, T = T, z = Z) + K_ad = eos.K_values_evaluator(c_ad) + K = result.K + x = result.liquid.mole_fractions + y = result.vapor.mole_fractions + + @. K = value(K_ad) + V = flash_2ph!(nothing, K, eos, c) + + pure_liquid = V <= 0.0 + pure_vapor = V >= 1.0 + if pure_vapor || pure_liquid + if pure_vapor + phase_state = MultiComponentFlash.single_phase_v + else + phase_state = MultiComponentFlash.single_phase_l + end + @. x = Z + @. y = Z + V = convert(Num_t, V) + else + phase_state = MultiComponentFlash.two_phase_lv + @. x = liquid_mole_fraction(Z, K, vapor_frac) + @. y = vapor_mole_fraction(x, K) + # TODO: Handle partial derivatives of RR here. + end + Z_L = Z_V = convert(Num_t, NaN) + new_result = FlashedMixture2Phase(phase_state, K, V, x, y, Z_L, Z_V) + return new_result::FR +end diff --git a/src/types.jl b/src/types.jl index 4803b6df..1077394c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -56,7 +56,7 @@ function MultiPhaseCompositionalSystemLV(equation_of_state, phases = (LiquidPhas end only(findall(isequal(LiquidPhase()), phases)) only(findall(isequal(VaporPhase()), phases)) - MultiPhaseCompositionalSystemLV{typeof(equation_of_state), T, O, typeof(reference_densities), N}(phases, c, equation_of_state, reference_densities) + return MultiPhaseCompositionalSystemLV{typeof(equation_of_state), T, O, typeof(reference_densities), N}(phases, c, equation_of_state, reference_densities) end function Base.show(io::IO, sys::MultiPhaseCompositionalSystemLV) @@ -94,7 +94,7 @@ end Set up a standard black-oil system. Keyword arguments `rs_max` and `rv_max` can either be nothing or callable objects / functions for the maximum Rs and Rv as a function of pressure. `phases` can be specified together with -`reference_densities` for each phase. +`reference_densities` for each phase. NOTE: For the black-oil model, the reference densities significantly impact many aspects of the PVT behavior. These should generally be set consistently with the From 74d0d8eaf9bf5b950f61ea5d49ea2b557ba77e8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 26 Jan 2024 14:52:46 +0100 Subject: [PATCH 19/23] Add partials to rachford rice via inverse function theorem --- src/multicomponent/variables/flash.jl | 36 ++++++++++++++++++++++++--- src/multicomponent/wells.jl | 20 ++++++++++++++- 2 files changed, 52 insertions(+), 4 deletions(-) diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index 74b26612..0cc929b6 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -137,6 +137,10 @@ function update_flash_buffer!(buf, eos, Pressure, Temperature, OverallMoleFracti end end +function update_flash_buffer!(buf, eos::KValuesEOS, Pressure, Temperature, OverallMoleFractions) + return nothing +end + function internal_flash!(f, S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) @inline function immiscible_sat(::Nothing, i) return 0.0 @@ -275,11 +279,37 @@ function k_value_flash!(result::FR, eos, P, T, Z, z) where FR V = convert(Num_t, V) else phase_state = MultiComponentFlash.two_phase_lv - @. x = liquid_mole_fraction(Z, K, vapor_frac) + V = add_derivatives_to_vapor_fraction_rachford_rice(V, K_ad, Z, K, z) + @. x = liquid_mole_fraction(Z, K, V) @. y = vapor_mole_fraction(x, K) - # TODO: Handle partial derivatives of RR here. end - Z_L = Z_V = convert(Num_t, NaN) + Z_L = Z_V = convert(Num_t, 1.0) new_result = FlashedMixture2Phase(phase_state, K, V, x, y, Z_L, Z_V) return new_result::FR end + +function add_derivatives_to_vapor_fraction_rachford_rice(V::Float64, K, z, K_val = value(K), z_val = value(z)) + Zt = eltype(z) + Kt = eltype(K) + T = Base.promote_type(Zt, Kt) + if T != Float64 + N = length(z) + V0 = V + V = convert(T, V) + ∂V = V.partials + if Kt == T + for i in 1:N + dK_i = MultiComponentFlash.objectiveRR_dK(V0, K_val, z_val, i) + ∂V += K[i].partials*dK_i + end + end + if Zt == T + for i in 1:N + dz_i = MultiComponentFlash.objectiveRR_dz(V0, K_val, z_val, i) + ∂V += z[i].partials*dz_i + end + end + V = T(V0, ∂V) + end + return V +end diff --git a/src/multicomponent/wells.jl b/src/multicomponent/wells.jl index 88fd0a38..c2debb09 100644 --- a/src/multicomponent/wells.jl +++ b/src/multicomponent/wells.jl @@ -3,11 +3,29 @@ function flash_wellstream_at_surface(var, well_model, system::S, state, rhoS, co nc = MultiComponentFlash.number_of_components(eos) z = SVector{nc}(state.OverallMoleFractions[:, 1]) flash, _, surface_moles = get_separator_intermediate_storage(var, system, z, 2) + S_l, S_v, rho_l, rho_v = separator_flash(flash, system, surface_moles, eos, cond, z) + return compositional_surface_densities(state, system, S_l, S_v, rho_l, rho_v) +end + +function separator_flash(flash, system, surface_moles, eos, cond, z) result = separator_flash!(flash, eos, cond, z) rho_l, rho_v = mass_densities(eos, cond.p, cond.T, result) S_l, S_v = phase_saturations(eos, cond.p, cond.T, result) @assert S_l + S_v ≈ 1.0 - return compositional_surface_densities(state, system, S_l, S_v, rho_l, rho_v) + return (S_l, S_v, rho_l, rho_v) +end + +function separator_flash(flash, system, surface_moles, eos::KValuesEOS, cond, z) + if has_other_phase(system) + _, rho_l, rho_v = reference_densities(system) + else + rho_l, rho_v = reference_densities(system) + end + # TODO Clean up, correct, derivatives? + V = flash_2ph(eos, (p = cond.p, T = cond.T, z = z)) + S_l = 1.0 - V + S_v = V + return (S_l, S_v, rho_l, rho_v) end Base.@propagate_inbounds function multisegment_well_perforation_flux!(out, sys::CompositionalSystem, state_res, state_well, rhoS, conn) From 37f3ac8e3602adfa8dbab3e748daf51ece04b109 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 28 Jan 2024 18:52:00 +0100 Subject: [PATCH 20/23] Bump version + Compat --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 469b1f77..0a08716d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "JutulDarcy" uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a" authors = ["Olav Møyner "] -version = "0.2.15" +version = "0.2.16" [deps] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" @@ -52,7 +52,7 @@ Krylov = "0.9.1" LinearAlgebra = "1" LoopVectorization = "0.12.115" MAT = "0.10.3" -MultiComponentFlash = "1.1.3" +MultiComponentFlash = "1.1.9" OrderedCollections = "1.6.2" Parsers = "2.7.1" PartitionedArrays = "0.3.2" From eaf15376f6d03b9c20f5cc734a6a74431f3d6636 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 28 Jan 2024 19:03:09 +0100 Subject: [PATCH 21/23] Remove duplicate keywords --- src/InputParser/deckinput/keywords/grid.jl | 10 ---------- src/InputParser/deckinput/keywords/props.jl | 4 ---- 2 files changed, 14 deletions(-) diff --git a/src/InputParser/deckinput/keywords/grid.jl b/src/InputParser/deckinput/keywords/grid.jl index 5613beac..d7d7fef3 100644 --- a/src/InputParser/deckinput/keywords/grid.jl +++ b/src/InputParser/deckinput/keywords/grid.jl @@ -49,16 +49,6 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COORDSYS}) parser_message(cfg, outer_data, "COORDSYS", PARSER_MISSING_SUPPORT) end -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MAPUNITS}) - # TODO: This needs to be handled - partial_parse!(data, outer_data, units, cfg, f, :MAPUNITS) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GRIDUNIT}) - # TODO: This needs to be handled - partial_parse!(data, outer_data, units, cfg, f, :GRIDUNIT) -end - function check_unit(unit_str, units, kw) ref = uppercase("$(deck_unit_system_label(units.from))") u = uppercase(unit_str) diff --git a/src/InputParser/deckinput/keywords/props.jl b/src/InputParser/deckinput/keywords/props.jl index 9dc0b3f3..1ae75700 100644 --- a/src/InputParser/deckinput/keywords/props.jl +++ b/src/InputParser/deckinput/keywords/props.jl @@ -316,7 +316,3 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUCT}) parser_message(cfg, outer_data, "AQUCT", PARSER_MISSING_SUPPORT) end -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUANCON}) - skip_record(f) - parser_message(cfg, outer_data, "AQUANCON", PARSER_MISSING_SUPPORT) -end From 21b2439d0c8bbe9aadcc81606fdd3b59b8f51623 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 28 Jan 2024 19:37:19 +0100 Subject: [PATCH 22/23] Fix to parser tests --- test/parser.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/parser.jl b/test/parser.jl index 9aed7c8a..6179fe81 100644 --- a/test/parser.jl +++ b/test/parser.jl @@ -6,8 +6,8 @@ import JutulDarcy.InputParser: clean_include_path, parse_defaulted_line @test clean_include_path("/some/path", " 'MYFILE'") == joinpath("/some/path", "MYFILE") @test clean_include_path("/some/path", " ./MYFILE") == joinpath("/some/path", "MYFILE") @test clean_include_path("/some/path", " './MYFILE'") == joinpath("/some/path", "MYFILE") - @test clean_include_path("/some/path", " 'INCLUDE/file.txt' / (Some comment)") == "/some/path/INCLUDE/file.txt" - @test clean_include_path("/some/path", " 'INCLUDE/file.txt'/(Some comment)") == "/some/path/INCLUDE/file.txt" + @test clean_include_path("/some/path", " 'file.txt' / (Some comment)") == joinpath("/some/path", "file.txt") + @test clean_include_path("/some/path", " 'file.txt'/(Some comment)") == joinpath("/some/path", "file.txt") @test parse_defaulted_line("3.0 2* 7", [1.0, 2, 3, 4]) == [3.0, 2, 3, 7] @test parse_defaulted_line("2.0", [1.0, 2, 3, 4]) == [2.0, 2, 3, 4] From c7c95b525fb49d47fcdd807a4464fb27cf5243b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 28 Jan 2024 20:32:24 +0100 Subject: [PATCH 23/23] Bump MCF compat for bugfix for Julia 1.9 and below --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0a08716d..1bcba205 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ Krylov = "0.9.1" LinearAlgebra = "1" LoopVectorization = "0.12.115" MAT = "0.10.3" -MultiComponentFlash = "1.1.9" +MultiComponentFlash = "1.1.10" OrderedCollections = "1.6.2" Parsers = "2.7.1" PartitionedArrays = "0.3.2"