Skip to content

Commit

Permalink
Refactor and improve pinch processing
Browse files Browse the repository at this point in the history
  • Loading branch information
@moyner
moyner committed Aug 29, 2024
1 parent b7e16a4 commit b6ace74
Showing 1 changed file with 77 additions and 45 deletions.
122 changes: 77 additions & 45 deletions src/CornerPointGrid/processing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -452,8 +452,11 @@ function grid_from_primitives(primitives; nnc = missing, pinch = missing)
insert_interior_face!(prev_cell, cell, nodes, is_vertical, is_idir, face_type)
end
end
# Create pinch maps
pinch_top_to_bottom = generate_pinch_map(pinch, primitives, lines, column_lines, columns)

# Horizontal faces (top/bottom and faces along column)
pinch_added = 0
node_buffer = Int[]
sizehint!(node_buffer, 10)
for (cl, col) in zip(column_lines, columns)
Expand All @@ -477,8 +480,6 @@ function grid_from_primitives(primitives; nnc = missing, pinch = missing)
bottom_is_boundary = cell_is_boundary(next)
cell_bnds = map(l -> find_cell_bounds(cell, l), current_column_lines)
for is_top in (true, false)
# TODO: c1/c2 definition will have to be modified to pick the
# right normal vector, check this later.
if is_top
if !top_is_boundary
# Avoid adding interior faces twice.
Expand All @@ -496,6 +497,13 @@ function grid_from_primitives(primitives; nnc = missing, pinch = missing)
c2 = next
ft = :bottom
end
# Pinch treatment
if haskey(pinch_top_to_bottom, c1)
# If the there is a mapping (going down) from c1 to some
# other cell due to pinch we should not add anything.
@assert cell_is_boundary(c2)
continue
end
# Index into pillars
node_in_pillar_indices = map(F, cell_bnds)
# Then find the global node indices
Expand All @@ -504,7 +512,30 @@ function grid_from_primitives(primitives; nnc = missing, pinch = missing)
end
end
end
# primitives.column_boundary
# We skipped a bunch faces that belonged to pinched layers. Time to add them
# back in by systematically going through the pinch list.
num_pinched = length(keys(pinch_top_to_bottom))
if num_pinched > 0
pinch_count = 0
for (cl, col) in zip(column_lines, columns)
number_of_cells_in_column = length(col.cells)
current_column_lines = map(l -> lines[l], cl)
for top_cell in col.cells
cell_bnds = map(l -> find_cell_bounds(top_cell, l), current_column_lines)
bottom_cell = get(pinch_top_to_bottom, top_cell, nothing)
if isnothing(bottom_cell)
continue
end
node_in_pillar_indices = map(last, cell_bnds)
# Then find the global node indices
node_indices = map((line, i) -> line.nodes[i], current_column_lines, node_in_pillar_indices)
insert_face!(top_cell, bottom_cell, node_indices, is_vertical = false, is_boundary = false, is_idir = false, face_type = :bottom)
pinch_count += 1
end
end
@assert num_pinched == pinch_count
end
# Vertical faces
for is_bnd in [true, false]
if is_bnd
col_neighbors = primitives.column_boundary
Expand Down Expand Up @@ -542,48 +573,6 @@ function grid_from_primitives(primitives; nnc = missing, pinch = missing)
end
end

if !ismissing(pinch)
# Loop over columns, look for gaps
(; pinch, minpv_removed) = pinch
@assert length(pinch) == 5
thres = pinch[1]
num_added = 0
for (cl, col) in zip(column_lines, columns)
number_of_cells_in_column = length(col.cells)
current_column_lines = map(l -> lines[l], cl)

start = 1
while start < number_of_cells_in_column
before_inactive, last_inactive, done = find_next_gap(col.cells, start)
if done || before_inactive == last_inactive
break
end
top_cell = col.cells[before_inactive]
bottom_cell = col.cells[last_inactive + 1]
@assert top_cell > 0
@assert bottom_cell > 0

# Indices of face on the bottom of upper cell
top_pos = map(l -> last(find_cell_bounds(top_cell, l)), current_column_lines)
node_indices_top = map((line, i) -> line.nodes[i], current_column_lines, top_pos)

bottom_pos = map(l -> first(find_cell_bounds(bottom_cell, l)), current_column_lines)
node_indices_bottom = map((line, i) -> line.nodes[i], current_column_lines, top_pos)
z(i) = primitives.nodes[i][3]
z_face(ix) = (z(ix[1]) + z(ix[2]) + z(ix[3]) + z(ix[4]))/4.0
depth_top = z_face(node_indices_top)
depth_bottom = z_face(node_indices_bottom)
if depth_bottom - depth_top > thres
continue
end
# This is a bit hackish, but we just add a new actual face on top of the other boundary face that is already there.
insert_face!(top_cell, bottom_cell, node_indices_top, is_vertical = false, is_boundary = false, is_idir = false, face_type = :bottom)
start = last_inactive + 1
num_added += 1
end
end
end

if !ismissing(nnc)
to_active_ix = zeros(Int, nx*ny*nz)
to_active_ix[active] = eachindex(active)
Expand Down Expand Up @@ -684,6 +673,49 @@ function grid_from_primitives(primitives; nnc = missing, pinch = missing)
return g
end

function generate_pinch_map(pinch, primitives, lines, column_lines, columns)
pinch_top_to_bottom = Dict{Int, Int}()
if !ismissing(pinch)
# Loop over columns, look for gaps
(; pinch, minpv_removed) = pinch
@assert length(pinch) == 5
thres = pinch[1]
num_added = 0
for (cl, col) in zip(column_lines, columns)
number_of_cells_in_column = length(col.cells)
current_column_lines = map(l -> lines[l], cl)

start = 1
while start < number_of_cells_in_column
before_inactive, last_inactive, done = find_next_gap(col.cells, start)
if done || before_inactive == last_inactive
break
end
top_cell = col.cells[before_inactive]
bottom_cell = col.cells[last_inactive + 1]
@assert top_cell > 0
@assert bottom_cell > 0

# Indices of face on the bottom of upper cell
top_pos = map(l -> last(find_cell_bounds(top_cell, l)), current_column_lines)
node_indices_top = map((line, i) -> line.nodes[i], current_column_lines, top_pos)

bottom_pos = map(l -> first(find_cell_bounds(bottom_cell, l)), current_column_lines)
node_indices_bottom = map((line, i) -> line.nodes[i], current_column_lines, top_pos)
z(i) = primitives.nodes[i][3]
z_face(ix) = (z(ix[1]) + z(ix[2]) + z(ix[3]) + z(ix[4]))/4.0
depth_top = z_face(node_indices_top)
depth_bottom = z_face(node_indices_bottom)
start = last_inactive + 1
if depth_bottom - depth_top < thres
pinch_top_to_bottom[top_cell] = bottom_cell
num_added += 1
end
end
end
end
return pinch_top_to_bottom
end

function find_next_gap(cells, start)
# For cells, starting at "start", find the next interval of negative values.
Expand Down

0 comments on commit b6ace74

Please sign in to comment.