diff --git a/gempy/API/compute_API.py b/gempy/API/compute_API.py index 25c1fcdd..9d70ad88 100644 --- a/gempy/API/compute_API.py +++ b/gempy/API/compute_API.py @@ -90,7 +90,7 @@ def compute_model_at(gempy_model: GeoModel, at: np.ndarray, reset=True ) - sol = compute_model(gempy_model, engine_config, validate_serialization=True) + sol = compute_model(gempy_model, engine_config, validate_serialization=False) return sol.raw_arrays.custom diff --git a/test/test_private/test_terranigma/test_2025_2.py b/test/test_private/test_terranigma/test_2025_2.py new file mode 100644 index 00000000..7860a6ae --- /dev/null +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -0,0 +1,159 @@ +import os +import dotenv + +import gempy as gp +from gempy.API.io_API import read_surface_points +from gempy.core.data.surface_points import SurfacePointsTable +import gempy_viewer as gpv + +dotenv.load_dotenv() + + +def test_2025_2(): + # * Here I just leave as variable both, the manual that + # * you did and "the proper" that injects it to gempy + # * before interpolation without modifying anything else + manual_rescale = 1 + proper_rescale = 20. + + range_ = 2.4 + orientation_loc = -690 * manual_rescale + path_to_data = os.getenv("TEST_DATA") + + data = { + "a": read_surface_points(f"{path_to_data}/a.dat"), + "b": read_surface_points(f"{path_to_data}/b.dat"), + "c": read_surface_points(f"{path_to_data}/c.dat"), + "d": read_surface_points(f"{path_to_data}/d.dat"), + "e": read_surface_points(f"{path_to_data}/e.dat"), + "f": read_surface_points(f"{path_to_data}/f.dat"), + } + + # rescale the Z values + data = { + k: SurfacePointsTable.from_arrays( + x=v.data["X"], + y=v.data["Y"], + z=manual_rescale * v.data["Z"], # rescaling the z values + names=[k] * len(v.data), + nugget=v.data["nugget"] + ) + for k, v in data.items() + } + + color_generator = gp.data.ColorsGenerator() + elements = [] + for event, pts in data.items(): + orientations = gp.data.OrientationsTable.initialize_empty() + element = gp.data.StructuralElement( + name=event, + color=next(color_generator), + surface_points=pts, + orientations=orientations, + ) + elements.append(element) + + group = gp.data.StructuralGroup( + name="Series1", + elements=elements, + structural_relation=gp.data.StackRelationType.ERODE, + fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_FORMATIONS, + ) + structural_frame = gp.data.StructuralFrame( + structural_groups=[group], color_gen=color_generator + ) + + xmin = 525816 + xmax = 543233 + ymin = 5652470 + ymax = 5657860 + zmin = -780 * manual_rescale + zmax = -636 * manual_rescale + + # * Add 20% to extent + xmin -= 0.2 * (xmax - xmin) + xmax += 0.2 * (xmax - xmin) + ymin -= 0.2 * (ymax - ymin) + ymax += 0.2 * (ymax - ymin) + zmin -= 0.2 * (zmax - zmin) + zmax += 1 * (zmax - zmin) + + geo_model = gp.create_geomodel( + project_name="test", + extent=[xmin, xmax, ymin, ymax, zmin, zmax], + refinement=5, + structural_frame=structural_frame, + ) + + # * Here it is the way of rescaling one of the axis. Input transform + # * is used (by default) to rescale data into a unit cube but it accepts any transformation matrix. + geo_model.input_transform.scale[2] *= proper_rescale + + if True: + gpv.plot_3d( + model=geo_model, + ve=1, + image=True, + kwargs_pyvista_bounds={ + 'show_xlabels': False, + 'show_ylabels': False, + }, + transformed_data=True # * This is interesting, transformed data shows the data as it goes to the interpolation (after applying the transform) + ) + + + geo_model.interpolation_options.evaluation_options.number_octree_levels_surface = 4 + geo_model.interpolation_options.kernel_options.range = range_ + gp.modify_surface_points(geo_model, nugget=1e-5) + gp.add_orientations( + geo_model=geo_model, + x=[525825], + y=[5651315], + z=[orientation_loc], # * Moving the orientation further + pole_vector=[[0, 0, 1]], + elements_names=["a"] + ) + solution = gp.compute_model( + geo_model, + engine_config=gp.data.GemPyEngineConfig( + backend=gp.data.AvailableBackends.numpy + ), + ) + + gpv.plot_3d( + model=geo_model, + ve=1., + show_lith=True, + image=True, + kwargs_pyvista_bounds={ + 'show_xlabels': False, + 'show_ylabels': False, + 'show_zlabels': False, + }, + transformed_data=True + ) + + + # region Exporting scalar field + gpv.plot_2d( + geo_model, + show_scalar=True, + series_n=0 + ) + + # * The scalar fields can be found for dense and octree grids: + print(geo_model.solutions.raw_arrays.scalar_field_matrix) + + # * For custom grids so far we do not have a property that gives it directly, but it can be accessed here + + octree_lvl = 0 # * All the grids that are not octree are computed on octree level 0 + stack_number = -1 # * Here we choose the stack that we need. At the moment boolean operations--for erosion-- are not calculated on the scalar field + gempy_output = geo_model.solutions.octrees_output[octree_lvl].outputs_centers[stack_number] + slice_ = gempy_output.grid.custom_grid_slice + scalar_field = gempy_output.scalar_fields.exported_fields.scalar_field[slice_] + print(scalar_field) + # endregion + + +if __name__ == "__main__": + test_2025_2() diff --git a/test/test_private/test_terranigma/test_2025_3.py b/test/test_private/test_terranigma/test_2025_3.py new file mode 100644 index 00000000..69dbbfa6 --- /dev/null +++ b/test/test_private/test_terranigma/test_2025_3.py @@ -0,0 +1,122 @@ +import os +import dotenv + +import time +import numpy as np +import gempy as gp +from gempy.API.io_API import read_surface_points + +dotenv.load_dotenv() + +# Load data +def test_2025_3(): + + path_to_data = os.getenv("TEST_DATA") + # path_to_data = r"C:/Users/Benjamink/OneDrive - Mira Geoscience Limited/Documents/projects/implicit modelling/Nutrien/demo_terranigma/from_miguel" + + data = { + "a": read_surface_points(f"{path_to_data}/a.dat"), + "b": read_surface_points(f"{path_to_data}/b.dat"), + "c": read_surface_points(f"{path_to_data}/c.dat"), + "d": read_surface_points(f"{path_to_data}/d.dat"), + "e": read_surface_points(f"{path_to_data}/e.dat"), + "f": read_surface_points(f"{path_to_data}/f.dat"), + } + + # Build structural frame + + color_generator = gp.core.data.ColorsGenerator() + elements = [] + for event, pts in data.items(): + orientations = gp.data.OrientationsTable.initialize_empty() + element = gp.core.data.StructuralElement( + name=event, + color=next(color_generator), + surface_points=pts, + orientations=orientations, + ) + elements.append(element) + + group = gp.core.data.StructuralGroup( + name="Series1", + elements=elements, + structural_relation=gp.core.data.StackRelationType.ERODE, + fault_relations=gp.core.data.FaultsRelationSpecialCase.OFFSET_FORMATIONS, + ) + structural_frame = gp.core.data.StructuralFrame( + structural_groups=[group], color_gen=color_generator + ) + + # create cell centers with similar size to the BlockMesh used for Nutrien modelling + + xmin = 525816 + xmax = 543233 + ymin = 5652470 + ymax = 5657860 + zmin = -780 - 40 + zmax = -636 + 40 + + x = np.arange(xmin, xmax, 50) + y = np.arange(ymin, ymax, 50) + z = np.arange(zmin, zmax, 1) + X, Y, Z = np.meshgrid(x, y, z) + centers = np.c_[X.flatten(), Y.flatten(), Z.flatten()] + + # Create geomodel + + geo_model = gp.create_geomodel( + project_name="test", + extent=[xmin, xmax, ymin, ymax, zmin, zmax], + refinement=6, + structural_frame=structural_frame, + ) + gp.add_orientations( + geo_model=geo_model, + x=[525825], + y=[5651315], + z=[-686], + pole_vector=[[0, 0, 1]], + elements_names=["a"] + ) + + # Ignore surface creation in timing lithology block creation + geo_model.interpolation_options.evaluation_options.mesh_extraction=False + + # Time interpolation into octree cells + + tic = time.perf_counter() + solution = gp.compute_model( + geo_model, + engine_config=gp.data.GemPyEngineConfig( + backend=gp.data.AvailableBackends.PYTORCH, + use_gpu=True, + ), + ) + toc = time.perf_counter() + elapsed = toc - tic + print(f"Octree interpolation runtime: {int(elapsed / 60)} minutes {int(elapsed % 60)} seconds.") + + octrees_outputs = solution.octrees_output + n_cells = 0 + for octree_output in octrees_outputs: + n_cells += octree_output.outputs_centers[0].exported_fields.scalar_field.size().numel() + if len(octree_output.outputs_corners)>0: + n_cells += octree_output.outputs_corners[0].exported_fields.scalar_field.size().numel() + print(f"Number of cells evaluated: {n_cells}") + + # Time extra interpolation on regular grid centers. I was expecting/hoping that this second step + # would just be an evaluation of the continuous scalar field solution from first step. + + tic = time.perf_counter() + model = gp.compute_model_at( + geo_model, + at=centers, + engine_config=gp.data.GemPyEngineConfig( + backend=gp.data.AvailableBackends.PYTORCH, + use_gpu=True, + ), + ) + toc = time.perf_counter() + elapsed = toc - tic + print(f"Evaluate model on regular grid centers: {int(elapsed / 60)} minutes {int(elapsed % 60)} seconds") + print(f"Number of cells evaluated: {centers.shape[0]}") diff --git a/test/test_private/test_terranigma/test_2025_I.py b/test/test_private/test_terranigma/test_2025_I.py new file mode 100644 index 00000000..afb5d1cb --- /dev/null +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -0,0 +1,116 @@ +import os + +import dotenv + +import gempy as gp +from gempy.API.io_API import read_surface_points +import gempy_viewer as gpv + +dotenv.load_dotenv() + + +def test_2025_1(): + range_ = 0.6 + orientation_loc = -286 + + path_to_data = os.getenv("TEST_DATA") + + data = { + "a": read_surface_points(f"{path_to_data}/a.dat"), + "b": read_surface_points(f"{path_to_data}/b.dat"), + "c": read_surface_points(f"{path_to_data}/c.dat"), + "d": read_surface_points(f"{path_to_data}/d.dat"), + "e": read_surface_points(f"{path_to_data}/e.dat"), + "f": read_surface_points(f"{path_to_data}/f.dat"), + } + + color_generator = gp.data.ColorsGenerator() + elements = [] + for event, pts in data.items(): + orientations = gp.data.OrientationsTable.initialize_empty() + element = gp.data.StructuralElement( + name=event, + color=next(color_generator), + surface_points=pts, + orientations=orientations, + ) + elements.append(element) + + group = gp.data.StructuralGroup( + name="Series1", + elements=elements, + structural_relation=gp.data.StackRelationType.ERODE, + fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_FORMATIONS, + ) + structural_frame = gp.data.StructuralFrame( + structural_groups=[group], color_gen=color_generator + ) + + xmin = 525816 + xmax = 543233 + ymin = 5652470 + ymax = 5657860 + zmin = -780 + zmax = -636 + + # * Add 20% to extent + xmin -= 0.2 * (xmax - xmin) + xmax += 0.2 * (xmax - xmin) + ymin -= 0.2 * (ymax - ymin) + ymax += 0.2 * (ymax - ymin) + zmin -= 0.2 * (zmax - zmin) + zmax += 0.2 * (zmax - zmin) + + geo_model = gp.create_geomodel( + project_name="test", + extent=[xmin, xmax, ymin, ymax, zmin, zmax], + refinement=5, + structural_frame=structural_frame, + ) + + geo_model.interpolation_options.evaluation_options.number_octree_levels_surface = 4 + geo_model.interpolation_options.kernel_options.range = range_ + + + if False: + gpv.plot_3d( + model=geo_model, + ve=10, + image=True, + kwargs_pyvista_bounds={ + 'show_xlabels': False, + 'show_ylabels': False, + } + ) + + geo_model.interpolation_options.evaluation_options.number_octree_levels_surface = 4 + geo_model.interpolation_options.kernel_options.range = range_ + + gp.add_orientations( + geo_model=geo_model, + x=[525825], + y=[5651315], + z=[orientation_loc], # * Moving the orientation further + pole_vector=[[0, 0, 1]], + elements_names=["a"] + ) + solution = gp.compute_model( + geo_model, + engine_config=gp.data.GemPyEngineConfig( + backend=gp.data.AvailableBackends.numpy + ), + ) + + + + gpv.plot_3d( + model=geo_model, + ve=10, + show_lith=False, + image=True, + kwargs_pyvista_bounds={ + 'show_xlabels': False, + 'show_ylabels': False, + 'show_zlabels': False, + } + )