Multiple surface creation in a gempy model with no stratigraphic sequences

[shakasaki]

Hello gempy folks,

I am trying to use gempy's cokriging abilities to mesh a complex fault network within a homogeneous rock (granite). I want to treat the faults as surfaces in gempy to allow for curvature. I do not care about chronological relationships, the faults can cross each other in any direction.

I am able to do this only by assigning one surface (structural element?) at a time. I can initialize the model and then run
gp_model.structural_frame.structural_elements[0].name = 'fault1'
to create the fault, and then populate its points and orientations. However, I would like to add several of these "independent) faults in the model before hitting the "compute" button. Below you see a minimal working example of the code that works (one fault at a time) and a comparison of what I would like to do (commented out). I'm pretty sure there is an easy workaround, but I haven't managed to find it.

Thank you in advance for any help!

import gempy as gp
import numpy as np
import pandas as pd
import pyvista as pv

# Dummy dataframe with structure points and orientations for 3 surfaces
structure_df = pd.DataFrame({
    'Easting (m)': [0, 10, 20],
    'Northing (m)': [0, 10, 20],
    'Elevation (m)': [0, 0, 0],
    'Nx': [0, 0, 0],
    'Ny': [0, 0, 0],
    'Nz': [1, 1, 1],
    'nugget': [100, 100, 100],
    'name': ['fault1', 'fault2', 'fault3']
})

extent = [ -10, 30, -10, 30, -10, 10]

gp_model = gp.create_geomodel(project_name='test surface from df',
                                structural_frame=gp.data.StructuralFrame.initialize_default_structure(),
                                extent=extent, )

for index, row in structure_df.iterrows():
    nugget = row['nugget']
    points = structure_df[['Easting (m)', 'Northing (m)', 'Elevation (m)']].to_numpy().tolist()
    orientations = structure_df[['Nx', 'Ny', 'Nz']].to_numpy().tolist()

    # this works but actually just overwrites the same structural element each time
    gp_model.structural_frame.structural_elements[0].name = 'fault1'
    surf_name = 'fault1'

    # READ HERE
    # this is what i would like to do, but it fails
    # gp_model.structural_frame.structural_elements[index].name = row['name']
    # surf_name = row['name']

    gp.add_surface_points(
        geo_model=gp_model,
        x=[point[0] for point in points],
        y=[point[1] for point in points],
        z=[point[2] for point in points],
        elements_names=[surf_name] * len(points),
    )
    gp.add_orientations(
        geo_model=gp_model,
        x=[point[0] for point in points],
        y=[point[1] for point in points],
        z=[point[2] for point in points],
        elements_names=[surf_name]*len(points),
        pole_vector=orientations,
        nugget=[nugget]*len(points),
    )
gp_model.update_transform(gp.data.GlobalAnisotropy.NONE)
gp.compute_model(gp_model, compute_mesh=True)

[NilsChudalla]

Hi @shakasaki ,

you are very close, I tweeked your code for this result (changed the orientations to see the result):

# Dummy dataframe with structure points and orientations for 3 surfaces
structure_df = pd.DataFrame({
    'Easting (m)': [0, 10, 20],
    'Northing (m)': [0, 10, 20],
    'Elevation (m)': [0, 0, 0],
    'Nx': [0, 0, 0],
    'Ny': [1, 1, 1],
    'Nz': [1, 1, 1],
    'nugget': [100, 100, 100],
    'name': ['fault1', 'fault2', 'fault3']
})

extent = [ -10, 30, -10, 30, -10, 10]

gp_model = gp.create_geomodel(project_name='test surface from df',
                                structural_frame=gp.data.StructuralFrame.initialize_default_structure(),
                                extent=extent,
                                resolution=[20,20,20])

for index, row in structure_df.iterrows():
    nugget = row['nugget']
    points = structure_df[['Easting (m)', 'Northing (m)', 'Elevation (m)']].to_numpy().tolist()
    orientations = structure_df[['Nx', 'Ny', 'Nz']].to_numpy()

    # this works but actually just overwrites the same structural element each time
    
    surf_name = structure_df.loc[index,"name"]

    # you got it almost right, but you have to create a new structural element first
    new_element = gp.data.StructuralElement(
        name = surf_name,
        color=next(gp_model.structural_frame.color_generator),
        surface_points=gp.data.SurfacePointsTable.from_arrays(
            x=structure_df["Easting (m)"].array[[index]],
            y=structure_df["Northing (m)"].array[[index]],
            z=structure_df["Elevation (m)"].array[[index]], 
            names=[surf_name]),
        orientations=gp.data.OrientationsTable.from_arrays(
            x=structure_df["Easting (m)"].array[[index]],
            y=structure_df["Northing (m)"].array[[index]],
            z=structure_df["Elevation (m)"].array[[index]], 
            G_x=structure_df["Nx"].array[[index]],
            G_y=structure_df["Ny"].array[[index]],
            G_z=structure_df["Nz"].array[[index]],
            names=[surf_name])
    )

    # new element has to be put into a group
    group_fault = gp.data.StructuralGroup(
        name=surf_name,
        elements=[new_element],
        structural_relation=gp.data.StackRelationType.ERODE,
    )

    # new group has to be inserted into structural frame

    gp_model.structural_frame.insert_group(index=0,group=group_fault)

gp.set_is_fault(gp_model, structure_df["name"].tolist())

# remove unused elements and groups
gp.remove_element_by_name(gp_model,element_name="surface1")
gp.remove_structural_group_by_name(gp_model,group_name="default_formations")

gp_model.update_transform(gp.data.GlobalAnisotropy.NONE)
gp.compute_model(gp_model, compute_mesh=True)

[NilsChudalla]

I uploaded a similar project I had a few months ago, feel free to have a look. The aim was to build a DFN with Gempy. It utilizes the scalarfields by adding parallel faults to single structural group, which reduces the computational load quite a bit. Check it out here

[shakasaki]

Hi @NilsChudalla and thank you for the quick response. Indeed, your example works great! The only thing I tried to change, although unsuccessfully, is to add several points to a given surface before computing the model. So, for example each surface would have 5 points. Not sure why this is not converging, but I basically get a "LinAlgError: Singular matrix" error when simply changing the original dataframe to

structure_df = pd.DataFrame({ 'Easting (m)': [0, 5, 10, 20], 'Northing (m)': [0, 5, 10, 20], 'Elevation (m)': [0, 5, 0, 0], 'Nx': [0, 0, 0, 0], 'Ny': [0, 0, 1, 1], 'Nz': [1, 1, 1, 1], 'nugget': [100, 100, 100, 100], 'name': ['fault1', 'fault1', 'fault2', 'fault3'] })

This example would add two points to fault1, and since the nugget is rather high, i assume it is not the orientation vector that causes a singular matrix. However, if I change the first occurrence of fault1 to fault0, I get a model. What am I doing wrong here?

ps. while we're at it, is it possible to suppress this message:

RuntimeWarning: overflow encountered in exp return 1.0 [/](https://file+.vscode-resource.vscode-cdn.net/) (1.0 + bt.t.exp(x))

[NilsChudalla]

These Singular matrices are stealing my nerves :D I will have a look next week, if I can come up with a way to prevent this. I found Pytorch actually often avoiding these, so maybe try that until then. Or add more points. In the repo I shared, I played around with some "air layer" from a meshgrid. This solved the singular matrix error in some other model I had. But this answer is quite unsatisfying to me.