cartesian dot product used in spherical coordate system

[MFraters]

@bangerth mentioned in #301 that the dot product as defined in the point class in only valid for cartesian coordiante systems. I agree that is correct, but I have been using it in 2d spherical system as well as an approximation. This approximation seems to be working well enough for the cases where I had tested it, but I think it will only hold if the point are not too far away from each other.

I am currently viewing fixing this issue as a way to gain some accuracy, but I am happy to be proven wrong.

[MFraters]

With 2D I mean that I omit the the radius, so the 2D point contains longitude and latitude, so unit wise it is fine (both angles in spherical coordinate system).

[bangerth]

Locally every smooth manifold is Cartesian, so if the two points are close to each other, then your approach works.

But take a point on the equator, x=(phi=0, psi=0), and any other point y. In your scheme, their dot product is zero suggesting that the two vectors are orthogonal when that is clearly not the case. In fact, x*x=0 even though the point is not the origin, just a point on the sphere. A similar paradoxon appears if you choose x=(phi=0,psi=pi/2), a point on the equator. Then choose y=(phi=pi/2,psi=0), the north pole. You have x*y=0 as expected because the two points are perpendicular to each other. But y=(phi=pi/2, psi=pi/2) is also the north pole and x*y=pi^2/4, suggesting that the two points are not perpendicular to each other when clearly they are.

[MFraters]

Thanks for the explanation.

Even though this is only used to define the details of the trench and fault lines, which are inherently small scale, the corner cases are worse that I expected. This approach is actually only in use in the old/simple system. The new continuous slab and faults use a completely different way of getting the information (actually computing the distance at all points which form the line of the line and three points in between, and then doing a search in the area with the smallest distance to get a better resolution). This new method also has the advantage of making the slabs and faults continuous.

I think everyone should switch to the new approach anyway, but it is slightly more expensive if I remember correctly and I am not sure how much that difference is affected by all the changes and optimalizations after that pull request. The following two benchmark results from the last merged pull request should be a fail comparison I think:

Benchmark	Master	Feature	Difference (99.9% CI)
Slab interpolation curved simple none	1.422 ± 0.039 (s=293)	1.421 ± 0.043 (s=343)	-0.8% .. +0.7%
Slab interpolation simple curved CMS	1.541 ± 0.051 (s=306)	1.545 ± 0.050 (s=280)	-0.7% .. +1.1%

For comparison, this is what they looked like at the time of Rene's optimalization (#289):

Benchmark	Master	Feature	Difference (99.9% CI)
Slab interpolation simple none	2.374 ± 0.074 (s=231)	1.512 ± 0.034 (s=236)	-37.1% .. -35.6%
Spherical slab interpolation simple none	2.874 ± 0.033 (s=168)	2.180 ± 0.028 (s=194)	-24.5% .. -23.8%

With this in mind I would personally be in favor of just getting rid of the old method. It would also allow me to clean up a lot of other code as well and make some new (potentially significant) optimizations, but it might be good if @alarshi first tried it to see if and how much of a difference in setup time it would make for her. The pull request (#237) has some pictures which explain the old options and the new option. I think the new option would be better in the long term anyway, because it is smoother and therefore better for the the solver.

[MFraters]

An other solution would be of course be to compute the correct dot product in spherical coordinates, but that would most likely take more computation time as well and I would need to investigate what that needed to be and how it would fit in the current computation of the fractions. So I think it would be easier to just only allow the dot product for Cartesian coordinates and get rid of the old interpolation option.