Tesseroid forward modeling about harmonica HOT 11 CLOSED

Ideally we should separate the discretization algorithms from the forward modeling. They would be a pre-processing step. @santisoler gave the idea of having the algorithms convert the tesseroids into a series of point masses. The GLQ weights can be integrated into the densities. Then we pass the point masses to a spherical point mass modeling code. This way we only maintain 1 set of forward modeling functions. The tesseroid code would be mostly a discretization code.

from harmonica.

One important feature I want for this Tesseroid implementation is to work under geodetic coordinates rather than geocentric spherical ones. It's not hard to solve, we just need to convert original geodetic coordinates to spherical ones, perform the forward modelling and then reconvert the computation points to geodetic coordinates.

Maybe we could have an argument in order to specify if we are passing geodetic or spherical coordinates.

from harmonica.

Also, we must define how Tesseroids mesher are going to be build.
@leouieda gave the idea to use xarray in order to build them. Maybe with a xr.Dataset containing the latitude and longitude coordinates and the top, bottom and density data values, where each latitude, longitude point represents the center of the Tesseroid (the mesh will be padded by half the horizontal size of each Tesseroid).

from harmonica.

Would be nice to be able to choose the GLQ degree and not to hard code it.
This could help to research on how lower or higher degrees affect the accuracy and the computation time.

I think the roots and weights for the Legendre polynomials can be computed through the
Numpy's Legendre Module (numpy.polynomial.legendre).

from harmonica.

One important feature I want for this Tesseroid implementation is to work under geodetic coordinates rather than geocentric spherical ones.

Is there a particular reason for this? If we're doing inversion, it's not good for the forward modeling to perform all of these calculations. We have a public function for doing the conversion which is very simple and not much effort to call before passing in the coordinates for modeling. It is done on top of the script once. There is also need to do the inverse transform since we have these values. Adding this requires maintaining several if in the function and another function argument. I would be hesitant to do this unless there is a very compelling reason to do so.

from harmonica.

Also, we must define how Tesseroids mesher are going to be build.

This should probably be a separate issue. It will also relate to how we represent prism and point mass meshes.

from harmonica.

Would be nice to be able to choose the GLQ degree and not to hard code it.

I completely agree (as you know :)

I think the roots and weights for the Legendre polynomials can be computed through the Numpy's Legendre Module (numpy.polynomial.legendre).

Ooooh that's new! OK, that sounds great!

from harmonica.

One important feature I want for this Tesseroid implementation is to work under geodetic coordinates rather than geocentric spherical ones.

Is there a particular reason for this? If we're doing inversion, it's not good for the forward modeling to perform all of these calculations. We have a public function for doing the conversion which is very simple and not much effort to call before passing in the coordinates for modeling. It is done on top of the script once. There is also need to do the inverse transform since we have these values. Adding this requires maintaining several if in the function and another function argument. I would be hesitant to do this unless there is a very compelling reason to do so.

I didn't thought about the repeated coordinate conversion that would take place if we use this forward modelling on an inversion.
I was thinking more on the user interface rather than computation performance: even tesseroids work on spherical coordinates (and computationally would be better to make inversions on that coordinate system), most data is given in geodetic coordinates. So, for example, if a user that wants to forward model a buried body, the would need to be aware to make the coordinate conversion, i.e. writing one more line (that's my "not-so-compelling" reason haha).

In conclusion, I agree we should not make this coordinates conversions inside the function. Moreover, this would make tests a lot easier.

from harmonica.

Also, we must define how Tesseroids mesher are going to be build.

This should probably be a separate issue. It will also relate to how we represent prism and point mass meshes.

I agree. We can leave this point to a separated issue too.

from harmonica.

Update on arbitrary GLQ degree

The numpy.polynomial.legendre.leggauss() function gets the degree of the quadrature as argument and returns two arrays: the first one contains the non-scaled coordinates of the nodes and the second one contains the corresponding weights.

For example:

from numpy.polynomial.legendre import leggauss

print(leggauss(2))

>> (array([-0.57735027,  0.57735027]), array([1., 1.]))

We get the nodes and weights for the second order degree, the ones Fatiando 0.5 uses.

But if we ask for a third order degree:

from numpy.polynomial.legendre import leggauss

print(leggauss(3))

>> (array([-0.77459667,  0.        ,  0.77459667]),
    array([0.55555556, 0.88888889, 0.55555556]))

With this function we can easily implement arbitrary degree GLQ. We can even choose GLQ degree for each direction (longitude, latitude or radial). I would like to experiment on how the radial GLQ degree changes the accuracy and computation time.

from harmonica.

👍 that's great! The GLQ degrees should definitely be an optional argument to the modeling function

from harmonica.