Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bug on safe_log and solves discontinuous magnetic fields #100

Merged
merged 5 commits into from
Oct 2, 2024

Conversation

santisoler
Copy link
Member

Fix how _safe_log checks if r == |x|. Add the other two shifted coordinates (y and z) as arguments for _safe_log so we can use them to assess that the observation point is inline with a pair of vertices or not. Make the check by evaluating if y == 0 and z == 0 instead of -x == r (with x < 0), since -x and r could be the same value (up to machine precision) while y and z are non-zero (specially when $|x| \gg |y|$ and/or $|x| \gg |z|$). Update calls to _safe_log, update its docstrings, and update docstrings of kernels that mentioned the safe log function. Add tests that catches the bug, and that pass with the bugfix.

Relevant issues/PRs:

Fixes #98

Add the other two shifted coordinates (`y` and `z`) as arguments for
`_safe_log` so we can use them to assess that the observation point is
inline with a pair of vertices or not. Make the check by evaluating if
`y == 0 and z == 0` instead of `-x == r` (with `x < 0`), since `-x` and `r`
could be the same value (up to machine precision) while `y` and `z` are
non-zero (specially when $|x| \gg |y|,|z|$). Update calls to
`_safe_log`, update its docstrings, and update docstrings of kernels
that mentioned the safe log function.
choclo/prism/_kernels.py Outdated Show resolved Hide resolved
choclo/prism/_kernels.py Outdated Show resolved Hide resolved
@santisoler
Copy link
Member Author

Running the snippet I pasted in #98 on this branch, I get the following image, showing that the fields are well evaluated with this bugfix:

Figure_2

@santisoler
Copy link
Member Author

@leouieda I'm merging this PR since the tests I wrote verify that the bugfix solves the issue without introducing new errors, and it doesn't impact the performance of the kernels. Feel free to comment if you have any thought!

@santisoler
Copy link
Member Author

Ok, one more check that this bugfixes improves the code.

Here is a snippet that evaluates the _safe_log on a set of observation points:

import numpy as np
import matplotlib.pyplot as plt

from choclo.prism._kernels import _safe_log


def eval_log(easting, northing, upward):
    radius = np.sqrt(easting**2 + northing**2 + upward**2)
    try:
        result = _safe_log(upward, easting, northing, radius)
    except TypeError:
        result = _safe_log(upward, radius)
    return result


spacing = 1e-9
vmax = 1e-5
size = int(2 * vmax / spacing) + 1
easting = np.linspace(-vmax, vmax, size)
northing, upward = 0.0, -50.0

results = np.array([eval_log(e, northing, upward) for e in easting])
results

plt.plot(easting, results)
plt.show()

With the code in main I get the following result:
Figure_4

While with this branch I get:
Figure_3

With this bugfix we solve the discontinuity issue, but we also improve the accuracy of the evaluations of the safe_log.

@santisoler santisoler merged commit 7bb37dc into main Oct 2, 2024
16 checks passed
@santisoler santisoler deleted the bugfix branch October 2, 2024 23:01
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Discontinuity in magnetic field on lines above horizontal edge of prism
1 participant