/tutorials/Tutorial_3_kernel.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Tutorial 3

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 constructed the class is never seen, however
to know how it works and what it does really helps to understand how COMET is
working.
"""
# comet packages
from comet import pyhed as ph
from comet import snmr
# numerics
import numpy as np
# plotting
import matplotlib.pyplot as plt

if __name__ == '__main__':
    # set log level to 20 == Info level
    ph.log.setLevel(20)

    # %%
    # Quickly initialize some default stuff covered by the first two tutorials
    # =========================================================================
    # three half overlapping loops, 1D earth + earth magnetic field
    # init all 9 combinations for FID's
    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)
    # %%
    # Now the kernel class
    # =========================================================================
    # get 1D kernel for the first fid
    kern = site.createKernel(fid=0, dimension=1)
    # that's basically all you need to generate a valid working kernel class
    # let's take a look
    print(kern)
    # The kernel class itself has only references to the survey, loops and FID
    # Changes of parameters like pulses, loop turns etc are made there and
    # piped to the kernel.

    # However the kernel class comes with some features too.
    # The z discretization for a 1D kernel for example.
    kern.createZVector(90, -60, min_thk=0.5)
    # 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 in 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.5
    # thick. That is nice for inversions later. The kernel values are summed up
    # for the small thicknesses to have a more accurate kernel function in the
    # swallow subsurface.
    # Bye the way kern.setZVector() is also possible.

    # With this a mesh can be generated.
    # The mesh generation tool has a lot of flexibility, however
    # we will explain the different values elsewhere, here we just build the
    # default mesh.
    kern.create1DKernelMesh()

    # with this discretization we can calculate the kernel function
    # we will allow the usage of up to 4 CPUs (depending on the hardware)
    kern.calculate(num_cpu=4)
    # if we do this right here, the kernel is requesting a magnetic field
    # first, however there is none. Following the principle of supply and
    # demand, the magnetic field calculation routine of the loop (in this case)
    # loop 0 tries to get a magnetic field. Next thing, there is no mesh yet,
    # so further traceback will result in the creation of a mesh for the first
    # loop in order to calculate a magnetic field.
    # So with calling "kern.calculate(num_cpu=4)" first a mesh of loop 0, then
    # a magnetic field and last a kernel for FID 0 is computed.
    # In fact even the manual call of "kern.create1DKernelMesh()"
    # or even "kern.createZVector(...)" is not needed, however internally only
    # the default input parameter are used, which in some cases may not b
    # enough.

    # overview in the console
    print(kern)

    # and show
    kern.show(toplot=['real', 'imag', '0D'])

    # Save and load
    # =========================================================================
    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.

    kern2 = snmr.kernel.Kernel(name='kern_fid_0')
    print(kern2)
    # you see an empty survey class and a None as FID!, no loops are known
    # plotting would be possible, inversions too, but a recalculation would
    # fail

    # Let's fix this by reconnecting with survey and FID
    kern2.setSurvey(site, fid=0)
    print(kern2)
    # now we are back online, "kern2" is now basically exact the same as "kern"
    plt.show()

# The End