Tutorials¶
The four tutorial files contain the most important commands and parameter options currently available in custEM. Users can comment or uncomment all arguments and keyword arguments to explore the usage of custEM.
NOTICE¶
As custEM requies MPI for all realistic simulations, the jupyter notebooks are currently not maintained because embedding MPI runs in jupyter notebooks does not work robustly. Please refer to the python files instead.
Only the meshgen tutorials work as jupyter notebook. For running simulations, please use the provided Python script equivalent to the jupyter notebook. Run simulations such as the run_tutorial.py python script by calling it from the command line with the mpirun syntax:
–> mpirun -n 4 python run_tutorial.py
Here, 4 parallel MPI processes are used.
Without modification, the provided “run” example requires about 7 GB of RAM and can be conducted in serial or parallel (please remember that even in serial, the “mpirun -n 1 python run_tutorial.py” syntax is required.
The html docs are also accessible from the main documentation of custEM The run_tutorial is independed from the meshgen_tutorial, but the plot_tutorial depends on the run_tutorial. Therefore, for using the plot_tutorial, it is necessary to run the FE computations before.
Mesh generation¶
Meshgen tutorial¶
This tutorial is a guide to create tetrahedral meshes for calculations whit custEM. All options are listed in the comment sections and can be uncomment for the required purposes. Please note that the documentation of custEM is still in development as is this tutorial. More details about incorporating digital-elevation models (DEM) will be added soon.
Import modules, utility functions and synthetic definitions¶
In [1]:
from custEM.meshgen import meshgen_utils as mu
from custEM.meshgen.meshgen_tools import BlankWorld
from custEM.misc.synthetic_definitions import example_3_line_tx
from custEM.misc.synthetic_definitions import surface_anomly_outcrop_1
# from custEM.misc.synthetic_definitions import topo_func1
Incorporating topography¶
- it is referred to the topography meshing tutorial
Create computational domain¶
In [ ]:
# The following keyword arguments can be uncommented as needed.
# #############################################################################
M = BlankWorld( # blank 3D domain
name='my_mesh', # name of the output mesh
# m_dir='./meshes', # main mesh directory
# s_dir # save dir., default: **m_dir**+ '/_mesh'
# w_dir # working dir., default:
# # # # # # # # # # # # # # # # # **m_dir**+ '/mesh_create_tmp'
# p_dir # para dir., default: **m_dir**+ '/para'
# t_dir # topo dir., **p_dir**+ '/topo'
# x_dim=[-1e4, 1e4], # x-dimension [m] of the mesh
# y_dim=[-1e4, 1e4], # y-dimension [m] of the mesh
# y_dim=[-1e4, 1e4], # z-dimension [m] of the mesh
# backup_script = True, # store mesh generation script as backup
#
# keyword arguments regarding building the 2D surface mesh.
#
# inner_area=None, # refinement of area at surface
# inner_area_size=[], # size of this area
# inner_area_shift=[0., 0.], # shift in x- and y-dir of area
# inner_area_cell_size=1e4, # max. triangle-size (MA) in area
# outer_area_cell_size=1e6, # MA outside of this area
# triangle_q=34., # quality of all 2D *triangle* meshes
#
# keyword arguments regarding building the 3D subsurface and air.
#
# interface_cell_sizes=None, # MA for subsurface interfaces
# airspace_cell_size=0., # max. tetrahedral volume (MV) for airspace
# subsurface_cell_size=0., # MV for halfspace subsurface
# layer_cell_sizes=None, # list of MV for layered subsurface
# tol=1e-2 # tolerance value for adding marker
#
# keyword arguments regarding incorporation of topography
#
# topo=None, # topography (from DEM or synthetic),
# # # # # # # # # # # # # # # # # e.g. topo=topo_func1
# subsurface_topos=None, # subsurface interface topo (DEM, synth.)
# easting_shift=None, # shift DEM east. relative to model domain
# northing_shift=None, # shift DEM north. relative to model domain
# centering=False, # center DEM to model domain
# rotation=None, # rotate DEM realative to model domain
#
# keyword arguments to preserve structures in the mesh
#
preserve_edges=True # important to preserve Tx paths
# preserve_nodes # preserve single points added to the mesh
)
Create surface mesh¶
In [3]:
# Build the 2D surface mesh and incorporate transmitter or observation lines.
# #############################################################################
# There are two different ways of adding and refining receiver locations
# option 1, recommended in general, see usage in topography meshing tutorial
# # the following function builds triangles around each Rx positions, so they
# # are located on face-centers of tetrahedra afterwards, leading to a better
# # interpolation accuracy in general
# rx = mu.line_x(-5e3, 5e3, 100)
# rx_tri = mu.refine_rx(rx, 1., 30)
# M.build_surface(
# insert_paths=rx_tri, # refine observation points
# insert_line_tx=[example_3_line_tx()] # add crooked dipole Tx
# )
# # the following command adds the Rx coordinates (note: above, we have
# # specified the refined triangles around the Rx locations)
# # to the parameters file for automated detection during the simulation
# M.add_rx(rx)
# option 2, only recommended for dense spacing on observation lines/slices
rx = mu.line_x(-5e3, 5e3, 1000)
M.build_surface(
insert_lines=[rx], # add observation line
insert_line_tx=[example_3_line_tx()] # add crooked dipole Tx
# insert_loop_tx=[mu.loop_r(...)], # add circ. or rect. loops
# insert_points=[mu.pointset(...)], # add points to the surface
)
In [4]:
# Define the intersecting area of anomalies reaching the surface.
# Here, we use a polygone as defined in *surface_anomly_outcrop_1*
# The outcrop is extended to depth later on.
# #############################################################################
M.add_surface_anomaly(insert_paths=[surface_anomly_outcrop_1()],
dips=[60.], # dip of anomaly
depths=[-1000.], # absolute depth of anomaly
dip_azimuths=[110.], # dip azimuth
cell_sizes=[1e5]) # max. cell volume in anomaly
Build 3D world based on 2D surface mesh¶
- Build either a fullspace, halfspace or layered-earth mesh
In [5]:
# The most simple world is a fullspace, in this case the **build_surface**
# method does not need to be executed or is ignored, respectively.
# #############################################################################
# M.build_fullspace_mesh()
In [ ]:
# Build a halfspace model or halfspace-like if topography is not *None*.
# #############################################################################
M.build_halfspace_mesh()
In [7]:
# Build a layered-earth mesh, topography can be set for subsurface layers
# when initializing the **BlankWorld** or in this function call as keyword arg.
# If subsurface topo is used, *layer_depths* will be ignored.
# #############################################################################
# M.build_layered_earth_mesh(n_layers=2, # number of layers
# layer_depths=[-300.], # n-1 layer depths
# )
Add anomalies to the subsurface¶
- Anomalies are not allowed to intersect layers (for now)!
In [ ]:
# Add brick anomaly
# #############################################################################
M.add_brick(start=[-1000., -300., -200.], # lower left back corner
stop=[-500.0, 700.0, -700.], # upper right front corner
cell_size=1e4) # max. cell volume constraint
In [9]:
# Add dipping plate anomaly
# #############################################################################
# M.add_plate(500., 500., 50., # length, width, height
# [1200.0, 0.0, -500.], # origin (center) of plate
# 45., 117., # dip, dip azimuth
# cell_size=1e4) # max. cell volume constraint
Add further structures¶
- For instance, flight observation lines
In [10]:
# Add three 6km observation flight lines from y=-500 to y=500 with 500m spacing
# in 50m height above surface with 20m observation point spacing on each line.
# #############################################################################
# M.add_paths([mu.line_y(-3e3, 3e3, n_segs=300, x=-500., z=50.,
# topo=M.topo, topo_dir=M.t_dir)])
# M.add_paths([mu.line_y(-3e3, 3e3, n_segs=300, x=0., z=50.,
# topo=M.topo, topo_dir=M.t_dir)])
# M.add_paths([mu.line_y(-3e3, 3e3, n_segs=300, x=500., z=50.,
# topo=M.topo, topo_dir=M.t_dir)])
Add tetrahedron-boundary¶
- Used to increase the domain size for reducing boundary effects (low freq.)
In [ ]:
# The computational domain size is increased in x-, y- z-dir. by the
# three given factors, respecively.
# #############################################################################
M.there_is_always_a_bigger_world(1e1, 1e1, 1e1)
Call TetGen for meshing¶
- Automatically export Omega to .poly file in w_dir
- Call TetGen and export in .vtk for Paraview and .mesh (Medit format)
- Automatically copy .mesh filed to s_dir for automated conversion
In [ ]:
# Call TetGen with different command line options and set further arguments.
# for the TetGen options, it is referred to the TetGen documentation
# #############################################################################
M.call_tetgen(tet_param='-pq1.3aA',
# tet_param='default, # '-pq1.2aA',
# tet_param='raw', # '-p'
# tet_param='faces', # '-pMA'
export_vtk=True # export *.vtk* for Paraview
# print_infos=True, # print mesh infos at the end
# export_before=True, # export *.poly* file
# # # # # # # # # # # # # # # # # # automatically
# copy=True # copy to **s_dir**
)