Clarify what is meant by "kernel"

[leouieda]

It would be great to clarify this in the Project Goals. For example, I noticed that the functions don't include universal constants or physical properties. For gravity, this is ok since getting the gravity field is a matter of simple multiplication. For magnetics, the combination is a bit more difficult since it involves a dot product with the vector magnetisation. So it would be worth having that in a function already instead of requiring users to code it. It also may be necessary since doing it out of the jit compiled function would mean losing parallel and low memory execution.

It would be worth clarifying the scope of the package with this in mind. A good way to go about it would be to try to add a magnetic dipole implementation to see how it fits with the current paradigm.

Or maybe the "kernel" could mean the version without the for loops? It may be worth making that available in the API for use in external jit-compiled code (like tesseroids).

[santisoler]

I agree that the use of "kernel" is a too loose at the moment. I still struggle to give a closed definition to it. I've been using this name for grouping a set of low-level functions used in the forward models for gravity and magnetics.

The name "kernel" is inherited from the integral kernels involved in the forward modelling functions of potential fields. For example, for any three dimensional body $\Gamma$, a potential field $\phi$ on an observation point $\mathbf{p}$ can be computed as:

$$ \phi(\mathbf{p}) = \int\limits_\Gamma g(\mathbf{p}) k(\mathbf{p}, \mathbf{q}) \text{d}\mathbf{q} $$

where $k(\mathbf{p}, \mathbf{q})$ is the integral kernel and $\mathbf{g(\mathbf{p})}$ is a function that doesn't depend on the integration variable and it's usually constant for gravity and magnetic forward modelling.

For the gravity potential it would be:

$$ V\text{grav}(\mathbf{p}) = G \int\limits\Gamma \frac{ \rho(\mathbf{q}) }{ \lVert \mathbf{p} - \mathbf{q} \rVert } \text{d}\mathbf{q} $$

And for magnetic potential (in absence of currents):

$$ V_\text{mag}(\mathbf{p}) = \frac{\mu0}{4\pi} \int\limits\Gamma \mathbf{M(\mathbf{q})} \cdot \nabla_\mathbf{q} \left( \frac{1}{\lVert \mathbf{p} - \mathbf{q} \rVert} \right) \text{d}\mathbf{q} $$

If we want to continue with this idea of a "kernel" function, then the physical properties (density and magnetization) must be included in them.

It would be worth clarifying the scope of the package with this in mind. A good way to go about it would be to try to add a magnetic dipole implementation to see how it fits with the current paradigm.

I agree. In my head the implementation of the magnetic dipole would ask the coordinates of the observation point, the coordinates of the dipole and its magnetization vector.

Or maybe the "kernel" could mean the version without the for loops? It may be worth making that available in the API for use in external jit-compiled code (like tesseroids).

For now I'm only adding functions without the for loops, so: single computation point and single source. There is a broad interest in having those functions without the for loops.

---

Nevertheless, there are some caveats. We perform the prism forward model by evaluating the analytical solution to these integrals, so we are not actually computing the integral kernel. Should we use a different name for that function then? Which one? Could we still include it in the "kernel" group?

[santisoler]

After a discussion we had during a meeting this week we decided to stick with the following design for Choclo.

Everything is JIT compiled

We could offer non-jit compiled versions of our functions, although there are some caveats. Some functions are written to support floats (like the safe_log and safe_atan ones). If the would have to support comparisons for arrays, we would need to introduce some checks there. We want these functions to run quick and dirty, no extra steps.

Besides, the forward modelling functions we will offer in Choclo (single source, single computation point) are meant to be run using Numba for loops, so it makes sense for everything to be JIT compiled.

No external loops

Choclo will not offer (at least for now) functions that compute forward models over sets of computations points and/or sets of sources. It will only offer forward modelling functions for single computation point and single source. The reason behind this is that the outer for loops can be optimized and written in many different combinations depending on the specific problem the user want to tackle. It's easier to leave that part to them rather than offering a collection of different functions that perform almost the same, but through different ways. Besides, with the functions that Choclo will offer, these outer loops only cost two to three extra lines, e.g.:

Full forward model

import choclo as ch

@numba.jit(parallel=True)
def gravity_potential(...):
    result = np.zeros(N)
    for point in prange(N):
        for source in range(M):
            result[point] += ch.point.gravity_pot(...)
    return result

Build Jacobian matrix

import choclo as ch

@numba.jit(parallel=True)
def jacobian(...):
    A = np.empty(shape)
    for i in prange(shape[0]):
        for j in range(shape[1]):
            A[i, j] = ch.point.gravity_pot(..., density=1)

Architecture

Choclo will be split in submodules, one for each source type plus a utils one (will include the distance functions).

Each submodule will have a two level system of functions.

1. Level 1: Low-level kernel functions that only perform the math. They don't call any other function within Choclo. For example, for point sources, they compute 1/r and its derivatives. They will use the following name convention: kernel_pot (for the kernel function used in the potential field), kernel_e, kernel_n, kernel_u (for the kernels used in the gradient of the potential) and kernel_ee, kernel_nn, kernel_uu, kernel_en, kernel_ez, kernel_nz (for the kernels used in the tensor components and magnetic field).
2. Level 2: Higher-level functions that perform the forward modelling of a single source on a single computation point. These will depend on the Level 1 functions and on the utils functions. These functions will take the physical property of the source (density, magnetization vector) and will make use of universal constants. They will return every result in SI units. Their naming conventions will be:

   gravity_pot(coordinates, source, density)
   gravity_e
   gravity_n
   gravity_u
   gravity_ee
   gravity_en
   gravity_eu
   gravity_nn
   gravity_nu
   gravity_uu
   magnetic_e(coordinates, source, magnetization_vector)
   magnetic_n
   magnetic_u
   magnetic_tf_anomaly(coordinates, source, magnetization_vector, igrf_inc, igrf_dec)

Docs

The docs will include instructions on how to use Choclo functions properly and in an optimized way, showing examples on how to build a full forward model and a Jacobian matrix.