Model of Ales, France: Plotting Sections and Maps

Explanation

This model is generally quite unstable and requires float64 precision to determine a solution. The lack of data in one corner for the TRIAS and LIAS series causes the model to bend unrealistically, eroding the CARBO layer, which disappears in that section.

import gempy as gp
import gempy_viewer as gpv
import os

Setting up paths

cwd = os.getcwd()
data_path = os.path.join(cwd, 'examples') if 'examples' not in cwd else os.path.join(cwd, '../..')

path_interf = os.path.join(data_path, "data/input_data/AlesModel/2018_interf.csv")
path_orient = os.path.join(data_path, "data/input_data/AlesModel/2018_orient_clust_n_init5_0.csv")
path_dem = os.path.join(data_path, "data/input_data/AlesModel/_cropped_DEM_coarse.tif")

Creating the geological model

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Claudius',
    extent=[729550.0, 751500.0, 1913500.0, 1923650.0, -1800.0, 800.0],
    resolution=None,
    refinement=6,
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=path_orient,
        path_to_surface_points=path_interf,
    )
)

Setting up the section grid

gp.set_section_grid(
    grid=geo_model.grid,
    section_dict={
            'section1': ([732000, 1916000], [745000, 1916000], [200, 150])
    }
)
Active grids: GridTypes.NONE|SECTIONS|OCTREE
start stop resolution dist
section1 [732000, 1916000] [745000, 1916000] [200, 150] 13000.0


Sorting lithologies

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={
            'fault_left': 'fault_left',
            'fault_right': 'fault_right',
            'fault_lr': 'fault_lr',
            'Trias_Series': ('TRIAS', 'LIAS'),
            'Carbon_Series': 'CARBO',
            'Basement_Series': 'basement'
    },
    remove_unused_series=True
)
Could not find element 'basement' in any group.
Structural Groups: StructuralGroup:
Name:fault_left
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:fault_left

StructuralGroup:
Name:fault_right
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:fault_right

StructuralGroup:
Name:fault_lr
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:fault_lr

StructuralGroup:
Name:Trias_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:TRIAS

StructuralElement:
Name:LIAS

StructuralGroup:
Name:Carbon_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:CARBO
Fault Relations:
fault_leftfault_righ...fault_lrTrias_Seri...Carbon_Ser...
fault_left
fault_right
fault_lr
Trias_Series
Carbon_Series
True
False


Changing colors

geo_model.structural_frame.get_element_by_name("LIAS").color = "#015482"
geo_model.structural_frame.get_element_by_name("TRIAS").color = "#9f0052"
geo_model.structural_frame.get_element_by_name("CARBO").color = "#ffbe00"

Plotting the 2D model

gpv.plot_2d(geo_model, direction='y')
Cell Number: mid Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7fbcb5ad3df0>

Plotting section traces

gpv.plot_section_traces(geo_model)
Cell Number: -1 Direction: z
<function plot_section_traces at 0x7fbd043bb010>

Setting faults

gp.set_is_fault(
    frame=geo_model.structural_frame,
    fault_groups=[
            geo_model.structural_frame.get_group_by_name('fault_left'),
            geo_model.structural_frame.get_group_by_name('fault_right'),
            geo_model.structural_frame.get_group_by_name('fault_lr')
    ],
    change_color=True
)
Structural Groups: StructuralGroup:
Name:fault_left
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_left

StructuralGroup:
Name:fault_right
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_right

StructuralGroup:
Name:fault_lr
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_lr

StructuralGroup:
Name:Trias_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:TRIAS

StructuralElement:
Name:LIAS

StructuralGroup:
Name:Carbon_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:CARBO
Fault Relations:
fault_leftfault_righ...fault_lrTrias_Seri...Carbon_Ser...
fault_left
fault_right
fault_lr
Trias_Series
Carbon_Series
True
False


Setting topography from file

gp.set_topography_from_file(
    grid=geo_model.grid,
    filepath=path_dem,
    crop_to_extent=[729550.0, 751500.0, 1913500.0, 1923650.0]
)

# Plotting 3D model with topography
gpv.plot_3d(geo_model, show_topography=True, ve=1, image=True)
AlesmodelAlesmodel
Active grids: GridTypes.NONE|SECTIONS|TOPOGRAPHY|OCTREE

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7fbcb6c18070>

Getting the Carbon Series

carbo = geo_model.structural_frame.get_group_by_name("Carbon_Series")

Modifying interpolation options for better model fitting

geo_model.interpolation_options.number_octree_levels_surface = 4
geo_model.interpolation_options.kernel_options.range = 0.8
gp.modify_surface_points(
    geo_model=geo_model,
    elements_names=["CARBO", "LIAS", "TRIAS"],
    nugget=0.005
)
Structural Groups: StructuralGroup:
Name:fault_left
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_left

StructuralGroup:
Name:fault_right
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_right

StructuralGroup:
Name:fault_lr
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_lr

StructuralGroup:
Name:Trias_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:TRIAS

StructuralElement:
Name:LIAS

StructuralGroup:
Name:Carbon_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:CARBO
Fault Relations:
fault_leftfault_righ...fault_lrTrias_Seri...Carbon_Ser...
fault_left
fault_right
fault_lr
Trias_Series
Carbon_Series
True
False


Displaying the structural frame

StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=fault_left,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault_left,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=fault_right,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault_right,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=fault_lr,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=fault_lr,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=Trias_Series,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=TRIAS,
        color=#9f0052,
        is_active=True
),
Element(
        name=LIAS,
        color=#015482,
        is_active=True
)
]
),
StructuralGroup(
        name=Carbon_Series,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=CARBO,
        color=#ffbe00,
        is_active=True
)
]
)
],
        fault_relations=
