simpeg / discretize

Discretization tools for finite volume and inverse problems.
http://discretize.simpeg.xyz/
MIT License
163 stars 33 forks source link

Cell state not preserved with negative cell sizes #345

Open domfournier opened 6 months ago

domfournier commented 6 months ago

This minimal example demonstrates the issue. The cell state appears not preserved if the cell size is negative.

image

from discretize import TreeMesh
import numpy as np

point = [0, 0, 0]

# All positive cell sizes
mesh = TreeMesh([[10] * 16, [10] * 4, [10] * 8], [0, 0, 0])
mesh.insert_cells(point, [3], finalize=True)

# New mesh using state of first
mesh2 = TreeMesh([[10] * 16, [10] * 4, [10] * 8], [0, 0, 0])
mesh2.__setstate__((mesh.cell_state['indexes'], mesh.cell_state['levels']))

assert np.all(mesh.cell_centers == mesh2.cell_centers)  # Both meshes are equal

# Repeat with negative z cell size
mesh = TreeMesh([[10] * 16, [10] * 4, [-10] * 8], [0, 0, 0])
mesh.insert_cells(point, [3], finalize=True)

# New mesh using state of first
mesh2 = TreeMesh([[10] * 16, [10] * 4, [-10] * 8], [0, 0, 0])
mesh2.__setstate__((mesh.cell_state['indexes'], mesh.cell_state['levels']))

assert np.all(mesh.cell_centers == mesh2.cell_centers)  # Not equal
jcapriot commented 6 months ago

Interesting, there must be some missing logic wrapping around negative levels when not all cell width arrays have the same length.

jcapriot commented 6 months ago

Also… why are you using negative cell widths?

We should probably make the meshes error on this case due to undefined behavior.

domfournier commented 6 months ago

Came about while testing the conversion between the octrees stored in geoh5 from and to discretize. GA can deal with either positive or negative cell sizes , so it made sense to check.

I can just block the conversion if it happens for now - no big deal. Just thought I would flat it.