Summary
extract_submesh(g::UnstructuredMesh, cells) sets the submesh's cell_map to
cells, i.e. to the parent's local cell numbers. When the parent itself has
a non-identity cell_map (for example because it was already produced by
extract_submesh, or has inactive cells), the parent's map is not composed
in. The submesh keeps the parent's structure, so cell_ijk, cell_index
and anything else that relies on cell_map return wrong logical indices.
Minimal reproducer
using Jutul
g0 = UnstructuredMesh(CartesianMesh((3, 2, 2))) # 12 cells, cell_map = 1:12
g1 = extract_submesh(g0, [2, 3, 5, 6, 8, 9, 11, 12]) # drop the i = 1 cells
g2 = extract_submesh(g1, [3, 4, 5, 6, 7, 8]) # drop g1's first two cells
collect(g1.cell_map) # [2, 3, 5, 6, 8, 9, 11, 12] — correct
collect(g2.cell_map) # [3, 4, 5, 6, 7, 8] — expected [5, 6, 8, 9, 11, 12]
cell_ijk(g2, 1) # (3, 1, 1) — expected (2, 2, 1)
Through ordinary JutulDarcy use
The same happens without calling extract_submesh directly. A mesh with
inactive cells goes through extract_submesh once in reservoir_mesh
(JutulDarcy/src/reservoir_mesh.jl:107), and reservoir_domain(...; min_porevolume) calls it again (JutulDarcy/src/utils.jl:217):
using Jutul, JutulDarcy
actnum = ones(Int, 3, 2, 2); actnum[1, 1, 1] = 0 # one inactive cell
m = reservoir_mesh((3, 2, 2), (300.0, 200.0, 20.0); actnum = actnum)
poro = fill(0.2, number_of_cells(m)); poro[1] = 1e-6 # tiny pore volume
d = reservoir_domain(m; porosity = poro, min_porevolume = 1.0)
g = physical_representation(d)
collect(g.cell_map) # [2, 3, ..., 11] — expected [3, 4, ..., 12]
cell_ijk(g, 1) # (2, 1, 1) — expected (3, 1, 1)
cell_index(g, (3, 1, 1)) # 2, but local cell 2 is really (1, 2, 1)
Every remaining cell reports a logical index shifted by one. Anything that
locates cells by (I, J, K) is affected, e.g. setup_well given IJK tuples
(facility/wells/wells.jl:84), which would silently perforate the wrong
cells.
Cause
Jutul/src/meshes/unstructured/utils.jl, end of extract_submesh:
return UnstructuredMesh(
...;
structure = g.structure,
cell_map = cells, # parent-local numbers, not composed with g.cell_map
kwarg...
)
Suggested fix
Compose with the parent's map:
cell_map = isnothing(g.cell_map) ? cells : g.cell_map[cells],
Checked by passing that value explicitly as a keyword argument (which
extract_submesh forwards, overriding the default): correct cell_map and
cell_ijk in the first reproducer, and — applying the same composition to
the extract_submesh call that reservoir_domain makes — in the second.
Summary
extract_submesh(g::UnstructuredMesh, cells)sets the submesh'scell_maptocells, i.e. to the parent's local cell numbers. When the parent itself hasa non-identity
cell_map(for example because it was already produced byextract_submesh, or has inactive cells), the parent's map is not composedin. The submesh keeps the parent's
structure, socell_ijk,cell_indexand anything else that relies on
cell_mapreturn wrong logical indices.Minimal reproducer
Through ordinary JutulDarcy use
The same happens without calling
extract_submeshdirectly. A mesh withinactive cells goes through
extract_submeshonce inreservoir_mesh(
JutulDarcy/src/reservoir_mesh.jl:107), andreservoir_domain(...; min_porevolume)calls it again (JutulDarcy/src/utils.jl:217):Every remaining cell reports a logical index shifted by one. Anything that
locates cells by
(I, J, K)is affected, e.g.setup_wellgiven IJK tuples(
facility/wells/wells.jl:84), which would silently perforate the wrongcells.
Cause
Jutul/src/meshes/unstructured/utils.jl, end ofextract_submesh:Suggested fix
Compose with the parent's map:
Checked by passing that value explicitly as a keyword argument (which
extract_submeshforwards, overriding the default): correctcell_mapandcell_ijkin the first reproducer, and — applying the same composition tothe
extract_submeshcall thatreservoir_domainmakes — in the second.