Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
hz = [(dz, nbcz)]
mesh = TreeMesh([hx, hy, hz], x0='CCN')
# Refine based on surface topography
mesh = refine_tree_xyz(
mesh, topo, octree_levels=[2, 2], method='surface', finalize=False
)
# Refine box base on region of interest
xp, yp, zp = np.meshgrid([-100., 100.], [-100., 100.], [-80., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[2, 2], method='box', finalize=False
)
mesh.finalize()
##########################################################
# Create Magnetic Vector Intensity Model (MVI)
# --------------------------------------------
#
# Magnetic vector models are defined by three-component effective
# susceptibilities. To create a magnetic vector
# model, we must
#
# 1) Define the magnetic susceptibility for each cell. Then multiply by the
# unit vector direction of the inducing field. (induced contribution)
# 2) Define the remanent magnetization vector for each cell and normalized
mesh = TreeMesh([hx, hy, hz], x0='CCN')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz, octree_levels=[0, 0, 0, 0, 1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
mesh = refine_tree_xyz(
mesh, electrode_locs, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, yp, zp = np.meshgrid([-600., 600.], [-300., 300.], [-500., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 2], method='box', finalize=False
)
mesh.finalize()
###############################################################
# Create Conductivity Model and Mapping for OcTree Mesh
# -----------------------------------------------------
#
# Here we define the conductivity model that will be used to predict DC
# resistivity data. The model consists of a conductive sphere and a
# resistive sphere within a moderately conductive background. Note that
# you can carry through this work flow with a resistivity model if desired.
#
# Define conductivity model in S/m (or resistivity model in Ohm m)
# Create OcTree Mesh
# ------------------
#
# Here we define the OcTree mesh that is used for this example.
#
dh = 20. # base cell width
dom_width = 3000. # domain width
nbc = 2**int(np.round(np.log(dom_width/dh)/np.log(2.))) # num. base cells
# Define the base mesh
h = [(dh, nbc)]
mesh = TreeMesh([h, h, h], x0='CCC')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz, octree_levels=[0, 0, 0, 1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
mesh = refine_tree_xyz(
mesh, rx_locs, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, yp, zp = np.meshgrid([-300., 300.], [-300., 300.], [-400., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 4], method='box', finalize=False
)
mesh.finalize()
# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
hz = [(dz, nbcz)]
mesh = TreeMesh([hx, hy, hz], x0='CCN')
# Refine based on surface topography
mesh = refine_tree_xyz(
mesh, topo, octree_levels=[2, 2], method='surface', finalize=False
)
# Refine box based on region of interest
xp, yp, zp = np.meshgrid([-100., 100.], [-100., 100.], [-80., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[2, 2], method='box', finalize=False
)
mesh.finalize()
#######################################################
# Density Contrast Model and Mapping on OcTree Mesh
# -------------------------------------------------
#
# Here, we create the density contrast model that will be used to predict gravity gradiometry
# data and the mapping from the model to the mesh. The model
# consists of a less dense block and a more dense sphere.
#
# Define density contrast values for each unit in g/cc
background_val = 0.
# resistivity and IP data.
#
dh = 10. # base cell width
dom_width_x = 2400. # domain width x # domain width y
dom_width_z = 1200. # domain width z
nbcx = 2**int(np.round(np.log(dom_width_x/dh)/np.log(2.))) # num. base cells x
nbcz = 2**int(np.round(np.log(dom_width_z/dh)/np.log(2.))) # num. base cells z
# Define the base mesh
hx = [(dh, nbcx)]
hz = [(dh, nbcz)]
mesh = TreeMesh([hx, hz], x0='CN')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz[:, [0, 2]], octree_levels=[1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
dc_survey.getABMN_locations()
electrode_locations = np.r_[
dc_survey.a_locations, dc_survey.b_locations,
dc_survey.m_locations, dc_survey.n_locations
]
unique_locations = np.unique(electrode_locations, axis=0)
mesh = refine_tree_xyz(
mesh, unique_locations, octree_levels=[2, 4], method='radial', finalize=False
)
# resistivity and IP data.
#
dh = 10. # base cell width
dom_width_x = 2400. # domain width x
dom_width_z = 1200. # domain width z
nbcx = 2**int(np.round(np.log(dom_width_x/dh)/np.log(2.))) # num. base cells x
nbcz = 2**int(np.round(np.log(dom_width_z/dh)/np.log(2.))) # num. base cells z
# Define the base mesh
hx = [(dh, nbcx)]
hz = [(dh, nbcz)]
mesh = TreeMesh([hx, hz], x0='CN')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, xyz_topo[:,[0, 2]], octree_levels=[1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers. First we need to obtain the
# set of unique electrode locations.
dc_survey.getABMN_locations()
electrode_locations = np.c_[
dc_survey.a_locations, dc_survey.b_locations,
dc_survey.m_locations, dc_survey.n_locations
]
unique_locations = np.unique(
np.reshape(electrode_locations, (4*dc_survey.nD, 2)), axis=0
)
mesh = refine_tree_xyz(
dh = 25. # base cell width
dom_width = 1600. # domain width
nbc = 2**int(np.round(np.log(dom_width/dh)/np.log(2.))) # num. base cells
# Define the base mesh
h = [(dh, nbc)]
mesh = TreeMesh([h, h, h], x0='CCC')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz, octree_levels=[0, 0, 0, 1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
mesh = refine_tree_xyz(
mesh, rx_locs, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, yp, zp = np.meshgrid([-300., 300.], [-300., 300.], [-300., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 4], method='box', finalize=False
)
mesh.finalize()
###############################################################
# Create Resistivity Model and Mapping for OcTree Mesh
# ----------------------------------------------------
#
# Create OcTree Mesh
# ------------------
#
# Here we define the OcTree mesh that is used for this example.
#
dh = 25. # base cell width
dom_width = 1600. # domain width
nbc = 2**int(np.round(np.log(dom_width/dh)/np.log(2.))) # num. base cells
# Define the base mesh
h = [(dh, nbc)]
mesh = TreeMesh([h, h, h], x0='CCC')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz, octree_levels=[0, 0, 0, 1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
mesh = refine_tree_xyz(
mesh, rx_locs, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, yp, zp = np.meshgrid([-300., 300.], [-300., 300.], [-300., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 4], method='box', finalize=False
)
mesh.finalize()
mesh = TreeMesh([h, h, h], x0='CCC')
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz, octree_levels=[0, 0, 0, 1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
mesh = refine_tree_xyz(
mesh, rx_locs, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, yp, zp = np.meshgrid([-300., 300.], [-300., 300.], [-300., 0.])
xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 4], method='box', finalize=False
)
mesh.finalize()
###############################################################
# Create Resistivity Model and Mapping for OcTree Mesh
# ----------------------------------------------------
#
# Here, we define the electrical properties of the Earth as a resistivity
# model. The model consists of a long vertical conductive pipe within a more
# resistive background.
#
# Log-Resistivity in log[Ohm m]
air_val = 1e8
# Mesh refinement based on topography
mesh = refine_tree_xyz(
mesh, topo_xyz[:, [0, 2]], octree_levels=[1], method='surface', finalize=False
)
# Mesh refinement near transmitters and receivers
dc_survey.getABMN_locations()
electrode_locations = np.r_[
dc_survey.a_locations, dc_survey.b_locations,
dc_survey.m_locations, dc_survey.n_locations
]
unique_locations = np.unique(electrode_locations, axis=0)
mesh = refine_tree_xyz(
mesh, unique_locations, octree_levels=[2, 4], method='radial', finalize=False
)
# Refine core mesh region
xp, zp = np.meshgrid([-800., 800.], [-800., 0.])
xyz = np.c_[mkvc(xp), mkvc(zp)]
mesh = refine_tree_xyz(
mesh, xyz, octree_levels=[0, 2, 2], method='box', finalize=False
)
mesh.finalize()
###############################################################
# Project Surveys to Discretized Topography
# -----------------------------------------