From 3cdfcac7fa4a5e2891709a285efabe359484cc79 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 18 Sep 2024 23:02:24 +0200 Subject: [PATCH] Support for numerical aquifers (#13) * Parse numerical aquifers * WIP on processing numerical aquifers * More aquifer setup work * Update nnc_and_aquifers.jl * Update nnc_and_aquifers.jl * Remove unused function * Improve matching logic for defaulted aquifer indices * Update Project.toml --- Project.toml | 2 +- src/CornerPointGrid/CornerPointGrid.jl | 4 + src/CornerPointGrid/faults.jl | 5 + src/CornerPointGrid/interface.jl | 45 +---- src/CornerPointGrid/nnc_and_aquifers.jl | 252 ++++++++++++++++++++++++ src/InputParser/InputParser.jl | 1 - src/InputParser/keywords/grid.jl | 72 ++++++- 7 files changed, 336 insertions(+), 45 deletions(-) create mode 100644 src/CornerPointGrid/nnc_and_aquifers.jl diff --git a/Project.toml b/Project.toml index 74b0d46..910ac27 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoEnergyIO" uuid = "3b1dd628-313a-45bb-9d8d-8f3c48dcb5d4" authors = ["Olav Møyner and contributors"] -version = "1.1.7" +version = "1.1.8" [deps] Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" diff --git a/src/CornerPointGrid/CornerPointGrid.jl b/src/CornerPointGrid/CornerPointGrid.jl index eb958bf..007eec6 100644 --- a/src/CornerPointGrid/CornerPointGrid.jl +++ b/src/CornerPointGrid/CornerPointGrid.jl @@ -4,8 +4,11 @@ module CornerPointGrid import Jutul: number_of_cells, number_of_faces, number_of_boundary_faces import Jutul: cell_ijk import Jutul: set_mesh_entity_tag! + import Jutul: mesh_entity_has_tag import Jutul: LinearInterpolant, BilinearInterpolant import Jutul: extract_submesh + import Jutul: cell_index + import Jutul: jutul_message import StaticArrays: SVector import LinearAlgebra: cross, dot, norm export mesh_from_grid_section @@ -13,6 +16,7 @@ module CornerPointGrid include("processing.jl") include("processing_utils.jl") include("faults.jl") + include("nnc_and_aquifers.jl") include("utils.jl") include("pinch.jl") end diff --git a/src/CornerPointGrid/faults.jl b/src/CornerPointGrid/faults.jl index 1df9098..e8ce868 100644 --- a/src/CornerPointGrid/faults.jl +++ b/src/CornerPointGrid/faults.jl @@ -44,6 +44,11 @@ function mesh_add_fault_tags!(G::UnstructuredMesh, faults) @debug "Fault $fault: Added $(length(fault_faces)) faces" set_mesh_entity_tag!(G, Faces(), :faults, Symbol(fault), fault_faces) end + return G +end + +function mesh_add_fault_tags!(G::UnstructuredMesh, faults::Missing) + return G end function get_sorted_face_pairs(N, IJK, ijk_ix) diff --git a/src/CornerPointGrid/interface.jl b/src/CornerPointGrid/interface.jl index 70fcf25..45eaf06 100644 --- a/src/CornerPointGrid/interface.jl +++ b/src/CornerPointGrid/interface.jl @@ -38,9 +38,9 @@ function mesh_from_grid_section(f, actnum = missing, repair_zcorn = true) else G = mesh_from_dxdydz_and_tops(grid, actnum = actnum) end - if haskey(grid, "FAULTS") - mesh_add_fault_tags!(G, grid["FAULTS"]) - end + # Handle faults + faults = get(grid, "FAULTS", missing) + mesh_add_fault_tags!(G, faults) return G end @@ -167,42 +167,3 @@ function cell_centers_from_deltas(dx, x0 = 0.0) end return (x, nx) end - -function insert_nnc_faces!(G::UnstructuredMesh, new_faces_neighbors, new_faces_nodes = fill(Int[], length(new_faces_neighbors))) - expand_indirection = x -> map(i -> copy(x[i]), 1:length(x)) - c2f = expand_indirection(G.faces.cells_to_faces) - f2n = expand_indirection(G.faces.faces_to_nodes) - - faceno = number_of_faces(G) - for neighbors in new_faces_neighbors - faceno += 1 - l, r = neighbors - push!(G.faces.neighbors, neighbors) - push!(c2f[l], faceno) - push!(c2f[r], faceno) - end - - for nodes in new_faces_nodes - push!(f2n, nodes) - end - replace_indirection!(G.faces.cells_to_faces, c2f) - replace_indirection!(G.faces.faces_to_nodes, f2n) - - @assert number_of_faces(G) == faceno - return G -end - -function replace_indirection!(x, expanded) - empty!(x.vals) - empty!(x.pos) - - push!(x.pos, 1) - for vals in expanded - n = length(vals) - for v in vals - push!(x.vals, v) - end - push!(x.pos, x.pos[end] + n) - end - return x -end diff --git a/src/CornerPointGrid/nnc_and_aquifers.jl b/src/CornerPointGrid/nnc_and_aquifers.jl new file mode 100644 index 0000000..1b5266e --- /dev/null +++ b/src/CornerPointGrid/nnc_and_aquifers.jl @@ -0,0 +1,252 @@ + +function insert_nnc_faces!(G::UnstructuredMesh, new_faces_neighbors, new_faces_nodes = fill(Int[], length(new_faces_neighbors))) + expand_indirection = x -> map(i -> copy(x[i]), 1:length(x)) + c2f = expand_indirection(G.faces.cells_to_faces) + f2n = expand_indirection(G.faces.faces_to_nodes) + + faceno = number_of_faces(G) + for neighbors in new_faces_neighbors + faceno += 1 + l, r = neighbors + push!(G.faces.neighbors, neighbors) + push!(c2f[l], faceno) + push!(c2f[r], faceno) + end + + for nodes in new_faces_nodes + push!(f2n, nodes) + end + replace_indirection!(G.faces.cells_to_faces, c2f) + replace_indirection!(G.faces.faces_to_nodes, f2n) + + @assert number_of_faces(G) == faceno + return G +end + +function replace_indirection!(x, expanded) + empty!(x.vals) + empty!(x.pos) + + push!(x.pos, 1) + for vals in expanded + n = length(vals) + for v in vals + push!(x.vals, v) + end + push!(x.pos, x.pos[end] + n) + end + return x +end + +function setup_numerical_aquifers!(data::AbstractDict) + aqunum = get(grid, "AQUNUM", missing) + aqucon = get(grid, "AQUCON", missing) + return mesh_add_numerical_aquifers!(G, aqunum, aqucon) +end + +function mesh_add_numerical_aquifers!(mesh, AQUNUM::Missing, AQUCON) + return nothing +end + +function mesh_add_numerical_aquifers!(mesh, AQUNUM, AQUCON) + if ismissing(AQUCON) + @warn "AQUNUM was defined by AQUCON was not found. Aquifer will not have any effect." + return nothing + end + # Remove repeats + AQUNUM = filter_aqunum(AQUNUM) + + dims = mesh.structure.I + actnum = fill(false, prod(dims)) + actnum[mesh.cell_map] .= true + # Create parameter list for aquifers + prm_T = @NamedTuple{ + cell::Int64, + area::Float64, + length::Float64, + porosity::Float64, + permeability::Float64, + depth::Float64, + pressure::Float64, + pvtnum::Int64, + satnum::Int64, + boundary_faces::Vector{Int}, + added_faces::Vector{Int}, + boundary_transmult::Vector{Float64}, + trans_option::Vector{Int} + } + aquifer_parameters = Dict{Int, prm_T}() + num_cells_start = number_of_cells(mesh) + new_cells_cell_map = Int[] + for aqunum in AQUNUM + id, I, J, K, A, L, phi, perm, D, p0, pvt, sat = aqunum + ix = ijk_to_linear(I, J, K, dims) + # "AQUNUM cannot declare aquifer in active cell ($I, $J, $K), cell must be inactive." + if actnum[ix] == false + # Add the cell + push!(new_cells_cell_map, ix) + cell = num_cells_start + length(new_cells_cell_map) + else + # Cell can be active and should then transition to being treated as + # aquifer. We find the matching cell. + cell = cell_index(mesh, (I, J, K)) + end + aquifer_parameters[id] = ( + cell = cell, + area = A, + length = L, + porosity = phi, + permeability = perm, + depth = D, + pressure = p0, + pvtnum = pvt, + satnum = sat, + boundary_faces = Int[], # Boundary faces that were connected + added_faces = Int[], # Corresponding fake faces that were aded + boundary_transmult = Float64[], # Trans mult of those fake faces + trans_option = Int[] # Option for how to compute that trans trans + ) + end + # Add all the faces + @assert length(keys(aquifer_parameters)) == length(AQUNUM) + IJK = map(i -> cell_ijk(mesh, i), 1:number_of_cells(mesh)) + nf0 = number_of_faces(mesh) + added_face_no = 0 + new_faces_neighbors = Tuple{Int, Int}[] + for (i, aqucon) in enumerate(AQUCON) + id, I_start, I_stop, J_start, J_stop, K_start, K_stop, dir, tranmult, opt, = aqucon + # Defaulted means start at 1 + if I_start < 1 + I_start = 1 + end + if J_start < 1 + J_start = 1 + end + if K_start < 1 + K_start = 1 + end + # Defaulted means no upper bound + if I_stop < 1 + I_stop = dims[1] + end + if J_stop < 1 + J_stop = dims[2] + end + if K_stop < 1 + K_stop = dims[3] + end + bfaces = find_faces_for_aquifer(mesh, I_start:I_stop, J_start:J_stop, K_start:K_stop, dir, IJK) + prm = aquifer_parameters[id] + + for bface in bfaces + added_face_no += 1 + push!(prm.boundary_faces, bface) + push!(prm.boundary_transmult, tranmult) + push!(prm.added_faces, added_face_no + nf0) + push!(prm.trans_option, opt) + # Add the new faces + c = mesh.boundary_faces.neighbors[bface] + push!(new_faces_neighbors, (c, prm.cell)) + end + end + fpos = mesh.faces.cells_to_faces.pos + bpos = mesh.boundary_faces.cells_to_faces.pos + for (i, c) in enumerate(new_cells_cell_map) + # Add cells without any face connections since these will be added + # afterwards as NNCs. + push!(fpos, fpos[end]) + push!(bpos, bpos[end]) + # Add the indices in the global enumeration of cells that have been + # added as active. + push!(mesh.cell_map, c) + end + insert_nnc_faces!(mesh, new_faces_neighbors) + @assert length(mesh.faces.neighbors) == nf0 + added_face_no + return aquifer_parameters +end + +function filter_aqunum(AQUNUM; warn = true) + indices = Dict{Int, Int}() + for (i, aqunum) in enumerate(AQUNUM) + id = first(aqunum) + if haskey(indices, id) && warn + # Can overwrite connections (assumed?) + jutul_message("AQUNUM", "Line $i: Replacing aquifer with id $id from line $(indices[id]) as new entry was provided.", color = :yellow) + end + indices[id] = i + end + # Do another loop to make sure this operation is stable and doesn't reorder + # the connections. Could use OrderedDict if we add OrderedCollections. + keep = Int[] + for (i, aqunum) in enumerate(AQUNUM) + id = first(aqunum) + if indices[id] == i + push!(keep, i) + end + end + return AQUNUM[keep] +end + +function find_faces_for_aquifer(mesh, I_range, J_range, K_range, dir, IJK) + bnd_faces = Int[] + if length(dir) == 1 || dir[2] == '+' + pos = true + else + @assert dir[2] == '-' + pos = false + end + + d = dir[1] + if d == 'X' || d == 'I' + ijk_orientation = :i + ijk_ix = 1 + if pos + dir_orientation = :right + else + dir_orientation = :left + end + elseif d == 'Y' || d == 'J' + ijk_orientation = :j + ijk_ix = 2 + if pos + dir_orientation = :upper + else + dir_orientation = :lower + end + elseif d == 'Z' || d == 'K' + ijk_orientation = :k + ijk_ix = 3 + if pos + # Positive direction is down in input + dir_orientation = :bottom + else + dir_orientation = :top + end + else + error("Bad direction for aquifer entry: $dir") + end + function boundary_face_check(faceno, e) + if !mesh_entity_has_tag(mesh, e, :ijk_orientation, ijk_orientation, faceno) + return false + end + if !mesh_entity_has_tag(mesh, e, :direction, dir_orientation, faceno) + return false + end + return true + end + + for cell in 1:number_of_cells(mesh) + I, J, K = IJK[cell] + if I in I_range && J in J_range && K in K_range + # TODO: We assume that aquifers are attached to the boundary here. + # Could maybe be generalized to handle interior faces too, if some + # tags were to be added. + for bfaceno in mesh.boundary_faces.cells_to_faces[cell] + if boundary_face_check(bfaceno, BoundaryFaces()) + push!(bnd_faces, bfaceno) + end + end + end + end + return bnd_faces +end diff --git a/src/InputParser/InputParser.jl b/src/InputParser/InputParser.jl index de04eec..892ad02 100644 --- a/src/InputParser/InputParser.jl +++ b/src/InputParser/InputParser.jl @@ -59,7 +59,6 @@ module InputParser skip_kw!(:PSPLITY, Inf, PARSER_MISSING_SUPPORT) skip_kw!(:PSPLITZ, Inf, PARSER_MISSING_SUPPORT) skip_kw!(:COMPLUMP, Inf, PARSER_MISSING_SUPPORT) - skip_kw!(:AQUNUM, Inf, PARSER_MISSING_SUPPORT) skip_kw!(:TRACER, Inf, PARSER_MISSING_SUPPORT) skip_kw!(:THPRES, Inf, PARSER_MISSING_SUPPORT) skip_kw!(:PIMULTAB, Inf, PARSER_MISSING_SUPPORT) diff --git a/src/InputParser/keywords/grid.jl b/src/InputParser/keywords/grid.jl index c7e1707..1eb5b20 100644 --- a/src/InputParser/keywords/grid.jl +++ b/src/InputParser/keywords/grid.jl @@ -335,7 +335,6 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:FAULTS}) end end - function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTFLT}) d = "*" tdims = [d, 1.0, 1.0] @@ -348,3 +347,74 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:MULTFLT}) data["MULTFLT"][name] = (mul_t, mul_d) end end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUNUM}) + def_and_u = [ + (1, :id), # Aquifer no + (-1, :id), # I + (-1, :id), # J + (-1, :id), # K + (NaN, :area), # surface area to aquifer + (NaN, :length), # length of aquifer + (NaN, :id), # poro of aquifer + (NaN, :permeability), # perm of aquifer + (NaN, :length), # depth of aquifer + (NaN, :pressure), # initial pressure aquifer + (1, :id), # PVTNUM for aquifer + (1, :id), # SATNUM for aquifer + ] + defaults = map(first, def_and_u) + utypes = map(last, def_and_u) + out = [] + while true + rec = read_record(f) + if length(rec) == 0 + break + end + l = parse_defaulted_line(rec, defaults) + swap_unit_system_axes!(l, units, utypes) + push!(out, l) + end + data["AQUNUM"] = out +end + +function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:AQUCON}) + def_and_u = [ + (1, :id), # Aquifer no + (-1, :id), # I start + (-1, :id), # I stop + (-1, :id), # J start + (-1, :id), # J stop + (-1, :id), # K start + (-1, :id), # K stop + ("Defaulted", :id), # Face orientation, I+, I-, ... + (1.0, :id), # Trans multiplier + (0, :id), # Type of trans calculator to use + ("NO", :id), # Allow internal aquifer connections + (1.0, :id), # Unsupported? + (1.0, :id) # Unsupported? + ] + defaults = map(first, def_and_u) + utypes = map(last, def_and_u) + out = [] + while true + rec = read_record(f) + if length(rec) == 0 + break + end + l = parse_defaulted_line(rec, defaults) + dir = l[8] + if length(dir) != 2 + ok = false + else + d1, d2 = dir + ok = d1 in ('I', 'J', 'K') && d2 in ('+', '-') + end + if !ok + throw(ArgumentError("Direction for AQUCON was $dir, must be on the format I/J/K and +-, i.e. I+")) + end + swap_unit_system_axes!(l, units, utypes) + push!(out, l) + end + data["AQUCON"] = out +end