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_4.py b/test/test_private/test_terranigma/test_2025_4.py
new file mode 100644
index 00000000..aecf6360
--- /dev/null
+++ b/test/test_private/test_terranigma/test_2025_4.py
@@ -0,0 +1,123 @@
+import dotenv
+import matplotlib.pyplot as plt
+
+import gempy as gp
+import numpy as np
+from gempy_engine.core.data.stack_relation_type import StackRelationType
+
+dotenv.load_dotenv()
+
+
+def stack_data(data, n_levels, min_level):
+    """
+    Dump cell centers/sizes and GemPy model for all octree levels.
+
+    :param data: gp.GeoModel containing GemPy Solution with interpolation data
+        on cell centers of an octree mesh
+    :param n_levels: Requested number of octree levels.
+    :param min_level: Requested min_level.
+    """
+
+    centers = []
+    model = []
+    h = []
+    for i in range(min_level, n_levels - min_level + 1):
+        output = data.solutions.octrees_output[i].outputs_centers[-1]
+        xyz = data.input_transform.apply_inverse(output.grid.octree_grid.values)
+        h = output.grid.octree_grid.dx / data.input_transform.scale[0]
+        centers.append(np.column_stack([xyz, np.repeat(h, len(xyz))]))
+        model.append(output.ids_block)
+
+    centers = np.vstack(centers)
+    model = np.hstack(model)
+
+    return centers, model
+
+
+def cells_inside_mask(centers):
+    """
+    Find cells that contain other cells.
+
+    :param centers: Cell centers (centers[:, :3]) and sizes (centers[:, -1]) of an
+        octree mesh
+    """
+
+    # Compute faces array where each row is a bounding box for a cell in the mesh.
+    h = centers[:, -1]
+    faces = np.column_stack([
+        centers[:, 0] - (h / 2),
+        centers[:, 0] + (h / 2),
+        centers[:, 1] - (h / 2),
+        centers[:, 1] + (h / 2),
+        centers[:, 2] - (h / 2),
+        centers[:, 2] + (h / 2),
+    ])
+
+    # Loop over bounding boxes and mask if any other cell centers
+    mask = np.ones(len(centers), dtype=bool)
+    for i, bbox in enumerate(faces):
+        contains_cells = (
+                (centers[:, 0] > bbox[0])
+                & (centers[:, 0] < bbox[1])
+                & (centers[:, 1] > bbox[2])
+                & (centers[:, 1] < bbox[3])
+                & (centers[:, 2] > bbox[4])
+                & (centers[:, 2] < bbox[5])
+        )
+
+        mask[i] = np.sum(contains_cells) == 1
+
+    return mask
+
+
+def test_2025_4():
+    data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'
+    path_to_data = data_path + "/data/input_data/jan_models/"
+    # Create a GeoModel instance
+
+    n_levels = 6
+    min_level = 2
+    extent = 1000
+
+    data = gp.create_geomodel(
+        project_name='fault',
+        extent=[0, extent, 0, extent, 0, extent],
+        refinement=n_levels,
+        importer_helper=gp.data.ImporterHelper(
+            path_to_orientations=path_to_data + "model5_orientations.csv",
+            path_to_surface_points=path_to_data + "model5_surface_points.csv"
+        )
+    )
+
+    # Map geological series to surfaces
+    gp.map_stack_to_surfaces(
+        gempy_model=data,
+        mapping_object={
+            "Fault_Series": 'fault',
+            "Strat_Series": ('rock2', 'rock1')
+        }
+    )
+
+    # Define fault groups
+    data.structural_frame.structural_groups[0].structural_relation = StackRelationType.FAULT
+    data.structural_frame.fault_relations = np.array([[0, 1], [0, 0]])
+
+    # Compute the geological model
+    data.interpolation_options.evaluation_options.octree_min_level = min_level
+    data.interpolation_options.evaluation_options.octree_error_threshold = 0.2
+    gp.compute_model(data)
+
+    cell_centers, model = stack_data(data, n_levels, min_level)
+    mask = cells_inside_mask(cell_centers)
+
+    filtered_centers = cell_centers[mask]
+    filtered_model = model[mask]
+
+    print(
+        "original n cells: ", len(cell_centers),
+        "filtered n cells: ", len(filtered_centers)
+    )
+
+
+if __name__ == "__main__":
+    test_2025_4()
\ No newline at end of file
diff --git a/test/test_private/test_terranigma/test_2025_5.py b/test/test_private/test_terranigma/test_2025_5.py
new file mode 100644
index 00000000..13562922
--- /dev/null
+++ b/test/test_private/test_terranigma/test_2025_5.py
@@ -0,0 +1,63 @@
+import dotenv
+import matplotlib.pyplot as plt
+
+import gempy as gp
+import numpy as np
+from gempy_engine.core.data.stack_relation_type import StackRelationType
+
+dotenv.load_dotenv()
+
+def test_2025_5():
+    data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'
+    path_to_data = data_path + "/data/input_data/jan_models/"
+    # Create a GeoModel instance
+
+    n_levels = 4
+    extent = 1000
+
+    levels = []
+    n_cells = []
+    for min_level in [0, 2]:
+        data = gp.create_geomodel(
+            project_name='fault',
+            extent=[0, extent, 0, extent, 0, extent],
+            refinement=n_levels,
+            importer_helper=gp.data.ImporterHelper(
+                path_to_orientations=path_to_data + "model5_orientations.csv",
+                path_to_surface_points=path_to_data + "model5_surface_points.csv"
+            )
+        )
+
+        # Map geological series to surfaces
+        gp.map_stack_to_surfaces(
+            gempy_model=data,
+            mapping_object={
+                "Fault_Series": 'fault',
+                "Strat_Series": ('rock2', 'rock1')
+            }
+        )
+
+        # Define fault groups
+        data.structural_frame.structural_groups[0].structural_relation = StackRelationType.FAULT
+        data.structural_frame.fault_relations = np.array([[0, 1], [0, 0]])
+
+        # Compute the geological model
+        data.interpolation_options.evaluation_options.octree_min_level = min_level
+        data.interpolation_options.evaluation_options.octree_error_threshold = 0.2
+        gp.compute_model(data)
+
+        levels.append([
+            k.outputs_centers[-1].grid.octree_grid.dx / data.input_transform.scale[0]
+            for k in data.solutions.octrees_output
+        ])
+        n_cells.append([
+            len(k.outputs_centers[-1].grid.octree_grid.values)
+            for k in data.solutions.octrees_output
+        ])
+
+    print(f"min level 0: {n_cells[0]} cells at {levels[0]}m")
+    print(f"min level 2: {n_cells[1]} cells at {levels[1]}m")
+
+
+if __name__ == "__main__":
+    test_2025_5()
\ No newline at end of file
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,
+        }
+    )