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

@@ -1,7 +1,7 @@
 name = "GeoEnergyIO"
 uuid = "3b1dd628-313a-45bb-9d8d-8f3c48dcb5d4"
 authors = ["Olav Møyner <olav.moyner@gmail.com> and contributors"]
-version = "1.1.7"
+version = "1.1.8"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

src/CornerPointGrid/CornerPointGrid.jl

@@ -4,15 +4,19 @@ 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
     include("interface.jl")
     include("processing.jl")
     include("processing_utils.jl")
     include("faults.jl")
+    include("nnc_and_aquifers.jl")
     include("utils.jl")
     include("pinch.jl")
 end

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)

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

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

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)