Tutorials¶
Tutorial 1 Loop¶
Tutorial 1 - The Loop class¶
This first tutorial aims to present the loop class to the user. The class itself represents a single closed or grounded loop of any shape in the field used for the calculation of magnetic or electric fields.
Only a few steps are necessary to initialize a loop instance, but many configurations are possible. In this sript the common methods of the loop class are presented and explained.
Start off with COMET¶
%matplotlib inline
# first we import the pyhed module from comet
from comet import pyhed as ph
# and the numeric library numpy
import numpy as np
The logger controls the console output (red background). The verbosity increases with lower values for the log level. default is 30 (warning) which means everything less important than a warning is prevented from appearing in the console. There are more verbose levels, including 20 (info), 16 (progress), 13 (trace) and 10 (debug) or more silent like 40 (error) and 50 (critical), however the last two are usually only printed when the code is to exit anyway.
ph.log.setLevel(20) # INFO
Defining the loop geometry¶
The loop is a general class that can have any shape and represents either a source or transmitter. The basic shape is defined additionally to the number of segments or maximum length for the dipole discretization. Predefined shapes utility functions are (see Loop API documentation):
buildCircle- perfect circular loopbuildRectangle/buildSquare- rectangle or squarebuildFig8/buildFig8Circle- figure-of-eight loopsbuildLoop- arbitrarily shaped loopbuildDipole/buildLine- single dipoles or a line of dipoles
We assume a 40m * 20m figure-of-eight loop, consisting of segments with a maximum length of 3m.
loop = ph.loop.buildFig8([[-20, -10], [20, 10]], max_length=2)
print(loop) # Show some information like positions & angles
# We can have a look at the loop by
ax = loop.show()
# and see the position of the segments
Define 1D resistivity and frequency¶
Those values are set and stored in the config of the loop class. The config is a small class that allows to share the values if multiple loops are handled. By default it contains a 1000 $\Omega$m halfspace with a Larmor frequency of 2000 Hz. Mostly, the B field is used which is defined by a TE mode for an ungrounded source.
print(loop.config)
# We set the model to a two-layer case (2nd layer infinite)
loop.setModel([1000, 10], d=[10])
# and can access those properties
print("rho =", loop.config.rho, ", d =", loop.config.d)
loop.setFrequency(2200) # loop.config.f
loop.setFType('B') # sets the field type ('B', 'H', or 'E')
print(loop.config)
Define mesh and export fields¶
The dipole positions of the loop don't have to be a part of any mesh. In order to calculate a field first coordinates the have to be given to the class object. The easiest way to do this is via a simple array. It is recommended to use pygimli Mesh for more sophisticated models or if values have to be interpolated afterwards.
first example: simple array¶
# If you want the field on the two coordinates (5, 5, -5) and
# (50, 50, -5), you would use the loop object like this:
coords = np.array([[5, 5, -5], [50, 50, -5]]).T
# shape: (3, n) for n points
loop.setLoopMesh(coords)
field = loop.calculate(num_cpu=1)
# the field values are stored in the variable loop.field
# we get another array of shape (3, n) for n points
# with complex valued vector field for each input coordinate
print('loop.field (shape =', loop.field.shape, '):')
print(loop.field)
# for more complex calculations we need to define more points or a mesh
# this can be done using the appropriate methods of the loops class:
# tetgen is needed for that.
loop.createLoopMesh()
# the mesh contains the source and is stored in loop.loopmesh
print(loop.loopmesh) # show mesh numbers
# now we calculate the new field for the mesh, this time with more cores
field = loop.calculate(num_cpu=4)
# for visualization the field can be exported as vtk.
# Note that this works only with a mesh, not with pure coordinates
loop.exportVTK('loop.vtk')
One can use Paraview (paraview.org) for visualization of the resulting magnetic field.
Tutorial 2 Survey and FID¶
Tutorial 2 - the Survey and FID classes¶
The tutorial for the Survey class is covering the definition and adjustments for simple FID experiments as well as the site parameters as magnetic field, tx and rx definition, pulse moments and many more.
Main Focus is on the import of a dataset using simple arrays. Thus tutorial explains all the steps needed and shows where COMET stores the different information.
Two complex classes are needed directly for this and presented here. Furthermore the loop class presented in the first example are used.
This script does not calculate anything extensive. It's only purpose is to give you a feeling on how to use the survey and fid class and where to put your input information.
%matplotlib inline
# import main modules from comet and other classes
from comet import pyhed as ph
from comet import snmr
import numpy as np # numerics
import matplotlib.pyplot as plt # plotting
ph.log.setLevel(20) # INFO, see Loop tutorial
Importing data 1¶
There are several ways to import data. Here we show you the manually way, convinience functions are available for MRSMatlab files. The basic function for importing raw data need at least four arrays:
- data (1D vector, complex),
- error (1D vector, complex),
- times (1D vector, float)
pulses (1D vector, float)
optionally phases (1D vector, float)
In this tutorial we take some previously saved vectors from a csv file (you have to modify the import with respect to your data file format). Numpy has some easy and very flexible options when it comes to importing ascii files, e.g.
data = np.genfromtxt('datafile', **kwargs)
See genfromtxt API doc for options.
# Here we simply import an npz (numpy compressed vector) file.
# importing those we get numpy arrays of various shapes (see console)
data_npz = np.load('tut2_data.npz')
# Data
data = data_npz['data']
print('data', data.shape, data.dtype)
# Error
error = data_npz['error']
print('error', error.shape, error.dtype)
# Pulses
pulses = data_npz['pulses']
print('pulses', pulses.shape, pulses.dtype)
# Phases (not needed for rotated amplitude inversion)
phases = data_npz['phases']
print('phases', phases.shape, phases.dtype)
# time vector
times = data_npz['times']
print('times', times.shape, times.dtype)
Importing data 2¶
Now we have several arrays and need to store them in the correct objects for COMET to use them. This is the main purpose of the survey class and the FID's. Every FID represents one measurement in the field, the survey represents the field and all the global configurations. The survey class also manages the loop classes introduced in the first tutorial.
# Now we create an empty survey class representing our site
site = snmr.survey.Survey()
print(site)
The survey class is empty. It holds no references to loops or measurements so far. In order to put our data in, we need to define our measurements. This needs all the loops used on this site so we can refer to them when initializing the FID's. We now quickly add some loops to the Survey (Loop classes are explained in detail in the first tutorial).
# We start defining three half-overlapping loops
# (see Tutorial 1 on Loop classes)
loop0 = ph.loop.buildRectangle([[-40, -20], [0, 20]], max_length=2)
loop1 = ph.loop.buildRectangle([[-20, -20], [20, 20]], max_length=2)
loop2 = ph.loop.buildRectangle([[0, -20], [40, 20]], max_length=2)
# Three 40 * 40 square loops from -40 to 40 meter in x
# let's add them to our manager survey class and have a look
site.addLoop([loop0, loop1, loop2])
print(site)
print(site.loops)
ph.plot.showLoopLayout(*site.loops);
The FID class¶
FID stands for free induction decay. See also its API Documentation. Finally we define which loop is used as transmitter and receiver to define our NMR experiments.
site.createSounding(tx=0, rx=0)
# A coincident sounding using the first (in Python 0!) loop as receiver and transmitter
# let's take all combinations:
for tx_i in range(3):
for rx_i in range(3):
site.createSounding(tx=tx_i, rx=rx_i)
print(site.fids)
# Note that the sounding using tx=0 and rx=0 is not created again.
As we are setting everthing up, we can might as well take a small detour to define our earth and resistivity.
For computing the NMR kernel, we require the strength and orientation of the magnetic field, given by the total field and its inclination and the declination of the profile. Note that this will set the frequency of the loops to the corresponding Larmor frequency.
# earth magnetic field at site
site.setEarth(mag=48000e-9, incl=60, decl=0)
print(site.earth)
# resistivtiy distribution at site
# 1 Layer case (1000 Ohmm, 10m thickness) + halfspace (10 Ohmm)
# take larmor frequency from earth magnetic field
config = ph.config(rho=[1000, 10], d=[10], f=site.earth.larmor)
# config is also a class, however it is little more than a container...
# ... but it can be piped to the survey class
site.setLoopConfig(config)
# This sets the resistivity configuration for all loops at once
# let's check the first loop right away
print(site.loops[0].config)
Finally the data¶
The FID's are the target for our array imported above. We show this at the first FID as example.
# we just grab the first one
fid = site.fids[0]
print(fid)
# Set data, error and times, phases (unrotated)
fid.setRawDataErrorAndTimes(data, error, times, phases=phases,
rotated=False)
# Pulse moments (in Ampere-seconds, NOT amperes!)
fid.setPulses(pulses)
# Note that if no pulse vector is given, a default one is available
# (e.g. for synthetic studies)
# For gating the data, define the number of desired gates
# (sometimesthe first gates are empty and therefore deleted,
# so you might end up a few gates short)
fid.gating(num_gates=42)
print(len(fid.gates), "gates")
# plot example for one fid
# make figure with one axis to plot stuff into
fig, ax = plt.subplots()
# fills a given plot axis with a gated fid plot
ax, cb = ph.plot.drawFid(ax, fid, to_plot='abs')
Some last things¶
There is more to the FID class. Here are a few examples on how to change the settings.
# There are a lot of setter functions, e.g.
fid.setFrequencyOffset(2.3) # in Hertz
# definition of time gates, if gating is not used
fid.setGates([0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0], midpoints=False)
# first gate between 0 and 0.1, second between 0.1 and 0.2 and so on...
# Note that instead of the gate centers its boundaries are set as
# the gate length is needed to compute individual noise levels.
# set pulse (s) and deadtime (s) of the measuring device
fid.setPulseDuration(0.04, deadtime_device=0.005)
# the effective deadtime is calculated internally summing up the deadtime
# of the device and half of the pulse duration
print('Effective deadtime: {:.3f} s'.format(fid.deadtime))
# Define the transmitter/receiver indices and turns of the loop
fid.setTx(0, turns=1) # stored in fid.tx_index and fid.tx_turns
# same for rx
fid.setRx(1, turns=1) # stored in fid.rx_index and fid.rx_turns
# Note: set tx and rx can be used if the loops of the survey class are imported after the FID data,
# however usually you define them at the creation of the fids via "createSounding(tx=...)"
# as shown above.
print(fid)
Tutorial 3 Kernel¶
Tutorial 3 - The Kernel class¶
After importing/defining everything in the Survey and FID classes the kernel class is used for the computation of the kernel functions.
If a simple synthetic case is conducted, the class remains hidden, however it is important to know how it works and what it does understand COMET.
# make sure graphics are plotted
%matplotlib inline
# import comet submodules
from comet import pyhed as ph
from comet import snmr
import numpy as np # numerics
import matplotlib.pyplot as plt # plotting
ph.log.setLevel(20)
Setting up the geometry¶
We first initialize some material that is covered by the first two tutorials. For the kernel, we assume three half overlapping 20x20m sized loops. All 9 combinations for FID's are considered.
loop0 = ph.loop.buildRectangle([[-40, -20], [0, 20]], max_length=2)
loop1 = ph.loop.buildRectangle([[-20, -20], [20, 20]], max_length=2)
loop2 = ph.loop.buildRectangle([[0, -20], [40, 20]], max_length=2)
site = snmr.survey.Survey()
site.addLoop([loop0, loop1, loop2])
site.setEarth(mag=48000e-9, incl=60, decl=0)
config = ph.config(rho=[1000, 10], d=[10], f=site.earth.larmor)
site.setLoopConfig(config)
for tx_i in range(3):
for rx_i in range(3):
site.createSounding(tx=tx_i, rx=rx_i)
print(site)
# Now we create a 1D kernel for the first FID
kern = site.createKernel(fid=0, dimension=1)
# This returns a Kernel class. let's take a look
print(kern)
# It contains references to the survey, loops and FIDs.
# Changes of survey parameters are made there and piped to the kernel.
# Let's have a look at the features of the kernel class,
# e.g. z discretization for a 1D kernel.
kern.createZVector(90, -60, min_thk=0.5) ## better keywords?
This creates a z-vector in 90 steps up to 60 meter depth. The layers are spaced using a sinus hyperbolicus function (more linear compared to log, which is nicer at larger depth).
The min_thk specifies the thickness of the final kernel discretization. Small layers around z = 0 are combined until each layer is at least 0.5m thick, which is better for inversions later. By summing up the values for the small thicknesses we obtain a more accurate kernel function in the shallow subsurface.
By kern.setZVector() the z vector can be set arbitrarily.
print(kern.zvec) # z always positive upwards
# With this a mesh can be generated.
kern.create1DKernelMesh()
# There are a lot of options, however we build a default mesh.
print(kern.kernelMesh1D)
# Now we calulate the kernel function using 4 CPUs
kern.calculate(num_cpu=4)
# First, it requests a magnetic field for all involved loops (here only 0).
# If it does not exist none already, it will trigger the computation
# following the supply-and-demand principle.
# This in turn requires a mesh, so it will be created.
# Last, the kernel for FID 0 is computed.
# Even the calls
# kern.create1DKernelMesh()
# kern.createZVector(...)
# are not needed, because they will be triggered using default parameters.
# Have a look at the kernels properties
print(kern)
# and at the kernel matrix itself
kern.show(toplot=['real', 'imag', '0D']);
# alternatively, amplitude and phase kernels can be plotted
kern.show(toplot=['amp', 'phase']);
kern.save(savename='kern_fid_0')
# saves kernel class including mesh and kernel matrix, but excluding survey
# and FID (just the references are saved == indices to the correct loops
# and the FID). A backup of the pulse moments is saved with the kernel but
# discarded, once the survey is loaded and connected back to the kernel.
# At a later stage, we can create a new kernel by loading
kern2 = snmr.kernel.Kernel(name='kern_fid_0')
print(kern2)
# Note the empty survey class and a None FID, no loops are known,
# plotting and inversion would be possible, but a recalculation fails
# We can fix this by reconnecting with survey and FID
kern2.setSurvey(site, fid=0)
print(kern2)
# now we are back at the previous stage (kern2==kern)