Skip to content

Commit f37523f

Browse files
committed
Azure CI deploy v0.25.0 from 5f060ec0a
1 parent 2f5d2b4 commit f37523f

File tree

18,609 files changed

+7159999
-1
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

18,609 files changed

+7159999
-1
lines changed

latest

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
v0.24.0
1+
v0.25.0

v0.25.0/.buildinfo

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
# Sphinx build info version 1
2+
# This file records the configuration used when building these files. When it is not found, a full rebuild will be done.
3+
config: b875a16033d1afb3fc33f43f228f2caa
4+
tags: 645f666f9bcd5a90fca523b33c5a78b7
Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
"""
2+
1D Forward Simulation for a Single Sounding
3+
===========================================
4+
5+
Here we use the module *simpeg.electromangetics.time_domain_1d* to predict
6+
the stepoff response for a single sounding over a 1D layered Earth.
7+
In this tutorial, we focus on the following:
8+
9+
- Defining receivers, waveforms, sources and the survey
10+
- The units of the model and resulting data
11+
- Defining and running the 1D simulation for a single sounding
12+
13+
Our survey geometry consists of a horizontal loop source with a radius of 6 m
14+
located 20 m above the Earth's surface. The receiver is located at the centre
15+
of the loop and measures the vertical component of the response.
16+
17+
18+
"""
19+
20+
#####################################################
21+
# Import Modules
22+
# --------------
23+
#
24+
25+
import numpy as np
26+
import os
27+
from matplotlib import pyplot as plt
28+
29+
from simpeg import maps
30+
import simpeg.electromagnetics.time_domain as tdem
31+
from simpeg.utils import plot_1d_layer_model
32+
33+
write_output = False
34+
plt.rcParams.update({"font.size": 16})
35+
36+
# sphinx_gallery_thumbnail_number = 2
37+
38+
#####################################################################
39+
# Create Survey
40+
# -------------
41+
#
42+
# Here we demonstrate a general way to define the receivers, sources, waveforms and survey.
43+
# For this tutorial, the source is a horizontal loop whose current waveform
44+
# is a unit step-off. The receiver measures the vertical component of the magnetic flux
45+
# density at the loop's center.
46+
#
47+
48+
# Source properties
49+
source_location = np.array([0.0, 0.0, 20.0])
50+
source_orientation = "z" # "x", "y" or "z"
51+
source_current = 1.0 # maximum on-time current
52+
source_radius = 6.0 # source loop radius
53+
54+
# Receiver properties
55+
receiver_location = np.array([0.0, 0.0, 20.0])
56+
receiver_orientation = "z" # "x", "y" or "z"
57+
times = np.logspace(-5, -2, 31) # time channels (s)
58+
59+
# Define receiver list. In our case, we have only a single receiver for each source.
60+
# When simulating the response for multiple component and/or field orientations,
61+
# the list consists of multiple receiver objects.
62+
receiver_list = []
63+
receiver_list.append(
64+
tdem.receivers.PointMagneticFluxDensity(
65+
receiver_location, times, orientation=receiver_orientation
66+
)
67+
)
68+
69+
# Define the source waveform. Here we define a unit step-off. The definition of
70+
# other waveform types is covered in a separate tutorial.
71+
waveform = tdem.sources.StepOffWaveform()
72+
73+
# Define source list. In our case, we have only a single source.
74+
source_list = [
75+
tdem.sources.CircularLoop(
76+
receiver_list=receiver_list,
77+
location=source_location,
78+
waveform=waveform,
79+
current=source_current,
80+
radius=source_radius,
81+
)
82+
]
83+
84+
# Define the survey
85+
survey = tdem.Survey(source_list)
86+
87+
88+
###############################################
89+
# Defining a 1D Layered Earth Model
90+
# ---------------------------------
91+
#
92+
# Here, we define the layer thicknesses and electrical conductivities for our
93+
# 1D simulation. If we have N layers, we define N electrical conductivity
94+
# values and N-1 layer thicknesses. The lowest layer is assumed to extend to
95+
# infinity. If the Earth is a halfspace, the thicknesses can be defined by
96+
# an empty array, and the physical property values by an array of length 1.
97+
#
98+
# In this case, we have a more conductive layer within a background halfspace.
99+
# This can be defined as a 3 layered Earth model.
100+
#
101+
102+
# Physical properties
103+
background_conductivity = 1e-1
104+
layer_conductivity = 1e0
105+
106+
# Layer thicknesses
107+
thicknesses = np.array([40.0, 40.0])
108+
n_layer = len(thicknesses) + 1
109+
110+
# physical property models
111+
model = background_conductivity * np.ones(n_layer)
112+
model[1] = layer_conductivity
113+
114+
# Define a mapping for conductivities
115+
model_mapping = maps.IdentityMap(nP=n_layer)
116+
117+
# Plot conductivity model
118+
thicknesses_for_plotting = np.r_[thicknesses, 40.0]
119+
120+
fig = plt.figure(figsize=(6, 5))
121+
ax = fig.add_axes([0.15, 0.15, 0.8, 0.8])
122+
plot_1d_layer_model(thicknesses_for_plotting, model, ax=ax, show_layers=False)
123+
plt.gca().invert_yaxis()
124+
125+
126+
#######################################################################
127+
# Define the Forward Simulation, Predict Data and Plot
128+
# ----------------------------------------------------
129+
#
130+
# Here we define the simulation and predict the 1D TDEM sounding data.
131+
# The simulation requires the user define the survey, the layer thicknesses
132+
# and a mapping from the model to the conductivities of the layers.
133+
#
134+
# When using the *simpeg.electromagnetics.time_domain_1d* module,
135+
# predicted data are organized by source, then by receiver, then by time channel.
136+
#
137+
138+
# Define the simulation
139+
simulation = tdem.Simulation1DLayered(
140+
survey=survey,
141+
thicknesses=thicknesses,
142+
sigmaMap=model_mapping,
143+
)
144+
145+
# Predict data for a given model
146+
dpred = simulation.dpred(model)
147+
148+
# Plot sounding
149+
fig = plt.figure(figsize=(6, 6))
150+
ax = fig.add_axes([0.2, 0.15, 0.75, 0.78])
151+
ax.loglog(times, dpred, "k-o", lw=2)
152+
ax.set_xlabel("Times (s)")
153+
ax.set_ylabel("|B| (T)")
154+
ax.set_title("Magnetic Flux")
155+
156+
157+
##################################################
158+
# Write Output (Optional)
159+
# -----------------------
160+
#
161+
162+
if write_output:
163+
dir_path = os.path.dirname(__file__).split(os.path.sep)
164+
dir_path.extend(["outputs"])
165+
dir_path = os.path.sep.join(dir_path) + os.path.sep
166+
167+
if not os.path.exists(dir_path):
168+
os.mkdir(dir_path)
169+
170+
np.random.seed(347)
171+
noise = 0.05 * np.abs(dpred) * np.random.randn(len(dpred))
172+
dpred += noise
173+
fname = dir_path + "em1dtm_data.txt"
174+
np.savetxt(fname, np.c_[times, dpred], fmt="%.4e", header="TIME BZ")
Lines changed: 158 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,158 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"\n# Forward Simulation of VRM Response on a Tree Mesh\n\nHere we use the module *simpeg.electromagnetics.viscous_remanent_magnetization*\nto predict the characteristic VRM response over magnetically viscous top soil.\nWe consider a small-loop, ground-based survey which uses a coincident loop\ngeometry. For this tutorial, we focus on the following:\n\n - How to define the transmitters and receivers\n - How to define the survey\n - How to define a diagnostic physical property\n - How to define the physics for the linear potential fields formulation\n - How to include surface topography (if desired)\n - Modeling on an OcTree mesh\n\n\nNote that for this tutorial, we are only modeling the VRM response. A separate\ntutorial have been developed for modeling both the inductive and VRM responses.\n"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"## Import modules\n\n\n"
15+
]
16+
},
17+
{
18+
"cell_type": "code",
19+
"execution_count": null,
20+
"metadata": {
21+
"collapsed": false
22+
},
23+
"outputs": [],
24+
"source": [
25+
"from simpeg.electromagnetics import viscous_remanent_magnetization as vrm\nfrom simpeg.utils import plot2Ddata\nfrom simpeg import maps\n\nfrom discretize import TreeMesh\nfrom discretize.utils import mkvc, refine_tree_xyz, active_from_xyz\n\nimport numpy as np\nfrom scipy.interpolate import LinearNDInterpolator\n\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\n\n# sphinx_gallery_thumbnail_number = 2"
26+
]
27+
},
28+
{
29+
"cell_type": "markdown",
30+
"metadata": {},
31+
"source": [
32+
"## Defining Topography\n\nSurface topography is defined as an (N, 3) numpy array. We create it here but\nthe topography could also be loaded from a file. To keep the example simple,\nwe set flat topography at z = 0 m.\n\n\n"
33+
]
34+
},
35+
{
36+
"cell_type": "code",
37+
"execution_count": null,
38+
"metadata": {
39+
"collapsed": false
40+
},
41+
"outputs": [],
42+
"source": [
43+
"[x_topo, y_topo, z_topo] = np.meshgrid(\n np.linspace(-100, 100, 41), np.linspace(-100, 100, 41), 0.0\n)\nx_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)\nxyz_topo = np.c_[x_topo, y_topo, z_topo]"
44+
]
45+
},
46+
{
47+
"cell_type": "markdown",
48+
"metadata": {},
49+
"source": [
50+
"## Survey\n\nHere we define the sources, the receivers and the survey. For this exercise,\na coincident loop-loop system measures the vertical component of the VRM\nresponse.\n\n\n"
51+
]
52+
},
53+
{
54+
"cell_type": "code",
55+
"execution_count": null,
56+
"metadata": {
57+
"collapsed": false
58+
},
59+
"outputs": [],
60+
"source": [
61+
"# Define the transmitter waveform. This strongly determines the behaviour of the\n# characteristic VRM response. Here we use a step-off. The off-time begins at\n# 0 s.\nwaveform = vrm.waveforms.StepOff(t0=0)\n\n# Define the time channels for the receivers. The time channels must ALL be\n# ALL the off-time defined by the waveform.\ntime_channels = np.logspace(-4, -1, 31)\n\n# Define the transmitter and receiver locations. This step will define the\n# receivers 0.5 m above the Earth in the even you use more general topography.\nx = np.linspace(-40.0, 40.0, 21)\ny = np.linspace(-40.0, 40.0, 21)\nx, y = np.meshgrid(x, y)\nx, y = mkvc(x.T), mkvc(y.T)\nfun_interp = LinearNDInterpolator(np.c_[x_topo, y_topo], z_topo)\nz = fun_interp(np.c_[x, y]) + 0.5 # sensor height 0.5 m above surface.\nlocations = np.c_[mkvc(x), y, z]\n\n# Define the source-receiver pairs\nsource_list = []\nfor pp in range(0, locations.shape[0]):\n # Define dbz/dt receiver\n loc_pp = np.reshape(locations[pp, :], (1, 3))\n receivers_list = [\n vrm.receivers.Point(\n loc_pp, times=time_channels, field_type=\"dbdt\", orientation=\"z\"\n )\n ]\n\n dipole_moment = [0.0, 0.0, 1.0]\n\n # Define the source\n source_list.append(\n vrm.sources.MagDipole(\n receivers_list, mkvc(locations[pp, :]), dipole_moment, waveform\n )\n )\n\n# Define the survey\nsurvey = vrm.Survey(source_list)"
62+
]
63+
},
64+
{
65+
"cell_type": "markdown",
66+
"metadata": {},
67+
"source": [
68+
"## Defining an OcTree Mesh\n\nHere, we create the OcTree mesh that will be used for the tutorial. Since only\nthe very near surface contributes significantly to the response, the dimensions\nof the domain in the z-direction can be small. Here, we are assuming the\nmagnetic viscosity is negligible below 8 metres.\n\n\n"
69+
]
70+
},
71+
{
72+
"cell_type": "code",
73+
"execution_count": null,
74+
"metadata": {
75+
"collapsed": false
76+
},
77+
"outputs": [],
78+
"source": [
79+
"dx = 2 # minimum cell width (base mesh cell width) in x\ndy = 2 # minimum cell width (base mesh cell width) in y\ndz = 1 # minimum cell width (base mesh cell width) in z\n\nx_length = 100.0 # domain width in x\ny_length = 100.0 # domain width in y\nz_length = 8.0 # domain width in y\n\n# Compute number of base mesh cells required in x and y\nnbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))\nnbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))\nnbcz = 2 ** int(np.round(np.log(z_length / dz) / np.log(2.0)))\n\n# Define the base mesh\nhx = [(dx, nbcx)]\nhy = [(dy, nbcy)]\nhz = [(dz, nbcz)]\nmesh = TreeMesh([hx, hy, hz], x0=\"CCN\")\n\n# Refine based on surface topography\nmesh = refine_tree_xyz(\n mesh, xyz_topo, octree_levels=[2, 2], method=\"surface\", finalize=False\n)\n\nmesh.finalize()"
80+
]
81+
},
82+
{
83+
"cell_type": "markdown",
84+
"metadata": {},
85+
"source": [
86+
"## Defining the Model\n\nFor the linear potential field formulation, the magnetic viscosity\ncharacterizing each cell can be defined by an \"amalgamated magnetic property\"\n(see Cowan, 2016). Here we define an amalgamated magnetic property model.\nThe model is made by summing a set of 3D Gaussian distributions.\n\nFor other formulations of the forward simulation, you may define the parameters\nassuming a log-uniform or log-normal distribution of time-relaxation constants.\n\n\n"
87+
]
88+
},
89+
{
90+
"cell_type": "code",
91+
"execution_count": null,
92+
"metadata": {
93+
"collapsed": false
94+
},
95+
"outputs": [],
96+
"source": [
97+
"# Find cells active in the forward simulation (cells below surface)\nind_active = active_from_xyz(mesh, xyz_topo)\n\n# Define 3D Gaussian distribution parameters\nxyzc = mesh.gridCC[ind_active, :]\nc = 3 * np.pi * 8**2\npc = np.r_[4e-4, 4e-4, 4e-4, 6e-4, 8e-4, 6e-4, 8e-4, 8e-4]\nx_0 = np.r_[50.0, -50.0, -40.0, -20.0, -15.0, 20.0, -10.0, 25.0]\ny_0 = np.r_[0.0, 0.0, 40.0, 10.0, -20.0, 15.0, 0.0, 0.0]\nz_0 = np.r_[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\nvar_x = c * np.r_[3.0, 3.0, 3.0, 1.0, 3.0, 0.5, 0.1, 0.1]\nvar_y = c * np.r_[20.0, 20.0, 1.0, 1.0, 0.4, 0.5, 0.1, 0.4]\nvar_z = c * np.r_[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]\n\n# Define model\nmodel = np.zeros(np.shape(xyzc[:, 0]))\nfor ii in range(0, 8):\n model += (\n pc[ii]\n * np.exp(-((xyzc[:, 0] - x_0[ii]) ** 2) / var_x[ii])\n * np.exp(-((xyzc[:, 1] - y_0[ii]) ** 2) / var_y[ii])\n * np.exp(-((xyzc[:, 2] - z_0[ii]) ** 2) / var_z[ii])\n )\n\n# Plot Model\nmpl.rcParams.update({\"font.size\": 12})\nfig = plt.figure(figsize=(7.5, 7))\n\nplotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan)\nax1 = fig.add_axes([0.09, 0.12, 0.72, 0.77])\nmesh.plot_slice(\n plotting_map * model,\n normal=\"Z\",\n ax=ax1,\n ind=0,\n grid=True,\n clim=(np.min(model), np.max(model)),\n pcolor_opts={\"cmap\": \"magma_r\"},\n)\nax1.set_title(\"Model slice at z = 0 m\")\nax1.set_xlabel(\"x (m)\")\nax1.set_ylabel(\"y (m)\")\n\nax2 = fig.add_axes([0.83, 0.12, 0.05, 0.77])\nnorm = mpl.colors.Normalize(vmin=np.min(model), vmax=np.max(model))\ncbar = mpl.colorbar.ColorbarBase(\n ax2, norm=norm, orientation=\"vertical\", cmap=mpl.cm.magma_r\n)\ncbar.set_label(\"Amalgamated Magnetic Property (SI)\", rotation=270, labelpad=15, size=12)\n\nplt.show()"
98+
]
99+
},
100+
{
101+
"cell_type": "markdown",
102+
"metadata": {},
103+
"source": [
104+
"## Define the Simulation\n\nHere we define the formulation for solving Maxwell's equations. We have chosen\nto model the off-time VRM response. There are two important keyword arguments,\n*refinement_factor* and *refinement_distance*. These are used to refine the\nsensitivities of the cells near the transmitters. This improves the accuracy\nof the forward simulation without having to refine the mesh near transmitters.\n\n\n"
105+
]
106+
},
107+
{
108+
"cell_type": "code",
109+
"execution_count": null,
110+
"metadata": {
111+
"collapsed": false
112+
},
113+
"outputs": [],
114+
"source": [
115+
"# For this example, cells lying within 2 m of a transmitter will be modeled\n# as if they are comprised of 4^3 equal smaller cells. Cells within 4 m of a\n# transmitter will be modeled as if they are comprised of 2^3 equal smaller\n# cells.\nsimulation = vrm.Simulation3DLinear(\n mesh,\n survey=survey,\n active_cells=ind_active,\n refinement_factor=2,\n refinement_distance=[2.0, 4.0],\n)"
116+
]
117+
},
118+
{
119+
"cell_type": "markdown",
120+
"metadata": {},
121+
"source": [
122+
"## Predict Data and Plot\n\n\n"
123+
]
124+
},
125+
{
126+
"cell_type": "code",
127+
"execution_count": null,
128+
"metadata": {
129+
"collapsed": false
130+
},
131+
"outputs": [],
132+
"source": [
133+
"# Predict VRM response\ndpred = simulation.dpred(model)\n\n# Reshape for plotting\nn_times = len(time_channels)\nn_loc = locations.shape[0]\ndpred = np.reshape(dpred, (n_loc, n_times))\n\n# Plot\nfig = plt.figure(figsize=(13, 5))\n\n# Index for what time channel you would like to see the data map.\ntime_index = 10\n\nv_max = np.max(np.abs(dpred[:, time_index]))\nv_min = np.min(np.abs(dpred[:, time_index]))\nax11 = fig.add_axes([0.12, 0.1, 0.33, 0.85])\nplot2Ddata(\n locations[:, 0:2],\n -dpred[:, time_index],\n ax=ax11,\n ncontour=30,\n clim=(v_min, v_max),\n contourOpts={\"cmap\": \"magma_r\"},\n)\nax11.set_xlabel(\"x (m)\")\nax11.set_ylabel(\"y (m)\")\ntitlestr = \"- dBz/dt at t=\" + \"{:.1e}\".format(time_channels[time_index]) + \" s\"\nax11.set_title(titlestr)\n\nax12 = fig.add_axes([0.46, 0.1, 0.02, 0.85])\nnorm1 = mpl.colors.Normalize(vmin=v_min, vmax=v_max)\ncbar1 = mpl.colorbar.ColorbarBase(\n ax12, norm=norm1, orientation=\"vertical\", cmap=mpl.cm.magma_r\n)\ncbar1.set_label(\"$T/s$\", rotation=270, labelpad=15, size=12)\n\n# Indicies for some locations you would like to see the decay\nlocation_indicies = [0, 65, 217]\ncolor_flags = [\"k\", \"r\", \"b\"]\nlegend_str = []\n\nax2 = fig.add_axes([0.6, 0.1, 0.35, 0.85])\nfor ii in range(0, len(location_indicies)):\n ax2.loglog(time_channels, -dpred[location_indicies[ii], :], color_flags[ii], lw=2)\n legend_str.append(\n \"(\"\n + \"{:.1f}\".format(locations[location_indicies[ii], 0])\n + \" m, \"\n + \"{:.1f}\".format(locations[location_indicies[ii], 1])\n + \" m)\"\n )\n\nax2.set_xlim((np.min(time_channels), np.max(time_channels)))\nax2.set_xlabel(\"time [s]\")\nax2.set_ylabel(\"-dBz/dt [T/s]\")\nax2.set_title(\"Decay Curve\")\nax2.legend(legend_str)"
134+
]
135+
}
136+
],
137+
"metadata": {
138+
"kernelspec": {
139+
"display_name": "Python 3",
140+
"language": "python",
141+
"name": "python3"
142+
},
143+
"language_info": {
144+
"codemirror_mode": {
145+
"name": "ipython",
146+
"version": 3
147+
},
148+
"file_extension": ".py",
149+
"mimetype": "text/x-python",
150+
"name": "python",
151+
"nbconvert_exporter": "python",
152+
"pygments_lexer": "ipython3",
153+
"version": "3.11.14"
154+
}
155+
},
156+
"nbformat": 4,
157+
"nbformat_minor": 0
158+
}

0 commit comments

Comments
 (0)