From 3ad37aafa0a3006e0b5ee1a7d3aa48b10f9dceaf Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 18 May 2025 08:42:39 +0100 Subject: [PATCH 01/14] [TEST] Added private test --- .../test_terranigma/test_2025_I.py | 73 +++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 test/test_private/test_terranigma/test_2025_I.py 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 000000000..2099e93fb --- /dev/null +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -0,0 +1,73 @@ +import os + +import dotenv + +import gempy as gp +from gempy.API.io_API import read_surface_points + + +dotenv.load_dotenv() + +def test_2025_1(): + + 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.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 + ) + + xmin = 525816 + xmax = 543233 + ymin = 5652470 + ymax = 5657860 + zmin = -780 + zmax = -636 + geo_model = gp.create_geomodel( + project_name="test", + extent=[xmin, xmax, ymin, ymax, zmin, zmax], + refinement=4, + 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"] + ) + solution = gp.compute_model( + geo_model, + engine_config=gp.data.GemPyEngineConfig( + backend=gp.data.AvailableBackends.PYTORCH, + use_gpu=True, + ), + ) + From 3c72e176cee1ff0fdb075b987511f9128c63cc91 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 18 May 2025 11:05:29 +0100 Subject: [PATCH 02/14] [TEST] Best configuration without big changes --- .../test_terranigma/test_2025_I.py | 66 ++++++++++++++----- 1 file changed, 51 insertions(+), 15 deletions(-) diff --git a/test/test_private/test_terranigma/test_2025_I.py b/test/test_private/test_terranigma/test_2025_I.py index 2099e93fb..4bd65bc21 100644 --- a/test/test_private/test_terranigma/test_2025_I.py +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -4,21 +4,24 @@ 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"), + "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"), + "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"), + "f": read_surface_points(f"{path_to_data}/f.dat"), } color_generator = gp.data.ColorsGenerator() @@ -36,10 +39,10 @@ def test_2025_1(): group = gp.data.StructuralGroup( name="Series1", elements=elements, - structural_relation=gp.core.data.StackRelationType.ERODE, - fault_relations=gp.core.data.FaultsRelationSpecialCase.OFFSET_FORMATIONS, + structural_relation=gp.data.StackRelationType.ERODE, + fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_FORMATIONS, ) - structural_frame = gp.core.data.StructuralFrame( + structural_frame = gp.data.StructuralFrame( structural_groups=[group], color_gen=color_generator ) @@ -49,25 +52,58 @@ def test_2025_1(): 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=4, + refinement=5, structural_frame=structural_frame, ) + + if False: + gpv.plot_3d( + model=geo_model, + ve=10, + 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=[-686], + 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.PYTORCH, - use_gpu=True, + backend=gp.data.AvailableBackends.PYTORCH ), ) - + + gpv.plot_3d( + model=geo_model, + ve=10, + show_lith=False, + image=False, + kwargs_pyvista_bounds={ + 'show_xlabels': False, + 'show_ylabels': False, + 'show_zlabels': False, + } + ) From e1d5e6d85d650e4fb88dac420f89c3fe87e41380 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Mon, 19 May 2025 12:26:09 +0100 Subject: [PATCH 03/14] [TEST] using numpy for ci --- test/test_private/test_terranigma/test_2025_I.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_private/test_terranigma/test_2025_I.py b/test/test_private/test_terranigma/test_2025_I.py index 4bd65bc21..2e2245bf7 100644 --- a/test/test_private/test_terranigma/test_2025_I.py +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -92,7 +92,7 @@ def test_2025_1(): solution = gp.compute_model( geo_model, engine_config=gp.data.GemPyEngineConfig( - backend=gp.data.AvailableBackends.PYTORCH + backend=gp.data.AvailableBackends.numpy ), ) From 824e8f8cb7c6f8183fd9998f727da8d793d5279d Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Mon, 19 May 2025 12:40:15 +0100 Subject: [PATCH 04/14] [TEST] Take 3d screenshots --- test/test_private/test_terranigma/test_2025_I.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_private/test_terranigma/test_2025_I.py b/test/test_private/test_terranigma/test_2025_I.py index 2e2245bf7..5ce55e025 100644 --- a/test/test_private/test_terranigma/test_2025_I.py +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -72,6 +72,7 @@ def test_2025_1(): gpv.plot_3d( model=geo_model, ve=10, + image=True, kwargs_pyvista_bounds={ 'show_xlabels': False, 'show_ylabels': False, @@ -100,7 +101,7 @@ def test_2025_1(): model=geo_model, ve=10, show_lith=False, - image=False, + image=True, kwargs_pyvista_bounds={ 'show_xlabels': False, 'show_ylabels': False, From b07c9b183b2be305f52e4359827db0ec86fafe7c Mon Sep 17 00:00:00 2001 From: benjamink Date: Tue, 27 May 2025 14:58:40 -0700 Subject: [PATCH 05/14] Add a test that shows improved contact point adherence for 20x scaled z data. --- .../test_terranigma/test_2025_2.py | 127 ++++++++++++++++++ 1 file changed, 127 insertions(+) create mode 100644 test/test_private/test_terranigma/test_2025_2.py 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 000000000..b6d4dc64f --- /dev/null +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -0,0 +1,127 @@ +import os +import pathlib +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(): + rescale = 20 + range_ = 2.4 + orientation_loc = -286 * rescale + # 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_cleaned.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=rescale * v.data["Z"], + 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 * rescale + zmax = -636 * 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, + ) + + if False: + gpv.plot_3d( + model=geo_model, + ve=40, + 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=2, + show_lith=True, + image=False, + kwargs_pyvista_bounds={ + 'show_xlabels': False, + 'show_ylabels': False, + 'show_zlabels': False, + }, + ) + +if __name__ == "__main__": + test_2025_2() From af71a18dd5358f25cddbf883b217582e7eed626c Mon Sep 17 00:00:00 2001 From: benjamink Date: Fri, 6 Jun 2025 15:00:23 -0700 Subject: [PATCH 06/14] Push new test demonstrating rescaling strategy and timing issue. --- .../test_terranigma/test_2025_2.py | 12 +- .../test_terranigma/test_2025_3.py | 113 ++++++++++++++++++ 2 files changed, 118 insertions(+), 7 deletions(-) create mode 100644 test/test_private/test_terranigma/test_2025_3.py diff --git a/test/test_private/test_terranigma/test_2025_2.py b/test/test_private/test_terranigma/test_2025_2.py index b6d4dc64f..8f557829d 100644 --- a/test/test_private/test_terranigma/test_2025_2.py +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -1,5 +1,4 @@ import os -import pathlib import dotenv import gempy as gp @@ -13,12 +12,11 @@ def test_2025_2(): rescale = 20 range_ = 2.4 - orientation_loc = -286 * rescale - # 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" + orientation_loc = -690 * rescale + path_to_data = os.getenv("TEST_DATA") data = { - "a": read_surface_points(f"{path_to_data}/a_cleaned.dat"), + "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"), @@ -31,7 +29,7 @@ def test_2025_2(): k: SurfacePointsTable.from_arrays( x=v.data["X"], y=v.data["Y"], - z=rescale * v.data["Z"], + z=rescale * v.data["Z"], # rescaling the z values names=[k] * len(v.data), nugget=v.data["nugget"] ) @@ -95,7 +93,7 @@ def test_2025_2(): 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], 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 000000000..457a729f1 --- /dev/null +++ b/test/test_private/test_terranigma/test_2025_3.py @@ -0,0 +1,113 @@ +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.") + + # 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") \ No newline at end of file From 4e5210d0a3165960c4135836a20b25d2b057e0bc Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 18 May 2025 08:42:39 +0100 Subject: [PATCH 07/14] [TEST] Added private test # Conflicts: # test/test_private/test_terranigma/test_2025_I.py --- test/test_private/test_terranigma/test_2025_I.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_private/test_terranigma/test_2025_I.py b/test/test_private/test_terranigma/test_2025_I.py index 5ce55e025..caae5f944 100644 --- a/test/test_private/test_terranigma/test_2025_I.py +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -6,6 +6,7 @@ from gempy.API.io_API import read_surface_points import gempy_viewer as gpv + dotenv.load_dotenv() From 27ee02489dfe679120411f4f524925bd77235564 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Sun, 18 May 2025 11:05:29 +0100 Subject: [PATCH 08/14] [TEST] Best configuration without big changes --- test/test_private/test_terranigma/test_2025_I.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/test/test_private/test_terranigma/test_2025_I.py b/test/test_private/test_terranigma/test_2025_I.py index caae5f944..afb5d1cb3 100644 --- a/test/test_private/test_terranigma/test_2025_I.py +++ b/test/test_private/test_terranigma/test_2025_I.py @@ -6,7 +6,6 @@ from gempy.API.io_API import read_surface_points import gempy_viewer as gpv - dotenv.load_dotenv() @@ -69,6 +68,10 @@ def test_2025_1(): 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, @@ -98,6 +101,8 @@ def test_2025_1(): ), ) + + gpv.plot_3d( model=geo_model, ve=10, From 5cd1d13541f9f19b7fb2a8292756d28f81f71276 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 16:22:32 +0200 Subject: [PATCH 09/14] [TEST] Test adapt test to include z anisotropy --- .../test_terranigma/test_2025_2.py | 31 +++++++++++++------ 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/test/test_private/test_terranigma/test_2025_2.py b/test/test_private/test_terranigma/test_2025_2.py index 8f557829d..17a3ec3b0 100644 --- a/test/test_private/test_terranigma/test_2025_2.py +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -10,9 +10,14 @@ def test_2025_2(): - rescale = 20 + # * 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 * rescale + orientation_loc = -690 * manual_rescale path_to_data = os.getenv("TEST_DATA") data = { @@ -29,7 +34,7 @@ def test_2025_2(): k: SurfacePointsTable.from_arrays( x=v.data["X"], y=v.data["Y"], - z=rescale * v.data["Z"], # rescaling the z values + z=manual_rescale * v.data["Z"], # rescaling the z values names=[k] * len(v.data), nugget=v.data["nugget"] ) @@ -62,8 +67,8 @@ def test_2025_2(): xmax = 543233 ymin = 5652470 ymax = 5657860 - zmin = -780 * rescale - zmax = -636 * rescale + zmin = -780 * manual_rescale + zmax = -636 * manual_rescale # * Add 20% to extent xmin -= 0.2 * (xmax - xmin) @@ -79,17 +84,23 @@ def test_2025_2(): 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 False: + if True: gpv.plot_3d( model=geo_model, - ve=40, - image=True, + ve=1, + image=False, 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_ @@ -111,7 +122,7 @@ def test_2025_2(): gpv.plot_3d( model=geo_model, - ve=2, + ve=proper_rescale, show_lith=True, image=False, kwargs_pyvista_bounds={ From 887e10c77a89f57dac3ec65a77244b9667838251 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Fri, 13 Jun 2025 10:03:10 +0200 Subject: [PATCH 10/14] [TEST] Add evaluation of octree cells in test Added logic to compute and print the number of cells evaluated during octree processing. Also included output for evaluated cell count on the regular grid centers. --- test/test_private/test_terranigma/test_2025_3.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/test/test_private/test_terranigma/test_2025_3.py b/test/test_private/test_terranigma/test_2025_3.py index 457a729f1..69dbbfa6f 100644 --- a/test/test_private/test_terranigma/test_2025_3.py +++ b/test/test_private/test_terranigma/test_2025_3.py @@ -12,7 +12,7 @@ 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" + # 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"), @@ -96,6 +96,14 @@ def test_2025_3(): 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. @@ -110,4 +118,5 @@ def test_2025_3(): ) toc = time.perf_counter() elapsed = toc - tic - print(f"Evaluate model on regular grid centers: {int(elapsed / 60)} minutes {int(elapsed % 60)} seconds") \ No newline at end of file + 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]}") From b39be6d27dd75b9323c7adeb173328aabf0693d6 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Fri, 13 Jun 2025 10:12:21 +0200 Subject: [PATCH 11/14] [WIP] Adding grabbing the scalar field --- .../test_terranigma/test_2025_2.py | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/test/test_private/test_terranigma/test_2025_2.py b/test/test_private/test_terranigma/test_2025_2.py index 17a3ec3b0..e44babd8f 100644 --- a/test/test_private/test_terranigma/test_2025_2.py +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -93,7 +93,7 @@ def test_2025_2(): gpv.plot_3d( model=geo_model, ve=1, - image=False, + image=True, kwargs_pyvista_bounds={ 'show_xlabels': False, 'show_ylabels': False, @@ -119,18 +119,32 @@ def test_2025_2(): backend=gp.data.AvailableBackends.numpy ), ) - + gpv.plot_3d( model=geo_model, ve=proper_rescale, show_lith=True, - image=False, + image=True, kwargs_pyvista_bounds={ 'show_xlabels': False, 'show_ylabels': False, 'show_zlabels': False, }, ) + + + # 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) + + # endregion + if __name__ == "__main__": test_2025_2() From 1908f2851e4f4d960ad6b861f302933d9d3017c9 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Fri, 13 Jun 2025 10:39:29 +0200 Subject: [PATCH 12/14] [TEST] Adding code to grab scalar field --- test/test_private/test_terranigma/test_2025_2.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_private/test_terranigma/test_2025_2.py b/test/test_private/test_terranigma/test_2025_2.py index e44babd8f..e9821e103 100644 --- a/test/test_private/test_terranigma/test_2025_2.py +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -143,6 +143,14 @@ def test_2025_2(): # * 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 From f7128cf9a6d533d31597e2520e654b7db353bd5d Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Fri, 13 Jun 2025 10:58:33 +0200 Subject: [PATCH 13/14] [TEST] Update 3D plot test with fixed ve and transformed data Modified the `ve` value to a fixed scale of 1. and enabled `transformed_data` in the 3D plot test for better coverage. --- test/test_private/test_terranigma/test_2025_2.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_private/test_terranigma/test_2025_2.py b/test/test_private/test_terranigma/test_2025_2.py index e9821e103..7860a6aec 100644 --- a/test/test_private/test_terranigma/test_2025_2.py +++ b/test/test_private/test_terranigma/test_2025_2.py @@ -122,7 +122,7 @@ def test_2025_2(): gpv.plot_3d( model=geo_model, - ve=proper_rescale, + ve=1., show_lith=True, image=True, kwargs_pyvista_bounds={ @@ -130,6 +130,7 @@ def test_2025_2(): 'show_ylabels': False, 'show_zlabels': False, }, + transformed_data=True ) From de6cc9bf1b5b30d0263733fec7f39ef64d0ff406 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Fri, 13 Jun 2025 11:25:50 +0200 Subject: [PATCH 14/14] [TEST] Not testing serialization for now when using compute_at --- gempy/API/compute_API.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gempy/API/compute_API.py b/gempy/API/compute_API.py index 25c1fcdd2..9d70ad88a 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