[[False, False, False,  True,  True],
 [False, False, False,  True,  True],
 [False, False, False,  True,  True],
 [False, False, False, False, False],
 [False, False, False, False, False]],
Structural Groups: StructuralGroup:
Name:fault_left
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_left

StructuralGroup:
Name:fault_right
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_right

StructuralGroup:
Name:fault_lr
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:fault_lr

StructuralGroup:
Name:Trias_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:TRIAS

StructuralElement:
Name:LIAS

StructuralGroup:
Name:Carbon_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:CARBO
Fault Relations:
fault_leftfault_righ...fault_lrTrias_Seri...Carbon_Ser...
fault_left
fault_right
fault_lr
Trias_Series
Carbon_Series
True
False


Explanation of model characteristics and adjustments This model has characteristics that make it difficult to get the right default values: - It is large, and we want high resolution - Some series have a large conditional number (i.e., the model input is not very stable) To address these issues: - Reduce the chunk size during evaluation to trade speed for memory - Reduce the std of the error parameter in octree refinement, which evaluates fewer voxels but may leave some without refinement Enable debugging options to help tune these parameters.

Setting verbose and condition number options for debugging

geo_model.interpolation_options.evaluation_options.verbose = True
geo_model.interpolation_options.kernel_options.compute_condition_number = True

Observations and parameter adjustments The octree refinement is making the octree grid almost dense, and smaller chunks are needed to avoid running out of memory. Adjusting parameters accordingly:

geo_model.interpolation_options.evaluation_options.octree_error_threshold = 0.5
geo_model.interpolation_options.evaluation_options.evaluation_chunk_size = 50_000

Computing the model with the adjusted settings

geo_model.interpolation_options.mesh_extraction = False
_ = gp.compute_model(
    geo_model,
    engine_config=gp.data.GemPyEngineConfig(
        backend=gp.data.AvailableBackends.PYTORCH,
        use_gpu=True,
        dtype="float64"
    )
)
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
  • Voxel Scalar Values with Refinement Status
Setting Backend To: AvailableBackends.PYTORCH
Condition number: 2917573.7407638184.
Chunking done: 19 chunks
Condition number: 1174481.1989654184.
Chunking done: 9 chunks
Condition number: 1642063.3117937546.
Chunking done: 11 chunks
Condition number: 1586673.9519253937.
Chunking done: 248 chunks
Condition number: 1122611.6721707566.
Chunking done: 198 chunks
Chunking done: 32 chunks
Chunking done: 26 chunks
Number of voxels marked by stats: 136 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 136
Dense Grid would be 512 voxels
Number of voxels marked by stats: 129 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 186
Dense Grid would be 512 voxels
Number of voxels marked by stats: 152 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 242
Dense Grid would be 512 voxels
Number of voxels marked by stats: 136 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 314
Dense Grid would be 512 voxels
Number of voxels marked by stats: 170 of torch.Size([512]).
 Number of voxels marked by corners : 328
Total voxels: 367
Dense Grid would be 512 voxels
Chunking done: 27 chunks
Chunking done: 22 chunks
Chunking done: 16 chunks
Chunking done: 8 chunks
Chunking done: 9 chunks
Chunking done: 205 chunks
Chunking done: 164 chunks
Number of voxels marked by stats: 1014 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 1014
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 976 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 1429
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 1135 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 1852
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 1051 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 2390
Dense Grid would be 4096 voxels
Number of voxels marked by stats: 1247 of torch.Size([3360]).
 Number of voxels marked by corners : 1656
Total voxels: 2741
Dense Grid would be 4096 voxels
Chunking done: 14 chunks
Chunking done: 7 chunks
Chunking done: 8 chunks
Chunking done: 184 chunks
Chunking done: 147 chunks
Chunking done: 108 chunks
Chunking done: 51 chunks
Chunking done: 62 chunks
Chunking done: 1459 chunks
Chunking done: 1167 chunks
Number of voxels marked by stats: 7902 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 7902
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 7629 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 11219
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 8784 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 14630
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 8070 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 18800
Dense Grid would be 32768 voxels
Number of voxels marked by stats: 9624 of torch.Size([24048]).
 Number of voxels marked by corners : 8203
Total voxels: 21430
Dense Grid would be 32768 voxels
Chunking done: 102 chunks
Chunking done: 48 chunks
Chunking done: 59 chunks
Chunking done: 1378 chunks
Chunking done: 1103 chunks

Plotting the 2D model with and without topography

gpv.plot_2d(geo_model, show_topography=False, section_names=['topography'], show_lith=True)
gpv.plot_2d(geo_model, cell_number=[4], direction=['y'], show_topography=True, show_data=True)
gpv.plot_2d(geo_model, cell_number=[-4], direction=['y'], show_topography=True, show_data=True)
  • Geological map
  • Cell Number: 4 Direction: y
  • Cell Number: -4 Direction: y
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_viewer/API/_plot_2d_sections_api.py:106: UserWarning: Section contacts not implemented yet. We need to pass scalar field for the sections grid
  warnings.warn(

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7fbc9ea0bf70>

Setting thumbnail number for Sphinx-Gallery sphinx_gallery_thumbnail_number = -1

gpv.plot_3d(geo_model, show_lith=True, show_topography=True, kwargs_plot_structured_grid={'opacity': 0.8})
Alesmodel
<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7fbcac3b7190>

Total running time of the script: (0 minutes 53.019 seconds)

Gallery generated by Sphinx-Gallery