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.

[shakasaki]

Thank you,i tried out the pytorch backend and indeed the error disappears! For documentation purposes I include the working script below.

Any ideas on how to include several points per surface (what I describe above)? Not sure if you meant you'll look into it next week (this week), that's also fine - I really appreciate the fast responses you provide :-)

`import os
os.environ["DEFAULT_BACKEND"] = "PYTORCH"

import gempy as gp
import pandas as pd
import gempy_viewer as gpv

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)

gpv.plot_3d(gp_model)`

[NilsChudalla]

Hi @shakasaki,

I hope to come up with some more appropriate tips to avoid singular matrices this week, but I do not want to make any promises :D

Regarding several points per surface, the code does not change much. You could either add the points after creating the new elements and groups or enter the points upon element creation. The important thing is, that the shape of objects should be same, e.g. ["fault1", "fault1"] if you add two points. The latter approach I did here:

import gempy as gp
import numpy as np
import pandas as pd
import pyvista as pv
import os
#os.environ["DEFAULT_BACKEND"] = "PYTORCH"

# Dummy dataframe with structure points and orientations for 3 surfaces
structure_df = pd.DataFrame({
    'Easting (m)': [0, 20, 10, 20],
    'Northing (m)': [0, 10, 10, 20],
    'Elevation (m)': [0, 0, 0, 0],
    'Nx': [0, 0, 0, 0],
    'Ny': [1, 1, 1, 1],
    'Nz': [1, 1, 1, 1],
    'nugget': [100, 100, 100, 100],
    'name': ['fault1', '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 formation in structure_df["name"].unique():

    subset_mask = structure_df["name"] == formation
    subset_df = structure_df.loc[subset_mask]



    # this works but actually just overwrites the same structural element each time
    
    surf_name = subset_df["name"].unique()[0] # it is ugly but works for this purpose

    # 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=subset_df["Easting (m)"].tolist(),
            y=subset_df["Northing (m)"].tolist(),
            z=subset_df["Elevation (m)"].tolist(), 
            names=subset_df["name"].tolist()),
        orientations=gp.data.OrientationsTable.from_arrays(
            x=subset_df["Easting (m)"].tolist(),
            y=subset_df["Northing (m)"].tolist(),
            z=subset_df["Elevation (m)"].tolist(), 
            G_x=subset_df["Nx"].tolist(),
            G_y=subset_df["Ny"].tolist(),
            G_z=subset_df["Nz"].tolist(),
            names=subset_df["name"].tolist())
    )

    # new element has to 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"].unique().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)


import gempy_viewer as gpv
gpv.plot_3d(gp_model)

[shakasaki]

Thank you @NilsChudalla that's fantastic, works like a charm. Solving the numpy warning would be nice, but I think you responded to the original question posted here fully!