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**
)
Mesh generation with topography data¶
This tutorial is a guide to create topography meshes for calculations whit custEM. It is possible to either describe a synthetic topography as f(x, y) or to incorporate real DEM (digital elevation model) data. Synthetic topographies can be definied in the synthetic_definitions.py or a custom file or directly in this script as will be explained following.
As real DEM data files are quite large, we present the usage and structure by exporting synthetic topography data to both supported file formats and use these data to simulate real-world topography interpolation.
In the following, we will handle two seperate cases, named I and II, to show different options or possibilites.
General information for DEM files¶
- supported file formats are either xyz or asc files
- for now, the DEM data points must be located on a regular grid
- the files must be located in the topography directory t_dir
- xyz files contain x-, y-, z- values for each single data point. The positions must be sorted either ascending in x-y coordinate order or vice versy, in y-x order
- asc files contain a header with lower left corner and grid spacing, followed by a matrix of z-values arranged corresponding to the grid
- for more information, generate the example files as coded below and have look at the them
Import modules, utility functions or synthetic definitions¶
In [ ]:
from custEM.meshgen import meshgen_utils as mu
from custEM.meshgen.meshgen_tools import BlankWorld
import numpy as np
# import synthetic example-topography description, here valley.
# #############################################################################
from custEM.misc.synthetic_definitions import valley
t_dir = '.' # directory for topography files
case I - Create example xyz file¶
Here we use the valley function from synthetic_definitions
In [ ]:
v = np.linspace(0., 2e4, 801) # vector with 25 m spacing
xy = np.zeros((len(v)**2, 2)) # empty 2D array for all nodes in grid
a = 0
for j in range(len(v)): # assign coordinates for rows
for k in range(len(v)): # assign coordinates for columns
xy[a, 0] = v[j]
xy[a, 1] = v[k]
a += 1
# use imported synthetic topo-function to calculate z-values
z_data = valley(xy[:, 0], xy[:, 1])
valley_data = np.hstack((xy, z_data.reshape(-1, 1)))
# write xyz file
mu.write_synth_topo_to_xyz(t_dir, 'valley_DEM', valley_data)
case II - Create example asc file¶
for this example, we will use a locally definied synthetic function
In [ ]:
# Define topography function
# #############################################################################
np.random.seed(42) # set seed for reproducibility
def slope1(x, y, h=300., z=200.):
return(1e-4 * h * x - 1e-4 * h * y + z)
v = np.linspace(-2e4, 2e4, 401) # vector with 100 m spacing
xy = np.zeros((len(v)**2, 2)) # empty 2D array for all nodes in grid
a = 0
for j in range(len(v)): # assign coordinates for rows
for k in range(len(v)): # assign coordinates for columns
xy[a, 0] = v[j]
xy[a, 1] = v[k]
a += 1
z_data = slope1(xy[:, 0], xy[:, 1]) # some basic slope
z_data += np.random.rand(len(z_data)) * 100. # add the "hills"
# write asc file and assume we have some UTM 32 N coordinates, that's why
# we shift our just generated 40 x 40 km outcrop of topo data with help of
# *x_min* and *y_min* to the assumed "real" coordinate world
mu.write_synth_topo_to_asc(t_dir,
'hilly_DEM',
z_data.reshape(len(v), len(v)),
spacing=100., # header infor
x_min=703000, # lower left x coordinate
y_min=5602000) # lower left y coordinate
Create computational domains for cases I and II¶
Incorporating topography¶
It is possible to define topography (otherwiese topo=None) either with functions or DEM files. The toolbox will automatically calculate z-values for corresponding nodes if a function is used or interpolate z-values from the DEM using the scipy regulargridinterpolator on the fly.
The examples aim to show the usage of the following important arguments:
centering: used to shift the computational domain Omega relative to the center of the coordinate frame. This is in particular useful to shift, for instance, a 20 x 20 km outcrop of DEM data to the modeling domain.
easting_shift & northing_shift: shift Omega relative to the extend of the DEM in x- and y- direction by these values. This can be used, for instance, to shift a start- or end-point of the transmitter in the model domain to the coordiantes measured in the field and hence, simplify the models setup.
rotation: rotate clockwise, starting with positive x-axis direction, the coordinate frame in the field relative to the model domain. This is useful if a survey setup with parallel and perpedicular lines should correspond to x- and y-axes in the model domain
In [ ]:
# Incorporate topography with *asc* file
# #############################################################################
M1 = BlankWorld( # create blank 3D domain Omega
name='hilly_mesh', # name of the output mesh
m_dir='./meshes', # main mesh directory
t_dir=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 = False, # store mesh generation script as backup
#
# keyword arguments regarding building the 2D surface mesh.
#
inner_area='box', # refinement of area at surface
inner_area_size=[2e3, 1e3], # size of this area
# inner_area_shift=[0., 0.], # shift in x- and y-dir of area
inner_area_cell_size=1e3, # max. triangle-size (MA) in area
outer_area_cell_size=1e5, # MA outside of this area
# keyword arguments regarding incorporation of topography
#
topo='hilly_DEM.asc', # topography (from DEM or synthetic),
# # # # # # # # # # # # # # # # # e.g. topo=valley
# subsurface_topos=None, # subsurface interface topo (DEM, synth.)
easting_shift=-723000, # shift DEM east. relative to model domain
northing_shift=-5622000, # shift DEM north. relative to model domain
# centering=False, # center DEM to model domain
rotation=127., # rotate DEM around relative to Omega in °
)
# Incorporate topography with *xyz* file
# #############################################################################
M2 = BlankWorld( # blank 3D domain
name='valley_mesh_from_xyz', # name of the output mesh
m_dir='./meshes', # main mesh directory
t_dir=t_dir, # topo dir., **p_dir**+ '/topo'
backup_script = False, # store mesh generation script as backup
#
# keyword arguments regarding building the 2D surface mesh.
#
# inner_area='ellipse', # refinement of area at surface
# inner_area_size=[2e3, 1e3], # size of this area
# inner_area_shift=[0., 0.], # shift in x- and y-dir of area
# inner_area_cell_size=1e3, # max. triangle-size (MA) in area
# outer_area_cell_size=1e5, # MA outside of this area
# keyword arguments regarding incorporation of topography
#
topo='valley_DEM.xyz', # topography (from DEM or synthetic),
# # # # # # # # # # # # # # # # # e.g. topo=valley
# subsurface_topos=None, # subsurface interface topo (DEM, synth.)
# easting_shift=-723000, # shift DEM east. relative to model domain
# northing_shift=-5622000, # shift DEM north. relative to model domain
centering=True, # center DEM to model domain
# rotation=None, # rotate DEM realative to model domain
)
# Incorporate topography with f(x,y) definition with a function
# #############################################################################
M3 = BlankWorld( # blank 3D domain
name='valley_mesh_from_func', # name of the output mesh
m_dir='./meshes', # main mesh directory
t_dir=t_dir, # topo dir., **p_dir**+ '/topo'
backup_script = False, # store mesh generation script as backup
#
# keyword arguments regarding building the 2D surface mesh.
#
# inner_area='ellipse', # refinement of area at surface
# inner_area_size=[2e3, 1e3], # size of this area
# inner_area_shift=[0., 0.], # shift in x- and y-dir of area
# inner_area_cell_size=1e3, # max. triangle-size (MA) in area
# outer_area_cell_size=1e5, # MA outside of this area
# keyword arguments regarding incorporation of topography
#
topo=valley, # topography (from DEM or synthetic),
# # # # # # # # # # # # # # # # # e.g. topo=valley
# subsurface_topos=None, # subsurface interface topo (DEM, synth.)
# easting_shift=-723000, # shift DEM east. relative to model domain
# northing_shift=-5622000, # shift DEM north. relative to model domain
# centering=False, # center DEM to model domain
# rotation=None, # rotate DEM realative to model domain
)
Create surface mesh¶
In [ ]:
# Build the 2D surface mesh and incorporate transmitter or observation lines.
# #############################################################################
M1.build_surface(
# insert_lines=[mu.line_x(-5e3, 5e3, 1000)], # add observation line
insert_loop_tx=[mu.loop_r(start=[-5e2, -5e2],
stop=[5e2, 5e2], n_segs=100)], # add Tx loop_r
)
M2.build_surface(
insert_line_tx=[mu.line_x(-1e3, 1e3, 1000)], # add line Tx
)
M3.build_surface(
insert_line_tx=[mu.line_x(-1e3, 1e3, 1000)], # add line Tx
)
Build 3D world based on 2D surface mesh¶
- Build either a fullspace, halfspace or layered-earth mesh
In [ ]:
# The most simple world is a fullspace, in this case the **build_surface**
# method does not need to be executed or is ignored, respectively.
# #############################################################################
# M1.build_fullspace_mesh()
In [ ]:
# Build a halfspace model or halfspace-like if topography is not *None*.
# #############################################################################
M1.build_halfspace_mesh()
In [ ]:
# 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.
# #############################################################################
M2.build_layered_earth_mesh(n_layers=2, # number of layers
layer_depths=[-1200.], # n-1 layer depths
)
M3.build_layered_earth_mesh(n_layers=2, # number of layers
layer_depths=[-1200.], # n-1 layer depths
)
Add anomalies to the subsurface¶
- Anomalies are not allowed to intersect layers (for now)!
In [ ]:
# Add brick anomaly
# #############################################################################
M1.add_brick(start=[-700., 300., -500.], # lower left back corner
stop=[-500., 700., -1000.], # upper right front corner
cell_size=1e4) # max. cell volume constraint
In [ ]:
# Add dipping plate anomaly
# #############################################################################
#M2.add_plate(500., 500., 50., # length, width, height
# [0.0, 0.0, -700.], # origin (center) of plate
# 45., 117., # dip, dip azimuth
# cell_size=1e4) # max. cell volume constraint
#M3.add_plate(500., 500., 50., # length, width, height
# [0.0, 0.0, -700.], # origin (center) of plate
# 45., 117., # dip, dip azimuth
# cell_size=1e4) # max. cell volume constraint
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.
# #############################################################################
M1.extend_world(1e1, 1e1, 1e1)
# M2.there_is_always_a_bigger_world(1e2, 1e2, 1e2)
# M3.there_is_always_a_bigger_world(1e2, 1e2, 1e2)
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
# #############################################################################
M1.call_tetgen(tet_param='-pq1.2aA',
# 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**
)
M2.call_tetgen(export_vtk=True)
M3.call_tetgen(export_vtk=True)
FEM computation
Run tutorial¶
This tutorial is a guide to compute finite-element solutions for CSEM setups with custEM. In the current state, we considered a small model which can be run without MPI in serial. It is similar to example '#'1 in the examples directory. In addition, a plate anomaly was incorporated in the model and the transmitter is a rectangular loop now. Please note that the documentation of custEM is still in development as is this tutorial.
In [ ]:
# imports
import dolfin as df
import numpy as np
from custEM.meshgen.meshgen_tools import BlankWorld
from custEM.meshgen import meshgen_utils as mu
from custEM.core import MOD
On-the-fly mesh generation¶
In [ ]:
# Alternatively to calling a mesh generation script like the *meshgen_tuorial*
# or the *topo_mesh_tutorial* before running a simulation script,
# we can also create the mesh on-thy-fly in the simulation script, which is
# called from the command line, e.g. *mpirun -n 4 python -u run_tutorial.py*
# #############################################################################
if df.MPI.rank(df.MPI.comm_world) == 0:
# create world with default dimensions of 1e4 m in each direction
O = BlankWorld(name='test_mesh',
m_dir='./meshes',
preserve_edges=True)
# specify 10 km x-directed surface observation line
rx1 = mu.line_x(start=-5e3, stop=5e3, n_segs=100)
rx1_tri = mu.refine_rx(rx1, 5., 30)
# specify 2x2 km airborne observation grid
n = 21
vec = np.linspace(-1000., 1000., n)
rx2 = np.zeros((n**2, 3))
counter = 0
for x in vec:
for y in vec:
rx2[counter, 0] = x
rx2[counter, 1] = y
counter += 1
rx2[:, 2] = 50. # set height of points to 50 m
rx2_tri = mu.refine_rx(rx2, 10., 30)
O.build_surface(
insert_line_tx=[mu.line_y(start=-5e2, stop=5e2,
n_segs=20)], # add first line Tx
insert_loop_tx=[mu.loop_r(start=[-5e2, -5e2], stop=[5e2, 5e2],
n_segs=40)], # add second loop Tx
insert_paths=rx1_tri) # add surface Rx positions
# build halfspace mesh and extend 2D surface mesh to 3D world
O.build_layered_earth_mesh(n_layers=3,
layer_depths=[-300., -1000.])
# add airborne Rx
O.add_paths(rx2_tri) # add refinement triangles
O.add_rx(rx1) # Rx info for parameter file
O.add_rx(rx2) # Rx info for parameter file
# Add plate anomaly
O.add_plate(500., 500., 100., # length, width, height
[100., 0., -700.], # origin (center) of plate
45., 117., # dip, dip azimuth
cell_size=1e4) # max. cell volume constraint
# Append boundary mesh
O.extend_world(1e1, 1e1, 1e1)
# call TetGen
O.call_tetgen(tet_param='-pq1.4aA', export_vtk=True)
else:
pass
Initialize MOD instance¶
In [ ]:
# Specify approach, meanwhile CSEM, TEM and MT approaches possible
# #############################################################################
# for controlled-source frequency-domain (CSEM)
approach = 'E_t' # total electric field formulation
# approach = 'E_s' # secondary electric field formulation
# approach = 'H_s' # secondary magnetic field formulation
# or one of the potential approaches 'Am_t', 'Am_s', 'An_t', 'An_s', 'Fm_s'
# for controlled-source time-domain (TEM)
# approach = 'E_IE' # implicit Euler time-stepping method
# approach = 'E_FT' # fourier-transform based method
# approach = 'E_RA' # rational Arnoldi method
# for natural-source frequency-domain (MT)
# approach = 'MT' # total electroc field formulation
# the following option is not implemented yet, but will be probably
# implemented in future as an alternative to the total-field formulation
# approach = 'MT_s' # secondary electric field formulation
M = MOD('test_mod', # mod name (required)
'test_mesh', # mesh name (required)
approach, # approach name (required)
p=1, # polynomial order
overwrite_mesh=True, # automatically overwrite existing mesh if a
# a newer version was created
overwrite_results=True, # overwrite existing results for this mesh
serial_ordering=True, # use serial ordering (for solver robustness)
debug_level=10, # define log level, 10 = info
# load_existing=False, # import existing results of this model,
# # requires *overwrite_results=False*
)
Set model parameters¶
- A single frequencies is specified as frequency or omega,
- Multiple frequencies are specified as frequencies or omegas
- The conductivity distribution needs to be set in the order layers or anoamlies were added to the mesh, e.g., this example requires a list of 4 entries for sigma_ground
- only secondary-field approaches require specification of sigma_0, here, we need 4 background resistivity values for sigma_0 to calculate delta_sigma on the fly
- sigma_air is 1e-8 by default and only relevant for total field approaches we found that a contrastto the subsurface resistivities of 4 orders of magnitude is totally sufficient to obtain accurate results
In [ ]:
# update model parameters
# #############################################################################
M.MP.update_model_parameters(
frequencies=[10, 100, 1000], # multiple frequencies (Hz)
# omega=1e1, # e.g., only one angular frequency
# current=1., # set equal current for all Tx
# currents=[1., -5.], # unequal source currents for each Tx
sigma_ground=[1e-3, 1e-2, # layer conductivities (top to bot.)
1e-4, 1e-1], # + one value for anomaly
# sigma_0=[1e-3, 1e-2, # background conducitivites, last one
# 1e-4, 1e-2], # for anomaly (in second layer)
# sigma_air=1e-8, # airspace conductivity
)
# - anisotropic parameters are specified by 3 (VTI) or 6 (general, upper
# triangle matrix of parameter tensor) values, e.g., to make the anomaly in
# our example general anistropic in an VTI anisotropic background layer, use
# sigma_ground=[1e-3,
# [2e-2, 2e-2, 1e-1],
# 1e-4,
# [1e-1, 2e-1, 3e-1, 4e-1, 5e-1, 6e-1]
# ] # note the general anisotropic values here are nonsense
# - user can also specifiy values for *mu_r*, *eps_r* in the same manner
# - users can also specify the IP parameters *ip_m*, *ip_c*, *ip_tau*, but
# this works only with isotropic conductivities for now (Cole-Cole model)
Build var form and solve main problem¶
In [ ]:
# set up variational formulation
# #############################################################################
# For the CSEM approaches, usually no additional keyword
# arguments are required, only for the full formulation or IP modeling
M.FE.build_var_form(# ip=True, # enable ip modeling
# quasi_static=True # enable full formulation
)
# time-domain approaches require specifications of the observation times,
# either provide an explicit *times* vector or let custEM choose appropriate
# *times* automatically based on the ranges (recommended for most applications)
# the following is only relevant for time-domain approaches
# M.FE.build_var_form(n_log_steps=41, # number of log time steps
# log_t_min=-4, # log-start of times
# log_t_may=0, # log-end of times
# n_lin_steps=30, # for IE approach only
# )
In [ ]:
# solve main problem and convert or export results
# #############################################################################
M.solve_main_problem(# convert_to_H=False, # turn off H-field conversion
# export_pvd=False, # turn off export for Paraview
# auto_interpolate=False # enable auto-interpolation
# # (only numpy array export)
)
Interpolation¶
In [ ]:
# create interpolation meshes and interpolate on them
# #############################################################################
M.IB.create_line_mesh(
'x', # x-directed line
start=-5000., # start
stop=5000., # stop
n_segs=100, # number of segments on line
z=-0.1, # go a bit into the earth
line_name='surf' # set custom line name,
)
M.IB.create_slice_mesh(
'z', # z-normal slice
dim=2e3, # total dimension,
n_segs=20, # numeber of segments in each direction
# dim=[2e3, 2e3], # use two values for defining a rectangle
# n_segs=[20, 20], # for different n_segs in each direction
slice_name='air', # set custom slice name
)
# lines in a specific direction get the ending *line_x/y/z*
# slices in a specific direction get the ending *slice_x/y/z*
M.IB.interpolate(
['surf_line_x', 'air_slice_z'], # interpolation meshes
# quantities=['E_t', 'H_t'] # specify fields to interpolate
)
# in secondary field approaches and if interested in interpolated secondary
# fields use *quantities=['E_t, 'H_t', 'E_s', 'H_s']*
# alternatively, users can use "auto_interpolate" to specify on all rx paths
# included during the mesh generation. In this case, no Paraview-files are
# exported and solutions for all Tx are written to a single numpy file
Visualization
Visualization tutorial¶
This plot tutorial introduces basic functions for importing interpolated data from CSEM or MT computations and depicting them latter in 1D, 2D or 3D. This tutorial is going to be extended to include even more useful commands for analyzing different models computed with custEM (mismatches, field variations) or compare ustEM models with external results for validation/cross-checking. Meanwhile, we refer to the API-documentation for many more options.
Import modules¶
In [ ]:
from custEM.post import make_plain_subfigure_box
from custEM.post import PlotFD as Plot
import matplotlib.pyplot as plt
Initialize Plot instance and import results¶
In [ ]:
# initialize Plot instance
# #############################################################################
P = Plot('test_mod', # model name
'test_mesh', # mesh name
'E_t', # approach name
# r_dir='results', # results directory
# s_dir='plots', # save directory for default plots
fig_size=[16, 12] # change figure size
)
line1 = 'surf_line_x'
slice1 = 'air_slice_z'
# If not overwritten by keywords, data will be always imported from the
# model specified before in the *Plot* class initialization. Overwriting
# the keyword arguemnts during the import allows adding furhter data from
# other models, meshes, or approaches for comparison puposes.
# Note that *tx* and *freq* keyword arguments only need to be specified in
# case of multiple Tx or Frequencies in the same model, otherwise results
# are only avaialble for one Tx or freq and identified automatically.
for ti in range(2):
P.import_line_data(line1, # line mesh name
# mod='mod_name', # specify alternative mod name
# mesh='mesh_name', # specify alternative mesh name
# mpproach='E_s', # specify alternative approach
tx=ti, # specify which Tx
freq=1, # load data of second freq 100 Hz
key='L' + str(ti) # give key to access data later on
)
P.import_slice_data(slice1, # line mesh name
tx=ti, # specify which Tx
freq=1, # load data of second freq 100 Hz
# stride=1, # downsample data to speed up plot
key='S' + str(ti) # give key to access data later on
)
Plot line and slice data¶
In [ ]:
# plot data and misfits
# #############################################################################
# if a title is specified, plots are saved in the *p_dir*, default 'plots'
# initialize a plot with all three components for E and H fields
P.plot_line_data(key='L0', title='Tx0')
# add further line data to the same plot
P.plot_line_data(key='L1', new=False, title='Tx1')
# make a more customizedline plot with amplitude and phase representation
# use custom coordinate and amplitude limits
# use specific matplotlib style keyword arguments
P.plot_line_data(key='L0', ap=True,
xlim=[-1., 1.], ylim=[1e-7, 1e-2],
EH='E', color='k', lw='0.5', ls='-',
marker='^', markersize='3', title='Tx0_amp_phase')
# plot slice data in the air
P.plot_slice_data(key='S1', EH='H', title='Tx0_air')
# plot misfits, rather nonsense to compare data from two different Tx
# in this case, but shows the functionality
P.plot_line_errors(key='L0', key2='L1', EH='E')
P.plot_slice_errors(key='S0', key2='S1', EH='H')
Plot 2D data in 3D¶
In [ ]:
# Will be updated soon with a proper example
# #############################################################################
# fig = plt.figure(figsize=(10, 7))
# P.add_3D_slice_plot(fig, 'K1_H_t', 'above_topo_slice_z', sp_pos=111,
# q_length=0.2, label=r'$\Re(\mathbf{|H|})$ $(V/ m)$',
# q_slicE='K2', q_name='K2_H_t', clim=[1e-8, 1e-4])
# plt.savefig('./plots/custom/H_r_3D_plot.png')
#
# fig = plt.figure(figsize=(10, 7))
# P.add_3D_slice_plot(fig, 'K1_H_t', 'above_topo_slice_z', sp_pos=111,
# q_length=0.2, label=r'$\Im(\mathbf{|H|})$ $(V/ m)$',
# q_slicE='K2', q_name='K2_H_t', clim=[1e-8, 1e-4],
# ri='imag')
# plt.savefig('./plots/custom/H_i_3D_plot.png')