Export a geological model from GemPy to use in MOOSE

import gempy as gp

Creating a geological model

The procedure of generating a geological model is presented in detail in Chapter 1-1 of the GemPy tutorials, so it will only be briefly presented here

Initiate a model

geo_model = gp.create_model('tutorial_moose_exp')
data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'

Import data from CSV-files with setting the resolution and model extent

gp.init_data(geo_model, [0, 2000., 0, 2000., 0, 2000.], [50, 50, 80],
             path_o=data_path + "/data/input_data/tut_chapter1/simple_fault_model_orientations.csv",
             path_i=data_path + "/data/input_data/tut_chapter1/simple_fault_model_points.csv",
             default_values=True)

Out:

Active grids: ['regular']

tutorial_moose_exp  2020-09-15 12:24

present the units and series

surface series order_surfaces color id
0 Shale Default series 1 #015482 1
1 Sandstone_1 Default series 2 #9f0052 2
2 Siltstone Default series 3 #ffbe00 3
3 Sandstone_2 Default series 4 #728f02 4
4 Main_Fault Default series 5 #443988 5
5 basement Basement 1 #ff3f20 6


combine units in series and make two series, as the fault needs its own

gp.map_stack_to_surfaces(geo_model,
                         {"Fault_Series": 'Main_Fault',
                          "Strat_Series": ('Sandstone_2', 'Siltstone', 'Shale', 'Sandstone_1', 'basement')},
                         remove_unused_series=True)

# set the fault series to be fault object
geo_model.set_is_fault(['Fault_Series'], change_color=False)
order_series BottomRelation isActive isFault isFinite
Fault_Series 1 Fault True True False
Strat_Series 2 Erosion True False False


check whether series were assigned correctly

surface series order_surfaces color id
4 Main_Fault Fault_Series 1 #443988 1
0 Shale Strat_Series 1 #015482 2
1 Sandstone_1 Strat_Series 2 #9f0052 3
2 Siltstone Strat_Series 3 #ffbe00 4
3 Sandstone_2 Strat_Series 4 #728f02 5
5 basement Strat_Series 5 #ff3f20 6


Model generation

After loading in the data, we set it up for interpolation and compute the model.

set up interpolator

gp.set_interpolator(geo_model,
                    compile_theano=True,
                    theano_optimizer='fast_compile',
                    verbose=[])

Out:

Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  1
Compilation Done!
Kriging values:
                   values
range            3.5e+03
$C_o$            2.9e+05
drift equations   [3, 3]

<gempy.core.interpolator.InterpolatorModel object at 0x7f859fd18690>

compute the model

gp.compute_model(geo_model, compute_mesh=False);

Out:

Lithology ids
  [6. 6. 6. ... 2. 2. 2.]

have a look at the data and computed model

gempy export MOOSE

Out:

<gempy.plot.vista.GemPyToVista object at 0x7f8564e97d10>

Exporting the Model to MOOSE

The voxel-model above already is the same as a model discretized in a hexahedral grid, so my immediately be used as input in a simulation tool, e.g. MOOSE. For this, we need to access to the unit IDs assigned to each voxel in GemPy. The array containing these IDs is called lith_block.

Out:

[6. 6. 6. ... 2. 2. 2.]

This array has the shape of (x,) and would be immediately useful, if GemPy and the chosen simulation code would populate a grid in the same way. Of course, however, that is not the case. This is why we have to restructure the lith_block array, so it can be read correctly by MOOSE.

model resolution

nx, ny, nz = geo_model.grid.regular_grid.resolution

# model extent
xmin, xmax, ymin, ymax, zmin, zmax = geo_model.grid.regular_grid.extent

These two parameters are important to, a) restructure lith_block, and b) write the input file for MOOSE correctly. For a), we need to reshape lith_block again to its three dimensions and re-flatten it in a MOOSE-conform way.

reshape to 3D array

units = ids.reshape((nx, ny, nz))
# flatten MOOSE conform
units = units.flatten('F')
The importance of nx, ny, nz is apparent from the cell above. But what about xmin, …, zmax?
A MOOSE input-file for mesh generation has the following syntax:
[MeshGenerators]
  [./gmg]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 50
    ny = 50
    nz = 80
    xmin = 0.0
    xmax = 2000.0
    yim = 0.0
    ymax = 2000.0
    zmin = 0.0
    zmax = 2000.0
    block_id = '1 2 3 4 5 6'
    block_name = 'Main_Fault Sandstone_2 Siltstone Shale Sandstone_1 basement'
  [../]

  [./subdomains]
    type = ElementSubdomainIDGenerator
    input = gmg
    subdomain_ids = ' ' # here you paste the transformed lith_block vector
  [../]
[]

[Mesh]
  type = MeshGeneratorMesh
[]

So these parameters are required inputs in the [MeshGenerators] object in the MOOSE input file. GemPy has a method to directly create such an input file, stored in gempy.utils.export.py.

The following cell shows how to call the method:

# sphinx_gallery_thumbnail_path = '_static/GemPy_model_combined.png'
import gempy.utils.export as export
export.export_moose_input(geo_model, path='')

Out:

Successfully exported geological model as moose input to ./

This method automatically stores a file geo_model_units_moose_input.i at the specified path. Either this input file could be extended with parameters to directly run a simulation, or it is used just for creating a mesh. In the latter case, the next step would be, to run the compiled MOOSE executable witch the optional flag --mesh-only.

E.g. with using the PorousFlow module:

$path_to_moose/moose/modules/porous_flow/porous_flow-opt -i pct_voxel_mesh.i --mesh-only

How to compile MOOSE is described in their documentation.

The now generated mesh with the name geo_model_units_moose_input_in.e can be used as input for another MOOSE input file, which contains the main simulation parameters. To call the file with the grid, the following part has to be added in the MOOSE simulation input file:

[Mesh]
  file = geo_model_units_moose_input_in.e
[]

The final output of the simulation may also be such an .e, which can, for instance, be opened with paraview. A simulated temperature field (purely conductive) of the created model would look like this:

gempy_temperature

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

Gallery generated by Sphinx-Gallery