Closed santisoler closed 3 weeks ago
I think I tracked down what's going on with this. The issue is coming when evaluating the kernel_en
function, and particularly the safe_log
, and due to machine precision.
The issue arises when evaluating the kernel on the two prism vertices that fall below the observation point (the two nodes marked in green circles):
When the observation point falls exactly above the nodes, the kernels are well evaluated and the result is correct.
But when the observation point is slightly away from it, then the evaluation on the nodes falls into a problem:
This leads to evaluate the kernels on this two nodes using different branches of the safe_log function and therefore creating this discontinuity.
I'll think how we can solve this and comment back with some ideas.
I think I arrived to a solution for it.
The issue is generated when two independent calls to _safe_log
interpret that the observation point is not above one of the nodes (the closest one to it), while it interprets that it is above the other one. This is because for the first one the r
is not equal to x
(x
could be equal to -50.0
and r
will be something like 50.00000001
-just a minor difference in in machine precision). But for the other node, r
and x
have the same absolute value (e.g. x
is -70.0
and r
is 70.0
).
This triggers that two different branches of the safe_log function are being executed: one that evaluates the actual log function, and the other one that evaluates one of the terms of the limits of the kernel function. This creates an inconsistency: their difference is not the value of the analytic solution on that point.
One way to solve this would be to make _safe_log
to not only take the x
and the r
, but all the shifted coordinates. This way we can ensure that $|x| = r$ by checking if y == 0 and z == 0
. This way, even if r
and |x|
have the same value up to machine precision, we will not neglect the non-zero values of the other components.
Moreover, since we are taking y
and z
as arguments, we can use the form of the log function that Fukushima suggests for cases in which $x < 0, r \ne |x|$:
$$ \ln(x + r) = \ln \left( \frac{y^2 + z^2}{r - x} \right) \quad x < 0, r \ne |x|$$
This form avoid some floating point errors that appear when trying to evaluate the log on x + r
.
I'll open a PR with the bugfix.
Description of the problem:
Today I found that the current forward modelling functions for magnetic fields experience some discontinuities when being evaluated on lines above one of the edges of the prism.
I found this misbehavior when exploring the values of
magnetic_e
along a N-S profile that falls right above one of the horizontal edges of the prism.Full code that generated the error
Output
System information
main
(eaf80e78329fabe52f45c047950d76d6fadb2b35)conda list
below:Output of conda list