diff --git a/Project.toml b/Project.toml index 35a34de8..6ba47496 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "JutulDarcy" uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a" authors = ["Olav Møyner "] -version = "0.2.7" +version = "0.2.8" [deps] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" @@ -14,6 +15,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" 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" @@ -23,44 +26,47 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" +[weakdeps] +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" +MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" +PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" + +[extensions] +JutulDarcyGLMakieExt = "GLMakie" +JutulDarcyHYPREExt = "HYPRE" +JutulDarcyPartitionedArraysExt = ["PartitionedArrays", "MPI", "HYPRE"] + [compat] AlgebraicMultigrid = "0.5.1" DataStructures = "0.18.13" +DelimitedFiles = "1.9.1" ForwardDiff = "0.10.30" -Jutul = "0.2.10" +HYPRE = "1.4.0" +Jutul = "0.2.12" +Krylov = "0.9.1" LoopVectorization = "0.12.115" MAT = "0.10.3" MultiComponentFlash = "1.1.3" +OrderedCollections = "1.6.2" +Parsers = "2.7.1" +PartitionedArrays = "0.3.2" Polyester = "0.6.11, 0.7.3" PrecompileTools = "1.0.1" -PartitionedArrays = "0.3.2" StaticArrays = "1.4.4" TestItems = "0.1.0" TimerOutputs = "0.5.19" Tullio = "0.3.4" -HYPRE = "1.4.0" julia = "1.7" -Krylov = "0.9.1" [extras] -TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" -TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" -PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" +PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" +TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" [targets] test = ["Test", "TestItemRunner", "TestItems", "HYPRE", "MPI", "PartitionedArrays"] - -[weakdeps] -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" -HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" -PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" -MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" - -[extensions] -JutulDarcyGLMakieExt = "GLMakie" -JutulDarcyHYPREExt = "HYPRE" -JutulDarcyPartitionedArraysExt = ["PartitionedArrays", "MPI", "HYPRE"] diff --git a/ext/JutulDarcyGLMakieExt/well_plots.jl b/ext/JutulDarcyGLMakieExt/well_plots.jl index 89c722ec..650ad22a 100644 --- a/ext/JutulDarcyGLMakieExt/well_plots.jl +++ b/ext/JutulDarcyGLMakieExt/well_plots.jl @@ -100,7 +100,7 @@ function JutulDarcy.plot_well_results(well_data::Vector, time = nothing; start_d else c_key = :Paired_10 end - cmap = cgrad(c_key, nw, categorical=true) + cmap = cgrad(c_key, max(nw, 2), categorical=true) end wellstr = [String(x) for x in wells] diff --git a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl index 566a538c..25396cd4 100644 --- a/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl +++ b/ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl @@ -40,7 +40,7 @@ module JutulDarcyPartitionedArraysExt function Jutul.parray_preconditioner_apply!(global_out, main_prec::CPRPreconditioner{<:BoomerAMGPreconditioner, <:Any}, X, preconditioners, simulator, arg...) global_cell_vector = simulator.storage.distributed_cell_buffer global_buf = simulator.storage.distributed_residual_buffer - map(local_values(X), preconditioners, ghost_values(X)) do x, prec, x_g + @tic "cpr first stage" map(local_values(X), preconditioners, ghost_values(X)) do x, prec, x_g @. x_g = 0.0 JutulDarcy.apply_cpr_first_stage!(prec, x, arg...) nothing @@ -49,7 +49,7 @@ module JutulDarcyPartitionedArraysExt # copy!(global_cell_vector, main_prec.p) p_h = main_prec.p @assert !isnothing(p_h) "CPR is not properly initialized." - map(own_values(global_cell_vector), preconditioners) do ov, prec + @tic "hypre GetValues" map(own_values(global_cell_vector), preconditioners) do ov, prec helper = prec.pressure_precond.data[:assembly_helper] indices = helper.indices indices::Vector{HYPRE.HYPRE_BigInt} @@ -59,29 +59,27 @@ module JutulDarcyPartitionedArraysExt # End unsafe shenanigans # consistent!(global_cell_vector) |> wait - map(own_values(global_buf), own_values(global_cell_vector), preconditioners) do dx, dp, prec + @tic "set dp" map(own_values(global_buf), own_values(global_cell_vector), preconditioners) do dx, dp, prec bz = prec.block_size for i in eachindex(dp) JutulDarcy.set_dp!(dx, bz, dp, i) end nothing end - if main_prec.full_system_correction - mul_ix = nothing - else - mul_ix = 1 + + @tic "correct residual" begin + mul!(X, main_prec.A_ps, global_buf, -1.0, true) + nothing end - full_op = Jutul.parray_linear_system_operator(simulator.storage.simulators, global_buf) - mul!(X, full_op, global_buf, -1.0, true) - map(local_values(global_out), local_values(X), preconditioners, local_values(global_cell_vector), ghost_values(X)) do y, x, prec, dp, x_g + @tic "increment dp" map(local_values(global_out), local_values(X), preconditioners, local_values(global_cell_vector), ghost_values(X)) do y, x, prec, dp, x_g @. x_g = 0.0 apply!(y, prec.system_precond, x, arg...) bz = prec.block_size JutulDarcy.increment_pressure!(y, dp, bz) nothing end - consistent!(global_out) |> wait + @tic "communication" consistent!(global_out) |> wait global_out end @@ -90,7 +88,7 @@ module JutulDarcyPartitionedArraysExt n = sim.storage.nc_process comm = sim.storage.comm if sim.storage[:number_of_processes] > 1 - @assert sim.backend isa Jutul.MPI_PArrayBackend "Cannot use HYPRE with emulated multiple processes." + @assert sim.backend isa Jutul.MPI_PArrayBackend "Cannot use HYPRE with emulated multiple processes. Backend was $(sim.backend)" end function create_hypre_vector() @@ -104,16 +102,25 @@ module JutulDarcyPartitionedArraysExt cpr.r_p = create_hypre_vector() cpr.p = create_hypre_vector() cpr.np = n + if cpr.full_system_correction + mul_ix = nothing + else + mul_ix = 1 + end + global_buf = sim.storage.distributed_residual_buffer + cpr.A_ps = Jutul.parray_linear_system_operator(sim.storage.simulators, global_buf) end A_p = cpr.A_p + A_ps = cpr.A_ps r_p = cpr.r_p x_p = cpr.p map(sim.storage.simulators, preconditioners) do sim, prec - sys = sim.storage.LinearizedSystem - model = sim.model - storage = sim.storage + storage = Jutul.get_simulator_storage(sim) + model = Jutul.get_simulator_model(sim) + sys = storage.LinearizedSystem prec.A_p = A_p + prec.A_ps = A_ps prec.p = x_p prec.r_p = r_p prec.np = n diff --git a/src/InputParser/InputParser.jl b/src/InputParser/InputParser.jl new file mode 100644 index 00000000..9c60e479 --- /dev/null +++ b/src/InputParser/InputParser.jl @@ -0,0 +1,6 @@ +module InputParser + using Parsers, DelimitedFiles, Jutul, OrderedCollections, Dates + 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 new file mode 100644 index 00000000..c07141a6 --- /dev/null +++ b/src/InputParser/deckinput/keywords/grid.jl @@ -0,0 +1,62 @@ + +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) + @warn "COORDSYS skipped." +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COORD}) + coord = parse_deck_matrix(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) + vals = parse_grid_vector(f, get_cartdims(outer_data), Float64) + vals = swap_unit_system!(vals, units, Val(:permeability)) + data["$k"] = vals +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PORO}) + data["PORO"] = parse_grid_vector(f, get_cartdims(outer_data), Float64) +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:NTG}) + data["NTG"] = parse_grid_vector(f, get_cartdims(outer_data), Float64) +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 + diff --git a/src/InputParser/deckinput/keywords/keywords.jl b/src/InputParser/deckinput/keywords/keywords.jl new file mode 100644 index 00000000..131e8512 --- /dev/null +++ b/src/InputParser/deckinput/keywords/keywords.jl @@ -0,0 +1,8 @@ +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 new file mode 100644 index 00000000..b0eba434 --- /dev/null +++ b/src/InputParser/deckinput/keywords/props.jl @@ -0,0 +1,88 @@ + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:RPTPROPS}) + read_record(f) +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, ::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{:PVTW}) + rec = read_record(f) + 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 + 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}) + rec = read_record(f) + 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 + 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{:ROCK}) + rec = read_record(f) + 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 + 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{:DENSITY}) + rec = read_record(f) + tdims = [NaN, NaN, NaN] + nreg = number_of_tables(outer_data, :pvt) + out = [] + for i = 1:nreg + t = parse_defaulted_line(rec, tdims) + swap_unit_system!(t, units, :density) + push!(out, t) + end + data["DENSITY"] = out +end diff --git a/src/InputParser/deckinput/keywords/regions.jl b/src/InputParser/deckinput/keywords/regions.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/InputParser/deckinput/keywords/runspec.jl b/src/InputParser/deckinput/keywords/runspec.jl new file mode 100644 index 00000000..8995e875 --- /dev/null +++ b/src/InputParser/deckinput/keywords/runspec.jl @@ -0,0 +1,108 @@ +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{: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{:NONNC}) + data["NONNC"] = true +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, ::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{:NUMRES}) + 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{: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{: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{: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 diff --git a/src/InputParser/deckinput/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl new file mode 100644 index 00000000..583d1355 --- /dev/null +++ b/src/InputParser/deckinput/keywords/schedule.jl @@ -0,0 +1,227 @@ +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"] + utypes = [:id, :id, :id, :id, :length, :id, :length, :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{:COMPDAT}) + d = "Default" + defaults = [d, -1, -1, -1, -1, "OPEN", -1, NaN, NaN, NaN, 0.0, -1, "Z", -1] + 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{:WCONPROD}) + d = "Default" + defaults = [d, "OPEN", d, Inf, Inf, Inf, Inf, NaN, 0.0, 0, 0.0] + wells = get_wells(outer_data) + wconprod = 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[7] = :liquid_rate_surface + utypes[8] = :liquid_rate_reservoir + utypes[9] = :pressure + utypes[10] = :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, Inf, Inf, NaN, Inf, 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 + 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, ::Union{Val{:WELLTARG}, Val{:WELTARG}}) + defaults = ["Default", "ORAT", NaN] + wells = get_wells(outer_data) + out = parse_defaulted_group_well(f, defaults, wells, 1) + if cfg.warn + @warn "WELLTARG not supported." + end + data["WELTARG"] = out +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WEFAC}) + defaults = ["Default", 1.0] + wells = get_wells(outer_data) + d = parse_defaulted_group_well(f, defaults, wells, 1) + if cfg.warn + @warn "WEFAC not fully handled." + end + data["WEFAC"] = d +end + +function convert_date_kw(t) + @assert length(t) == 4 + function get_month(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" + 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" + 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 + 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 + diff --git a/src/InputParser/deckinput/keywords/solution.jl b/src/InputParser/deckinput/keywords/solution.jl new file mode 100644 index 00000000..cbd5dbb3 --- /dev/null +++ b/src/InputParser/deckinput/keywords/solution.jl @@ -0,0 +1,67 @@ +# 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 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 + gdata["CURRENT_BOX"] = (lower = (1, 1, 1), upper = dim) +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{: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{: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}) + data["ACTNUM"] = parse_grid_vector(f, get_cartdims(outer_data), 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{:EQUIL}) + n = number_of_tables(outer_data, :equil) + def = [0.0, NaN, 0.0, 0.0, 0.0, 0.0, 0, 0, 0] + eunits = (:length, :pressure, :length, :pressure, :length, :pressure, :id, :id, :id) + 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 diff --git a/src/InputParser/deckinput/keywords/special.jl b/src/InputParser/deckinput/keywords/special.jl new file mode 100644 index 00000000..56b1acc6 --- /dev/null +++ b/src/InputParser/deckinput/keywords/special.jl @@ -0,0 +1,94 @@ + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:COPY}) + 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, d, il, iu, jl, ju, kl, ku]) + src = parsed[1] + dst = parsed[2] + @assert src != "Default" "Source was defaulted? rec = $rec" + @assert dst != "Default" "Destination was defaulted? rec = $rec" + + # Box can be kept. + il = parsed[3] + iu = parsed[4] + jl = parsed[5] + ju = parsed[6] + kl = parsed[7] + ku = parsed[8] + apply_copy!(data, dst, data[src], (il, iu), (jl, ju), (kl, ku), dims) + rec = read_record(f) + end +end + +function apply_copy!(data, dst, src, I, J, K, dims) + if haskey(data, dst) + error("Not implemented") + else + @assert I[1] == J[1] == K[1] == 1 + @assert I[2] == dims[1] + @assert J[2] == dims[2] + @assert K[2] == dims[3] + + data[dst] = copy(src) + end +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTIPLY}) + # TODO: Merge shared code with COPY + rec = read_record(f) + l, u = outer_data["GRID"]["CURRENT_BOX"] + dims = outer_data["GRID"]["cartDims"] + + 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] + @assert haskey(data, dst) "Unable to apply MULTIPLY to non-declared field $dst" + apply_multiply!(data[dst], factor, (il, iu), (jl, ju), (kl, ku), dims) + 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 diff --git a/src/InputParser/deckinput/keywords/summary.jl b/src/InputParser/deckinput/keywords/summary.jl new file mode 100644 index 00000000..36b34b11 --- /dev/null +++ b/src/InputParser/deckinput/keywords/summary.jl @@ -0,0 +1,102 @@ + +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 new file mode 100644 index 00000000..d8cd7bcf --- /dev/null +++ b/src/InputParser/deckinput/parser.jl @@ -0,0 +1,112 @@ +include("units.jl") +include("utils.jl") +include("keywords/keywords.jl") + +""" + parse_data_file(filename; unit = :si) + +Parse a .DATA file (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. + +NOTE: This function is experimental and only covers a small portion of the +keywords that exist for various simulators. +""" +function parse_data_file(filename; kwarg...) + outer_data = Dict{String, Any}() + parse_data_file!(outer_data, filename; kwarg...) + 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], + skip = [:SUMMARY], + units::Union{Symbol, Nothing} = :si, + 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) + basedir = dirname(filename) + function msg(s) + if verbose + println("PARSER: $s") + end + end + msg("Starting parse of $filename...") + f = open(filename, "r") + + cfg = (warn = true, ) + try + allsections = vcat(sections, skip) + while !eof(f) + m = next_keyword!(f) + if isnothing(m) + # Do nothing + elseif m in allsections + msg("Starting section $m") + data = new_section(outer_data, m) + skip_mode = m in skip + # 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 + elseif m == :INCLUDE + next = strip(readline(f)) + if endswith(next, '/') + next = rstrip(next, '/') + else + readline(f) + end + include_path = clean_include_path(basedir, next) + msg("Including file $include_path...") + parse_data_file!( + outer_data, include_path, data, + verbose = verbose, + sections = sections, + skip = skip, + skip_mode = skip_mode, + units = units, + unit_systems = unit_systems + ) + elseif m in (:DATES, :TIME, :TSTEP) + 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! + break + elseif !skip_mode + parse_keyword!(data, outer_data, unit_systems, cfg, f, Val(m)) + end + end + finally + close(f) + end + return outer_data +end + diff --git a/src/InputParser/deckinput/units.jl b/src/InputParser/deckinput/units.jl new file mode 100644 index 00000000..49ee4f32 --- /dev/null +++ b/src/InputParser/deckinput/units.jl @@ -0,0 +1,270 @@ +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 + 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] + 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 + 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 + 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, + 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 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) + 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) + 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) + return swap_unit_system(val, systems, Val(k)) +end + +function swap_unit_system(v, systems::Nothing, k) + # 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}) where k + (; to, from) = systems + 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{: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 diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl new file mode 100644 index 00000000..5bcceb3d --- /dev/null +++ b/src/InputParser/deckinput/utils.jl @@ -0,0 +1,409 @@ +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 = strip(readline(f)) + if !startswith(line, "--") + if endswith(line, '/') + line = strip(rstrip(line, '/')) + active = false + end + 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 + parsed = parse_defaulted_line(line, defaults) + 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) + return parse_defaulted_line([lines], defaults) +end + +function parse_defaulted_line(lines, defaults) + out = similar(defaults, 0) + sizehint!(out, length(defaults)) + pos = 1 + for line in lines + lsplit = split(strip(line), DECK_SPLIT_REGEX) + for s in lsplit + if length(s) == 0 + continue + end + if occursin('*', s) && !startswith(s, '\'') # Could be inside a string for wildcard matching + if s == "*" + num_defaulted = 1 + else + num_defaulted = Parsers.parse(Int, match(r"\d+\*", s).match[1:end-1]) + end + for i in 1:num_defaulted + push!(out, defaults[pos]) + pos += 1 + end + else + default = defaults[pos] + if default isa String + converted = strip(s, [' ', '\'']) + else + T = typeof(default) + converted = Parsers.parse(T, s) + end + push!(out, converted) + pos += 1 + end + end + end + n = length(defaults) + pos = length(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 + m = length(seg) + if m == 0 + continue + elseif n == -1 + n = m + else + @assert m == n "Expected $n was $m" + end + for d in seg + push!(data, parse(T, d)) + end + end + return reshape(data, n, length(data) ÷ n)' +end + +function preprocess_delim_records(split_lines) + # Strip end whitespace + split_lines = map(strip, split_lines) + # Remove comments + filter!(x -> !startswith(x, "--"), split_lines) + # Split into entries (could be whitespace + a comma anywhere in between) + split_rec = map(x -> split(x, r"\s*,?\s+"), split_lines) + # Remove entries + for recs in split_rec + filter!(x -> length(x) > 0, recs) + 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 + for i in 1:n_rep + parsed = parse(T, el)::T + 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 + push!(out, parse_deck_matrix(f)) + end + return out +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{T}) where T + # Do nothing + error("Unhandled keyword $T encountered.") +end + +function next_keyword!(f) + m = nothing + while isnothing(m) && !eof(f) + line = readline(f) + m = keyword_start(line) + end + return m +end + +function number_of_tables(outer_data, t::Symbol) + rs = outer_data["RUNSPEC"] + td = rs["TABDIMS"] + 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 + else + error(":$t is not known") + end +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"] + if t == :saturation + d = reg["SATNUM"] + elseif t == :pvt + d = reg["PVTNUM"] + elseif t == :equil + d = reg["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) + include_file_name = strip(include_file_name) + include_file_name = replace(include_file_name, "./" => "") + include_file_name = replace(include_file_name, "'" => "") + include_path = joinpath(basedir, include_file_name) + return include_path +end + +function get_section(outer_data, name::Symbol) + s = "$name" + is_sched = name == :SCHEDULE + T = OrderedDict{String, Any} + if is_sched + if !haskey(outer_data, s) + outer_data[s] = Dict("STEPS" => [T()], "WELSPECS" => 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 diff --git a/src/JutulDarcy.jl b/src/JutulDarcy.jl index fdef38c5..1e059813 100644 --- a/src/JutulDarcy.jl +++ b/src/JutulDarcy.jl @@ -74,7 +74,12 @@ module JutulDarcy include("porousmedia.jl") # MRST inputs and test cases that use MRST input - include("mrst_input.jl") + # and .DATA file simulation + 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") @@ -91,6 +96,11 @@ module JutulDarcy include("ext.jl") + include("InputParser/InputParser.jl") + using .InputParser + import .InputParser: parse_data_file + export 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/blackoil/blackoil.jl b/src/blackoil/blackoil.jl index 79044ce6..505bc297 100644 --- a/src/blackoil/blackoil.jl +++ b/src/blackoil/blackoil.jl @@ -103,3 +103,53 @@ function cnv_mb_errors_bo(r, Φ, b, dt, rhoS, ::Val{N}) where N end return (Tuple(cnv), Tuple(mb)) end + +function handle_alternate_primary_variable_spec!(init, found, sys::StandardBlackOilSystem) + # Internal utility to handle non-trivial specification of primary variables + nph = number_of_phases(sys) + @assert haskey(init, :Pressure) + @assert haskey(init, :Saturations) || haskey(init, :BlackOilUnknown) + + if nph == 3 && !haskey(init, :ImmiscibleSaturation) + S = init[:Saturations] + a, l, v = phase_indices(sys) + sw = S[a, :] + init[:ImmiscibleSaturation] = sw + push!(found, :ImmiscibleSaturation) + end + + if !haskey(init, :BlackOilUnknown) + S = init[:Saturations] + pressure = init[:Pressure] + nc = length(pressure) + if nph == 2 + sw = zeros(nc) + l, v = phase_indices(sys) + else + a, l, v = phase_indices(sys) + sw = S[a, :] + end + so = S[l, :] + sg = S[v, :] + + F_rs = sys.rs_max + F_rv = sys.rv_max + if has_disgas(sys) + rs = init[:Rs] + else + rs = zeros(nc) + end + if has_vapoil(sys) + rv = init[:Rs] + else + rv = zeros(nc) + end + so = @. 1.0 - so - sg + bo = map( + (w, o, g, r, v, p) -> blackoil_unknown_init(F_rs, F_rv, w, o, g, r, v, p), + sw, so, sg, rs, rv, pressure) + init[:BlackOilUnknown] = bo + push!(found, :BlackOilUnknown) + end + return init +end diff --git a/src/cpgrid/cpgrid.jl b/src/cpgrid/cpgrid.jl new file mode 100644 index 00000000..9825cd1c --- /dev/null +++ b/src/cpgrid/cpgrid.jl @@ -0,0 +1,708 @@ +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::Tuple, z) + p0, p1 = line + return interp_coord(p0, p1, z) +end + +function interp_coord(line::NamedTuple, z) + return interp_coord(line.x1, line.x2, z) +end + +function interp_coord(p0::SVector{3, T}, p1::SVector{3, T}, z::T) where T<:Real + z0 = p0[3] + z1 = p1[3] + + weight = (z - z0)/(z1 - z0) + + interp_pt = p0 .+ weight.*(p1 .- p0) + @assert interp_pt[3] ≈ z "expected $z was $(interp_pt[3])" + 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 + # offset = 2*prod(dims) + # near/far, left/right, up/down + if corner == (0, 0, 0) + major = 0 + minor = 1 + elseif corner == (0, 0, 1) + major = 2 + minor = 1 + elseif corner == (0, 1, 1) + major = 2 + minor = 2 + elseif corner == (0, 1, 0) + major = 0 + minor = 2 + elseif corner == (1, 0, 0) + major = 1 + minor = 1 + elseif corner == (1, 0, 1) + major = 3 + minor = 1 + elseif corner == (1, 1, 1) + major = 3 + minor = 2 + elseif corner == (1, 1, 0) + major = 1 + minor = 2 + else + error("Unsupported $corner_type") + end + planar_offset = 8*nx*ny*(k-1) + top_bottom_offset = major*2*nx*ny + y_offset = 2*nx*(j-1) + x_offset = 2*(i-1) + return planar_offset + top_bottom_offset + y_offset + x_offset + minor +end + + +function cpgrid_primitives(coord, zcorn, cartdims; actnum = missing) + # Add all lines that have at least one active neighbor + nx, ny, nz = cartdims + if ismissing(actnum) + actnum = Array{Bool, 3}(undef, nx, ny, nz) + @. actnum = true + end + # 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) + ) + end + # active_lines = BitArray(undef, nlinex, nliney) + x1, x2 = get_line(coord, 1, 1, nlinex, nliney) + line0 = generate_line(x1, x2) + function cell_index(i, j, k) + ix = ijk_to_linear(i, j, k, cartdims) + if actnum[i, j, k] + cell = remapped_indices[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) + get_corner(pos) = zcorn[corner_index(ix, pos, cartdims)] + for I1 in (0, 1) + for I2 in (0, 1) + L = lines[i + I2, j + I1] + for I3 in (0, 1) + c = get_corner((I1, I2, I3)) + # Note reversed indices, this is a bit of a mess + push!(L.z, c) + push!(L.cells, cell_index(i, j, k)) + end + 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}}() + # for i = 1:nx + # for j = 1:ny + # 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) + + # ll_act = lines_active[ll] + # rl_act = lines_active[rl] + # lr_act = lines_active[lr] + # rr_act = lines_active[rr] + + # if ll_act && rl_act && lr_act && rr_act + # push!(column_lines, (ll, rl, rr, lr)) + # end + # end + # end + # 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[typemin(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 = 0 + 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 + if col.cells[end] > 0 + push!(col.cells, -(nx*ny*nz+1)) + end + 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) + if i == 1 + e = get_boundary_edge(column_indices[i, j], i, j, :left) + push!(column_boundary, e) + end + else + # Add right edge to boundary + e = get_boundary_edge(self, i, j, :right) + 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) + if j == 1 + e = get_boundary_edge(column_indices[i, j], i, j, :lower) + push!(column_boundary, e) + end + else + e = get_boundary_edge(self, i, j, :upper) + 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 + # 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}() + + nx, ny, nz = cartdims + + nlinex = nx+1 + nliney = ny+1 + + function add_face_from_nodes!(V, Vpos, nodes) + for n in nodes + push!(V, n) + end + push!(Vpos, length(nodes) + Vpos[end]) + end + + function insert_boundary_face!(prev_cell, cell, nodes) + orient = cell_is_boundary(prev_cell) && !cell_is_boundary(cell) + @assert orient || (cell_is_boundary(cell) && !cell_is_boundary(prev_cell)) + 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) + boundary_faceno += 1 + end + + function insert_interior_face!(prev_cell, cell, nodes) + @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) + faceno += 1 + end + + function insert_face!(prev_cell, cell, nodes; is_boundary) + if is_boundary + insert_boundary_face!(prev_cell, cell, nodes) + else + insert_interior_face!(prev_cell, cell, nodes) + end + end + + # Horizontal faces (top/bottom and faces along column) + for (cl, col) in zip(column_lines, columns) + # Traverse each column one by one and then figure out what cells are + # connected, generate faces and push to the respective arrays + ncorners = length(cl) + line_ptr = ones(Int, ncorners) + self_lines = map(i -> lines[i], cl) + + + cells = col.cells + @assert cell_is_boundary(cells[1]) + @assert cell_is_boundary(cells[end]) + cell = cells[2] + # Deal with the top first + nodes_first = map(x -> first(x.nodes), self_lines) + n = length(cells) + for i in 2:(n-1) + prev_cell = cells[i-1] + current_cell = cells[i] + next_cell = cells[i+1] + if cell_is_boundary(current_cell) + # Keep going until we find actual cells + continue + end + next_cell_line_ptr!(line_ptr, self_lines, current_cell) + bottom_nodes = map((l, ix) -> l.nodes[ix+1], self_lines, line_ptr) + if cell_is_boundary(prev_cell) + # There was a gap or we were at the top, add boundary + top_nodes = map((l, ix) -> l.nodes[ix], self_lines, line_ptr) + insert_face!(prev_cell, current_cell, top_nodes; is_boundary = true) + end + # Next cell down in column might be a gap or at the bottom + is_bnd = cell_is_boundary(next_cell) + insert_face!(current_cell, next_cell, bottom_nodes; is_boundary = is_bnd) + end + # error() + end + # Vertical faces (between active columns) + function initialize_cpair(pillars, cols, i) + c = columns[cols[i]] + @assert cell_is_boundary(c.cells[1]) + return (lines[pillars[i]], c, c.cells[2]) + end + + + function seek_range(cells, cpos, p, A, B) + foundA = false + foundB = false + for i in cpos[p]:(cpos[p+1]-1) + c = cells[i] + foundA = foundA || c == A + foundB = foundB || c == B + if foundA && foundB + break + end + end + return (foundA, foundB) + end + + function seek_line(line, start, A, B) + cpos = line.cellpos + cells = line.cells + # @info "Looking for a node missing either of A=$A B=$B from $start"# cells + n = length(cpos)-1 + for p in start:n + foundA, foundB = seek_range(cells, cpos, p, A, B) + done = !foundA || !foundB + # @info "?!" p done + if done + # @info "Did find" (p, foundA, foundB) + # We reached the end of our range. + return (p, foundA, foundB) + end + end + # @info "Did not find" (n+1, false, false) + return (n+1, false, false) + end + + node_buf = zeros(Int, 10) + for (cols, pillars) in column_neighbors + lineA, colA, current_cellA = initialize_cpair(pillars, cols, 1) + lineB, colB, current_cellB = initialize_cpair(pillars, cols, 2) + ptrA = ptrB = 1 + colposA = colposB = 2 + # @info "?!" colA lineA + + nA = length(lineA.z) + nB = length(lineB.z) + it = 1 + # @warn "Starting" current_cellA current_cellB nA nB + # for i in 1:nA + # p = lineA.cellpos + # @error "Node $i:" lineA.cells[p[i]:(p[i+1]-1)] + # end + while ptrA < nA && ptrB < nB + + # Both functions should get both cell pairs + next_ptrA, foundA1, foundB1 = seek_line(lineA, ptrA, current_cellA, current_cellB) + next_ptrB, foundA2, foundB2 = seek_line(lineB, ptrB, current_cellA, current_cellB) + + @assert !(foundA1 && foundA2 && foundB1 && foundB2) + @assert foundA1 == foundA2 + @assert foundB1 == foundB2 + # @info "Found?" foundA1 foundB1 ptrA ptrB next_ptrA next_ptrB current_cellA current_cellB + + # Insert new face + bndA = cell_is_boundary(current_cellA) + bndB = cell_is_boundary(current_cellB) + if !bndA || !bndB + Arange = ptrA:(next_ptrA-1) + Brange = ptrB:(next_ptrB-1) + resize!(node_buf, 0) + for i in Arange + push!(node_buf, lineA.nodes[i]) + end + for i in Iterators.reverse(Brange) + push!(node_buf, lineB.nodes[i]) + end + # @info "Inserting face" primitives.nodes[node_buf] Arange Brange + is_bnd = bndA || bndB + insert_face!(current_cellA, current_cellB, node_buf; is_boundary = is_bnd) + end + + # @warn "Iteration $it:" zA zB + if !foundA1 + colposA += 1 + # @warn "Replacing A = $current_cellA with $next_cellA" colA.cells + current_cellA = colA.cells[min(colposA, length(colA.cells))] + end + if !foundB1 + colposB += 1 + # @warn "Replacing B = $current_cellB with $next_cellB" colA.cells + current_cellB = colB.cells[min(colposB, length(colB.cells))] + end + # ptrA:(next_ptrA-1) ptrB:(next_ptrB-1) + # @info "It $it" ptrA next_ptrA foundA1 current_cellA + # @info "It $it" ptrB next_ptrB foundB1 current_cellB + + ptrA = next_ptrA-1 + ptrB = next_ptrB-1 + it += 1 + end + + # error() + end + # Treat boundary + for (col, pillars) in primitives.column_boundary + # lineA, lineB = pillars + lineA, colA, current_cell = initialize_cpair(pillars, (col, col), 1) + lineB, colB, current_cellB = initialize_cpair(pillars, (col, col), 2) + @assert current_cell == current_cellB + ptrA = ptrB = 1 + colpos = 2 + # @info "?!" colA lineA + + nA = length(lineA.z) + nB = length(lineB.z) + it = 1 + # @warn "Starting" current_cellA current_cellB nA nB + # for i in 1:nA + # p = lineA.cellpos + # @error "Node $i:" lineA.cells[p[i]:(p[i+1]-1)] + # end + while ptrA < nA && ptrB < nB + next_ptrA, foundA1, foundB1 = seek_line(lineA, ptrA, current_cell, current_cell) + next_ptrB, foundA2, foundB2 = seek_line(lineB, ptrB, current_cell, current_cell) + + @assert !(foundA1 && foundA2 && foundB1 && foundB2) + @assert foundA1 == foundA2 + @assert foundB1 == foundB2 + + if !cell_is_boundary(current_cell) + Arange = ptrA:(next_ptrA-1) + Brange = ptrB:(next_ptrB-1) + resize!(node_buf, 0) + for i in Arange + push!(node_buf, lineA.nodes[i]) + end + for i in Iterators.reverse(Brange) + push!(node_buf, lineB.nodes[i]) + end + insert_face!(current_cell, 0, node_buf; is_boundary = true) + end + if !foundA1 + colpos = min(colpos+1, length(colA.cells)) + current_cell = colA.cells[colpos] + end + ptrA = next_ptrA-1 + ptrB = next_ptrB-1 + it += 1 + end + + # error() + 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) + + return 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 + ) +end + +function next_cell_line_ptr!(line_ptr, self_lines, next_cell) + for (i, line_i) in enumerate(self_lines) + cells_i = line_i.cells + ptr_i = line_ptr[i] + cellpos = line_i.cellpos + current_i = cellpos[ptr_i] + nc = length(cells_i) + while current_i <= nc && cells_i[current_i] != next_cell + current_i += 1 + end + next_ptr = -1 + for j in ptr_i:(length(cellpos)-1) + if current_i >= cellpos[j] && current_i < cellpos[j+1] + next_ptr = j + break + end + end + if next_ptr == -1 + @info "?!" line_i cells_i next_cell + error("This should not occur.") + end + pos = cellpos[next_ptr]:(cellpos[next_ptr+1]-1) + @assert next_cell in cells_i[pos] + line_ptr[i] = next_ptr + end +end \ No newline at end of file diff --git a/src/deck_types.jl b/src/deck_types.jl index 940d9228..4dd66dfe 100644 --- a/src/deck_types.jl +++ b/src/deck_types.jl @@ -82,6 +82,9 @@ struct ConstMuBTable{R} end function ConstMuBTable(pvtw::M) where M<:AbstractVector + pvtw = flat_region_expand(pvtw) + # Only one region supported atm + pvtw = first(pvtw) return ConstMuBTable(pvtw[1], 1.0/pvtw[2], pvtw[3], pvtw[4], pvtw[5]) end diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl index f50b0bcf..b3ac3aac 100644 --- a/src/facility/wells/wells.jl +++ b/src/facility/wells/wells.jl @@ -46,7 +46,7 @@ function common_well_setup(nr; dz = nothing, WI = nothing, gravity = gravity_con @warn "No well indices provided. Using 1e-12." WI = repeat(1e-12, nr) else - @assert length(WI) == nr "Must have one well index per perforated cell" + @assert length(WI) == nr "Must have one well index per perforated cell ($(length(WI)) well indices, $nr reservoir cells))" end return (WI, gdz) end @@ -89,19 +89,27 @@ function setup_well(g, K, reservoir_cells::AbstractVector; volumes = zeros(n) WI_computed = zeros(n) dz = zeros(n) + + function get_entry(x::AbstractVector, i) + return x[i] + end + function get_entry(x, i) + return x + end for (i, c) in enumerate(reservoir_cells) if K isa AbstractVector - k = K[c] + k_i = K[c] else - k = K[:, c] + k_i = K[:, c] end - if ismissing(WI) - WI_i = compute_peaceman_index(g, k, radius, c, skin = skin, Kh = Kh) - elseif WI isa AbstractFloat - WI_i = WI - else - WI_i = WI[i] + WI_i = get_entry(WI, i) + Kh_i = get_entry(Kh, i) + r_i = get_entry(radius, i) + s_i = get_entry(skin, i) + if ismissing(WI_i) || isnan(WI_i) + WI_i = compute_peaceman_index(g, k_i, r_i, c, skin = s_i, Kh = Kh_i) end + WI_i::AbstractFloat WI_computed[i] = WI_i center = vec(centers[:, i]) dz[i] = center[3] - reference_depth @@ -113,11 +121,13 @@ function setup_well(g, K, reservoir_cells::AbstractVector; Δ = cell_dims(g, c) d_index = findfirst(isequal(d), [:x, :y, :z]) h = Δ[d_index] - volumes[i] = h*π*radius^2 + volumes[i] = h*π*r_i^2 end if simple_well W = SimpleWell(reservoir_cells, WI = WI_computed, dz = dz, reference_depth = reference_depth, kwarg...) else + # Depth differences are taken care of via centers. + dz *= 0.0 W = MultiSegmentWell(reservoir_cells, volumes, centers; WI = WI_computed, dz = dz, reference_depth = reference_depth, kwarg...) end return W diff --git a/src/flux.jl b/src/flux.jl index acae6eb4..6e3beb94 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -43,24 +43,34 @@ end end @inline function darcy_phase_kgrad_potential(face, phase, state, model, flux_type, tpfa::TPFA, upw, common = flux_primitives(face, state, model, flux_type, upw, tpfa)) - ρ = state.PhaseMassDensities pc, ref_index = capillary_pressure(model, state) ∇p, T_f, gΔz = common l = tpfa.left r = tpfa.right Δpc = capillary_gradient(pc, l, r, phase, ref_index) - @inbounds ρ_c = ρ[phase, l] - @inbounds ρ_i = ρ[phase, r] - ρ_avg = 0.5*(ρ_i + ρ_c) + ρ_avg = face_average_density(model, state, tpfa, phase) q = -T_f*(∇p + Δpc + gΔz*ρ_avg) return q end +function face_average_density(model, state, tpfa, phase) + ρ = state.PhaseMassDensities + l = tpfa.left + r = tpfa.right + @inbounds ρ_c = ρ[phase, l] + @inbounds ρ_i = ρ[phase, r] + return 0.5*(ρ_i + ρ_c) +end + @inline function gradient(X, tpfa::TPFA) return @inbounds X[tpfa.right] - X[tpfa.left] end +@inline function gradient(X::AbstractMatrix, i, tpfa::TPFA) + return @inbounds X[i, tpfa.right] - X[i, tpfa.left] +end + pressure_gradient(state, disc) = gradient(state.Pressure, disc) @inline function upwind(upw::SPU, F, q) diff --git a/src/forces/bc.jl b/src/forces/bc.jl index 7e3935d7..1b4bef46 100644 --- a/src/forces/bc.jl +++ b/src/forces/bc.jl @@ -128,3 +128,90 @@ function apply_flow_bc!(acc, q, bc, model::SimulationModel{<:Any, T}, state, tim end end end + +function Jutul.vectorization_length(bc::FlowBoundaryCondition, variant) + if variant == :all + n = 4 # pressure, temperature, T_flow, T_thermal + f = bc.fractional_flow + if !isnothing(f) + n += length(f) + end + if !isnothing(bc.density) + n += 1 + end + return n + elseif variant == :control + return 1 + else + error("Variant $variant not supported") + end +end + +function Jutul.vectorize_force!(v, bc::FlowBoundaryCondition, variant) + names = [] + if variant == :all + v[1] = bc.pressure + push!(names, :pressure) + v[2] = bc.temperature + push!(names, :temperature) + v[3] = bc.trans_flow + push!(names, :trans_flow) + v[4] = bc.trans_thermal + push!(names, :trans_thermal) + offset = 4 + f = bc.fractional_flow + if !isnothing(f) + for (i, f_i) in enumerate(f) + offset += 1 + v[offset] = f_i + push!(names, Symbol("fractional_flow$i")) + end + end + if !isnothing(bc.density) + offset += 1 + v[offset] = bc.density + push!(names, :density) + end + elseif variant == :control + v[1] = bc.pressure + push!(names, :pressure) + else + error("Variant $variant not supported") + end + return (names = names, ) +end + +function Jutul.devectorize_force(bc::FlowBoundaryCondition, X, meta, variant) + p = X[1] + T = bc.temperature + trans_flow = bc.trans_flow + trans_thermal = bc.trans_thermal + f = bc.fractional_flow + ρ = bc.density + if variant == :all + T = X[2] + trans_flow = X[3] + trans_thermal = X[4] + + offset = 4 + if !isnothing(f) + f = X[(offset+1):(offset+length(f))] + end + if !isnothing(ρ) + ρ = X[end] + end + elseif variant == :control + # DO nothing + else + error("Variant $variant not supported") + end + return FlowBoundaryCondition( + bc.cell, + p, + T, + trans_flow, + trans_thermal, + f, + ρ + ) +end diff --git a/src/forces/sources.jl b/src/forces/sources.jl index 721c84d1..69e88bdd 100644 --- a/src/forces/sources.jl +++ b/src/forces/sources.jl @@ -39,3 +39,53 @@ function Jutul.subforce(s::AbstractVector{S}, model) where S<:SourceTerm return s[keep] end +function Jutul.vectorization_length(src::SourceTerm, variant) + n = 1 + if variant == :all + f = src.fractional_flow + if !isnothing(f) + n += length(f) + end + elseif variant == :control + # Do nothing + else + error("Variant $variant not supported") + end + return n +end + +function Jutul.vectorize_force!(v, src::SourceTerm, variant) + v[1] = src.value + names = [:value] + if variant == :all + offset = 1 + f = src.fractional_flow + if !isnothing(f) + for (i, f_i) in enumerate(f) + offset += 1 + v[offset] = f_i + push!(names, Symbol("fractional_flow$i")) + end + end + elseif variant == :control + # Do nothing + else + error("Variant $variant not supported") + end + return (names = names, ) +end + +function Jutul.devectorize_force(src::SourceTerm, X, meta, variant) + val = X[1] + f = src.fractional_flow + if variant == :all + if !isnothing(f) + f = tuple(X[2:end]...) + end + elseif variant == :control + # Do nothing + else + error("Variant $variant not supported") + end + return SourceTerm(src.cell, val, f, src.type) +end diff --git a/src/init/init.jl b/src/init/init.jl new file mode 100644 index 00000000..72a32028 --- /dev/null +++ b/src/init/init.jl @@ -0,0 +1,353 @@ + +function equilibriate_state(model, contacts, datum_depth = missing, datum_pressure = JutulDarcy.DEFAULT_MINIMUM_PRESSURE; + cells = missing, + rs = missing, + kwarg...) + D = model.data_domain + G = physical_representation(D) + cc = D[:cell_centroids][3, :] + if ismissing(cells) + cells = 1:number_of_cells(G) + pts = cc + else + pts = view(cc, cells) + end + + if ismissing(datum_depth) + datum_depth = pmin + end + sys = model.system + + init = Dict{Symbol, Any}() + init = equilibriate_state!(init, pts, model, sys, contacts, datum_depth, datum_pressure; cells = cells, rs = rs, kwarg...) + + is_blackoil = sys isa StandardBlackOilSystem + if is_blackoil + sat = init[:Saturations] + pressure = init[:Pressure] + nph, nc = size(sat) + if ismissing(rs) + Rs = zeros(nc) + else + Rs = rs.(pts) + end + # TODO: Handle Rv. + Rv = zeros(nc) + init[:Rs] = Rs + end + return init +end + +function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pressure; + cells = 1:length(depths), + rs = missing, + contacts_pc = missing, + kwarg... + ) + if ismissing(contacts_pc) + contacts_pc = zeros(number_of_phases(sys)-1) + end + zmin = minimum(depths) + zmax = maximum(depths) + + nph = number_of_phases(sys) + + @assert length(contacts) == nph-1 + + rho_s = JutulDarcy.reference_densities(sys) + phases = JutulDarcy.get_phases(sys) + disgas = JutulDarcy.has_disgas(sys) + if disgas + if JutulDarcy.has_other_phase(sys) + _, rhoOS, rhoGS = rho_s + else + rhoOS, rhoGS = rho_s + end + end + pvt = model.secondary_variables[:PhaseMassDensities].pvt + relperm = model.secondary_variables[:RelativePermeabilities] + function density_f(p, z, ph) + if phases[ph] == LiquidPhase() && disgas + Rs = min(rs(z), sys.rs_max(p)) + b = JutulDarcy.shrinkage(pvt[ph], 1, p, Rs, 1) + rho = b*(rhoOS + Rs*rhoGS) + else + rho = rho_s[ph]*JutulDarcy.shrinkage(only(pvt[ph].tab), p) + end + return rho + end + pressures = determine_hydrostatic_pressures(depths, depth, zmin, zmax, contacts, datum_pressure, density_f, contacts_pc) + s, pc = determine_saturations(depths, contacts, pressures; kwarg...) + + kr = similar(s) + JutulDarcy.update_kr!(kr, relperm, model, s, cells) + + init[:Saturations] = s + init[:Pressure] = init_reference_pressure(pressures, contacts, pc, kr, 2) + return init +end + + +function parse_state0_equil(model, datafile) + sys = model.system + init = Dict{Symbol, Any}() + sol = datafile["SOLUTION"] + G = physical_representation(model.data_domain) + nc = number_of_cells(G) + nph = number_of_phases(model.system) + ix = G.cell_map + is_blackoil = sys isa StandardBlackOilSystem + has_water = haskey(datafile["RUNSPEC"], "WATER") + has_oil = haskey(datafile["RUNSPEC"], "OIL") + has_gas = haskey(datafile["RUNSPEC"], "GAS") + disgas = JutulDarcy.has_disgas(model.system) + + equil = sol["EQUIL"] + nequil = JutulDarcy.InputParser.number_of_tables(datafile, :equil) + @assert length(equil) == nequil + inits = [] + for ereg in 1:nequil + eq = equil[ereg] + cells = findall(isequal(ereg), model.data_domain[:eqlnum]) + datum_depth = eq[1] + datum_pressure = eq[2] + + woc = eq[3] + woc_pc = eq[4] + goc = eq[5] + goc_pc = eq[6] + # Contact depths + s_max = 1.0 + s_min = 0.0 + + non_connate = 1.0 + s_max = Float64[] + s_min = Float64[] + # @assert !haskey(model.secondary_variables, :CapillaryPressure) "Capillary initialization not yet implemented." + has_pc = haskey(model.secondary_variables, :CapillaryPressure) + if has_pc + pc_f = model.secondary_variables[:CapillaryPressure].pc + pc = [] + for (i, f) in enumerate(pc_f) + s = only(f).X + cap = only(f).F + if s[1] < 0 + s = s[2:end] + cap = cap[2:end] + end + ix = unique(i -> cap[i], 1:length(cap)) + + if nph == 3 && i == 1 + @. cap *= -1 + end + push!(pc, (s = s[ix], pc = cap[ix])) + end + else + pc = nothing + end + if has_water + krw = only(model.secondary_variables[:RelativePermeabilities].krw) + swcon = krw.connate + push!(s_min, swcon) + push!(s_max, non_connate) + non_connate = 1.0 - swcon + end + if has_oil + push!(s_min, 0.0) + push!(s_max, non_connate) + end + if has_gas + push!(s_min, 0.0) + push!(s_max, non_connate) + end + + if nph == 1 + error("Not implemented.") + elseif nph == 2 + if has_oil && has_gas + contacts = (goc, ) + contacts_pc = (goc_pc, ) + else + contacts = (woc, ) + contacts_pc = (woc_pc, ) + end + else + contacts = (woc, goc) + contacts_pc = (woc_pc, goc_pc) + end + + if disgas + @assert haskey(sol, "RSVD") + rsvd = sol["RSVD"][ereg] + z = rsvd[:, 1] + Rs = rsvd[:, 2] + rs = Jutul.LinearInterpolant(z, Rs) + else + rs = missing + end + + subinit = equilibriate_state(model, contacts, datum_depth, datum_pressure, cells = cells, contacts_pc = contacts_pc, s_min = s_min, s_max = s_max, rs = rs, pc = pc) + push!(inits, subinit) + end + # TODO: Handle multiple regions by merging each init + return only(inits) +end + +function init_reference_pressure(pressures, contacts, kr, pc, ref_ix = 2) + nph, nc = size(kr) + p = zeros(nc) + for i in eachindex(p) + p[i] = pressures[ref_ix, i] + if kr[ref_ix, i] <= 0 + for ph in 1:nph + if kr[ph, i] > 0 + p[i] = pressures[ph, i]# - pc[ph, i] + end + end + end + end + return p +end + +function determine_hydrostatic_pressures(depths, depth, zmin, zmax, contacts, datum_pressure, density_f, contacts_pc) + nc = length(depths) + nph = length(contacts) + 1 + ref_ix = 2 + I_ref = phase_pressure_depth_table(depth, zmin, zmax, datum_pressure, density_f, ref_ix) + pressures = zeros(nph, nc) + pos = 1 + for ph in 1:nph + if ph == ref_ix + I = I_ref + else + contact = contacts[pos] + datum_pressure_ph = I_ref(contact) + I = phase_pressure_depth_table(contact, zmin, zmax, datum_pressure_ph, density_f, ph) + pos += 1 + end + for (i, z) in enumerate(depths) + pressures[ph, i] = I(z) + end + end + return pressures +end + + + +function phase_pressure_depth_table(depth, zmin, zmax, datum_pressure, density_f, phase) + if zmin > depth + z_up = Float64[] + p_up = Float64[] + else + z_up, p_up = integrate_phase_density(depth, zmin, datum_pressure, density_f, phase) + # Remove top point + z_up = reverse!(z_up[2:end]) + p_up = reverse!(p_up[2:end]) + end + if zmax < depth + z_down = Float64[] + p_down = Float64[] + else + z_down, p_down = integrate_phase_density(depth, zmax, datum_pressure, density_f, phase) + end + z = vcat(z_up, z_down) + p = vcat(p_up, p_down) + return Jutul.LinearInterpolant(z, p) +end + +function integrate_phase_density(z_datum, z_end, p0, density_f, phase; n = 1000, g = Jutul.gravity_constant) + dz = (z_end - z_datum)/n + pressure = zeros(n+1) + z = zeros(n+1) + pressure[1] = p0 + z[1] = z_datum + + for i in 2:(n+1) + p = pressure[i-1] + depth = z[i-1] + dz + pressure[i] = p + dz*density_f(p, depth, phase)*g + z[i] = depth + end + return (z, pressure) +end + +function determine_saturations(depths, contacts, pressures; s_min = missing, s_max = missing, pc = nothing) + nc = length(depths) + nph = length(contacts) + 1 + if ismissing(s_max) + s_max = zeros(nph) + end + if ismissing(s_max) + s_max = zeros(nph) + end + sat = zeros(nph, nc) + sat_pc = similar(sat) + if isnothing(pc) + for i in eachindex(depths) + z = depths[i] + ph = current_phase_index(z, contacts) + for j in axes(sat, 1) + is_main = ph == j + s = is_main*s_max[j] + !is_main*s_min[j] + sat[j, i] = s + end + end + else + ref_ix = 2 + offset = 1 + for ph in 1:nph + if ph != ref_ix + s, pc_pair = pc[offset] + pc_max = maximum(pc_pair) + pc_min = minimum(pc_pair) + + I = Jutul.LinearInterpolant(pc_pair, s) + I_pc = Jutul.LinearInterpolant(s, pc_pair) + for i in eachindex(depths) + z = depths[i] + + dp = pressures[ph, i] - pressures[ref_ix, i] + if dp > pc_max + s_eff = s_max[ph] + elseif dp < pc_min + s_eff = s_min[ph] + else + s_eff = I(dp) + end + s_eff = clamp(s_eff, s_min[ph], s_max[ph]) + sat[ph, i] = s_eff + sat_pc[ph, i] = I_pc(s_eff) + end + offset += 1 + end + end + for i in eachindex(depths) + s_fill = 1 - sum(view(sat, :, i)) + if s_fill < 0 + @warn "Negative saturation in cell $i: $s_fill" + end + sat[ref_ix, i] = s_fill + end + end + return (sat, sat_pc) +end + +function current_phase_index(z, depths; reverse = true) + n = length(depths)+1 + if reverse + i = current_phase_index(z, Base.reverse(depths), reverse = false) + return n - i + 1 + else + if z < depths[1] + return 1 + elseif z > depths[end] + return n + else + for (i, d) in enumerate(depths) + if d > z + return i + end + end + end + end +end diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl new file mode 100644 index 00000000..68b5ae4b --- /dev/null +++ b/src/input_simulation/data_input.jl @@ -0,0 +1,574 @@ +export simulate_data_file, case_from_data_input + +""" + simulate_data_file(inp; parse_arg = NamedTuple(), kwarg...) + +Simulate standard input file (with extension .DATA). `inp` can either be the +output from `case_from_data_input` or a String for the path of an input file. + +Additional arguments are passed onto `simulate_reservoir`. Extra inputs to the +parser can be sent as a `parse_arg` `NamedTuple`. +""" +function simulate_data_file(fn::String; parse_arg = NamedTuple(), kwarg...) + data = parse_data_file(fn; parse_arg...) + return simulate_data_file(data; kwarg...) +end + +function simulate_data_file(data::AbstractDict; setup_arg = NamedTuple(), kwarg...) + case = case_from_data_input(data; setup_arg...) + return simulate_reservoir(case; kwarg...) +end + +function case_from_data_input(datafile; kwarg...) + sys, pvt = parse_physics_types(datafile) + domain = parse_reservoir(datafile) + wells, controls, limits, cstep, dt, well_mul = parse_schedule(domain, sys, datafile) + + model, parameters0 = setup_reservoir_model(domain, sys; wells = wells, kwarg...) + + for (k, submodel) in pairs(model.models) + if submodel.system isa MultiPhaseSystem + # Modify secondary variables + svar = submodel.secondary_variables + # PVT + pvt = tuple(pvt...) + rho = DeckDensity(pvt) + if sys isa StandardBlackOilSystem + set_secondary_variables!(submodel, ShrinkageFactors = JutulDarcy.DeckShrinkageFactors(pvt)) + end + mu = DeckViscosity(pvt) + set_secondary_variables!(submodel, PhaseViscosities = mu, PhaseMassDensities = rho) + if k == :Reservoir + rs = datafile["RUNSPEC"] + oil = haskey(rs, "OIL") + water = haskey(rs, "WATER") + gas = haskey(rs, "GAS") + + JutulDarcy.set_deck_specialization!(submodel, datafile["PROPS"], domain[:satnum], oil, water, gas) + end + end + end + parameters = setup_parameters(model) + forces = parse_forces(model, wells, controls, limits, cstep, dt, well_mul) + state0 = parse_state0(model, datafile) + return JutulCase(model, dt, forces, state0 = state0, parameters = parameters) +end + +function parse_well_from_compdat(domain, wname, v, wspecs) + wc, WI, open = compdat_to_connection_factors(domain, v) + rd = wspecs.ref_depth + if isnan(rd) + rd = nothing + end + W = setup_well(domain, wc, name = Symbol(wname), WI = WI, reference_depth = rd) + return (W, WI, open) +end + +function compdat_to_connection_factors(domain, v) + G = physical_representation(domain) + K = domain[:permeability] + + ij_ix = collect(keys(v)) + wc = map(i -> cell_index(G, i), ij_ix) + + getf(k) = map(x -> v[x][k], ij_ix) + open = getf(:open) + d = getf(:diameter) + Kh = getf(:Kh) + WI = getf(:WI) + skin = getf(:skin) + dir = getf(:dir) + + for i in eachindex(WI) + W_i = WI[i] + if isnan(W_i) || W_i < 0.0 + c = wc[i] + if K isa AbstractVector + k_i = K[c] + else + k_i = K[:, c] + end + WI[i] = compute_peaceman_index(G, k_i, d[i]/2, c, skin = skin[i], Kh = Kh[i], dir = Symbol(dir[i])) + end + end + return (wc, WI, open) +end + +function parse_schedule(domain, runspec, schedule, sys) + dt, cstep, controls, completions, limits = parse_control_steps(runspec, schedule, sys) + ncomp = length(completions) + wells = [] + well_mul = [] + for (k, v) in pairs(completions[end]) + W, WI, open = parse_well_from_compdat(domain, k, v, schedule["WELSPECS"][k]) + wi_mul = zeros(length(WI), ncomp) + for (i, c) in enumerate(completions) + compdat = c[k] + _, WI, open = compdat_to_connection_factors(domain, compdat) + @. wi_mul[:, i] = (WI*open)/WI + end + push!(wells, W) + if all(isapprox(1.0), wi_mul) + push!(well_mul, nothing) + else + push!(well_mul, wi_mul) + end + end + return (wells, controls, limits, cstep, dt, well_mul) +end + + +function parse_forces(model, wells, controls, limits, cstep, dt, well_mul) + forces = [] + for (ctrl, lim) in zip(controls, limits) + ctrl_s = Dict{Symbol, Any}() + for (k, v) in pairs(ctrl) + ctrl_s[Symbol(k)] = v + end + lim_s = Dict{Symbol, Any}() + for (k, v) in pairs(lim) + lim_s[Symbol(k)] = v + end + f = setup_reservoir_forces(model, control = ctrl_s, limits = lim_s) + push!(forces, f) + end + return forces[cstep] +end + +function parse_state0(model, datafile) + rmodel = JutulDarcy.reservoir_model(model) + init = Dict{Symbol, Any}() + sol = datafile["SOLUTION"] + + if haskey(sol, "EQUIL") + init = parse_state0_equil(rmodel, datafile) + else + init = parse_state0_direct_assignment(rmodel, datafile) + end + state0 = setup_reservoir_state(model, init) +end + + +function parse_state0_direct_assignment(model, datafile) + sys = model.system + init = Dict{Symbol, Any}() + sol = datafile["SOLUTION"] + G = physical_representation(model.data_domain) + nc = number_of_cells(G) + ix = G.cell_map + is_blackoil = sys isa StandardBlackOilSystem + + function get_active(k; to_zero = false) + if haskey(sol, k) + x = zeros(nc) + for (i, c) in enumerate(G.cell_map) + x[i] = sol[k][c] + end + elseif to_zero + x = zeros(nc) + else + x = missing + end + return x + end + + nph = number_of_phases(sys) + sat = zeros(nph, nc) + if is_blackoil + rs = get_active("RS", to_zero = true) + rv = get_active("RV", to_zero = true) + end + + pressure = get_active("PRESSURE") + s_rem = ones(nc) + + sw = get_active("SWAT") + if !ismissing(sw) + s_rem -= sw + end + + sg = get_active("SGAS") + if !ismissing(sg) + s_rem -= sg + end + + @assert !ismissing(pressure) + init[:Pressure] = pressure + + if is_blackoil + if ismissing(sw) + sw = zeros(nc) + end + @assert !ismissing(sg) + if nph == 3 + init[:ImmiscibleSaturation] = sw + end + F_rs = sys.rs_max + F_rv = sys.rv_max + + so = s_rem + init[:BlackOilUnknown] = map( + (w, o, g, r, v, p) -> JutulDarcy.blackoil_unknown_init(F_rs, F_rv, w, o, g, r, v, p), + sw, so, sg, rs, rv, pressure) + else + error("Not implemented yet.") + end + return init +end + +function parse_reservoir(data_file) + grid = data_file["GRID"] + coord = grid["COORD"] + zcorn = grid["ZCORN"] + actnum = grid["ACTNUM"] + cartdims = grid["cartDims"] + primitives = JutulDarcy.cpgrid_primitives(coord, zcorn, cartdims, actnum = actnum) + G = JutulDarcy.grid_from_primitives(primitives) + + nc = number_of_cells(G) + ix = G.cell_map + # TODO: PERMYY etc for full tensor perm + perm = zeros(3, nc) + poro = zeros(nc) + for (i, c) in enumerate(G.cell_map) + perm[1, i] = grid["PERMX"][c] + perm[2, i] = grid["PERMY"][c] + perm[3, i] = grid["PERMZ"][c] + poro[i] = grid["PORO"][c] + end + satnum = JutulDarcy.InputParser.table_region(data_file, :saturation, active = G.cell_map) + eqlnum = JutulDarcy.InputParser.table_region(data_file, :equil, active = G.cell_map) + + domain = reservoir_domain(G, permeability = perm, porosity = poro, satnum = satnum, eqlnum = eqlnum) +end + + +function parse_physics_types(datafile) + runspec = datafile["RUNSPEC"] + props = datafile["PROPS"] + has(name) = haskey(runspec, name) && runspec[name] + has_wat = has("WATER") + has_oil = has("OIL") + has_gas = has("GAS") + has_disgas = has("DISGAS") + has_vapoil = has("VAPOIL") + + is_immiscible = !has_disgas && !has_vapoil + # TODO: compositional detection + # is_compositional = haskey(mrst_data, "mixture") + + phases = [] + rhoS = Vector{Float64}() + pvt = [] + if haskey(props, "DENSITY") + deck_density = props["DENSITY"] + if deck_density isa Matrix + deck_density = JutulDarcy.flat_region_expand(deck_density, 3) + end + if length(deck_density) > 1 + @warn "Multiple PVT regions found. Picking first one." deck_density + end + deck_density = deck_density[1] + rhoOS = deck_density[1] + rhoWS = deck_density[2] + rhoGS = deck_density[3] + else + @assert is_compositional + rhoOS = rhoWS = rhoGS = 1.0 + has_oil = true + has_gas = true + end + + if has_wat + push!(pvt, JutulDarcy.deck_pvt_water(props)) + push!(phases, AqueousPhase()) + push!(rhoS, rhoWS) + end + + if has_oil + push!(pvt, JutulDarcy.deck_pvt_oil(props)) + push!(phases, LiquidPhase()) + push!(rhoS, rhoOS) + end + + if has_gas + push!(pvt, JutulDarcy.deck_pvt_gas(props)) + push!(phases, VaporPhase()) + push!(rhoS, rhoGS) + end + if is_immiscible + sys = ImmiscibleSystem(phases, reference_densities = rhoS) + else + has_water = length(pvt) == 3 + oil_pvt = pvt[1 + has_water] + if oil_pvt isa JutulDarcy.PVTO + rs_max = JutulDarcy.saturated_table(oil_pvt) + else + rs_max = nothing + end + gas_pvt = pvt[2 + has_water] + if gas_pvt isa JutulDarcy.PVTG + rv_max = JutulDarcy.saturated_table(gas_pvt) + else + rv_max = nothing + end + sys = JutulDarcy.StandardBlackOilSystem(rs_max = rs_max, rv_max = rv_max, phases = phases, + reference_densities = rhoS) + end + + return (system = sys, pvt = pvt) +end + +function parse_schedule(domain, sys, datafile) + schedule = datafile["SCHEDULE"] + return parse_schedule(domain, datafile["RUNSPEC"], schedule, sys) +end + +function parse_control_steps(runspec, schedule, sys) + rho_s = JutulDarcy.reference_densities(sys) + phases = JutulDarcy.get_phases(sys) + + wells = schedule["WELSPECS"] + steps = schedule["STEPS"] + if length(keys(steps[end])) == 0 + # Prune empty final record + steps = steps[1:end-1] + end + + + tstep = Vector{Float64}() + cstep = Vector{Int}() + compdat = Dict{String, OrderedDict}() + controls = Dict{String, Any}() + limits = Dict{String, Any}() + for k in keys(wells) + compdat[k] = OrderedDict{NTuple{3, Int}, Any}() + controls[k] = nothing + limits[k] = nothing + end + all_compdat = [] + all_controls = [] + all_limits = [] + + if haskey(runspec, "START") + start_date = runspec["START"] + else + start_date = missing + end + current_time = 0.0 + function add_dt!(dt, ctrl_ix) + @assert dt > 0.0 + push!(tstep, dt) + push!(cstep, ctrl_ix) + end + + for (ctrl_ix, step) in enumerate(steps) + found_time = false + for (key, kword) in pairs(step) + if key == "DATES" + @assert !ismissing(start_date) + @assert !found_time + found_time = true + for date in kword + cdate = start_date + Second(current_time) + dt = Float64(Second(date - cdate).value) + add_dt!(dt, ctrl_ix) + + current_time += dt + end + elseif key == "TIME" + @assert !found_time + found_time = true + + for time in kword + dt = time - current_time + add_dt!(dt, ctrl_ix) + current_time = time + end + elseif key == "TSTEP" + @assert !found_time + found_time = true + for dt in kword + add_dt!(dt, ctrl_ix) + current_time = current_time + dt + end + elseif key == "COMPDAT" + for cd in kword + wname, I, J, K1, K2, flag, satnum, WI, diam, Kh, skin, Dfac, dir = cd + @assert haskey(wells, wname) + head = wells[wname].head + if I < 1 + I = head[1] + end + if J < 1 + J = head[2] + end + entry = compdat[wname] + for K in K1:K2 + entry[(I, J, K)] = ( + open = flag == "OPEN", + satnum = satnum, + WI = WI, + diameter = diam, + Kh = Kh, + skin = skin, + dir = dir, + ctrl = ctrl_ix) + end + end + elseif key in ("WCONINJE", "WCONPROD", "WCONHIST", "WCONINJ") + for wk in kword + name = wk[1] + controls[name], limits[name] = keyword_to_control(sys, wk, key) + end + elseif key in ("WEFAC", "WELTARG") + @warn "$key Not supported properly." + else + error("Unhandled keyword $key") + end + end + push!(all_compdat, deepcopy(compdat)) + push!(all_controls, deepcopy(controls)) + push!(all_limits, deepcopy(limits)) + if !found_time + error("Did not find supported time kw in step $ctrl_ix: Keys were $(keys(step))") + end + end + return (dt = tstep, control_step = cstep, controls = all_controls, completions = all_compdat, limits = all_limits) +end + + +function keyword_to_control(sys, kw, k::String) + return keyword_to_control(sys, kw, Val(Symbol(k))) +end + +function keyword_to_control(sys, kw, ::Val{:WCONPROD}) + rho_s = JutulDarcy.reference_densities(sys) + phases = JutulDarcy.get_phases(sys) + + flag = kw[2] + ctrl = kw[3] + orat = kw[4] + wrat = kw[5] + grat = kw[6] + lrat = kw[7] + bhp = kw[9] + return producer_control(sys, flag, ctrl, orat, wrat, grat, lrat, bhp) +end + +function keyword_to_control(sys, kw, ::Val{:WCONHIST}) + rho_s = JutulDarcy.reference_densities(sys) + phases = JutulDarcy.get_phases(sys) + # 1 name + error("Not implemented yet.") +end + +function producer_limits(; bhp = Inf, lrat = Inf, orat = Inf, wrat = Inf, grat = Inf) + lims = Dict{Symbol, Any}() + if isfinite(bhp) + lims[:bhp] = bhp + end + if isfinite(lrat) + lims[:lrat] = -lrat + end + if isfinite(orat) + lims[:orat] = -orat + end + if isfinite(wrat) + lims[:wrat] = -wrat + end + if isfinite(grat) + lims[:grat] = -grat + end + return NamedTuple(pairs(lims)) +end + +function producer_control(sys::Union{ImmiscibleSystem, StandardBlackOilSystem}, flag, ctrl, orat, wrat, grat, lrat, bhp) + rho_s = JutulDarcy.reference_densities(sys) + phases = JutulDarcy.get_phases(sys) + + if flag == "SHUT" || flag == "STOP" + ctrl = DisabledControl() + lims = nothing + else + @assert flag == "OPEN" + if ctrl == "LRAT" + t = SurfaceLiquidRateTarget(-lrat) + elseif ctrl == "WRAT" + t = SurfaceWaterRateTarget(-wrat) + elseif ctrl == "ORAT" + t = SurfaceOilRateTarget(-orat) + elseif ctrl == "GRAT" + t = SurfaceGasRateTarget(-grat) + elseif ctrl == "BHP" + t = BottomHolePressureTarget(bhp) + else + error("$ctype control not supported") + end + ctrl = ProducerControl(t) + lims = producer_limits(bhp = bhp, orat = orat, wrat = wrat, grat = grat, lrat = lrat) + end + return (ctrl, lims) +end + +function injector_limits(; bhp = Inf, surface_rate = Inf, reservoir_rate = Inf) + lims = Dict{Symbol, Any}() + if isfinite(bhp) + lims[:bhp] = bhp + end + if isfinite(surface_rate) + lims[:rate] = surface_rate + end + if !isinf(reservoir_rate) + @warn "Non-defaulted reservoir rate limit not supported: $reservoir_rate" + end + return NamedTuple(pairs(lims)) +end + +function injector_control(sys::Union{ImmiscibleSystem, StandardBlackOilSystem}, flag, type, ctype, surf_rate, res_rate, bhp) + rho_s = JutulDarcy.reference_densities(sys) + phases = JutulDarcy.get_phases(sys) + + if flag == "SHUT" || flag == "STOP" + ctrl = DisabledControl() + lims = nothing + else + @assert flag == "OPEN" + if ctype == "RATE" + t = TotalRateTarget(surf_rate) + elseif ctype == "BHP" + t = BottomHolePressureTarget(bhp) + else + # RESV, GRUP, THP + error("$ctype control not supported") + end + mix = Float64[] + rho = 0.0 + for (phase, rho_ph) in zip(phases, rho_s) + if phase == LiquidPhase() + v = Float64(type == "OIL") + elseif phase == AqueousPhase() + v = Float64(type == "WATER") + else + @assert phase isa VaporPhase + v = Float64(type == "GAS") + end + rho += rho_ph*v + push!(mix, v) + end + @assert sum(mix) ≈ 1.0 + ctrl = InjectorControl(t, mix, density = rho) + lims = injector_limits(bhp = bhp, surface_rate = surf_rate, reservoir_rate = res_rate) + end + return (ctrl, lims) +end + +function keyword_to_control(sys, kw, ::Val{:WCONINJE}) + type = kw[2] + flag = kw[3] + ctype = kw[4] + surf_rate = kw[5] + res_rate = kw[6] + bhp = kw[7] + return injector_control(sys, flag, type, ctype, surf_rate, res_rate, bhp) +end diff --git a/src/input_simulation/input_simulation.jl b/src/input_simulation/input_simulation.jl new file mode 100644 index 00000000..0631329e --- /dev/null +++ b/src/input_simulation/input_simulation.jl @@ -0,0 +1,2 @@ +include("mrst_input.jl") +include("data_input.jl") diff --git a/src/mrst_input.jl b/src/input_simulation/mrst_input.jl similarity index 97% rename from src/mrst_input.jl rename to src/input_simulation/mrst_input.jl index c646435a..3aa42cd7 100644 --- a/src/mrst_input.jl +++ b/src/input_simulation/mrst_input.jl @@ -52,11 +52,6 @@ function reservoir_domain_from_mrst(name::String; extraout = false) @debug "File read complete. Unpacking data..." g = MRSTWrapMesh(exported["G"]) has_trans = haskey(exported, "T") && length(exported["T"]) > 0 - if haskey(exported, "N") - N = Int64.(exported["N"]') - else - N = nothing - end function get_vec(d) if isa(d, AbstractArray) @@ -69,6 +64,19 @@ function reservoir_domain_from_mrst(name::String; extraout = false) perm = copy((exported["rock"]["perm"])') domain = reservoir_domain(g, porosity = poro, permeability = perm) + nf = number_of_faces(domain) + if haskey(exported, "N") + N = Int64.(exported["N"]') + nf = size(N, 2) + if nf != number_of_faces(domain) + d = Jutul.dim(g) + domain.entities[Faces()] = nf + domain[:areas, Faces()] = fill!(zeros(nf), NaN) + domain[:normals, Faces()] = fill!(zeros(d, nf), NaN) + domain[:face_centroids, Faces()] = fill!(zeros(d, nf), NaN) + end + domain[:neighbors, Faces()] = N + end # Deal with face data if has_trans @debug "Found precomputed transmissibilities, reusing" @@ -171,6 +179,13 @@ function get_well_from_mrst_data( L = vec(segs["length"]) D = vec(segs["diameter"]) rough = vec(segs["roughness"]) + n_segs = size(well_topo, 2) + if length(L) != n_segs + @warn "Inconsistent segments. Adding averages." + L = repeat([mean(L)], n_segs) + D = repeat([mean(D)], n_segs) + rough = repeat([mean(rough)], n_segs) + end @assert size(well_topo, 2) == length(L) == length(D) == length(rough) segment_models = map(SegmentWellBoreFrictionHB, L, rough, D) else @@ -345,6 +360,25 @@ function deck_relperm(props; oil, water, gas, satnum = nothing) return ReservoirRelativePermeability(; krarg..., regions = satnum, scaling = scaling) end +function flat_region_expand(x::AbstractMatrix, n = nothing) + # Utility to handle mismatch between MRST and Jutul parsers in simple PVT + # table format. + if !isnothing(n) + @assert size(x, 2) == n + end + x = map(i -> x[i, :], axes(x, 1)) + return x +end + +function flat_region_expand(x::Vector{Float64}, n = nothing) + return [x] +end + + +function flat_region_expand(x::Vector, n = nothing) + return x +end + function deck_pc(props; oil, water, gas, satnum = nothing) function get_pc(T, pc_ix, reverse = false) found = false @@ -428,7 +462,6 @@ function model_from_mat_deck(G, data_domain, mrst_data, res_context) is_immiscible = !has_disgas && !has_vapoil is_compositional = haskey(mrst_data, "mixture") - pvt = [] phases = [] rhoS = Vector{Float64}() if haskey(props, "DENSITY") @@ -447,7 +480,7 @@ function model_from_mat_deck(G, data_domain, mrst_data, res_context) has_oil = true has_gas = true end - + pvt = [] if is_compositional if has_wat push!(rhoS, rhoWS) @@ -617,11 +650,11 @@ end function set_deck_pvmult!(vars, param, props) # Rock compressibility (if present) if haskey(props, "ROCK") - rock = props["ROCK"] - if size(rock, 1) > 1 + rock = JutulDarcy.flat_region_expand(props["ROCK"]) + if length(rock) > 1 @warn "Rock has multiple regions, taking the first..." rock - rock = rock[1, :] end + rock = first(rock) if rock[2] > 0 static = param[:FluidVolume] delete!(param, :FluidVolume) diff --git a/src/multicomponent/flux.jl b/src/multicomponent/flux.jl index b6b28ae2..e8cb21aa 100644 --- a/src/multicomponent/flux.jl +++ b/src/multicomponent/flux.jl @@ -42,3 +42,15 @@ end end return q end + +function face_average_density(model::CompositionalModel, state, tpfa, phase) + ρ = state.PhaseMassDensities + s = state.Saturations + l = tpfa.left + r = tpfa.right + @inbounds s_l = s[phase, l] + @inbounds s_r = s[phase, r] + @inbounds ρ_l = ρ[phase, l] + @inbounds ρ_r = ρ[phase, r] + return (s_l*ρ_r + s_r*ρ_l)/max(s_l + s_r, 1e-8) +end diff --git a/src/multiphase.jl b/src/multiphase.jl index b26c5dfa..59934a3d 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -151,6 +151,10 @@ function Jutul.default_parameter_values(data_domain, model, param::Transmissibil U = data_domain[:permeability] g = physical_representation(data_domain) T = compute_face_trans(g, U) + if any(x -> x < 0, T) + c = count(x -> x < 0, T) + @warn "$c negative transmissibilities detected." + end else error(":permeability or :transmissibilities symbol must be present in DataDomain to initialize parameter $symb, had keys: $(keys(data_domain))") end diff --git a/src/porousmedia.jl b/src/porousmedia.jl index d34ad033..99874591 100644 --- a/src/porousmedia.jl +++ b/src/porousmedia.jl @@ -9,17 +9,20 @@ end function compute_peaceman_index(Δ, K, radius; dir::Symbol = :z, constant = 0.14, Kh = nothing, skin = 0, check = true) K_d = diag(K) - if dir == :x + if dir == :x || dir == :X L, d1, d2 = Δ i, j = 2, 3 - elseif dir == :y + elseif dir == :y || dir == :Y d1, L, d2 = Δ i, j = 1, 3 else d1, d2, L = Δ i, j = 1, 2 - @assert dir == :z "dir must be either :x, :y or :z (was :$dir)" + @assert dir == :z || dir == :Z "dir must be either :x, :y or :z (was :$dir)" end + @assert L > 0 + @assert d1 > 0 + @assert d2 > 0 k1, k2 = K_d[i], K_d[j] function kratio(l, v) @@ -39,7 +42,7 @@ function compute_peaceman_index(Δ, K, radius; dir::Symbol = :z, constant = 0.14 re = kratio(re1, re2) ke = sqrt(k1*k2) - if isnothing(Kh) + if isnothing(Kh) || isnan(Kh) Kh = L*ke end WI = 2 * π * Kh / (log(re / radius) + skin) diff --git a/src/utils.jl b/src/utils.jl index ca2c8f82..0b7082c6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -397,30 +397,18 @@ export setup_reservoir_state # Ex: For immiscible two-phase setup_reservoir_state(model, Pressure = 1e5, Saturations = [0.2, 0.8]) -Convenience constructor that initializes a state for a `MultiModel` set up using [`setup_reservoir_model`](@ref). -The main convenience over [`setup_state`](@ref) is only the reservoir initialization values need be provided: wells -are automatically initialized from the connected reservoir cells. +Convenience constructor that initializes a state for a `MultiModel` set up using +[`setup_reservoir_model`](@ref). The main convenience over [`setup_state`](@ref) +is only the reservoir initialization values need be provided: wells are +automatically initialized from the connected reservoir cells. + +As an alternative to passing keyword arguments, a `Dict{Symbol, Any}` instance +can be sent in as a second, non-keyword argument. """ -function setup_reservoir_state(model; kwarg...) +function setup_reservoir_state(model::MultiModel; kwarg...) rmodel = reservoir_model(model) - pvars = [k for k in keys(Jutul.get_primary_variables(rmodel))] - np = length(pvars) - ok = repeat([false], np) - res_init = Dict{Symbol, Any}() - for (k, v) in kwarg - I = findfirst(isequal(k), pvars) - if isnothing(I) - @warn "Recieved primary variable $k, but this is not known to reservoir model... Adding anyway." - else - ok[I] = true - end - res_init[k] = v - end - if !all(ok) - missing_primary_variables = pvars[.!ok] - @warn "Not all primary variables were initialized for reservoir model." missing_primary_variables - end - res_state = setup_state(rmodel, res_init) + pvars = collect(keys(Jutul.get_primary_variables(rmodel))) + res_state = setup_reservoir_state(rmodel; kwarg...) # Next, we initialize the wells. init = Dict(:Reservoir => res_state) perf_subset(v::AbstractVector, i) = v[i] @@ -459,6 +447,40 @@ function setup_reservoir_state(model; kwarg...) return state end +function setup_reservoir_state(model, init) + return setup_reservoir_state(model; pairs(init)...) +end + +function setup_reservoir_state(rmodel::SimulationModel; kwarg...) + pvars = collect(keys(Jutul.get_primary_variables(rmodel))) + svars = collect(keys(Jutul.get_secondary_variables(rmodel))) + np = length(pvars) + found = Symbol[] + res_init = Dict{Symbol, Any}() + for (k, v) in kwarg + I = findfirst(isequal(k), pvars) + if isnothing(I) + if !(k in svars) + @warn "Recieved primary variable $k, but this is not known to reservoir model... Adding anyway." + end + else + push!(found, k) + end + res_init[k] = v + end + handle_alternate_primary_variable_spec!(res_init, found, rmodel.system) + if length(found) != length(pvars) + missing_primary_variables = setdiff(pvars, found) + @warn "Not all primary variables were initialized for reservoir model." missing_primary_variables + end + return setup_state(rmodel, res_init) +end + +function handle_alternate_primary_variable_spec!(res_init, found, system) + # Internal utility to handle non-trivial specification of primary variables + return res_init +end + export setup_reservoir_forces """ setup_reservoir_forces(model; control = nothing, limits = nothing, set_default_limits = true, ) diff --git a/src/variables/pvt.jl b/src/variables/pvt.jl index 8cd7f5ae..81d263a3 100644 --- a/src/variables/pvt.jl +++ b/src/variables/pvt.jl @@ -123,7 +123,15 @@ function Jutul.default_parameter_values(data_domain, model, param::FluidVolume, return copy(vol) end -struct Temperature <: ScalarVariable end +Base.@kwdef struct Temperature{T} <: ScalarVariable + min::T = 0.0 + max::T = Inf + max_rel::Union{T, Nothing} = nothing + max_abs::Union{T, Nothing} = nothing +end -Jutul.default_value(model, ::Temperature) = 303.15 # 30.15 C° -Jutul.minimum_value(::Temperature) = 0.0 +Jutul.default_value(model, T::Temperature) = 303.15 # 30.15 C° +Jutul.minimum_value(T::Temperature) = T.min +Jutul.maximum_value(T::Temperature) = T.max +Jutul.absolute_increment_limit(T::Temperature) = T.max_abs +Jutul.relative_increment_limit(T::Temperature) = T.max_rel diff --git a/test/parser.jl b/test/parser.jl new file mode 100644 index 00000000..86ababab --- /dev/null +++ b/test/parser.jl @@ -0,0 +1,17 @@ +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 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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index afae40d3..8d07c5b1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,14 +28,18 @@ end include("sens_multimodel.jl") end -@testset "MRST input cases" begin +@testitem "MRST input cases" begin include("mrst_cases.jl") end -@testset "PArray solve" begin +@testitem "PArray solve" begin include("parray.jl") end +@testitem "Parser" begin + include("parser.jl") +end + @run_package_tests nothing # include("gpu.jl") diff --git a/test/utils.jl b/test/utils.jl index b12162b6..2a93c653 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -24,3 +24,18 @@ r = 0.2; end end + +import JutulDarcy: current_phase_index +@testset "current_phase_index" begin + depths = (1.0, 2.0) # W O G + @test current_phase_index(0.1, depths, reverse = false) == 1 + @test current_phase_index(1.01, depths, reverse = false) == 2 + @test current_phase_index(1.5, depths, reverse = false) == 2 + @test current_phase_index(2.5, depths, reverse = false) == 3 + + depths = (2.0, 1.0) # WO OG + @test current_phase_index(0.1, depths, reverse = true) == 3 # Gas + @test current_phase_index(1.01, depths, reverse = true) == 2 # Oil + @test current_phase_index(1.5, depths, reverse = true) == 2 # Oil + @test current_phase_index(2.5, depths, reverse = true) == 1 # Water +end