diff --git a/Project.toml b/Project.toml index 1bcba205..f07a8445 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +GeoEnergyIO = "3b1dd628-313a-45bb-9d8d-8f3c48dcb5d4" Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -17,7 +18,6 @@ LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MAT = "23992714-dd62-5051-b70f-ba57cb901cac" MultiComponentFlash = "35e5bd01-9722-4017-9deb-64a5d32478ff" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -46,6 +46,7 @@ DelimitedFiles = "1.9.1" DocStringExtensions = "0.9.3" ForwardDiff = "0.10.30" GLMakie = "0.9.2" +GeoEnergyIO = "1" HYPRE = "1.4.0" Jutul = "0.2.20" Krylov = "0.9.1" @@ -54,7 +55,6 @@ LoopVectorization = "0.12.115" MAT = "0.10.3" MultiComponentFlash = "1.1.10" OrderedCollections = "1.6.2" -Parsers = "2.7.1" PartitionedArrays = "0.3.2" Polyester = "0.6.11, 0.7.3" PrecompileTools = "1.0.1" diff --git a/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl b/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl index 707efe1e..3bafa294 100644 --- a/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl +++ b/ext/JutulDarcyHYPREExt/JutulDarcyHYPREExt.jl @@ -4,6 +4,10 @@ module JutulDarcyHYPREExt using Jutul using SparseArrays using PrecompileTools + using TimerOutputs + + timeit_debug_enabled() = Jutul.timeit_debug_enabled() + include("cpr.jl") # @compile_workload begin diff --git a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl index 46b688c4..593d156e 100644 --- a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl +++ b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl @@ -1,6 +1,6 @@ module JutulDarcyPartitionedArraysExt using Jutul, JutulDarcy - using PrecompileTools + using PrecompileTools, TimerOutputs # Specific dependencies using PartitionedArrays, MPI, HYPRE using LinearAlgebra @@ -15,6 +15,8 @@ module JutulDarcyPartitionedArraysExt increment_pressure!, set_dp!, correct_residual!, CPRStorage, CPRPreconditioner, number_of_components + timeit_debug_enabled() = Jutul.timeit_debug_enabled() + function setup_reservoir_simulator_parray( case::JutulCase, backend::PArrayBackend; diff --git a/src/InputParser/InputParser.jl b/src/InputParser/InputParser.jl deleted file mode 100644 index 2c14fed7..00000000 --- a/src/InputParser/InputParser.jl +++ /dev/null @@ -1,6 +0,0 @@ -module InputParser - using Parsers, DelimitedFiles, Jutul, OrderedCollections, Dates, LinearAlgebra - export parse_data_file - - include("deckinput/parser.jl") -end diff --git a/src/InputParser/deckinput/keywords/grid.jl b/src/InputParser/deckinput/keywords/grid.jl deleted file mode 100644 index d7d7fef3..00000000 --- a/src/InputParser/deckinput/keywords/grid.jl +++ /dev/null @@ -1,238 +0,0 @@ -function parse_and_set_grid_data!(data, outer_data, units, cfg, f, k; unit = :id, T = Float64, default = zero(T)) - bdims = get_boxdims(outer_data) - cdims = get_cartdims(outer_data) - vals = parse_grid_vector(f, bdims, T) - if unit != :id - vals = swap_unit_system!(vals, units, Val(unit)) - end - skey = "$k" - if bdims == cdims - data[skey] = vals - else - if !haskey(data, skey) - data[skey] = fill(default, cdims) - end - d = data[skey] - @assert size(d) == cdims - I, J, K = get_box_indices(outer_data) - d[I, J, K] = vals - end -end - -function finish_current_section!(data, units, cfg, outer_data, ::Val{:GRID}) - if !haskey(data, "MINPV") - io = IOBuffer("1e-6\n/\n") - parse_keyword!(data, outer_data, units, cfg, io, Val(:MINPV)) - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GRIDFILE}) - rec = read_record(f) - tdims = [0, 1]; - data["GRIDFILE"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Union{Val{:MINPORV}, Val{:MINPV}}) - rec = read_record(f) - tdims = [1e-6]; - rec = parse_defaulted_line(rec, tdims) - zcorn = swap_unit_system!(rec, units, :volume) - data["MINPV"] = only(rec) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:INIT}) - data["INIT"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COORDSYS}) - read_record(f) - parser_message(cfg, outer_data, "COORDSYS", PARSER_MISSING_SUPPORT) -end - -function check_unit(unit_str, units, kw) - ref = uppercase("$(deck_unit_system_label(units.from))") - u = uppercase(unit_str) - if u != 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 - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{: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}) - # TODO: This needs to be handled - partial_parse!(data, outer_data, units, cfg, f, :GDORIENT) -end - -function partial_parse!(data, outer_data, units, cfg, f, k::Symbol) - rec = read_record(f) - parser_message(cfg, outer_data, "$k", PARSER_MISSING_SUPPORT) - data["$k"] = rec -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MAPAXES}) - rec = parse_deck_vector(f, Float64) - data["MAPAXES"] = rec -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COORD}) - coord = parse_deck_vector(f, Float64) - coord = swap_unit_system!(coord, units, Val(:length)) - data["COORD"] = coord -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ZCORN}) - zcorn = parse_deck_vector(f, Float64) - zcorn = swap_unit_system!(zcorn, units, Val(:length)) - data["ZCORN"] = zcorn -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:PERMX}, Val{:PERMY}, Val{:PERMZ}}) - k = unpack_val(v) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, k, unit = :permeability) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Val{:MULTREGT}) - parser_message(cfg, outer_data, "MULTREGT", PARSER_MISSING_SUPPORT) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{ - Val{:MULTX}, Val{:MULTY}, Val{:MULTZ}, - Val{Symbol("MULTX-")}, Val{Symbol("MULTY-")}, Val{Symbol("MULTZ-")} - } - ) - k = unpack_val(v) - parser_message(cfg, outer_data, "$k", PARSER_JUTULDARCY_MISSING_SUPPORT) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, k, unit = :id) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Val{:MULTPV}) - k = unpack_val(v) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, k) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:PRATIO}, Val{:BIOTCOEF}}) - k = unpack_val(v) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, k) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:YMODULE}}) - k = unpack_val(v) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, k, unit = :gigapascal) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:POELCOEF}, Val{:THELCOEF}, Val{:THERMEXR}, Val{:THCONR}}) - k = unpack_val(v) - vals = parse_grid_vector(f, get_cartdims(outer_data), Float64) - parser_message(cfg, outer_data, "$k", PARSER_PARTIAL_SUPPORT) - data["$k"] = vals -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:FIPNUM}, Val{:PVTNUM}, Val{:SATNUM}, Val{:EQLNUM}, Val{:ROCKNUM}}) - k = unpack_val(v) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, k, T = Int) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PORO}) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, :PORO) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NTG}) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, :NTG) -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:DX}, Val{:DY}, Val{:DZ}}) - k = unpack_val(v) - data["$k"] = parse_grid_vector(f, get_cartdims(outer_data), Float64) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TOPS}) - tops = parse_deck_vector(f, Float64) - data["TOPS"] = swap_unit_system!(tops, units, Val(:length)) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DIMENS}) - rec = read_record(f) - to_int = x -> Parsers.parse(Int, x) - d = to_int.(filter!(x -> length(x) > 0, split(only(rec), DECK_SPLIT_REGEX))) - data["DIMENS"] = d - set_cartdims!(outer_data, d) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SPECGRID}) - rec = read_record(f) - tdims = [1, 1, 1, 1, "F"] - data["SPECGRID"] = parse_defaulted_line(rec, tdims) - set_cartdims!(outer_data, data["SPECGRID"][1:3]) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PINCH}) - rec = read_record(f) - tdims = [0.001, "GAP", Inf, "TOPBOT", "TOP"] - parser_message(cfg, outer_data, "PINCH", PARSER_JUTULDARCY_PARTIAL_SUPPORT) - data["PINCH"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FAULTS}) - read_record - tdims = ["NAME", -1, -1, -1, -1, -1, -1, "XYZ_IJK"] - if !haskey(data, "FAULTS") - data["FAULTS"] = Dict{String, Any}() - end - faults = data["FAULTS"] - while true - rec = read_record(f) - if length(rec) == 0 - break - end - parsed = parse_defaulted_line(rec, tdims, required_num = length(tdims), keyword = "FAULTS") - name = parsed[1] - flt = ( - i = parsed[2]:parsed[3], - j = parsed[4]:parsed[5], - k = parsed[6]:parsed[7], - direction = parsed[8] - ) - if haskey(faults, name) - push!(faults[name], flt) - else - faults[name] = [flt] - end - end -end - - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTFLT}) - d = "*" - tdims = [d, 1.0, 1.0] - faults = outer_data["GRID"]["FAULTS"] - parser_message(cfg, outer_data, "MULTFLT", PARSER_JUTULDARCY_MISSING_SUPPORT) - out = parse_defaulted_group_well(f, tdims, faults); - if !haskey(data, "MULTFLT") - data["MULTFLT"] = Dict{String, Tuple{Float64, Float64}}() - end - for (name, mul_t, mul_d) in out - data["MULTFLT"][name] = (mul_t, mul_d) - end -end diff --git a/src/InputParser/deckinput/keywords/keywords.jl b/src/InputParser/deckinput/keywords/keywords.jl deleted file mode 100644 index 131e8512..00000000 --- a/src/InputParser/deckinput/keywords/keywords.jl +++ /dev/null @@ -1,8 +0,0 @@ -include("runspec.jl") -include("grid.jl") -include("props.jl") -include("regions.jl") -include("schedule.jl") -include("solution.jl") -include("summary.jl") -include("special.jl") \ No newline at end of file diff --git a/src/InputParser/deckinput/keywords/props.jl b/src/InputParser/deckinput/keywords/props.jl deleted file mode 100644 index 1ae75700..00000000 --- a/src/InputParser/deckinput/keywords/props.jl +++ /dev/null @@ -1,318 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:PROPS}) - -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTPROPS}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FULLIMP}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:CNAMES}) - n = compositional_number_of_components(outer_data) - templ = fill("", n) - rec = read_record(f) - data["CNAMES"] = parse_defaulted_line(rec, templ) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EOS}) - rec = read_record(f) - data["EOS"] = only(parse_defaulted_line(rec, ["PR"])) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NCOMPS}) - rec = read_record(f) - data["NCOMPS"] = parse(Int, only(rec)) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:OMEGAA}) - parser_message(cfg, outer_data, "OMEGAA", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["OMEGAA"] = parse_deck_vector(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:OMEGAB}) - parser_message(cfg, outer_data, "OMEGAB", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["OMEGAB"] = parse_deck_vector(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:LBCCOEF}) - parser_message(cfg, outer_data, "LBCCOEF", PARSER_JUTULDARCY_MISSING_SUPPORT) - rec = read_record(f) - defaults = [0.1023, 0.023364, 0.058533, -0.040758, 0.0093324] - data["LBCCOEF"] = parse_defaulted_line(rec, defaults) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:STCOND}) - std = parse_deck_vector(f) - @assert length(std) == 2 - swap_unit_system_axes!(std, units, [:relative_temperature, :pressure]) - data["STCOND"] = std -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BIC}) - n = compositional_number_of_components(outer_data) - bic = parse_deck_vector(f) - @assert length(bic) == n*(n-1)÷2 "Bad length for BIC input." - m = zeros(n, n) - ix = 1 - for i in 1:n - for j in 1:(i-1) - m[i, j] = bic[ix] - ix += 1 - end - end - data["BIC"] = Symmetric(collect(m')) -end - -function parse_compositional_helper!(f, outer_data, data, k) - n = compositional_number_of_components(outer_data) - val = parse_deck_vector(f) - @assert length(val) == n "One $k should be provided per component (expected $n, was $(length(val)))." - return val -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ACF}) - data["ACF"] = parse_compositional_helper!(f, outer_data, data, "ACF") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ZCRIT}) - parser_message(cfg, outer_data, "ZCRIT", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["ZCRIT"] = parse_compositional_helper!(f, outer_data, data, "ZCRIT") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SSHIFT}) - parser_message(cfg, outer_data, "SSHIFT", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["SSHIFT"] = parse_compositional_helper!(f, outer_data, data, "SSHIFT") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PARACHOR}) - parser_message(cfg, outer_data, "PARACHOR", PARSER_MISSING_SUPPORT) - # TODO: Units. - data["PARACHOR"] = parse_compositional_helper!(f, outer_data, data, "PARACHOR") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:VCRITVIS}) - parser_message(cfg, outer_data, "VCRITVIS", PARSER_MISSING_SUPPORT) - # TODO: Units. - data["VCRITVIS"] = parse_compositional_helper!(f, outer_data, data, "VCRITVIS") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ZI}) - parser_message(cfg, outer_data, "ZI", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["ZI"] = parse_compositional_helper!(f, outer_data, data, "ZI") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PCRIT}) - p_c = parse_compositional_helper!(f, outer_data, data, "PCRIT") - swap_unit_system!(p_c, units, :pressure) - data["PCRIT"] = p_c -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TCRIT}) - p_c = parse_compositional_helper!(f, outer_data, data, "TCRIT") - swap_unit_system!(p_c, units, :absolute_temperature) - data["TCRIT"] = p_c -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MW}) - mw = parse_compositional_helper!(f, outer_data, data, "MW") - swap_unit_system!(mw, units, :molar_mass) - data["MW"] = mw -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:VCRIT}) - V = parse_compositional_helper!(f, outer_data, data, "VCRIT") - swap_unit_system!(V, units, :critical_volume) - data["VCRIT"] = V -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:SWOF}, Val{:SGOF}}) - k = unpack_val(v) - sat_tab = parse_saturation_table(f, outer_data) - for tab in sat_tab - swap_unit_system_axes!(tab, units, (:identity, :identity, :identity, :pressure)) - end - data["$k"] = sat_tab -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:SWFN}, Val{:SGFN}}) - k = unpack_val(v) - sat_tab = parse_saturation_table(f, outer_data) - for tab in sat_tab - swap_unit_system_axes!(tab, units, (:identity, :identity, :pressure)) - end - data["$k"] = sat_tab -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:SOF2}, Val{:SOF3}}) - k = unpack_val(v) - sat_tab = parse_saturation_table(f, outer_data) - data["$k"] = sat_tab -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVDG}) - pvdg = parse_dead_pvt_table(f, outer_data) - for tab in pvdg - swap_unit_system_axes!(tab, units, (:pressure, :gas_formation_volume_factor, :viscosity)) - end - data["PVDG"] = pvdg -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTO}) - pvto = parse_live_pvt_table(f, outer_data) - for tab in pvto - swap_unit_system_axes!(tab["data"], units, (:pressure, :liquid_formation_volume_factor, :viscosity)) - swap_unit_system!(tab["key"], units, :u_rs) - end - data["PVTO"] = pvto -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTG}) - pvtg = parse_live_pvt_table(f, outer_data) - for tab in pvtg - swap_unit_system_axes!(tab["data"], units, (:u_rv, :gas_formation_volume_factor, :viscosity)) - swap_unit_system!(tab["key"], units, :pressure) - end - data["PVTG"] = pvtg -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTW}) - tdims = [NaN, NaN, NaN, NaN, 0.0] - utypes = (:pressure, :liquid_formation_volume_factor, :compressibility, :viscosity, :compressibility) - nreg = number_of_tables(outer_data, :pvt) - out = [] - for i = 1:nreg - rec = read_record(f) - t = parse_defaulted_line(rec, tdims) - swap_unit_system_axes!(t, units, utypes) - @assert all(isfinite, t) "PVTW cannot be defaulted, found defaulted record in region $i" - push!(out, t) - end - data["PVTW"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVCDO}) - tdims = [NaN, NaN, NaN, NaN, NaN] - utypes = (:pressure, :liquid_formation_volume_factor, :compressibility, :viscosity, :compressibility) - nreg = number_of_tables(outer_data, :pvt) - out = [] - for i = 1:nreg - rec = read_record(f) - t = parse_defaulted_line(rec, tdims) - swap_unit_system_axes!(t, units, utypes) - @assert all(isfinite, t) "PVCDO cannot be defaulted, found defaulted record in region $i" - push!(out, t) - end - data["PVCDO"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVDO}) - pvdo = parse_dead_pvt_table(f, outer_data) - for tab in pvdo - swap_unit_system_axes!(tab, units, (:pressure, :liquid_formation_volume_factor, :viscosity)) - end - data["PVDO"] = pvdo -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ROCK}) - tdims = [NaN, NaN, NaN, NaN, NaN, NaN] - utypes = [:pressure, :compressibility, :compressibility, :compressibility, :id, :id] - out = [] - nreg = number_of_tables(outer_data, :pvt) - for i = 1:nreg - rec = read_record(f) - l = parse_defaulted_line(rec, tdims) - swap_unit_system_axes!(l, units, utypes) - push!(out, l) - end - data["ROCK"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DIFFC}) - parser_message(cfg, outer_data, "DIFFC", PARSER_MISSING_SUPPORT) - nreg = number_of_tables(outer_data, :pvt) - for i = 1:nreg - rec = read_record(f) - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SPECROCK}) - parser_message(cfg, outer_data, "SPECROCK", PARSER_MISSING_SUPPORT) - nreg = number_of_tables(outer_data, :saturation) - for i = 1:nreg - rec = read_record(f) - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COMPS}) - rec = read_record(f) - ncomp = only(parse_defaulted_line(rec, [0])) - data["COMPS"] = ncomp -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SCALECRS}) - rec = read_record(f) - scale = only(parse_defaulted_line(rec, ["NO"])) - scale = uppercase(scale) - if scale == "Y" - scale = "YES" - elseif scale == "N" - scale = "NO" - end - if scale == "YES" || scale == "NO" - data["SCALECRS"] = scale - else - error("SCALECRS must be one of the following: Y, N, YES, NO") - end -end - - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DENSITY}) - tdims = [NaN, NaN, NaN] - nreg = number_of_tables(outer_data, :pvt) - out = [] - for i = 1:nreg - rec = read_record(f) - t = parse_defaulted_line(rec, tdims) - swap_unit_system!(t, units, :density) - push!(out, t) - end - data["DENSITY"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Union{Val{:RSCONST}, Val{:RSCONSTT}}) - rec = read_record(f) - # TODO: This is missing units. - tdims = [NaN, NaN] - parsed = parse_defaulted_line(rec, tdims, required_num = length(tdims), keyword = "RSCONST") - parser_message(cfg, outer_data, "RSCONST", PARSER_PARTIAL_SUPPORT) - data["RSCONST"] = parsed -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FILLEPS}) - data["FILLEPS"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ROCKOPTS}) - rec = read_record(f) - tdims = ["PRESSURE", "NOSTORE", "PVTNUM", "DEFLATION"] - parsed = parse_defaulted_line(rec, tdims) - 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 - diff --git a/src/InputParser/deckinput/keywords/regions.jl b/src/InputParser/deckinput/keywords/regions.jl deleted file mode 100644 index 253fd45f..00000000 --- a/src/InputParser/deckinput/keywords/regions.jl +++ /dev/null @@ -1,9 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:REGIONS}) - -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTREGS}) - parser_message(cfg, outer_data, "RPTREGS", PARSER_MISSING_SUPPORT) - rec = read_record(f) - # TODO: Process the record. -end diff --git a/src/InputParser/deckinput/keywords/runspec.jl b/src/InputParser/deckinput/keywords/runspec.jl deleted file mode 100644 index e9648928..00000000 --- a/src/InputParser/deckinput/keywords/runspec.jl +++ /dev/null @@ -1,227 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:RUNSPEC}) - -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NOECHO}) - # Do nothing -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ECHO}) - # Do nothing -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MESSAGES}) - parser_message(cfg, outer_data, "MESSAGES", PARSER_MISSING_SUPPORT) - rec = read_record(f) - # TODO: Process the record. -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:START}) - rec = read_record(f) - tdims = [1, "JAN", 1970, "00:00:00"]; - start = parse_defaulted_line(rec, tdims) - data["START"] = convert_date_kw(start) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TITLE}) - m = next_keyword!(f) - data["TITLE"] = m -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:METRIC}) - data["METRIC"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FIELD}) - data["FIELD"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WATER}) - data["WATER"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:OIL}) - data["OIL"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GAS}) - data["GAS"] = true -end - -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 - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:HWELLS}) - data["HWELLS"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PRCORR}) - data["PRCORR"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MONITOR}) - data["MONITOR"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:THERMAL}) - data["THERMAL"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MECH}) - data["MECH"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:VAPOIL}) - data["VAPOIL"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:UNIFOUT}) - data["UNIFOUT"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:UNIFIN}) - data["UNIFIN"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:CPR}) - data["CPR"] = true -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NUMRES}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTSAVE}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TABDIMS}) - rec = read_record(f) - tdims = [1, 1, 20, 20, 1, 20, 20, 1, - 1, -1, 10, 1, -1, 0, 0, -1, - 10, 10, 10, -1, 5, 5, 5, 0, -1]; - # TODO: Special logic for -1 entries - 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]; - data["FAULTDIM"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GRIDOPTS}) - rec = read_record(f) - tdims = ["NO", 0, 0]; - data["GRIDOPTS"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ENDSCALE}) - rec = read_record(f) - parser_message(cfg, outer_data, "ENDSCALE", PARSER_JUTULDARCY_MISSING_SUPPORT) - tdims = ["NODIR", "REVERS", 1, 20, 0]; - data["ENDSCALE"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EQLDIMS}) - rec = read_record(f) - tdims = [1, 100, 50, 1, 50]; - data["EQLDIMS"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MSGFILE}) - rec = read_record(f) - tdims = [1]; - data["MSGFILE"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:REGDIMS}) - rec = read_record(f) - tdims = [1, 1, 0, 0, 0, 1, 0, 0, 0]; - data["REGDIMS"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELLDIMS}) - rec = read_record(f) - tdims = [0, 0, 0, 0, 5, 10, 5, 4, 3, 0, 1, 1, 10, 201] - 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] - data["VFPPDIMS"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:VFPIDIMS}) - rec = read_record(f) - tdims = [0, 0, 0] - data["VFPIDIMS"] = parse_defaulted_line(rec, tdims) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUDIMS}) - rec = read_record(f) - 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/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl deleted file mode 100644 index 97375053..00000000 --- a/src/InputParser/deckinput/keywords/schedule.jl +++ /dev/null @@ -1,454 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:SCHEDULE}) - compord = outer_data["SCHEDULE"]["COMPORD"] - for k in keys(get_wells(outer_data)) - if !haskey(compord, k) - compord[k] = "TRACK" - end - end -end - -function welspecs_to_well(ws) - name = ws[1]::AbstractString - # TODO: Parse more fields. - w = ( - head = (ws[3], ws[4]), - ref_depth = ws[5], - preferred_phase = ws[6], - drainage_radius = ws[7] - ) - return (name, w) -end - -function get_wells(outer_data) - return outer_data["SCHEDULE"]["WELSPECS"] -end - -function get_well(outer_data, name) - return get_wells(outer_data)[name] -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELSPECS}) - d = "Default" - defaults = [d, d, -1, -1, NaN, d, 0.0, "STD", "SHUT", "YES", 0, "SEG", 0, d, d, "STD", 0] - utypes = [:id, :id, :id, :id, :length, :id, :length, :id, :id, :id, :id, :id, :id, :id, :id, :id, :id] - wspecs = parse_defaulted_group(f, defaults) - wells = get_wells(outer_data) - for ws in wspecs - swap_unit_system_axes!(ws, units, utypes) - name, well = welspecs_to_well(ws) - # @assert !haskey(wells, name) "Well $name is already defined." - wells[name] = well - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COMPORD}) - d = "Default" - defaults = [d, "TRACK"] - compord = parse_defaulted_group(f, defaults) - out = outer_data["SCHEDULE"]["COMPORD"] - for co in compord - well, val = co - out[well] = val - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COMPDAT}) - d = "Default" - defaults = [d, -1, -1, -1, -1, "OPEN", -1, NaN, NaN, NaN, 0.0, -1.0, "Z", -1.0] - wells = get_wells(outer_data) - compdat = parse_defaulted_group_well(f, defaults, wells, 1) - # Unit conversion - utypes = identity_unit_vector(defaults) - utypes[8] = :transmissibility - utypes[9] = :length - utypes[10] = :Kh - utypes[12] = :time_over_volume - utypes[14] = :length - for cd in compdat - swap_unit_system_axes!(cd, units, utypes) - end - data["COMPDAT"] = compdat -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELOPEN}) - d = "Default" - defaults = [d, "OPEN", -1, -1, -1, -1, -1] - wells = get_wells(outer_data) - welopen = parse_defaulted_group_well(f, defaults, wells, 1) - parser_message(cfg, outer_data, "WELTARG", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["WELOPEN"] = welopen -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WCONPROD}) - d = "Default" - defaults = [ - d, "OPEN", d, -1.0, -1.0, - -1.0, -1.0, -1.0, -1.0, -1.0, - 0, 0.0, NaN, NaN, NaN, - NaN, NaN, NaN, NaN, NaN - ] - wells = get_wells(outer_data) - wconprod = parse_defaulted_group_well(f, defaults, wells, 1) - for kw in wconprod - for i in 4:10 - # Handle defaults - if kw[i] < 0 - kw[i] = Inf - end - end - end - utypes = identity_unit_vector(defaults) - utypes[4] = :liquid_rate_surface - utypes[5] = :liquid_rate_surface - utypes[6] = :gas_rate_surface - utypes[7] = :liquid_rate_surface - utypes[8] = :liquid_rate_reservoir - utypes[9] = :pressure - for wp in wconprod - swap_unit_system_axes!(wp, units, utypes) - end - data["WCONPROD"] = wconprod -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WCONHIST}) - d = "Default" - defaults = [d, "OPEN", d, 0.0, 0.0, 0.0, 0, d, 0.0, 0.0, 0.0] - wells = get_wells(outer_data) - out = parse_defaulted_group_well(f, defaults, wells, 1) - utypes = identity_unit_vector(defaults) - utypes[4] = :liquid_rate_surface - utypes[5] = :liquid_rate_surface - utypes[6] = :gas_rate_surface - utypes[9] = :pressure - utypes[10] = :pressure - utypes[11] = :gas_rate_surface - for o in out - swap_unit_system_axes!(o, units, utypes) - end - data["WCONHIST"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WCONINJ}) - d = "Default" - defaults = [d, d, "OPEN", d, Inf, Inf, 0.0, "NONE", -1, Inf, 0.0, 0.0] - wells = get_wells(outer_data) - out = parse_defaulted_group_well(f, defaults, wells, 1) - utypes = identity_unit_vector(defaults) - utypes[6] = :liquid_rate_surface - utypes[9] = :pressure - utypes[10] = :pressure - for o in out - w = get_well(outer_data, o[1]) - if w.preferred_phase == "GAS" - utypes[5] = :gas_rate_surface - utypes[12] = :u_rv - else - utypes[5] = :liquid_rate_surface - if w.preferred_phase == "OIL" - utypes[12] = :u_rs - end - end - swap_unit_system_axes!(o, units, utypes) - end - data["WCONINJ"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WCONINJE}) - d = "Default" - defaults = [d, d, "OPEN", d, -1.0, -1.0, NaN, -1.0, 0, 0.0, 0.0, 0, 0, 0] - wells = get_wells(outer_data) - out = parse_defaulted_group_well(f, defaults, wells, 1) - utypes = identity_unit_vector(defaults) - utypes[6] = :liquid_rate_reservoir - utypes[7] = :pressure - utypes[8] = :pressure - for o in out - for i in [5, 6, 8] - if o[i] <= 0 - o[i] = Inf - end - end - w = get_well(outer_data, o[1]) - if w.preferred_phase == "GAS" - utypes[5] = :gas_rate_surface - utypes[10] = :u_rv - else - utypes[5] = :liquid_rate_surface - if w.preferred_phase == "OIL" - utypes[12] = :u_rs - end - end - swap_unit_system_axes!(o, units, utypes) - end - data["WCONINJE"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WCONINJH}) - d = "Default" - defaults = [d, d, "OPEN", 0.0, 0.0, 0.0, 0, 0.0, 0.0, 0.0, 0.0, "RATE"] - wells = get_wells(outer_data) - out = parse_defaulted_group_well(f, defaults, wells, 1) - utypes = identity_unit_vector(defaults) - utypes[5] = :pressure - utypes[6] = :pressure - for o in out - w = get_well(outer_data, o[1]) - if w.preferred_phase == "GAS" - utypes[4] = :gas_rate_surface - utypes[8] = :u_rv - else - utypes[4] = :liquid_rate_surface - if w.preferred_phase == "OIL" - utypes[8] = :u_rs - end - end - swap_unit_system_axes!(o, units, utypes) - end - data["WCONINJH"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Union{Val{:WELLTARG}, Val{:WELTARG}}) - defaults = ["Default", "ORAT", NaN] - wells = get_wells(outer_data) - out = parse_defaulted_group_well(f, defaults, wells, 1) - parser_message(cfg, outer_data, "WELTARG", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["WELTARG"] = out -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, "YES"] - wells = get_wells(outer_data) - parsed = parse_defaulted_group_well(f, defaults, wells, 1) - push_and_create!(data, "WEFAC", parsed) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WPIMULT}) - parser_message(cfg, outer_data, "WPIMULT", PARSER_JUTULDARCY_MISSING_SUPPORT) - - defaults = ["Default", 1.0, -1, -1, -1, -1, -1, -1] - wells = get_wells(outer_data) - parsed = parse_defaulted_group_well(f, defaults, wells, 1) - push_and_create!(data, "WPIMULT", parsed) -end - -function convert_date_kw(t) - @assert length(t) == 4 - function get_month(s) - s = uppercase(s) - if s == "JAN" - return 1 - elseif s == "FEB" - return 2 - elseif s == "MAR" - return 3 - elseif s == "APR" - return 4 - elseif s == "MAY" - return 5 - elseif s == "JUN" - return 6 - elseif s == "JLY" || s == "JUL" - return 7 - elseif s == "AUG" - return 8 - elseif s == "SEP" - return 9 - elseif s == "OCT" - return 10 - elseif s == "NOV" - return 11 - else - @assert s == "DEC" "Did not understand month format $s" - return 12 - end - end - yr = t[3] - m = get_month(t[2]) - d = t[1] - - return parse(DateTime, "$yr-$d-$m $(t[4])", dateformat"y-d-m H:M:S") -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DATES}) - defaults = [1, "JAN", 1970, "00:00:00"] - dt = parse_defaulted_group(f, defaults) - out = Vector{DateTime}() - sizehint!(out, length(dt)) - for t in dt - if length(t[end]) > 8 - # TODO: Fix dirty hack for skipping ms - t[end] = t[end][1:8] - end - push!(out, convert_date_kw(t)) - end - data["DATES"] = out -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TSTEP}) - dt = parse_deck_vector(f) - swap_unit_system!(dt, units, :time) - data["TSTEP"] = dt -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TUNING}) - skip_records(f, 3) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTSCHED}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NUPCOL}) - rec = read_record(f) - tdims = [3]; - data["NUPCOL"] = only(parse_defaulted_line(rec, tdims)) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WELLSTRE}) - n = compositional_number_of_components(outer_data) - defaults = ["DEFAULT", 0.0] - for i in 1:(n-1) - push!(defaults, 0.0) - end - rec = parse_defaulted_group(f, defaults) - for r in rec - sum = 0.0 - for i in 1:n - sum += r[i+1] - end - @assert sum ≈ 1.0 "Compositions in well stream $(r[1]) defined by WELLSTRE must sum up to one (was $sum)" - end - data["WELLSTRE"] = rec -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WINJGAS}) - defaults = ["Default", "GRUP", "Default", "Default", 0] - wells = get_wells(outer_data) - d = parse_defaulted_group_well(f, defaults, wells, 1) - data["WINJGAS"] = d -end - - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WRFTPLT}) - parser_message(cfg, outer_data, "WRFTPLT", PARSER_MISSING_SUPPORT) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:GRUPTREE}) - parser_message(cfg, outer_data, "GRUPTREE", PARSER_MISSING_SUPPORT) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUFETP}) - parser_message(cfg, outer_data, "AQUFETP", PARSER_MISSING_SUPPORT) - skip_record(f) -end - -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] - @assert wname != d - wrec = wheader[2:end] - - # Now follows the segments - 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 - -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 - -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 deleted file mode 100644 index 546cdb20..00000000 --- a/src/InputParser/deckinput/keywords/solution.jl +++ /dev/null @@ -1,213 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:SOLUTION}) - -end - -# Utilities - -function get_cartdims(outer_data) - g = get_section(outer_data, :GRID) - @assert haskey(g, "cartDims") "Cannot access cartDims, has not been set." - return g["cartDims"] -end - -function get_boxdims(outer_data) - gdata = get_section(outer_data, :GRID) - box = gdata["CURRENT_BOX"] - dim = @. box.upper - box.lower + 1 - return dim -end - -function get_box_indices(outer_data) - gdata = get_section(outer_data, :GRID) - box = gdata["CURRENT_BOX"] - l = box.lower - u = box.upper - return (l[1]:u[1], l[2]:u[2], l[3]:u[3]) -end - -function get_box_indices(outer_data, il, iu, jl, ju, kl, ku) - return (il:iu, jl:ju, kl:ku) -end - -function set_cartdims!(outer_data, dim) - @assert length(dim) == 3 - g = get_section(outer_data, :GRID) - dim = tuple(dim...) - gdata = get_section(outer_data, :GRID) - gdata["cartDims"] = dim - reset_current_box!(outer_data) -end - -function reset_current_box!(outer_data) - gdata = get_section(outer_data, :GRID) - gdata["CURRENT_BOX"] = (lower = (1, 1, 1), upper = gdata["cartDims"]) -end - -# Keywords follow - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SGAS}) - data["SGAS"] = parse_grid_vector(f, outer_data["GRID"]["cartDims"], Float64) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SWAT}) - data["SWAT"] = parse_grid_vector(f, outer_data["GRID"]["cartDims"], Float64) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TEMPI}) - T_i = parse_grid_vector(f, outer_data["GRID"]["cartDims"], Float64) - swap_unit_system!(T_i, units, :relative_temperature) - data["TEMPI"] = T_i -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DATUM}) - rec = read_record(f) - l = parse_defaulted_line(rec, [NaN]) - swap_unit_system!(l, units, :length) - data["DATUM"] = only(l) -end - -function parse_mole_fractions!(f, outer_data) - d = outer_data["GRID"]["cartDims"] - nc = compositional_number_of_components(outer_data) - return parse_grid_vector(f, (d[1], d[2], d[3], nc), Float64) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ZMF}) - data["ZMF"] = parse_mole_fractions!(f, outer_data) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:XMF}) - data["XMF"] = parse_mole_fractions!(f, outer_data) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:YMF}) - data["YMF"] = parse_mole_fractions!(f, outer_data) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PRESSURE}) - p = parse_grid_vector(f, outer_data["GRID"]["cartDims"], Float64) - swap_unit_system!(p, units, :pressure) - data["PRESSURE"] = p -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RTEMP}) - 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}) - rs = parse_grid_vector(f, outer_data["GRID"]["cartDims"], Float64) - swap_unit_system!(rs, units, :u_rs) - data["RS"] = rs -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ACTNUM}) - parse_and_set_grid_data!(data, outer_data, units, cfg, f, :ACTNUM, T = Bool) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RSVD}) - n = number_of_tables(outer_data, :equil) - out = [] - for i = 1:n - rs = parse_deck_matrix(f) - swap_unit_system_axes!(rs, units, (:length, :u_rs)) - push!(out, rs) - end - 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 = [] - 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 - tvd = parse_deck_matrix(f) - swap_unit_system_axes!(tvd, units, (:length, :relative_temperature)) - push!(out, tvd) - end - data["RTEMPVD"] = out -end - - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RVVD}) - n = number_of_tables(outer_data, :equil) - out = [] - for i = 1:n - rs = parse_deck_matrix(f) - swap_unit_system_axes!(rs, units, (:length, :u_rv)) - push!(out, rs) - end - data["RVVD"] = out -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, 1, 0] - eunits = (:length, :pressure, :length, :pressure, :length, :pressure, :id, :id, :id, :id, :id) - out = [] - for i = 1:n - rec = read_record(f) - result = parse_defaulted_line(rec, def) - swap_unit_system_axes!(result, units, eunits) - push!(out, result) - 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/keywords/special.jl b/src/InputParser/deckinput/keywords/special.jl deleted file mode 100644 index be8388ff..00000000 --- a/src/InputParser/deckinput/keywords/special.jl +++ /dev/null @@ -1,353 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:EDIT}) - -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BOX}) - rec = read_record(f) - tdims = [1]; - gdata = get_section(outer_data, :GRID) - l, u = gdata["CURRENT_BOX"] - il, jl, kl = l - iu, ju, ku = u - - il, iu, jl, ju, kl, ku = parse_defaulted_line(rec, [il, iu, jl, ju, kl, ku]) - gdata["CURRENT_BOX"] = (lower = (il, jl, kl), upper = (iu, ju, ku)) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ENDBOX}) - reset_current_box!(outer_data) -end - -function parse_keyword!(data, outer_data, units, cfg, f, kval::Union{Val{:COPY}, Val{:ADD}, Val{:MULTIPLY}}) - k = unpack_val(kval) - is_copy = k == :COPY - rec = read_record(f) - gdata = get_section(outer_data, :GRID) - l, u = gdata["CURRENT_BOX"] - dims = get_cartdims(outer_data) - d = "Default" - - il = l[1] - iu = u[1] - jl = l[2] - ju = u[2] - kl = l[3] - ku = u[3] - if is_copy - d_op = d - else - d_op = NaN - end - - while length(rec) > 0 - parsed = parse_defaulted_line(rec, [d, d_op, il, iu, jl, ju, kl, ku]) - if is_copy - dst = parsed[2] - src = parsed[1] - op = NaN - @assert src != d "Source was defaulted for $k. rec = $rec" - else - dst = parsed[1] - op = parsed[2] - src = missing - @assert op != d_op "Operator was defaulted for $k. rec = $rec" - end - @assert dst != d "Destination was defaulted for $k. rec = $rec" - - # Box can be kept. - il = parsed[3] - iu = parsed[4] - jl = parsed[5] - ju = parsed[6] - kl = parsed[7] - ku = parsed[8] - IJK = get_box_indices(outer_data, il, iu, jl, ju, kl, ku) - if is_copy - if !haskey(data, dst) - T = eltype(data[src]) - data[dst] = zeros(T, dims) - end - apply_copy!(data[dst], data[src], IJK, dims) - else - if !haskey(data, dst) - data[dst] = zeros(Float64, dims) - end - if k == :ADD - # add is a const - apply_add!(data[dst], op, IJK, dims) - else - # multiply is a const - apply_multiply!(data[dst], op, IJK, dims) - end - end - rec = read_record(f) - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:OPERATE}) - rec = read_record(f) - gdata = get_section(outer_data, :GRID) - l, u = gdata["CURRENT_BOX"] - dims = get_cartdims(outer_data) - - 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, il, iu, jl, ju, kl, ku, d, d, NaN, NaN]) - target = parsed[1] - op = parsed[8] - source = parsed[9] - @assert target != d "Target was defaulted? rec = $rec" - @assert op != d "Operator was defaulted? rec = $rec" - @assert source != d "Source was defaulted? rec = $rec" - - # Box can be kept. - il = parsed[2] - iu = parsed[3] - jl = parsed[4] - ju = parsed[5] - kl = parsed[6] - ku = parsed[7] - - op_prm1 = parsed[10] - op_prm2 = parsed[11] - - IJK = get_box_indices(outer_data, il, iu, jl, ju, kl, ku) - 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!(operation_target, operation_source, IJK, op, op_prm1, op_prm2) - # On to the next one. - rec = read_record(f) - end -end - -function apply_operate!(target, source, IJK, op, prm1, prm2) - I, J, K = IJK - op = lowercase(op) - if op == "multx" - F = (t, s) -> prm1*s - elseif op == "addx" - F = (t, s) -> s + prm1 - elseif op == "multa" - F = (t, s) -> prm1*s + prm2 - elseif op == "abs" - F = (t, s) -> abs(s) - elseif op == "loge" - F = (t, s) -> ln(s) - elseif op == "log10" - F = (t, s) -> log10(s) - elseif op == "slog" - F = (t, s) -> 10^(prm1 + prm2*s) - elseif op == "poly" - F = (t, s) -> t + prm1*s^prm2 - elseif op == "inv" - F = (t, s) -> 1.0/s - elseif op == "multiply" - F = (t, s) -> t*s - elseif op == "multp" - F = (t, s) -> prm1*s^prm2 - elseif op == "minlim" - F = (t, s) -> max(prm1, s) - elseif op == "maxlim" - F = (t, s) -> min(prm1, s) - elseif op == "copy" - F = (t, s) -> s - else - error("OPERATE option $(uppercase(op)) is not implemented.") - end - for i in I - for j in J - for k in K - target[i, j, k] = F(target[i, j, k], source[i, j, k]) - end - end - end -end - -function apply_copy!(vals, src, IX, dims) - I, J, K = IX - vals[I, J, K] = src -end - -function apply_add!(vals, src, IX, dims) - I, J, K = IX - vals[I, J, K] .+= src -end - -function apply_multiply!(vals, src, IX, dims) - I, J, K = IX - 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) - grid = outer_data["GRID"] - l, u = grid["CURRENT_BOX"] - dims = 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] - factor = parsed[2] - @assert dst != "Default" - - # Box can be kept. - il = parsed[3] - iu = parsed[4] - jl = parsed[5] - ju = parsed[6] - kl = parsed[7] - ku = parsed[8] - 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 - - -function apply_multiply!(target::AbstractVector, factor, I, J, K, dims) - apply_multiply!(reshape(target, dims), factor, I, J, K, dims) -end - - -function apply_multiply!(target, factor, I, J, K, dims) - 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] *= factor - end - end - end -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EQUALS}) - # TODO: Merge shared code with COPY - 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] - constval = parsed[2] - @assert dst != "Default" - 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] - jl = parsed[5] - ju = parsed[6] - kl = parsed[7] - ku = parsed[8] - apply_equals!(target, constval, (il, iu), (jl, ju), (kl, ku), dims) - end - rec = read_record(f) - end -end - -function apply_equals!(target, constval, I, J, K, dims) - 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] = constval - end - 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 diff --git a/src/InputParser/deckinput/keywords/summary.jl b/src/InputParser/deckinput/keywords/summary.jl deleted file mode 100644 index 6e0b12ee..00000000 --- a/src/InputParser/deckinput/keywords/summary.jl +++ /dev/null @@ -1,105 +0,0 @@ -function finish_current_section!(data, units, cfg, outer_data, ::Val{:SUMMARY}) - -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NSTACK}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTSOL}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:EXCEL}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RUNSUM}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:SEPARATE}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NEWTON}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:TCPU}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:ELAPSED}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MAXDPR}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FOPR}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FWPR}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FWIR}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FOPT}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FWPT}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FWIT}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTONLY}) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WMCTL}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTRST}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WLPR}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WGOR}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WWIR}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WWPR}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WOPR}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WBHP}) - read_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BGSAT}) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BOSAT}) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BOKR}) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BWSAT}) - skip_record(f) -end - -function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:BPR}) - skip_record(f) -end diff --git a/src/InputParser/deckinput/parser.jl b/src/InputParser/deckinput/parser.jl deleted file mode 100644 index 7bb5ddf8..00000000 --- a/src/InputParser/deckinput/parser.jl +++ /dev/null @@ -1,259 +0,0 @@ -include("units.jl") -include("utils.jl") -include("keywords/keywords.jl") - -Base.@kwdef struct ParserVerbosityConfig - silent::Bool = false - verbose::Bool = true - warn_feature::Bool = true - warn_parsing::Bool = true - warn_limit::Int = 1 - warn_count::Dict{String, Int} = Dict{String, Int}() -end - -@enum PARSER_WARNING PARSER_MISSING_SUPPORT PARSER_JUTULDARCY_MISSING_SUPPORT PARSER_JUTULDARCY_PARTIAL_SUPPORT PARSER_PARTIAL_SUPPORT - -function keyword_header(current_data, keyword) - if haskey(current_data, "CURRENT_SECTION") - section = current_data["CURRENT_SECTION"] - out = "$section/$keyword" - else - out = "$keyword" - end - return out -end - -function parser_message(cfg::ParserVerbosityConfig, outer_data, keyword, msg::String; important = false) - if !cfg.silent - if important || cfg.verbose - println("$(keyword_header(outer_data, keyword)): $msg") - end - end -end - -function parser_message(cfg::ParserVerbosityConfig, outer_data, keyword, msg::PARSER_WARNING, reason = missing) - if cfg.silent - return - end - if msg == PARSER_MISSING_SUPPORT - text_msg = "Unsupported keyword. It will be ignored." - do_print = cfg.warn_parsing - elseif msg == PARSER_JUTULDARCY_MISSING_SUPPORT - text_msg = "$keyword is not supported by JutulDarcy solvers. It will be ignored in simulations." - do_print = cfg.warn_feature - elseif msg == PARSER_JUTULDARCY_PARTIAL_SUPPORT - text_msg = "$keyword is only partially supported by JutulDarcy solvers." - do_print = cfg.warn_feature - else - @assert msg == PARSER_PARTIAL_SUPPORT - text_msg = "Parser has only partial support. $keyword may have missing or wrong entries." - do_print = cfg.warn_parsing - end - if !ismissing(reason) - text_msg = "$text_msg\n\n$reason" - end - if haskey(cfg.warn_count, keyword) - cfg.warn_count[keyword] += 1 - else - cfg.warn_count[keyword] = 1 - end - if cfg.warn_count[keyword] <= cfg.warn_limit - jutul_message("Parser", "$(keyword_header(outer_data, keyword)) - $text_msg", color = :yellow) - end -end - -""" - parse_data_file(filename; unit = :si) - data = parse_data_file("MY_MODEL.DATA") - -Parse a .DATA file given by the `filename` (industry standard input file) into a -Dict. Units will be converted to strict SI unless you pass something else like -`units = :field`. Setting `units = nothing` will skip unit conversion. Note that -JutulDarcy assumes that the unit system is internally consistent. It is highly -recommended to parse to the SI units if you want to perform simulations. - -The best publicly available documentation on this format is available from the -Open Porous Media (OPM) project's webpages: [OPM Flow manual -](https://opm-project.org/?page_id=955). - -# Keyword arguments -- `warn_parsing=true`: Produce a warning when keywords are not supported (or - partially supported) by the parser. -- `warn_feature`=true`: Produce a warning when keywords are supported, but have - limited or missing support in the numerical solvers. -- `units=:si`: Symbol that indicates the unit system to be used in the output. - Setting this to `nothing` will return values without conversion, i.e. exactly - what is in the input files. `:si` will use strict SI. Other alternatives are - `:field` and `:metric`. `:lab` is currently unsupported. - -# Note -This function is experimental and only covers a small portion of the keywords -that exist for various simulators. You will get warnings that indicate the level -of support for keywords in both the parser and the numerical solvers when known -keywords with limited support. Pull requests for new keywords are welcome! -""" -function parse_data_file(filename; kwarg...) - outer_data = Dict{String, Any}() - parse_data_file!(outer_data, filename; kwarg...) - delete!(outer_data, "CURRENT_SECTION") - return outer_data -end - -function parse_data_file!(outer_data, filename, data = outer_data; - skip_mode = false, - verbose = false, - sections = [:RUNSPEC, :GRID, :PROPS, :REGIONS, :SOLUTION, :SCHEDULE, :EDIT], - skip = [:SUMMARY], - units::Union{Symbol, Nothing} = :si, - warn_parsing = true, - warn_feature = true, - basedir = missing, - silent = false, - is_outer = true, - input_units::Union{Symbol, Nothing} = nothing, - unit_systems = missing - ) - function get_unit_system_pair(from::Symbol, target::Symbol) - from_sys = DeckUnitSystem(from) - target_sys = DeckUnitSystem(target) - # We have a pair of unit systems to convert between - return (from = from_sys, to = target_sys) - end - - if isnothing(units) - # nothing for units means no conversion - @assert ismissing(unit_systems) - unit_systems = nothing - end - if !isnothing(units) && !isnothing(input_units) - # Input and output units are both specified, potentially overriding - # what's in RUNSPEC. This will also check that they are valid symbols - # for us. - unit_systems = get_unit_system_pair(input_units, units) - end - filename = realpath(filename) - if ismissing(basedir) - basedir = dirname(filename) - end - f = open(filename, "r") - - cfg = ParserVerbosityConfig( - verbose = verbose, - warn_parsing = warn_parsing, - warn_feature = warn_feature, - silent = silent - ) - parser_message(cfg, outer_data, "PARSER", "Starting parse of $filename...") - try - allsections = vcat(sections, skip) - while !eof(f) - m = next_keyword!(f) - if isnothing(m) - continue - end - parsed_keyword = false - parser_message(cfg, outer_data, "$m", "Starting parse...") - t_p = @elapsed if m == :SKIP || m == :ENDSKIP - if m == :SKIP - skip_mode = true - else - @assert skip_mode - skip_mode = false - end - elseif m in allsections - # Check if we have passed RUNSPEC so that units can be set - runspec_passed = m != :RUNSPEC && haskey(outer_data, "RUNSPEC") - unit_system_not_initialized = ismissing(unit_systems) - if runspec_passed && unit_system_not_initialized - unit_systems = get_unit_system_pair(current_unit_system(outer_data), units) - end - finish_current_section!(data, unit_systems, cfg, outer_data) - parser_message(cfg, outer_data, "$m", "Starting new section.") - data = new_section(outer_data, m) - skip_mode = m in skip - elseif m == :INCLUDE - next = strip(readline(f)) - if endswith(next, '/') - next = rstrip(next, '/') - else - readline(f) - end - include_path = clean_include_path(basedir, next) - 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, - warn_parsing = warn_parsing, - warn_feature = warn_feature, - silent = silent, - sections = sections, - skip = skip, - basedir = basedir, - skip_mode = skip_mode, - is_outer = false, - units = units, - unit_systems = unit_systems - ) - elseif m in (:DATES, :TIME, :TSTEP) - parsed_keyword = true - parse_keyword!(data, outer_data, unit_systems, cfg, f, Val(m)) - # New control step starts after this - data = OrderedDict{String, Any}() - push!(outer_data["SCHEDULE"]["STEPS"], data) - elseif m == :END - # All done! - finish_current_section!(data, unit_systems, cfg, outer_data) - break - elseif skip_mode - parser_message(cfg, outer_data, "$m", "Keyword skipped.") - else - parse_keyword!(data, outer_data, unit_systems, cfg, f, Val(m)) - parsed_keyword = true - end - if parsed_keyword - parser_message(cfg, outer_data, "$m", "Keyword parsed in $(t_p)s.") - end - end - finally - if is_outer - finish_current_section!(data, unit_systems, cfg, outer_data) - end - close(f) - end - return outer_data -end - -""" - parse_grdecl_file("mygrid.grdecl"; actnum_path = missing, kwarg...) - -Parse a GRDECL file separately from the full input file. Note that the GRID -section does not contain units - passing the `input_units` keyword is therefore -highly recommended. - -# Keyword arguments - - `units=:si`: Units to use for return values. Requires `input_units` to be set. - - `input_units=nothing`: The units the file is given in. - - `verbose=false`: Toggle verbosity. - -""" -function parse_grdecl_file(filename; actnum_path = missing, kwarg...) - outer_data = Dict{String, Any}() - data = new_section(outer_data, :GRID) - parse_data_file!(outer_data, filename, data; kwarg...) - if !ismissing(actnum_path) - parse_data_file!(outer_data, actnum_path, data; kwarg...) - end - if !haskey(data, "ACTNUM") - data["ACTNUM"] = fill(true, data["cartDims"]) - end - delete!(data, "CURRENT_BOX") - return data -end - -function finish_current_section!(data, units, cfg, outer_data) - if haskey(outer_data, "CURRENT_SECTION") - v = outer_data["CURRENT_SECTION"] - v::Symbol - finish_current_section!(data, units, cfg, outer_data, Val(v)) - end -end diff --git a/src/InputParser/deckinput/units.jl b/src/InputParser/deckinput/units.jl deleted file mode 100644 index 564141d8..00000000 --- a/src/InputParser/deckinput/units.jl +++ /dev/null @@ -1,298 +0,0 @@ -function current_unit_system(deck) - rs = deck["RUNSPEC"] - choices = ["METRIC", "SI", "LAB", "FIELD"] - found = false - sys = missing - for choice in choices - if haskey(rs, choice) - @assert ismissing(sys) "Cannot have multiple unit keywords (encountered $choice and $sys)" - sys = choice - end - end - if ismissing(sys) - @warn "Units not found. Assuming field units." - out = :field - else - out = Symbol(lowercase(sys)) - end - return out -end - -Base.@kwdef struct DeckUnitSystem{S, T} - length::T = 1.0 - area::T = 1.0 - time::T = 1.0 - density::T = 1.0 - pressure::T = 1.0 - mol::T = 1.0 - mass::T = 1.0 - u_rs::T = 1.0 - u_rv::T = 1.0 - concentration::T = 1.0 - compressibility::T = 1.0 - viscosity::T = 1.0 - surface_tension::T = 1.0 - jsurface_tension::T = 1.0 - permeability::T = 1.0 - liquid_volume_surface::T = 1.0 - liquid_volume_reservoir::T = 1.0 - liquid_formation_volume_factor::T = 1.0 - gas_volume_surface::T = 1.0 - gas_volume_reservoir::T = 1.0 - gas_formation_volume_factor::T = 1.0 - volume::T = 1.0 - transmissibility::T = 1.0 - rock_conductivity::T = 1.0 - volume_heat_capacity::T = 1.0 - mass_heat_capacity::T = 1.0 - molar_mass::T = 1.0 - relative_temperature::Symbol = :Celsius - absolute_temperature::Symbol = :Kelvin -end - -function DeckUnitSystem(sys::Symbol, T = Float64) - #meter, day, kilogram, bar = si_units(:meter, :day, :kilogram, :bar) - u = Jutul.all_units() - m = u[:meter] - K = u[:kelvin] - day = u[:day] - centi = u[:centi] - kilogram = u[:kilogram] - ft = u[:feet] - psi = u[:psi] - pound = u[:pound] - kilo = u[:kilo] - stb = u[:stb] - rankine = u[:rankine] - btu = u[:btu] - - # Commons - cP = u[:centi]*u[:poise] - mD = u[:milli]*u[:darcy] - if sys == :metric - len = m - kJ = u[:kilo]*u[:joule] - volume = m^3 - time = day - pressure = u[:bar] - mol = u[:kilo] - molar_mass = 1/(mol) - mass = kilogram - viscosity = cP - surface_tension = u[:newton]/m - jsurface_tension = u[:dyne]/(centi*m) - permeability = mD - liquid_volume_surface = volume - liquid_volume_reservoir = volume - gas_volume_surface = volume - gas_volume_reservoir = volume - rock_conductivity = kJ/(m*day*K) - volume_heat_capacity = kJ/(volume*K) - mass_heat_capacity = kJ/(mass*K) - relative_temperature = :Celsius - absolute_temperature = :Kelvin - elseif sys == :field - len = ft - time = day - pressure = psi - mol = pound*kilo - molar_mass = 1/(mol) - mass = pound - viscosity = cP - surface_tension = u[:lbf]/u[:inch] - jsurface_tension = u[:dyne]/(centi*m) - permeability = mD - liquid_volume_surface = stb - liquid_volume_reservoir = stb - gas_volume_surface = kilo*ft^3 - gas_volume_reservoir = stb - rock_conductivity = btu / (ft*day*rankine) - volume_heat_capacity = btu / (ft^3*rankine) - mass_heat_capacity = btu / (pound*rankine) - relative_temperature = :Fahrenheit - absolute_temperature = :Rankine - elseif sys == :lab - error("Not implemented") - else - @assert sys == :si - len = 1.0 - time = 1.0 - pressure = 1.0 - mol = 1.0 - mass = 1.0 - viscosity = 1.0 - surface_tension = 1.0 - jsurface_tension = 1.0 - permeability = 1.0 - liquid_volume_surface = 1.0 - liquid_volume_reservoir = 1.0 - gas_volume_surface = 1.0 - gas_volume_reservoir = 1.0 - rock_conductivity = 1.0 - volume_heat_capacity = 1.0 - mass_heat_capacity = 1.0 - molar_mass = 1.0 - relative_temperature = :Celsius - absolute_temperature = :Kelvin - end - area = len^2 - volume = len^3 - density = mass/volume - concentration = mass/volume - compressibility = 1.0/pressure - transmissibility = viscosity * liquid_volume_reservoir / (time * pressure) - return DeckUnitSystem{sys, T}( - length = len, - area = area, - time = time, - u_rs = gas_volume_surface/liquid_volume_surface, - u_rv = liquid_volume_surface/gas_volume_surface, - density = density, - pressure = pressure, - mol = mol, - mass = mass, - concentration = concentration, - compressibility = compressibility, - viscosity = viscosity, - surface_tension = surface_tension, - jsurface_tension = jsurface_tension, - permeability = permeability, - liquid_volume_surface = liquid_volume_surface, - liquid_volume_reservoir = liquid_volume_reservoir, - liquid_formation_volume_factor = liquid_volume_reservoir/liquid_volume_surface, - gas_volume_surface = gas_volume_surface, - gas_volume_reservoir = gas_volume_reservoir, - gas_formation_volume_factor = gas_volume_reservoir/gas_volume_surface, - volume = volume, - molar_mass = molar_mass, - transmissibility = transmissibility, - rock_conductivity = rock_conductivity, - volume_heat_capacity = volume_heat_capacity, - mass_heat_capacity = mass_heat_capacity, - relative_temperature = relative_temperature, - absolute_temperature = absolute_temperature - ) -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) - if dim == 1 - x_t = x' - else - x_t = x - end - for i in axes(x_t, 2) - x_i = view(x_t, :, i) - swap_unit_system!(x_i, systems, eachunit[i]) - end - return x -end - -function swap_unit_system_axes!(x::AbstractVector, systems, eachunit) - @assert eltype(eachunit)<:Symbol - @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 - return x -end - -function swap_unit_system!(x::AbstractArray, systems, k) - return swap_unit_system!(x, systems, Val(k)) -end - -function swap_unit_system!(x::AbstractArray, systems, k::Val) - for i in eachindex(x) - x[i] = swap_unit_system(x[i], systems, k) - end - return x -end - -function swap_unit_system(val, systems, k::Symbol) - return swap_unit_system(val, systems, Val(k)) -end - -function swap_unit_system(v, systems::Union{Nothing, Missing}, k::Val) - # No systems - trivial conversion - return v -end - -function swap_unit_system(v, systems::Union{Nothing, Missing}, k::Symbol) - # No systems - trivial conversion - return v -end - -function swap_unit_system(val, systems::NamedTuple, ::Union{Val{:identity}, Val{:id}}) - # Identity specifically means no unit. - return val -end - -function identity_unit_vector(x) - return identity_unit_vector(length(x)) -end - -function identity_unit_vector(n::Int) - utypes = Vector{Symbol}(undef, n) - fill!(utypes, :id) - return utypes -end - -function swap_unit_system(val, systems::NamedTuple, ::Val{k}; reverse = false) where k - (; to, from) = systems - if reverse - to, from = from, to - end - to_unit = deck_unit(to, k) - from_unit = deck_unit(from, k) - - val_si = convert_to_si(val, from_unit) - val_final = convert_from_si(val_si, to_unit) - return val_final -end - -function deck_unit(sys::DeckUnitSystem, s::Symbol) - return deck_unit(sys, Val(s)) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{k}) where k - return getproperty(sys, k) -end - -# Magic type overloads - -function deck_unit(sys::DeckUnitSystem, ::Val{:Kh}) - return deck_unit(sys, :permeability)*deck_unit(sys, :length) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:gigapascal}) - return si_unit(:Pa)*si_unit(:giga) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:time_over_volume}) - return deck_unit(sys, :time)/deck_unit(sys, :volume) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:liquid_rate_surface}) - return deck_unit(sys, :liquid_volume_surface)/deck_unit(sys, :time) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:gas_rate_surface}) - return deck_unit(sys, :gas_volume_surface)/deck_unit(sys, :time) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:liquid_rate_reservoir}) - return deck_unit(sys, :liquid_volume_reservoir)/deck_unit(sys, :time) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:gas_rate_reservoir}) - return deck_unit(sys, :gas_volume_reservoir)/deck_unit(sys, :time) -end - -function deck_unit(sys::DeckUnitSystem, ::Val{:critical_volume}) - return deck_unit(sys, :volume)/deck_unit(sys, :mol) -end diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl deleted file mode 100644 index 9ee7ac9a..00000000 --- a/src/InputParser/deckinput/utils.jl +++ /dev/null @@ -1,634 +0,0 @@ -function unpack_val(::Val{X}) where X - return X -end - -const DECK_SPLIT_REGEX = r"[ \t,]+" - -function read_record(f; fix = true) - split_lines = Vector{String}() - active = true - while !eof(f) && active - line = readline(f) - cpos = findfirst("--", line) - if !isnothing(cpos) - line = line[1:(first(cpos)-1)] - end - line = strip(line) - if !startswith(line, "--") - line = String(line) - if contains(line, '/') - # TODO: Think this is OK for parsing ASCII. - ix = findfirst('/', line) - line = line[1:ix-1] - active = false - end - line::String - if length(line) > 0 - push!(split_lines, line) - end - end - end - return split_lines -end - -function keyword_start(line) - if isnothing(match(r"^\s*--", line)) - m = match(r"\w+", line) - if m === nothing - return nothing - else - return Symbol(uppercase(m.match)) - end - else - return nothing - end -end - -function parse_defaulted_group_well(f, defaults, wells, namepos = 1) - out = [] - line = read_record(f) - while length(line) > 0 - allow_wildcard = fill(true, length(defaults)) - allow_wildcard[1] = false - parsed = parse_defaulted_line(line, defaults, allow_wildcard = allow_wildcard) - name = parsed[namepos] - if occursin('*', name) || occursin('?', name) - re = Regex(replace(name, "*" => ".*", "?" => ".")) - for wname in keys(wells) - if occursin(re, wname) - replaced_parsed = copy(parsed) - replaced_parsed[namepos] = wname - push!(out, replaced_parsed) - end - end - else - push!(out, parsed) - end - line = read_record(f) - end - return out -end - -function parse_defaulted_group(f, defaults) - out = [] - line = read_record(f) - while length(line) > 0 - parsed = parse_defaulted_line(line, defaults) - push!(out, parsed) - line = read_record(f) - end - return out -end - -function parse_defaulted_line(lines::String, defaults; kwarg...) - return parse_defaulted_line([lines], defaults; kwarg...) -end - -function parse_defaulted_line(lines, defaults; required_num = 0, keyword = "", allow_wildcard = missing) - out = similar(defaults, 0) - sizehint!(out, length(defaults)) - pos = 1 - for line in lines - line = replace_quotes(line) - lsplit = split(strip(line), DECK_SPLIT_REGEX) - for s in lsplit - if length(s) == 0 - continue - end - default = defaults[pos] - if ismissing(allow_wildcard) - allow_star = true - else - allow_star = allow_wildcard[pos] - end - if allow_star && occursin('*', s) && !startswith(s, '\'') # Could be inside a string for wildcard matching - if s == "*" - num_defaulted = 1 - else - parse_wildcard = match(r"\d+\*", s) - if isnothing(parse_wildcard) - error("Unable to parse string for * expansion: $s") - end - num_defaulted = Parsers.parse(Int, parse_wildcard.match[1:end-1]) - end - for i in 1:num_defaulted - push!(out, defaults[pos]) - pos += 1 - end - else - if default isa String - converted = strip(s, [' ', '\'']) - else - T = typeof(default) - converted = Parsers.tryparse(T, s) - if isnothing(converted) - converted = T.(Parsers.tryparse(Float64, s)) - end - end - push!(out, converted) - pos += 1 - end - end - end - n = length(defaults) - n_out = length(out) - if required_num > n - error("Bad record: $required_num entries required for keyword $keyword, but only $n records were present.") - end - pos = n_out + 1 - if pos < n + 1 - for i in pos:n - push!(out, defaults[i]) - end - end - return out -end - -## - -function parse_deck_matrix(f, T = Float64) - # TODO: This is probably a bad way to do large numerical datasets. - rec = read_record(f) - split_lines = preprocess_delim_records(rec) - data = Vector{T}() - n = -1 - for seg in split_lines - n = parse_deck_matrix_line!(data, seg, n) - end - if length(data) == 0 - out = missing - else - out = reshape(data, n, length(data) ÷ n)' - end - return out -end - -function parse_deck_matrix_line!(data::Vector{T}, seg, n) where T - m = length(seg) - if m == 0 - return n - elseif n == -1 - n = m - else - @assert m == n "Expected $n was $m" - end - for d in seg - if d == raw"1*" - # Defaulted... - @assert T == Float64 - push!(data, NaN) - else - push!(data, Parsers.parse(T, d)) - end - end - return n -end - -function preprocess_delim_records(split_lines) - # Old slow code: - # split_lines = map(strip, split_lines) - # filter!(x -> !startswith(x, "--"), split_lines) - # split_rec = map(x -> split(x, r"\s*,?\s+"), split_lines) - split_rec = Vector{Vector{String}}() - sizehint!(split_rec, length(split_lines)) - for line in split_lines - # Strip end whitespace - line = strip(line) - # Remove comments - if !startswith(line, "--") - # Split into entries (could be whitespace + a comma anywhere in between) - sub_rec = Vector{String}() - sizehint!(sub_rec, 10) - for rec in eachsplit(line, r"\s*,?\s+") - push!(sub_rec, String(rec)) - end - if length(sub_rec) > 0 - push!(split_rec, sub_rec) - end - end - end - return split_rec -end - -function parse_deck_vector(f, T = Float64) - # TODO: Speed up. - rec = read_record(f) - record_lines = preprocess_delim_records(rec) - n = length(record_lines) - out = Vector{T}() - sizehint!(out, n) - for split_rec in record_lines - for el in split_rec - if occursin('*', el) - n_rep, el = split(el, '*') - n_rep = parse(Int, n_rep) - else - n_rep = 1 - end - parsed = parse(T, el)::T - for i in 1:n_rep - push!(out, parsed) - end - end - end - return out -end - -function skip_record(f) - rec = read_record(f) - while length(rec) > 0 - rec = read_record(f) - end -end - -function skip_records(f, n) - for i = 1:n - rec = read_record(f) - end -end - -function parse_grid_vector(f, dims, T = Float64) - v = parse_deck_vector(f, T) - return reshape(v, dims) -end - -function parse_saturation_table(f, outer_data) - ns = number_of_tables(outer_data, :saturation) - return parse_region_matrix_table(f, ns) -end - -function parse_dead_pvt_table(f, outer_data) - np = number_of_tables(outer_data, :pvt) - return parse_region_matrix_table(f, np) -end - -function parse_live_pvt_table(f, outer_data) - nreg = number_of_tables(outer_data, :pvt) - out = [] - for i = 1:nreg - current = Vector{Vector{Float64}}() - while true - next = parse_deck_vector(f) - if length(next) == 0 - break - end - push!(current, next) - end - push!(out, restructure_pvt_table(current)) - end - return out -end - -function restructure_pvt_table(tab) - nvals_per_rec = 3 - function record_length(x) - # Actual number of records: 1 key value + nrec*N entries. Return N. - return (length(x) - 1) ÷ nvals_per_rec - end - @assert record_length(last(tab)) > 1 - nrecords = length(tab) - keys = map(first, tab) - current = 1 - for tab_ix in eachindex(tab) - rec = tab[tab_ix] - interpolate_missing_usat!(tab, tab_ix, record_length, nvals_per_rec) - end - # Generate final table - ntab = sum(record_length, tab) - data = zeros(ntab, nvals_per_rec) - for tab_ix in eachindex(tab) - rec = tab[tab_ix] - for i in 1:record_length(rec) - for j in 1:nvals_per_rec - linear_ix = (i-1)*nvals_per_rec + j + 1 - data[current, j] = rec[linear_ix] - end - current += 1 - end - end - - # Generate pos - pos = Int[1] - sizehint!(pos, nrecords+1) - for rec in tab - push!(pos, pos[end] + record_length(rec)) - end - return Dict("data" => data, "key" => keys, "pos" => pos) -end - -function interpolate_missing_usat!(tab, tab_ix, record_length, nvals_per_rec) - rec = tab[tab_ix] - if record_length(rec) == 1 - @assert nvals_per_rec == 3 - next_rec = missing - for j in (tab_ix):length(tab) - if record_length(tab[j]) > 1 - next_rec = tab[j] - break - end - end - @assert record_length(rec) == 1 - next_rec_length = record_length(next_rec) - sizehint!(rec, 1 + nvals_per_rec*next_rec_length) - - get_index(major, minor) = nvals_per_rec*(major-1) + minor + 1 - pressure(x, idx) = x[get_index(idx, 1)] - B(x, idx) = x[get_index(idx, 2)] - viscosity(x, idx) = x[get_index(idx, 3)] - - function constant_comp_interp(F, F_r, F_l) - # So that dF/dp * F = constant over the pair of points extrapolated from F - w = 2.0*(F_l - F_r)/(F_l + F_r) - return F*(1.0 + w/2.0)/(1.0 - w/2.0) - end - @assert !ismissing(next_rec) "Final table must be saturated." - - for idx in 2:next_rec_length - # Each of these gets added as new unsaturated points - p_0 = pressure(rec, idx - 1) - p_l = pressure(next_rec, idx - 1) - p_r = pressure(next_rec, idx) - - mu_0 = viscosity(rec, idx - 1) - mu_l = viscosity(next_rec, idx - 1) - mu_r = viscosity(next_rec, idx) - - B_0 = B(rec, idx - 1) - B_l = B(next_rec, idx - 1) - B_r = B(next_rec, idx) - - p_next = p_0 + p_r - p_l - B_next = constant_comp_interp(B_0, B_l, B_r) - mu_next = constant_comp_interp(mu_0, mu_l, mu_r) - - push!(rec, p_next) - push!(rec, B_next) - push!(rec, mu_next) - end - end - return tab -end - -function parse_region_matrix_table(f, nreg) - out = [] - for i = 1:nreg - next = parse_deck_matrix(f) - if ismissing(next) - if length(out) == 0 - error("First region table cannot be defaulted.") - end - next = copy(out[end]) - end - push!(out, next) - end - return out -end - -function parse_keyword!(data, outer_data, units, cfg, f, v::Val{T}) where T - # Keywords where we read a single record and don't do anything proper - - skip_kw_with_warn = Symbol[ - :SATOPTS, - :EQLOPTS, - :TRACERS, - :PIMTDIMS, - :FLUXNUM, - :OPTIONS - ] - - skip_list = [] - function skip_kw!(kw, num, msg = nothing) - push!(skip_list, (kw, num, msg)) - end - skip_kw!(:PETOPTS, 1) - skip_kw!(:PARALLEL, 1) - 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) - skip_kw!(:NONNC, 0) - skip_kw!(:NEWTRAN, 0) - 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) - skip_kw!(:TRACERS, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:PIMTDIMS, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:FLUXNUM, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:OPTIONS, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:EHYSTR, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:SWATINIT, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:ZIPPY2, 1, PARSER_MISSING_SUPPORT) - skip_kw!(:DRSDT, 1, PARSER_MISSING_SUPPORT) - 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) - skip_kw!(:PIMULTAB, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:VFPPROD, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:VFPINJ, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:WTRACER, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:GCONINJE, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:WTEST, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:WTEMP, Inf, PARSER_MISSING_SUPPORT) - - found = false - for (kw, num, msg) in skip_list - if kw != T - continue - end - if !isnothing(msg) - parser_message(cfg, outer_data, "$kw", msg) - end - if num == 0 - # Single word keywords are trivial to parse, just set a true flag. - data["$T"] = true - elseif num == 1 - data["$T"] = read_record(f) - else - skip_record(f) - end - found = true - break - end - - if !found - if startswith("$T", "TVDP") - parser_message(cfg, outer_data, "$T", PARSER_MISSING_SUPPORT) - read_record(f) - else - error("Unhandled keyword $T encountered.") - end - end -end - -function next_keyword!(f) - m = nothing - while isnothing(m) && !eof(f) - line = strip(readline(f)) - if line == "/" - continue - end - m = keyword_start(line) - end - return m -end - -function number_of_tables(outer_data, t::Symbol) - rs = outer_data["RUNSPEC"] - if haskey(rs, "TABDIMS") - td = rs["TABDIMS"] - else - td = [1 1] - end - if t == :saturation - return td[1] - elseif t == :pvt - return td[2] - elseif t == :equil - if haskey(rs, "EQLDIMS") - return rs["EQLDIMS"][1] - else - return 1 - end - 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) - return outer_data["RUNSPEC"]["COMPS"] -end - -function table_region(outer_data, t::Symbol; active = nothing) - num = number_of_tables(outer_data, t) - if num == 1 - dim = outer_data["GRID"]["cartDims"] - D = ones(Int, prod(dim)) - else - reg = outer_data["REGIONS"] - - function get_or_default(k) - if haskey(reg, k) - return vec(reg[k]) - else - dim = outer_data["GRID"]["cartDims"] - return ones(Int, prod(dim)) - end - end - - if t == :saturation - d = get_or_default("SATNUM") - elseif t == :pvt - d = get_or_default("PVTNUM") - elseif t == :equil - d = get_or_default("EQLNUM") - else - error(":$t is not known") - end - D = vec(d) - end - if !isnothing(active) - D = D[active] - end - return D -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] - end - include_file_name = replace(include_file_name, "'" => "") - # Do this one more time in case we have nested string and ./ - 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 - -function get_section(outer_data, name::Symbol) - s = "$name" - is_sched = name == :SCHEDULE - outer_data["CURRENT_SECTION"] = name - T = OrderedDict{String, Any} - if is_sched - if !haskey(outer_data, s) - outer_data[s] = Dict( - "STEPS" => [T()], - "WELSPECS" => T(), - "COMPORD" => T() - ) - end - out = outer_data[s]["STEPS"][end] - else - if !haskey(outer_data, s) - outer_data[s] = T() - end - out = outer_data[s] - end - return out -end - -function new_section(outer_data, name::Symbol) - data = get_section(outer_data, name) - return data -end - -function replace_quotes(str::String) - if '\'' in str - v = collect(str) - in_quote = false - new_char = Char[] - for i in eachindex(v) - v_i = v[i] - if v_i == '\'' - in_quote = !in_quote - elseif in_quote && v_i == ' ' - # TODO: Is this a safe replacement? - push!(new_char, '-') - else - push!(new_char, v_i) - end - end - str = String(new_char) - end - return str -end - -function push_and_create!(data, k, vals, T = Any) - if !haskey(data, k) - data[k] = T[] - end - out = data[k] - for v in vals - v::T - push!(out, v) - end - return data -end diff --git a/src/JutulDarcy.jl b/src/JutulDarcy.jl index 04db83fd..d18531ac 100644 --- a/src/JutulDarcy.jl +++ b/src/JutulDarcy.jl @@ -11,7 +11,6 @@ $(EXPORTS) module JutulDarcy export MultiPhaseSystem, ImmiscibleSystem, SinglePhaseSystem export reservoir_linsolve - export parse_data_file, parse_grdecl_file export get_1d_reservoir export DeckPhaseViscosities export DeckShrinkageFactors @@ -142,9 +141,11 @@ module JutulDarcy using TimerOutputs using PrecompileTools using Dates + using GeoEnergyIO import DataStructures: OrderedDict using DocStringExtensions + timeit_debug_enabled() = Jutul.timeit_debug_enabled() include("types.jl") include("deck_types.jl") @@ -172,8 +173,6 @@ module JutulDarcy include("input_simulation/input_simulation.jl") # Initialization by equilibriation include("init/init.jl") - # Corner point grids - include("cpgrid/cpgrid.jl") # Gradients, objective functions, etc include("gradients/gradients.jl") @@ -190,11 +189,6 @@ module JutulDarcy include("ext.jl") - include("InputParser/InputParser.jl") - using .InputParser - import .InputParser: parse_data_file - - @compile_workload begin precompile_darcy_multimodels() # We run a tiny MRST case to precompile the .MAT file loading diff --git a/src/cpgrid/cpgrid.jl b/src/cpgrid/cpgrid.jl deleted file mode 100644 index 6f700935..00000000 --- a/src/cpgrid/cpgrid.jl +++ /dev/null @@ -1,694 +0,0 @@ -include("processing.jl") -function ijk_to_linear(i, j, k, dims) - nx, ny, nz = dims - return (k-1)*nx*ny + (j-1)*nx + i -end - -function ij_to_linear(i, j, dims) - nx, ny = dims - return (j-1)*nx + i -end - -function linear_to_ijk(ix, dims) - nx, ny, nz = dims - linear_index = ix - x = mod(linear_index - 1, nx) + 1 - y = mod((linear_index - x) ÷ nx, ny) + 1 - leftover = (linear_index - x - (y-1)*nx) - z = (leftover ÷ (nx*ny)) + 1 - return (x, y, z) -end - -function get_line(coord, i, j, nx, ny) - ix = ijk_to_linear(i, j, 1, (nx, ny, 1)) - T = SVector{3, Float64} - x1 = T(coord[ix, 1:3]) - x2 = T(coord[ix, 4:end]) - - return (x1, x2) -end - -function interp_coord(line::NamedTuple, z) - if line.equal_points - x = line.x1 - T = eltype(x) - pt = SVector{3, T}(x[1], x[2], z) - else - pt = interp_coord(line.x1, line.x2, z) - end - return pt -end - -function interp_coord(p0::SVector{3, T}, p1::SVector{3, T}, z::T) where T<:Real - z0 = p0[3] - z1 = p1[3] - if z0 ≈ z1 - # Coinciding corner points. Just return the point and hope the pillar is - # inactive. - interp_pt = p0 - else - weight = (z - z0)/(z1 - z0) - interp_pt = p0 .+ weight.*(p1 .- p0) - @assert interp_pt[3] ≈ z "expected $z was $(interp_pt[3]) != $z" - end - return interp_pt -end - -function corner_index(cell, corner::NTuple, dims) - # Cell 1 [corner1, corner2], cell 2 [corner1, corner2], ..., cell n [corner1, corner2] repeated for nx in top plane - # Cell 1 [corner3, corner4], cell 2 [corner3, corner4], ..., cell n [corner3, corner4] - # Cell 1 [corner5, corner6], cell 2 [corner5, corner6], ..., cell n [corner5, corner6] - # Cell 1 [corner7, corner8], cell 2 [corner7, corner8], ..., cell n [corner7, corner8] - - if cell isa Int - i, j, k = linear_to_ijk(cell, dims) - else - i, j, k = cell - cell = ijk_to_linear(i, j, k, cartdims) - end - nx, ny, nz = dims - i_is_upper, j_is_upper, k_is_upper = corner - - j_is_upper, i_is_upper, k_is_upper = corner - cell_i_offset = 2*(i-1) - cell_j_offset = 4*nx*(j-1) - cell_k_offset = 8*nx*ny*(k-1) - - cell_offset = cell_i_offset + cell_j_offset + cell_k_offset - - i_offset = i_is_upper+1 - j_offset = j_is_upper*2*nx - k_offset = k_is_upper*4*nx*ny - - ijk = i_offset + j_offset + k_offset - return cell_offset + ijk -end - - -function cpgrid_primitives(coord, zcorn, cartdims; actnum = missing, pinch = missing) - # Add all lines that have at least one active neighbor - coord = reshape(coord, 6, :)' - nx, ny, nz = cartdims - if ismissing(actnum) - actnum = Array{Bool, 3}(undef, nx, ny, nz) - @. actnum = true - end - # Create bounding box used to get boundaries correct - z_mean = max(sum(zcorn)/length(zcorn), 1.0) - z_max = maximum(zcorn) + z_mean - z_min = minimum(zcorn) - z_mean - # actnum, pinch_map = handle_pinch!(actnum, zcorn, cartdims, pinch) - # remapped_indices = findall(vec(actnum)) - nactive = sum(vec(actnum)) - remapped_indices = Vector{Int}(undef, nx*ny*nz) - tmp = vec(actnum) - active_cell_indices = findall(isequal(1), tmp) - @. remapped_indices[tmp] = 1:nactive - - nlinex = nx+1 - nliney = ny+1 - @assert nliney*nlinex == size(coord, 1) - - function generate_line(p1, p2) - return ( - z = Vector{Float64}(), - cells = Vector{Int}(), - cellpos = Vector{Int}(), - nodes = Vector{Int}(), - x1 = SVector{3, Float64}(p1), - x2 = SVector{3, Float64}(p2), - equal_points = p1 ≈ p2 - ) - end - # active_lines = BitArray(undef, nlinex, nliney) - x1, x2 = get_line(coord, 1, 1, nlinex, nliney) - line0 = generate_line(x1, x2) - function boundary_index(i, j, is_top) - if is_top - layer_offset = -(2*nx*ny*nz) - else - layer_offset = -(nx*ny*nz) - end - return layer_offset - ij_to_linear(i, j, cartdims[1:2]) - end - - function cell_index(i, j, k) - ix = ijk_to_linear(i, j, k, cartdims) - if actnum[i, j, k] - cell = remapped_indices[ix] - # elseif haskey(pinch_map, ix) - # cell = remapped_indices[pinch_map[ix]] - else - cell = -ix - @assert cell <= 0 - end - return cell - end - - linear_line_ix(i, j) = ij_to_linear(i, j, (nlinex, nliney)) - lines = Matrix{typeof(line0)}(undef, nlinex, nliney) - for i in 1:nlinex - for j in 1:nliney - p1, p2 = get_line(coord, i, j, nlinex, nliney) - lines[i, j] = generate_line(p1, p2) - end - end - for i = 1:nx - for j = 1:ny - for k = 1:nz - ix = ijk_to_linear(i, j, k, cartdims) - active_cell_index = cell_index(i, j, k) - for I1 in (0, 1) - for I2 in (0, 1) - L = lines[i + I2, j + I1] - for I3 in (0, 1) - zcorn_ix = corner_index(ix, (I1, I2, I3), cartdims) - c = zcorn[zcorn_ix] - # Note reversed indices, this is a bit of a mess - push!(L.z, c) - push!(L.cells, active_cell_index) - end - end - end - end - end - end - # Add fake boundary nodes with corresponding cells - for i in 1:(nx+1) - for j in 1:(ny+1) - L = lines[i, j] - z_top = minimum(L.z) - z_bottom = maximum(L.z) - - for i_offset in (-1, 0) - for j_offset in (-1, 0) - I = i+i_offset - J = j+j_offset - if I > 0 && J > 0 - t = boundary_index(I, J, true) - b = boundary_index(I, J, false) - # Top layer - push!(L.z, z_min, z_top) - push!(L.cells, t, t) - # Bottom layer - push!(L.z, z_bottom, z_max) - push!(L.cells, b, b) - end - end - end - end - end - - # Process lines and merge similar nodes - nodes, lines_active = process_lines!(lines) - - # The four lines making up each column - column_lines = Vector{NTuple{4, Int64}}() - - # Tag columns as active or inactive - active_columns = Matrix{Bool}(undef, nx, ny) - for i in 1:nx - for j in 1:ny - is_active = false - for k in 1:nz - is_active = is_active || actnum[i, j, k] - end - active_columns[i, j] = is_active - end - end - # Generate the columns with cell lists - make_column(i, j) = (cells = Int[], i = i, j = j) - cT = typeof(make_column(1, 1)) - col_counter = 1 - columns = Vector{cT}() - column_indices = zeros(Int, nx, ny) - for i in 1:nx - for j in 1:ny - if active_columns[i, j] - col = make_column(i, j) - prev = boundary_index(i, j, true) - push!(col.cells, prev) - for k in 1:nz - cell = cell_index(i, j, k) - if cell != prev - push!(col.cells, cell) - end - prev = cell - end - # Put a boundary at the end - push!(col.cells, boundary_index(i, j, false)) - push!(columns, col) - column_indices[i, j] = col_counter - col_counter += 1 - - ll = linear_line_ix(i, j) - rl = linear_line_ix(i+1, j) - lr = linear_line_ix(i, j+1) - rr = linear_line_ix(i+1, j+1) - push!(column_lines, (ll, rl, rr, lr)) - end - end - end - ncol = length(columns) - ncoll = length(column_lines) - @assert ncol == ncoll "Mismatch in columns ($ncol) and column lines ($ncoll)" - - function get_edge(i, j, t) - if t == :right - p1 = linear_line_ix(i+1, j+1) - p2 = linear_line_ix(i+1, j) - elseif t == :left - p1 = linear_line_ix(i, j) - p2 = linear_line_ix(i, j+1) - elseif t == :upper - p1 = linear_line_ix(i, j+1) - p2 = linear_line_ix(i+1, j+1) - else - @assert t == :lower - p1 = linear_line_ix(i, j) - p2 = linear_line_ix(i+1, j) - end - return (p1, p2) - end - - function get_boundary_edge(self, i, j, t) - return (column = self, pillars = get_edge(i, j, t)) - end - - function get_interior_edge(c1, c2, i, j, t) - return (columns = (c1, c2), pillars = get_edge(i, j, t)) - end - - tmp = get_boundary_edge(1, 1, 1, :left) - column_boundary = Vector{typeof(tmp)}() - - tmp = get_interior_edge(1, 1, 1, 1, :left) - column_neighbors = Vector{typeof(tmp)}() - - for i in 1:nx - for j in 1:ny - if active_columns[i, j] - self = column_indices[i, j] - if i < nx && active_columns[i+1, j] - other = column_indices[i+1, j] - e = get_interior_edge(self, other, i, j, :right) - push!(column_neighbors, e) - else - # Add right edge to boundary - e = get_boundary_edge(self, i, j, :right) - push!(column_boundary, e) - end - if i == 1 - e = get_boundary_edge(column_indices[i, j], i, j, :left) - push!(column_boundary, e) - end - if j < ny && active_columns[i, j+1] - other = column_indices[i, j+1] - e = get_interior_edge(self, other, i, j, :upper) - push!(column_neighbors, e) - else - e = get_boundary_edge(self, i, j, :upper) - push!(column_boundary, e) - end - if j == 1 - e = get_boundary_edge(column_indices[i, j], i, j, :lower) - push!(column_boundary, e) - end - else - if i < nx && active_columns[i+1, j] - e = get_boundary_edge(column_indices[i+1, j], i+1, j, :left) - push!(column_boundary, e) - end - if j < ny && active_columns[i, j+1] - e = get_boundary_edge(column_indices[i, j+1], i, j+1, :lower) - push!(column_boundary, e) - end - end - end - end - - return ( - lines = lines, - lines_active = lines_active, - column_neighbors = column_neighbors, - column_boundary = column_boundary, - column_lines = column_lines, - columns = columns, - active = active_cell_indices, - nodes = nodes, - cartdims = cartdims - ) -end - -function process_lines!(lines) - nodes = Vector{SVector{3, Float64}}() - active_lines = BitVector(undef, length(lines)) - node_counter = 1 - for (line_ix, line) in enumerate(lines) - z = line.z - active = length(z) > 0 && !all(x -> x <= 0, line.cells) - active_lines[line_ix] = active - if active - p = sortperm(z) - @. line.z = z[p] - @. line.cells = line.cells[p] - pos = line.cellpos - push!(pos, 1) - - counter = 1 - for i in 2:length(z) - if z[i] ≈ z[i-1] - counter += 1 - else - push!(pos, pos[end] + counter) - counter = 1 - end - end - push!(pos, pos[end] + counter) - # Sort each set of cells - for i in 1:(length(pos)-1) - ci = view(line.cells, pos[i]:(pos[i+1]-1)) - sort!(ci) - end - ix = pos[1:end-1] - unique_z = line.z[ix] - # Put back the unique points only - resize!(line.z, 0) - for z_i in unique_z - push!(line.z, z_i) - push!(line.nodes, node_counter) - node_counter += 1 - new_node = interp_coord(line, z_i) - push!(nodes, new_node) - end - end - end - - return (nodes, active_lines) -end - -function grid_from_primitives(primitives) - (; - lines, - lines_active, - # column_neighbors, - column_lines, - columns, - active, - nodes, - cartdims - ) = primitives - if sum(active) == 0 - error("Grid has no active cells.") - end - # Faces mapping to nodes - faces = Vector{Int}() - face_pos = [1] - faceno = 1 - - cell_is_boundary(x) = x < 1 - # Boundary faces mapping to nodes - boundary_faces = Vector{Int}() - boundary_face_pos = [1] - boundary_faceno = 1 - - # Mapping from cell to faces - cell_faces = Vector{Vector{Int}}() - # Mapping from cell to boundary faces - cell_boundary_faces = Vector{Vector{Int}}() - - for c in eachindex(active) - push!(cell_faces, Vector{Int}()) - push!(cell_boundary_faces, Vector{Int}()) - end - face_neighbors = Vector{Tuple{Int, Int}}() - boundary_cells = Vector{Int}() - - vertical_face_tag = Vector{Int}() - horizontal_face_tag = Vector{Int}() - - vertical_bnd_face_tag = Vector{Int}() - horizontal_bnd_face_tag = Vector{Int}() - - nx, ny, nz = cartdims - - nlinex = nx+1 - nliney = ny+1 - - extra_node_lookup = Dict() - - function add_face_from_nodes!(V, Vpos, nodes) - n_global_nodes = length(primitives.nodes) - n_local_nodes = length(nodes) - @assert n_local_nodes > 2 - for n in nodes - @assert n <= n_global_nodes - @assert n > 0 - push!(V, n) - end - push!(Vpos, length(nodes) + Vpos[end]) - end - - function insert_boundary_face!(prev_cell, cell, nodes, is_vertical) - orient = cell_is_boundary(prev_cell) && !cell_is_boundary(cell) - @assert orient || (cell_is_boundary(cell) && !cell_is_boundary(prev_cell)) "cell pair $((cell, prev_cell)) is not on boundary" - if orient - self = cell - else - self = prev_cell - nodes = reverse(nodes) - end - add_face_from_nodes!(boundary_faces, boundary_face_pos, nodes) - push!(cell_boundary_faces[self], boundary_faceno) - push!(boundary_cells, self) - if is_vertical - push!(vertical_bnd_face_tag, boundary_faceno) - else - push!(horizontal_bnd_face_tag, boundary_faceno) - end - - boundary_faceno += 1 - end - - function insert_interior_face!(prev_cell, cell, nodes, is_vertical) - @assert cell > 0 - @assert prev_cell > 0 - @assert prev_cell != cell - add_face_from_nodes!(faces, face_pos, nodes) - # Note order here. - push!(face_neighbors, (prev_cell, cell)) - push!(cell_faces[cell], faceno) - push!(cell_faces[prev_cell], faceno) - if is_vertical - push!(vertical_face_tag, faceno) - else - push!(horizontal_face_tag, faceno) - end - faceno += 1 - end - - function insert_face!(prev_cell, cell, nodes; is_boundary, is_vertical) - if is_boundary - insert_boundary_face!(prev_cell, cell, nodes, is_vertical) - else - insert_interior_face!(prev_cell, cell, nodes, is_vertical) - end - end - - # Horizontal faces (top/bottom and faces along column) - node_buffer = Int[] - sizehint!(node_buffer, 10) - for (cl, col) in zip(column_lines, columns) - number_of_cells_in_column = length(col.cells) - current_column_lines = map(l -> lines[l], cl) - for (i, cell) in enumerate(col.cells) - if cell_is_boundary(cell) - continue - end - if i == 1 - prev = 0 - else - prev = col.cells[i-1] - end - if i == number_of_cells_in_column - next = 0 - else - next = col.cells[i+1] - end - top_is_boundary = cell_is_boundary(prev) - bottom_is_boundary = cell_is_boundary(next) - cell_bnds = map(l -> find_cell_bounds(cell, l), current_column_lines) - for is_top in (true, false) - # TODO: c1/c2 definition will have to be modified to pick the - # right normal vector, check this later. - if is_top - if !top_is_boundary - # Avoid adding interior faces twice. - continue - end - is_bnd = top_is_boundary - F = first - c1 = prev - c2 = cell - else - is_bnd = bottom_is_boundary - F = last - c1 = cell - c2 = next - end - # Index into pillars - node_in_pillar_indices = map(F, cell_bnds) - # Then find the global node indices - node_indices = map((line, i) -> line.nodes[i], current_column_lines, node_in_pillar_indices) - insert_face!(c1, c2, node_indices, is_vertical = false, is_boundary = is_bnd) - end - end - end - # primitives.column_boundary - for is_bnd in [true, false] - if is_bnd - col_neighbors = primitives.column_boundary - else - col_neighbors = primitives.column_neighbors - end - for (cols, pillars) in col_neighbors - # Get the pair of lines we are processing - p1, p2 = pillars - l1 = lines[p1] - l2 = lines[p2] - if length(cols) == 1 - a = b = only(cols) - else - a, b = cols - end - - col_a = columns[a] - col_b = columns[b] - - cell_pairs, overlaps = JutulDarcy.traverse_column_pair(col_a, col_b, l1, l2) - int_pairs, int_overlaps, bnd_pairs, bnd_overlaps = split_overlaps_into_interior_and_boundary(cell_pairs, overlaps) - - F_interior = (l, r, node_indices) -> insert_face!(l, r, node_indices, is_boundary = false, is_vertical = true) - F_bnd = (l, r, node_indices) -> insert_face!(l, r, node_indices, is_boundary = true, is_vertical = true) - - if is_bnd - # We are dealing with a boundary column, everything is boundary - add_vertical_cells_from_overlaps!(extra_node_lookup, F_bnd, nodes, int_pairs, int_overlaps, l1, l2) - else - add_vertical_cells_from_overlaps!(extra_node_lookup, F_interior, nodes, int_pairs, int_overlaps, l1, l2) - add_vertical_cells_from_overlaps!(extra_node_lookup, F_bnd, nodes, bnd_pairs, bnd_overlaps, l1, l2) - end - end - end - - function convert_to_flat(v) - flat_vals = Int[] - flat_pos = Int[1] - for cf in v - for face in cf - push!(flat_vals, face) - end - push!(flat_pos, flat_pos[end]+length(cf)) - end - return (flat_vals, flat_pos) - end - - c2f, c2f_pos = convert_to_flat(cell_faces) - b2f, b2f_pos = convert_to_flat(cell_boundary_faces) - - g = UnstructuredMesh( - c2f, - c2f_pos, - b2f, - b2f_pos, - faces, - face_pos, - boundary_faces, - boundary_face_pos, - primitives.nodes, - face_neighbors, - boundary_cells; - structure = CartesianIndex(cartdims[1], cartdims[2], cartdims[3]), - cell_map = primitives.active, - z_is_depth = true - ) - set_mesh_entity_tag!(g, Faces(), :orientation, :horizontal, horizontal_face_tag) - set_mesh_entity_tag!(g, Faces(), :orientation, :vertical, vertical_face_tag) - set_mesh_entity_tag!(g, BoundaryFaces(), :orientation, :horizontal, horizontal_bnd_face_tag) - set_mesh_entity_tag!(g, BoundaryFaces(), :orientation, :vertical, vertical_bnd_face_tag) - return g -end - - -function handle_pinch!(actnum, zcorn, cartdims, pinch) - if ismissing(pinch) - pinch = 0.001 - end - nx, ny, nz = cartdims - # Note: Remapped is from and to logical indices since this can potentially - # modify ACTNUM. - remapped = Dict{Int, Int}() - pinched = fill(false, nz) - pinch_count = 0 - for i in 1:nx - for j in 1:ny - pinched .= false - prev_active = 0 - pinch_found_in_col = false - for k in 1:nz - linear_ix = ijk_to_linear(i, j, k, cartdims) - active = actnum[i, j, k] - if active - get_zcorn(I1, I2, I3) = zcorn[corner_index(linear_ix, (I1, I2, I3), cartdims)] - get_pair(I, J) = (get_zcorn(I, J, 0), get_zcorn(I, J, 1)) - - # Check if cross part is small - # Check if maximum width is small - l_11, t_11 = get_pair(0, 0) - l_12, t_12 = get_pair(0, 1) - l_21, t_21 = get_pair(1, 0) - l_22, t_22 = get_pair(1, 1) - - d_11 = t_11 - l_11 - d_12 = t_12 - l_12 - d_21 = t_21 - l_21 - d_22 = t_22 - l_22 - - pinched_simple = max(d_11, d_12, d_21, d_22) <= pinch - cell_is_pinched = pinched_simple - - pinched[k] = cell_is_pinched - pinch_count += cell_is_pinched - - pinch_found_in_col = pinch_found_in_col || cell_is_pinched - end - end - if !pinch_found_in_col - continue - end - - for k in 1:nz - if !pinched[k] - continue - end - actnum[i, j, k] = false - - if k > 1 - prev_ix = ijk_to_linear(i, j, k-1, cartdims) - linear_ix = ijk_to_linear(i, j, k, cartdims) - # Could have a section of inactive cells, so look up just in case. - if haskey(remapped, prev_ix) - new_index = remapped[prev_ix] - else - if !actnum[i, j, k-1] - continue - end - new_index = prev_ix - end - remapped[linear_ix] = prev_ix - end - end - end - end - return (actnum, remapped) -end diff --git a/src/cpgrid/processing.jl b/src/cpgrid/processing.jl deleted file mode 100644 index 7da8c66e..00000000 --- a/src/cpgrid/processing.jl +++ /dev/null @@ -1,517 +0,0 @@ - -@enum CPGRID_PILLAR_AB_INTERSECTION begin - AB_RANGES_MATCH # A and B have the same top and bottom points on this pillar - AB_OVERLAP_A_FIRST # A and B overlap, with A on top and then B starting after A but continuing further along - AB_OVERLAP_B_FIRST - A_CONTAINS_B # A fully contains B and is larger in both upwards and downwards direction - B_CONTAINS_A - TOP_MATCHES_A_LONG # The top point matches and A goes further down than B - TOP_MATCHES_B_LONG - BOTTOM_MATCHES_A_LONG # The bottom point matches and A goes further up than B - BOTTOM_MATCHES_B_LONG - DISTINCT_A_ABOVE # The ranges do not overlap and A is above B - DISTINCT_B_ABOVE -end - - -function find_cell_bounds(cell, line) - get_pos(i) = line.cellpos[i]:(line.cellpos[i+1]-1) - get_cells(i) = @view line.cells[get_pos(i)] - - n = length(line.cellpos)-1 - start = 0 - for i in 1:n - if cell in get_cells(i) - start = i - break - end - end - @assert start > 0 "$cell start point not found in line $line" - stop = 0 - for i in (n:-1:start) - if cell in get_cells(i) - stop = i - break - end - end - @assert stop > 0 "$cell end point not found in line $line" - return (start, stop) -end - -function cell_top_bottom(cells, line1, line2; check = true) - # out = Vector{Tuple{Tuple{Int, Int}, Tuple{Int, Int}}}() - T = @NamedTuple{cell::Int64, line1::Tuple{Int64, Int64}, line2::Tuple{Int64, Int64}} - prev1 = prev2 = 1 - out = T[] - for cell in cells - pos1 = find_cell_bounds(cell, line1) - pos2 = find_cell_bounds(cell, line2) - line1_skip = prev1 != pos1[1] - line2_skip = prev2 != pos2[1] - if line1_skip || line2_skip - # Add boundary ghost cell - ghost1 = (prev1, pos1[1]) - ghost2 = (prev2, pos2[1]) - push!(out, (cell = 0, line1 = ghost1, line2 = ghost2)) - end - prev1 = pos1[end] - prev2 = pos2[end] - - # Positions in line1, positions in line 2 - push!(out, (cell = cell, line1 = pos1, line2 = pos2)) - end - if check - start1 = start2 = 1 - for (i, o) in enumerate(out) - @assert o.line1[1] == start1 - @assert o.line2[1] == start2 - - start1 = o.line1[2] - start2 = o.line2[2] - end - @assert start1 == out[end].line1[2] - @assert start2 == out[end].line2[2] - end - return out -end - -function add_vertical_cells_from_overlaps!(extra_node_lookup, F, nodes, cell_pairs, overlaps, l1, l2) - complicated_overlap_categories = (DISTINCT_A_ABOVE, DISTINCT_B_ABOVE) - node_pos = Int[] - function global_node_index(l, local_node) - return l.nodes[local_node] - end - function global_node_point(l, local_node::Int) - return nodes[global_node_index(l, local_node)] - end - - function global_node_point_and_index(l, local_node) - ix = global_node_index(l, local_node) - return (nodes[ix], ix) - end - - sizehint!(node_pos, 6) - for (overlap, cell_pair) in zip(overlaps, cell_pairs) - cell_a, cell_b = cell_pair - - # @info "Starting" overlap cell_pair - cat1 = overlap.line1.category - cat2 = overlap.line2.category - - # TODO: Figure out sign - edge1 = overlap.line1.overlap - edge2 = overlap.line2.overlap - n1 = length(edge1) - n2 = length(edge2) - - @assert n1 == 0 || maximum(edge1) <= length(l1.nodes) - @assert n2 == 0 || maximum(edge2) <= length(l2.nodes) - - line1_distinct = cat1 == DISTINCT_A_ABOVE || cat1 == DISTINCT_B_ABOVE - line2_distinct = cat2 == DISTINCT_A_ABOVE || cat2 == DISTINCT_B_ABOVE - - empty!(node_pos) - if cat1 == cat2 == AB_RANGES_MATCH - handle_matching_overlap!(node_pos, edge1, edge2, l1, l2) - else - handle_generic_intersections!(node_pos, extra_node_lookup, nodes, cell_a, cell_b, l1, l2, overlap, global_node_point) - end - if cell_a == cell_b - # If both cells are the same, we are dealing with a mirrored - # column pair at the boundary. - cell_a = 0 - end - if length(node_pos) > 2 - F(cell_a, cell_b, node_pos) - end - end - return nothing -end - -function handle_matching_overlap!(node_pos, edge1, edge2, l1, l2) - for local_node_ix in edge1 - push!(node_pos, l1.nodes[local_node_ix]) - end - for local_node_ix in reverse(edge2) - push!(node_pos, l2.nodes[local_node_ix]) - end -end - -function handle_matching_overlap_with_intersections!(node_pos, extra_node_lookup, nodes, cell_a, cell_b, l1, l2, overlap, global_node_point) - # Full ranges for line 1 - edge1 = overlap.line1.overlap - range1_a = overlap.line1.range_a - range1_b = overlap.line1.range_b - # Full ranges for line 2 - edge2 = overlap.line2.overlap - range2_a = overlap.line2.range_a - range2_b = overlap.line2.range_b - - # Fist traverse down left vertical edge, then find potential crossing point - # at lower horizontal edge, then up right vertical edge and finally - # potential crossing point at top horizontal edge. - - # 1. Left vertical edge overlap - for local_node_ix in edge1 - push!(node_pos, l1.nodes[local_node_ix]) - end - - # 2. Potential crossing point at lower horizontal edge - a1_l = range1_a[end] - b1_l = range1_b[end] - a2_l = range2_a[end] - b2_l = range2_b[end] - a_above_left_b_above_right = a1_l < b1_l && a2_l > b2_l - b_above_left_a_above_right = a1_l > b1_l && a2_l < b2_l - if a_above_left_b_above_right || b_above_left_a_above_right - # Crossing between two top lines - l1_p1 = global_node_point(l1, a1_l) - l1_p2 = global_node_point(l2, a2_l) - l2_p1 = global_node_point(l1, b1_l) - l2_p2 = global_node_point(l2, b2_l) - bottom_edge_node = handle_crossing_node!(extra_node_lookup, nodes, l1_p1, l1_p2, l2_p1, l2_p2) - push!(node_pos, bottom_edge_node) - end - - # 3. Reverse traversal of right edge - for local_node_ix in reverse(overlap.line2.overlap) - push!(node_pos, l2.nodes[local_node_ix]) - end - - # 4. Finally potential crossing point at top horizontal edge - a1_t = range1_a[1] - b1_t = range1_b[1] - a2_t = range2_a[1] - b2_t = range2_b[1] - - a_above_left_b_above_right = a1_t < b1_t && a2_t > b2_t - b_above_left_a_above_right = a1_t > b1_t && a2_t < b2_t - - if a_above_left_b_above_right || b_above_left_a_above_right - l1_p1 = global_node_point(l1, a1_t) - l1_p2 = global_node_point(l2, a2_t) - l2_p1 = global_node_point(l1, b1_t) - l2_p2 = global_node_point(l2, b2_t) - bottom_edge_node = handle_crossing_node!(extra_node_lookup, nodes, l1_p1, l1_p2, l2_p1, l2_p2) - push!(node_pos, bottom_edge_node) - end -end - -function handle_no_overlapping_edges(node_pos, extra_node_lookup, nodes, cell_a, cell_b, l1, l2, overlap, global_node_point) - # Full ranges for line 1 - edge1 = overlap.line1.overlap - range1_a = overlap.line1.range_a - range1_b = overlap.line1.range_b - # Full ranges for line 2 - edge2 = overlap.line2.overlap - range2_a = overlap.line2.range_a - range2_b = overlap.line2.range_b - - first_and_last(x) = (first(x), last(x)) - a1_t, a1_l = first_and_last(range1_a) - a2_t, a2_l = first_and_last(range2_a) - b1_t, b1_l = first_and_last(range1_b) - b2_t, b2_l = first_and_last(range2_b) - - # Lower nodes - l1_p1_l = global_node_point(l1, a1_l) - l1_p2_l = global_node_point(l2, a2_l) - l2_p1_l = global_node_point(l1, b1_l) - l2_p2_l = global_node_point(l2, b2_l) - # Top nodes - l1_p1_t = global_node_point(l1, a1_t) - l1_p2_t = global_node_point(l2, a2_t) - l2_p1_t = global_node_point(l1, b1_t) - l2_p2_t = global_node_point(l2, b2_t) - - # lower to lower - n_ll = handle_crossing_node!(extra_node_lookup, nodes, l1_p1_l, l1_p2_l, l2_p1_l, l2_p2_l) - push!(node_pos, n_ll) - - # lower [1] to top [2] - n_lt = handle_crossing_node!(extra_node_lookup, nodes, l1_p1_l, l1_p2_l, l2_p1_t, l2_p2_t) - push!(node_pos, n_lt) - - # top to top - n_tt = handle_crossing_node!(extra_node_lookup, nodes, l1_p1_t, l1_p2_t, l2_p1_t, l2_p2_t) - push!(node_pos, n_tt) - - # top [1] to lower [2] - n_tl = handle_crossing_node!(extra_node_lookup, nodes, l1_p1_t, l1_p2_t, l2_p1_l, l2_p2_l) - push!(node_pos, n_tl) -end - -function handle_generic_intersections!(node_pos, extra_node_lookup, nodes, cell_a, cell_b, l1, l2, overlap, global_node_point) - # Full ranges for line 1 - edge1 = overlap.line1.overlap - range1_a = overlap.line1.range_a - range1_b = overlap.line1.range_b - # Full ranges for line 2 - edge2 = overlap.line2.overlap - range2_a = overlap.line2.range_a - range2_b = overlap.line2.range_b - - first_and_last(x) = (first(x), last(x)) - a1_t, a1_l = first_and_last(range1_a) - a2_t, a2_l = first_and_last(range2_a) - b1_t, b1_l = first_and_last(range1_b) - b2_t, b2_l = first_and_last(range2_b) - - function pos_diff(a, b) - return (a < b, a == b) - end - - # Top-top matching - at_before_bt_1, matching_tt_1 = pos_diff(a1_t, b1_t) - at_before_bt_2, matching_tt_2 = pos_diff(a2_t, b2_t) - # Low-low matching - al_before_bl_1, matching_ll_1 = pos_diff(a1_l, b1_l) - al_before_bl_2, matching_ll_2 = pos_diff(a2_l, b2_l) - - # Top nodes - a1_t_pt = global_node_point(l1, a1_t) - a2_t_pt = global_node_point(l2, a2_t) - b1_t_pt = global_node_point(l1, b1_t) - b2_t_pt = global_node_point(l2, b2_t) - - line_top_a = (a1_t_pt, a2_t_pt) - line_top_b = (b1_t_pt, b2_t_pt) - - # Lower nodes - a1_l_pt = global_node_point(l1, a1_l) - a2_l_pt = global_node_point(l2, a2_l) - b1_l_pt = global_node_point(l1, b1_l) - b2_l_pt = global_node_point(l2, b2_l) - - line_low_a = (a1_l_pt, a2_l_pt) - line_low_b = (b1_l_pt, b2_l_pt) - - if at_before_bt_1 != at_before_bt_2 && !(matching_tt_1 || matching_tt_2) - @assert a1_t != b1_t - @assert a2_t != b2_t - - n_tt = handle_crossing_node!(extra_node_lookup, nodes, line_top_a, line_top_b) - push!(node_pos, n_tt) - end - - if length(edge1) == 0 - # Check which one is the lowest - # 1_top crossing 2_low (reversal ab_top_1 vs ab_low_1 over pair) - if at_before_bt_1 - # A is above B at left edge. This extra node is formed from A lower - # cutting B upper. - n_tl = handle_crossing_node!(extra_node_lookup, nodes, line_low_a, line_top_b) - else - n_tl = handle_crossing_node!(extra_node_lookup, nodes, line_top_a, line_low_b) - end - push!(node_pos, n_tl) - else - for local_node_ix in edge1 - push!(node_pos, l1.nodes[local_node_ix]) - end - end - - if al_before_bl_1 != al_before_bl_2 && !(matching_ll_1 || matching_ll_2) - @assert a1_l != b1_l - @assert a2_l != b2_l - - # 1_low crossing 2_low (reversal a/b low over pair) - n_ll = handle_crossing_node!(extra_node_lookup, nodes, line_low_a, line_low_b) - push!(node_pos, n_ll) - end - - if length(edge2) == 0 - if at_before_bt_2 - # A is above B at right edge. This extra node is formed from A lower - # cutting B upper. - n_lt = handle_crossing_node!(extra_node_lookup, nodes, line_low_a, line_top_b) - else - n_lt = handle_crossing_node!(extra_node_lookup, nodes, line_top_a, line_low_b) - end - push!(node_pos, n_lt) - else - for local_node_ix in reverse(edge2) - push!(node_pos, l2.nodes[local_node_ix]) - end - end - return node_pos -end - - - -function handle_crossing_node!(extra_node_lookup, nodes, line_a, line_b) - pt = find_crossing_node(line_a, line_b) - return cpgrid_get_or_add_crossing_node!(extra_node_lookup, nodes, pt) -end - -function find_crossing_node(l1, l2) - @assert length(l1) == 2 == length(l2) - return find_crossing_node(l1[1], l1[2], l2[1], l2[2]) -end - -function find_crossing_node(x1, x2, x3, x4) - # Line (x1 -> x2) crossing line (x3 -> x4) - a = x2 - x1 - b = x4 - x3 - c = x3 - x1 - axb = cross(a, b) - cxb = cross(c, b) - nrm = norm(axb, 2)^2 - if nrm == 0.0 - @warn "Bad intercept: $((x1, x2)) to $((x3, x4))" a b c - error() - end - x_intercept = x1 + a*(dot(cxb, axb)/nrm) - return x_intercept -end - -function cpgrid_get_or_add_crossing_node!(extra_node_lookup, nodes, pt) - if haskey(extra_node_lookup, pt) - ix = extra_node_lookup[pt] - else - push!(nodes, pt) - ix = length(nodes) - extra_node_lookup[pt] = ix - end - return ix -end - -function traverse_column_pair(col_a, col_b, l1, l2) - cell_pairs = Tuple{Int, Int}[] - # TODO: Deal with horrible type. We could just compute this on the fly. - overlaps = @NamedTuple{line1::@NamedTuple{category::JutulDarcy.CPGRID_PILLAR_AB_INTERSECTION, overlap::UnitRange{Int64}, range_a::UnitRange{Int64}, range_b::UnitRange{Int64}}, line2::@NamedTuple{category::JutulDarcy.CPGRID_PILLAR_AB_INTERSECTION, overlap::UnitRange{Int64}, range_a::UnitRange{Int64}, range_b::UnitRange{Int64}}}[] - - function find_end(a, b, s::Symbol) - end_a = a[s][2] - end_b = b[s][2] - if end_a == end_b - # Reached end of both. - return (end_a, true, true) - elseif end_a < end_b - # Reached end of a but not b - return (end_a, true, false) - else - # Reached end of b but not a - return (end_b, false, true) - end - end - - ord_a = cell_top_bottom(col_a.cells, l1, l2) - ord_b = cell_top_bottom(col_b.cells, l1, l2) - function get_local_line(pos, is_line1::Bool) - if is_line1 - d = pos.line1 - else - d = pos.line2 - end - start, stop = d - return (start, stop) - end - - function gen_category(t::CPGRID_PILLAR_AB_INTERSECTION, rng, range_a, range_b) - return ( - category = t, - overlap = rng, - range_a = range_a, - range_b = range_b - ) -end - - function determine_overlap(pos_a, pos_b, is_line1) - a_start, a_stop = get_local_line(pos_a, is_line1) - b_start, b_stop = get_local_line(pos_b, is_line1) - return determine_cell_overlap_inside_line(a_start, a_stop, b_start, b_stop) - end - for pos_a in ord_a - for pos_b in ord_b - # @info "Cell pair: $((pos_a.cell, pos_b.cell))" - - t1, overlap_1, range_1a, range_1b = determine_overlap(pos_a, pos_b, true) - t2, overlap_2, range_2a, range_2b = determine_overlap(pos_a, pos_b, false) - - t_equal = t1 == t2 - if t_equal && t1 in (DISTINCT_A_ABOVE, DISTINCT_B_ABOVE) - continue - end - # Unless A > B for both pillars, or B > A for both pillars, - # we have found a face. - push!(overlaps, - ( - line1 = gen_category(t1, overlap_1, range_1a, range_1b), - line2 = gen_category(t2, overlap_2, range_2a, range_2b) - ) - ) - push!(cell_pairs, (pos_a.cell, pos_b.cell)) - end - end - return (cell_pairs, overlaps) -end - -function split_overlaps_into_interior_and_boundary(cell_pairs, overlaps) - cell_is_boundary(x) = x < 1 - interior_pairs = similar(cell_pairs, 0) - interior_overlaps = similar(overlaps, 0) - - boundary_pairs = similar(interior_pairs) - boundary_overlaps = similar(interior_overlaps) - - for (cell_pair, overlap) in zip(cell_pairs, overlaps) - l, r = cell_pair - l_bnd = cell_is_boundary(l) - r_bnd = cell_is_boundary(r) - if l_bnd && r_bnd - # Two inactive cells, can be skipped. - continue - elseif l_bnd || r_bnd - push!(boundary_pairs, cell_pair) - push!(boundary_overlaps, overlap) - else - push!(interior_pairs, cell_pair) - push!(interior_overlaps, overlap) - end - end - return (interior_pairs, interior_overlaps, boundary_pairs, boundary_overlaps) - # return (interior = (interior_pairs, interior_overlaps), boundary = (boundary_pairs, boundary_overlaps)) -end -## -function determine_cell_overlap_inside_line(a_start, a_stop, b_start, b_stop) - a_range = a_start:a_stop - b_range = b_start:b_stop - # @info "Checking overlap" a_range b_range a_start a_stop b_start b_stop - a_starts_before_b = a_start < b_start - b_starts_before_a = b_start < a_start - - if a_range == b_range - out = (AB_RANGES_MATCH, a_range) - elseif a_starts_before_b && b_stop < a_stop - out = (A_CONTAINS_B, b_range) - elseif b_starts_before_a && a_stop < b_stop - out = (B_CONTAINS_A, a_range) - elseif a_start == b_start - # Take the shortest of the two ranges. - if a_stop > b_stop - out = (TOP_MATCHES_A_LONG, b_range) - else - out = (TOP_MATCHES_B_LONG, a_range) - end - elseif a_stop == b_stop - if a_starts_before_b - # A starts earlier and is therefore longer - out = (BOTTOM_MATCHES_A_LONG, b_range) - else - out = (BOTTOM_MATCHES_B_LONG, a_range) - end - elseif b_start > a_stop - out = (DISTINCT_A_ABOVE, 0:-1) - elseif a_start > b_stop - out = (DISTINCT_B_ABOVE, 0:-1) - elseif a_starts_before_b && b_stop > a_stop - out = (AB_OVERLAP_A_FIRST, b_start:a_stop) - elseif b_starts_before_a && a_stop > b_stop - out = (AB_OVERLAP_B_FIRST, a_start:b_stop) - else - error("Programming error") - end - # @info "Result:" out - return (out[1], out[2], a_range, b_range) -end \ No newline at end of file diff --git a/src/cpr.jl b/src/cpr.jl index 855c98a0..40a37b19 100644 --- a/src/cpr.jl +++ b/src/cpr.jl @@ -254,14 +254,14 @@ function apply!(x, cpr::CPRPreconditioner, r, arg...) # Zero out buffer, just in case @. x = 0.0 # presmooth - apply_cpr_smoother!(x, r, buf, smoother, A_ps, cpr.npre) - apply_cpr_pressure_stage!(cpr, cpr_s, r, arg...) + @tic "cpr smoother" apply_cpr_smoother!(x, r, buf, smoother, A_ps, cpr.npre) + @tic "cpr pressure stage" apply_cpr_pressure_stage!(cpr, cpr_s, r, arg...) # postsmooth if cpr.npost > 0 correct_residual_and_increment_pressure!(r, x, cpr_s.p, bz, buf, cpr_s.A_ps) - apply_cpr_smoother!(x, r, buf, smoother, A_ps, cpr.npost, skip_last = true) + @tic "cpr smoother" apply_cpr_smoother!(x, r, buf, smoother, A_ps, cpr.npost, skip_last = true) else - increment_pressure!(x, cpr_s.p, bz) + @tic "p increment" increment_pressure!(x, cpr_s.p, bz) end end @@ -276,7 +276,7 @@ function apply_cpr_smoother!(x, r, buf, smoother, A_ps, n; skip_last = false) end function correct_residual!(r, A, x) - mul!(r, A, x, -1.0, true) + @tic "residual correction" mul!(r, A, x, -1.0, true) end function apply_cpr_pressure_stage!(cpr::CPRPreconditioner, cpr_s::CPRStorage, r, arg...) @@ -291,15 +291,6 @@ function apply_cpr_pressure_stage!(cpr::CPRPreconditioner, cpr_s::CPRStorage, r, end end -function apply_cpr_second_stage!(x, cpr::CPRPreconditioner, cpr_s::CPRStorage, r, arg...) - bz, Δp = cpr_s.block_size, cpr_s.p - # We currently mutate r and it seems ok. Could copy here if needed. - y = r - @tic "r update" correct_residual_for_dp!(y, x, Δp, bz, cpr_s.x_ps, cpr_s.A_ps) - @tic "s apply" apply!(x, cpr.system_precond, y) - @tic "Δp" increment_pressure!(x, Δp, bz) -end - function cpr_p_apply!(Δp, cpr, p_precond, r_p, p_rtol) apply!(Δp, p_precond, r_p) if !isnothing(p_rtol) @@ -501,7 +492,7 @@ function correct_residual_and_increment_pressure!(r, x, Δp, bz, buf, A_ps) # y' = y - A*Δx # x = A \ y' + Δx @. buf = 0 - for i in 1:length(Δp) + @tic "p increment" @inbounds for i in eachindex(Δp) ix = (i-1)*bz + 1 p_i = Δp[i] buf[ix] = p_i diff --git a/src/facility/types.jl b/src/facility/types.jl index e13e65e3..41febe3b 100644 --- a/src/facility/types.jl +++ b/src/facility/types.jl @@ -126,6 +126,10 @@ bottom-hole-pressure constraint. """ struct SurfaceLiquidRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T + function SurfaceLiquidRateTarget(v::T) where T + isfinite(v) || throw(ArgumentError("Rate must be finite, was $v"))' + return new{T}(v) + end end lumped_phases(::SurfaceLiquidRateTarget) = (AqueousPhase(), LiquidPhase()) @@ -139,6 +143,10 @@ rarely injected into the subsurface. Abbreviated as ORAT in some settings. """ struct SurfaceOilRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T + function SurfaceOilRateTarget(v::T) where T + isfinite(v) || throw(ArgumentError("Rate must be finite, was $v"))' + return new{T}(v) + end end lumped_phases(::SurfaceOilRateTarget) = (LiquidPhase(), ) @@ -155,6 +163,10 @@ present. """ struct SurfaceGasRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T + function SurfaceGasRateTarget(v::T) where T + isfinite(v) || throw(ArgumentError("Rate must be finite, was $v"))' + return new{T}(v) + end end lumped_phases(::SurfaceGasRateTarget) = (VaporPhase(), ) @@ -170,6 +182,10 @@ become very high if there is little water present. """ struct SurfaceWaterRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T + function SurfaceWaterRateTarget(v::T) where T + isfinite(v) || throw(ArgumentError("Rate must be finite, was $v"))' + return new{T}(v) + end end lumped_phases(::SurfaceWaterRateTarget) = (AqueousPhase(), ) @@ -184,6 +200,10 @@ Often used for both [`InjectorControl`](@ref) [`ProducerControl`](@ref). """ struct TotalRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T + function TotalRateTarget(v::T) where T + isfinite(v) || throw(ArgumentError("Rate must be finite, was $v"))' + return new{T}(v) + end end Base.show(io::IO, t::TotalRateTarget) = print(io, "TotalRateTarget with value $(t.value) [m^3/s]") diff --git a/src/init/init.jl b/src/init/init.jl index a945c4d9..ae8e5f9c 100644 --- a/src/init/init.jl +++ b/src/init/init.jl @@ -1,9 +1,12 @@ -function equilibriate_state(model, contacts, datum_depth = missing, datum_pressure = JutulDarcy.DEFAULT_MINIMUM_PRESSURE; - cells = missing, - rs = missing, - rv = missing, - kwarg...) +function equilibriate_state(model, contacts, + datum_depth = missing, + datum_pressure = JutulDarcy.DEFAULT_MINIMUM_PRESSURE; + cells = missing, + rs = missing, + rv = missing, + kwarg... + ) model = reservoir_model(model) D = model.data_domain G = physical_representation(D) @@ -170,7 +173,7 @@ function parse_state0_equil(model, datafile) vapoil = JutulDarcy.has_vapoil(model.system) equil = sol["EQUIL"] - nequil = JutulDarcy.InputParser.number_of_tables(datafile, :equil) + nequil = GeoEnergyIO.InputParser.number_of_tables(datafile, :eqlnum) @assert length(equil) == nequil inits = [] inits_cells = [] @@ -215,7 +218,7 @@ function parse_state0_equil(model, datafile) end ix = unique(i -> cap[i], 1:length(cap)) - if nph == 3 && i == 1 + if i == 1 && get_phases(model.system)[1] isa AqueousPhase @. cap *= -1 end push!(pc, (s = s[ix], pc = cap[ix])) @@ -232,7 +235,7 @@ function parse_state0_equil(model, datafile) swcon = fill(krw.connate, ncells_reg) end push!(s_min, swcon) - push!(s_max, non_connate) + push!(s_max, ones(ncells_reg)) @. non_connate -= swcon end if has_oil @@ -252,11 +255,11 @@ function parse_state0_equil(model, datafile) contacts_pc = (goc_pc, ) else contacts = (woc, ) - contacts_pc = (woc_pc, ) + contacts_pc = (-woc_pc, ) end else contacts = (woc, goc) - contacts_pc = (woc_pc, goc_pc) + contacts_pc = (-woc_pc, goc_pc) end if disgas @@ -323,13 +326,7 @@ function parse_state0_equil(model, datafile) end for (k, v) in subinit - for (i, c) in enumerate(cells) - if v isa AbstractVector - init[k][cells] .= v - else - init[k][:, cells] .= v - end - end + fill_subinit!(init[k], cells, v) end end @assert all(touched) "Some cells are not initialized by equil: $(findall(!, touched))" @@ -337,6 +334,23 @@ function parse_state0_equil(model, datafile) return init end +function fill_subinit!(x::Vector, cells, v::Vector) + @assert length(v) == length(cells) + for (i, c) in enumerate(cells) + x[c] = v[i] + end +end + +function fill_subinit!(x::Matrix, cells, v::Matrix) + @assert size(x, 1) == size(v, 1) + @assert size(v, 2) == length(cells) + for (i, c) in enumerate(cells) + for j in axes(x, 1) + x[j, c] = v[j, i] + end + end +end + function init_reference_pressure(pressures, contacts, kr, pc, ref_ix = 2) nph, nc = size(kr) p = zeros(nc) @@ -368,7 +382,7 @@ function determine_hydrostatic_pressures(depths, depth, zmin, zmax, contacts, da I = I_ref else contact = contacts[pos] - datum_pressure_ph = I_ref(contact) + datum_pressure_ph = I_ref(contact) + contacts_pc[pos] I = phase_pressure_depth_table(contact, zmin, zmax, datum_pressure_ph, density_f, ph) pos += 1 end diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl index b65fec1d..0e1fc4d1 100644 --- a/src/input_simulation/data_input.jl +++ b/src/input_simulation/data_input.jl @@ -428,46 +428,6 @@ function parse_state0_direct_assignment(model, datafile) return init end -function mesh_from_grid_section(f, actnum = missing) - if f isa String - f = InputParser.parse_grdecl_file(f) - end - f::AbstractDict - if haskey(f, "GRID") - grid = f["GRID"] - else - grid = f - end - if ismissing(actnum) - actnum = get_effective_actnum(grid) - end - cartdims = grid["cartDims"] - if haskey(grid, "COORD") - coord = grid["COORD"] - zcorn = grid["ZCORN"] - primitives = JutulDarcy.cpgrid_primitives(coord, zcorn, cartdims, actnum = actnum) - G = JutulDarcy.grid_from_primitives(primitives) - else - @assert haskey(grid, "DX") - @assert haskey(grid, "DY") - @assert haskey(grid, "DZ") - @assert haskey(grid, "TOPS") - @warn "DX+DY+DZ+TOPS format is only supported if all cells are equally sized. If you get an error, this is the cause." - @assert all(actnum) - dx = only(unique(grid["DX"])) - dy = only(unique(grid["DY"])) - dz = only(unique(grid["DZ"])) - tops = only(unique(grid["TOPS"])) - G = CartesianMesh(cartdims, cartdims.*(dx, dy, dz)) - # We always want to return an unstructured mesh. - G = UnstructuredMesh(G) - end - if haskey(grid, "FAULTS") - mesh_add_fault_tags!(G, grid["FAULTS"]) - end - return G -end - function parse_reservoir(data_file) grid = data_file["GRID"] cartdims = grid["cartDims"] @@ -527,8 +487,8 @@ function parse_reservoir(data_file) end extra_data_arg[:temperature] = temperature end - satnum = JutulDarcy.InputParser.table_region(data_file, :saturation, active = active_ix) - eqlnum = JutulDarcy.InputParser.table_region(data_file, :equil, active = active_ix) + satnum = GeoEnergyIO.InputParser.get_data_file_cell_region(data_file, :satnum, active = active_ix) + eqlnum = GeoEnergyIO.InputParser.get_data_file_cell_region(data_file, :eqlnum, active = active_ix) domain = reservoir_domain(G; permeability = perm, @@ -543,111 +503,6 @@ function parse_reservoir(data_file) return domain end -function get_effective_actnum(g) - if haskey(g, "ACTNUM") - actnum = copy(g["ACTNUM"]) - else - actnum = fill(true, g["cartDims"]) - end - handle_zero_effective_porosity!(actnum, g) - return actnum -end - -function handle_zero_effective_porosity!(actnum, g) - if haskey(g, "MINPV") - minpv = g["MINPV"] - else - minpv = 1e-6 - end - added = 0 - active = 0 - - if haskey(g, "PORV") - porv = G["PORV"] - for i in eachindex(actnum) - if actnum[i] - pv = porv[i] - active += active - if pv < minpv - added += 1 - actnum[i] = false - end - end - end - elseif haskey(g, "PORO") - if haskey(g, "ZCORN") - zcorn = g["ZCORN"] - coord = reshape(g["COORD"], 6, :)' - cartdims = g["cartDims"] - else - zcorn = coord = cartdims = missing - end - # Have to handle zero or negligble porosity. - if haskey(g, "NTG") - ntg = g["NTG"] - else - ntg = ones(size(actnum)) - end - poro = g["PORO"] - for i in eachindex(actnum) - if actnum[i] - vol = zcorn_volume(g, zcorn, coord, cartdims, i) - pv = poro[i]*ntg[i]*vol - active += active - if pv < minpv - added += 1 - actnum[i] = false - end - end - end - end - @debug "$added disabled cells out of $(length(actnum)) due to low effective pore-volume." - return actnum -end - -function zcorn_volume(g, zcorn, coord, dims, linear_ix) - if ismissing(zcorn) - return 1.0 - end - nx, ny, nz = dims - i, j, k = linear_to_ijk(linear_ix, dims) - - get_zcorn(I1, I2, I3) = zcorn[corner_index(linear_ix, (I1, I2, I3), dims)] - get_pair(I, J) = (get_zcorn(I, J, 0), get_zcorn(I, J, 1)) - function pillar_line(I, J) - x1, x2 = get_line(coord, i+I, j+J, nx+1, ny+1) - return (x1 = x1, x2 = x2, equal_points = false) - end - - function interpolate_line(I, J, L) - pl = pillar_line(I, J) - return interp_coord(pl, L) - end - - l_11, t_11 = get_pair(0, 0) - l_12, t_12 = get_pair(0, 1) - l_21, t_21 = get_pair(1, 0) - l_22, t_22 = get_pair(1, 1) - - pt_11 = interpolate_line(0, 0, l_11) - pt_12 = interpolate_line(0, 1, l_12) - pt_21 = interpolate_line(1, 0, l_21) - pt_22 = interpolate_line(1, 1, l_22) - - - A_1 = norm(cross(pt_21 - pt_11, pt_12 - pt_11), 2) - A_2 = norm(cross(pt_21 - pt_22, pt_12 - pt_22), 2) - area = (A_1 + A_2)/2.0 - - d_11 = t_11 - l_11 - d_12 = t_12 - l_12 - d_21 = t_21 - l_21 - d_22 = t_22 - l_22 - - d_avg = 0.25*(d_11 + d_12 + d_21 + d_22) - return d_avg*area -end - function parse_physics_types(datafile) runspec = datafile["RUNSPEC"] props = datafile["PROPS"] @@ -1064,17 +919,23 @@ function injector_control(sys, streams, name, flag, type, ctype, surf_rate, res_ # This is a bit of a hack. flag = "OPEN" end - if flag == "SHUT" || flag == "STOP" - ctrl = DisabledControl() - lims = nothing - else + well_is_disabled = flag == "SHUT" || flag == "STOP" + if !well_is_disabled @assert flag == "OPEN" "Unsupported well flag: $flag" if ctype == "RATE" is_rate = true - t = TotalRateTarget(surf_rate) + if isfinite(surf_rate) + t = TotalRateTarget(surf_rate) + else + well_is_disabled = true + end elseif ctype == "BHP" is_rate = false - t = BottomHolePressureTarget(bhp) + if isfinite(bhp) + t = BottomHolePressureTarget(bhp) + else + well_is_disabled = true + end else # RESV, GRUP, THP error("$ctype control not supported") @@ -1082,10 +943,14 @@ function injector_control(sys, streams, name, flag, type, ctype, surf_rate, res_ rho, mix = select_injector_mixture_spec(sys, name, streams, type) if is_rate && surf_rate < MIN_ACTIVE_WELL_RATE @debug "Disabling injector $name with $ctype ctrl due to zero rate" surf_rate - ctrl = DisabledControl() - else - ctrl = InjectorControl(t, mix, density = rho) + well_is_disabled = true end + end + if well_is_disabled + ctrl = DisabledControl() + lims = nothing + else + ctrl = InjectorControl(t, mix, density = rho) if is_hist # TODO: This magic number comes from MRST. bhp_lim = convert_to_si(6895.0, :bar) @@ -1312,92 +1177,3 @@ function well_completion_sortperm(domain, wspec, order_t0, wc, dir) @assert length(sorted) == n return sorted end - -function mesh_add_fault_tags!(G::UnstructuredMesh, faults) - ijk = map(x -> cell_ijk(G, x), 1:number_of_cells(G)) - N = G.faces.neighbors - for (fault, specs) in faults - fault_faces = Int[] - for (I, J, K, dir) in specs - IJK = (I, J, K) - @assert length(dir) == 1 || length(dir) == 2 - d = dir[1] - if d == 'X' || d == 'I' - ix_self = 1 - ix_1 = 2 - ix_2 = 3 - elseif d == 'Y' || d == 'J' - ix_self = 2 - ix_1 = 1 - ix_2 = 3 - elseif d == 'Z' || d == 'K' - ix_self = 3 - ix_1 = 1 - ix_2 = 2 - else - error("Bad direction for fault $fault entry: $dir") - end - if length(dir) == 1 || dir[2] == '+' - inc = 1 - else - @assert dir[2] == '-' - inc = -1 - end - range_self = IJK[ix_self] - range_1 = IJK[ix_1] - range_2 = IJK[ix_2] - @assert length(range_self) == 1 - self = range_self[1] - - match_fault_to_faces!(fault_faces, N, ijk, range_1, ix_1, range_2, ix_2, self, ix_self, inc) - end - @debug "Fault $fault: Added $(length(fault_faces)) faces" - set_mesh_entity_tag!(G, Faces(), :faults, Symbol(fault), fault_faces) - end -end - -function match_fault_to_faces!(fault_faces, N, ijk_cells, range_1, ix_1, range_2, ix_2, self, ix_self, inc) - function sorted_tuple(a, b) - if a < b - pair = (a, b) - else - pair = (b, a) - end - end - pair = sorted_tuple(self, self+inc) - face = 0 - for (l, r) in N - face += 1 - - ijk_l = ijk_cells[l] - ijk_r = ijk_cells[r] - - self_l = ijk_l[ix_self] - self_r = ijk_r[ix_self] - - fpair = sorted_tuple(self_l, self_r) - if fpair != pair - continue - end - - # Keep going unless fixed indices match - cr_1 = ijk_r[ix_1] - if !(cr_1 in range_1) - continue - end - cr_2 = ijk_r[ix_2] - if !(cr_2 in range_2) - continue - end - cl_1 = ijk_l[ix_1] - if !(cl_1 in range_1) - continue - end - cl_2 = ijk_l[ix_2] - if !(cl_2 in range_2) - continue - end - - push!(fault_faces, face) - end -end diff --git a/src/multiphase.jl b/src/multiphase.jl index ce30f266..effd337f 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -9,16 +9,8 @@ get_phases(sys::MultiPhaseSystem) = sys.phases flow_system(sys::MultiPhaseSystem) = sys flow_system(sys::CompositeSystem) = sys.systems.flow - - - - number_of_components(sys::ImmiscibleSystem) = number_of_phases(sys) -# function ImmiscibleSystem(phases) -# @assert length(phases) > 1 "System should have at least two phases. For single-phase, use SinglePhaseSystem instead." -# end - # Single-phase phase_names(system) = get_name.(get_phases(system)) diff --git a/test/cpgrid.jl b/test/cpgrid.jl deleted file mode 100644 index 9326203a..00000000 --- a/test/cpgrid.jl +++ /dev/null @@ -1,51 +0,0 @@ -using Test -import JutulDarcy: determine_cell_overlap_inside_line -@testset "corner point pillar point overlap" begin - for start in -10:25 - for increment in 1:25 - top = start - mid = start + increment - bottom = mid + increment - far = bottom + increment - # A_CONTAINS_B and B_CONTAINS_A - # A B - # | | | - # | | = | - # | | | - @test determine_cell_overlap_inside_line(top, bottom, top, bottom) == (JutulDarcy.AB_RANGES_MATCH, top:bottom, top:bottom, top:bottom) - # AB_OVERLAP_A_FIRST and AB_OVERLAP_B_FIRST - # A B - # | x - # | | = | - # | | x - @test determine_cell_overlap_inside_line(mid, bottom, top, mid) == (JutulDarcy.AB_OVERLAP_B_FIRST, mid:mid, mid:bottom, top:mid) - # Reversed case. - @test determine_cell_overlap_inside_line(top, mid, mid, bottom) == (JutulDarcy.AB_OVERLAP_A_FIRST, mid:mid, top:mid, mid:bottom) - # TOP_MATCHES_A_LONG and TOP_MATCHES_B_LONG - # A B - # | | | - # | | = | - # | x - @test determine_cell_overlap_inside_line(top, bottom, top, mid) == (JutulDarcy.TOP_MATCHES_A_LONG, top:mid, top:bottom, top:mid) - # Reversed case. - @test determine_cell_overlap_inside_line(top, mid, top, bottom) == (JutulDarcy.TOP_MATCHES_B_LONG, top:mid, top:mid, top:bottom) - # BOTTOM_MATCHES_A_LONG and BOTTOM_MATCHES_B_LONG - # A B - # | x x - # | | = | - # | | | - @test determine_cell_overlap_inside_line(top, bottom, mid, bottom) == (JutulDarcy.BOTTOM_MATCHES_A_LONG, mid:bottom, top:bottom, mid:bottom) - # Reversed case. - @test determine_cell_overlap_inside_line(mid, bottom, top, bottom) == (JutulDarcy.BOTTOM_MATCHES_B_LONG, mid:bottom, mid:bottom, top:bottom) - # DISTINCT_A_ABOVE and DISTINCT_B_ABOVE - # A B - # | x - # | x - # | = x - # | x - @test determine_cell_overlap_inside_line(top, mid, bottom, far) == (JutulDarcy.DISTINCT_A_ABOVE, 0:-1, top:mid, bottom:far) - # Reversed case. - @test determine_cell_overlap_inside_line(bottom, far, top, mid) == (JutulDarcy.DISTINCT_B_ABOVE, 0:-1, bottom:far, top:mid) - end - end -end diff --git a/test/parser.jl b/test/parser.jl deleted file mode 100644 index 6179fe81..00000000 --- a/test/parser.jl +++ /dev/null @@ -1,19 +0,0 @@ -using Test -import JutulDarcy.InputParser: clean_include_path, parse_defaulted_line -@testset "InputParser" begin - @test clean_include_path("", " MYFILE") == "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", " ./MYFILE") == joinpath("/some/path", "MYFILE") - @test clean_include_path("/some/path", " './MYFILE'") == joinpath("/some/path", "MYFILE") - @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] - @test parse_defaulted_line("5 *", [1, 2]) == [5, 2] - @test parse_defaulted_line("5 2*", [1, 2, 3]) == [5, 2, 3] - @test parse_defaulted_line("5, 2", [1, 2]) == [5, 2] - @test parse_defaulted_line("5 HEI", [1, "Hallo"]) == [5, "HEI"] - @test parse_defaulted_line("2*", [1, "Hallo"]) == [1, "Hallo"] -end diff --git a/test/runtests.jl b/test/runtests.jl index 913ade1e..8f8be6ee 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,18 +36,10 @@ end include("parray.jl") end -@testitem "Parser" begin - include("parser.jl") -end - @testitem "Scalarization" begin include("scalarization.jl") end -@testitem "Corner point grids" begin - include("cpgrid.jl") -end - @testitem "Multi-phase systems basics" begin include("systems.jl") end diff --git a/test/sens_multimodel.jl b/test/sens_multimodel.jl index 9a4fc1c3..1e7f2948 100644 --- a/test/sens_multimodel.jl +++ b/test/sens_multimodel.jl @@ -5,10 +5,10 @@ function well_test_objective(model, state) return 2.0*q[1] + 0.5*q[2] end -function solve_adjoint_forward_test_system(; block_backend = true, kwarg...) +function solve_adjoint_forward_test_system(casename; block_backend = true, kwarg...) Darcy = si_unit(:darcy) states, reports, setup = JutulDarcy.simulate_mini_wellcase( - :immiscible_2ph; + casename; kwarg..., block_backend = block_backend, permeability = 0.65*Darcy @@ -16,8 +16,8 @@ function solve_adjoint_forward_test_system(; block_backend = true, kwarg...) return (setup[:model], setup[:state0], states, reports, setup[:parameters], setup[:forces]) end -function test_optimization_gradient(; use_scaling = true, use_log = false, kwarg...) - model, state0, states, reports, param, forces = solve_adjoint_forward_test_system(; kwarg...) +function test_optimization_gradient(casename = :immiscible_2ph; use_scaling = true, use_log = false, kwarg...) + model, state0, states, reports, param, forces = solve_adjoint_forward_test_system(casename; kwarg...) dt = report_timesteps(reports) ϵ = 1e-3 num_tol = 1e-2