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

Prevent frontogenesis from returning nans #3696

Merged
merged 2 commits into from
Nov 21, 2024

Conversation

sgdecker
Copy link
Contributor

@sgdecker sgdecker commented Nov 19, 2024

Description Of Changes

Previously, frontogenesis could return nans when there was a constant theta field (division by zero would occur) or if round-off error resulted in the argument to arcsin being slightly outside the valid domain of the function (-1 to 1). In this commit, edits are made to set points to zero where nans occur due to division by zero (the frontogenesis is zero when the magnitude of the theta gradient is zero anyway) and to use np.clip to ensure the argument to arcsin is valid.

I could not come up with a simplified test case that triggers the round-off error issue with arcsin, but I do include a test case for the constant theta situation. Because the test case results in a division by zero by design, it is currently failing since that triggers a RuntimeWarning.

Checklist

@sgdecker sgdecker requested a review from a team as a code owner November 19, 2024 15:49
@sgdecker sgdecker requested review from dopplershift and removed request for a team November 19, 2024 15:49
@sgdecker
Copy link
Contributor Author

Ideas on how to make my test pass would be greatly appreciated!

@sgdecker sgdecker force-pushed the front_invest branch 2 times, most recently from a2f713c to 98d632d Compare November 19, 2024 15:59
Fixes #3768 by making sure the argument to the arcsin function is
valid.

Previously, frontogenesis could return nans when there was a constant
theta field (division by zero would occur) or if round-off error
resulted in the argument to arcsin being slightly outside the valid
domain of the function (-1 to 1).  In this commit, edits are made to
set points to zero where nans occur due to division by zero (the
frontogenesis is zero when the magnitude of the theta gradient is zero
anyway) and to use np.clip to ensure the argument to arcsin is valid.

I could not come up with a simplified test case that triggers the
round-off error issue with arcsin, but I do include a test case for
the constant theta situation.  Because the test case results in a
division by zero by design, it is currently failing since that
triggers a RuntimeWarning.
src/metpy/calc/kinematics.py Outdated Show resolved Hide resolved
src/metpy/calc/kinematics.py Outdated Show resolved Hide resolved
Cleaning up numerical issues like division by zero and round-off error
after the fact are kludges, so this commit replaces the previous
approach with more careful numerics.
@dopplershift dopplershift added Type: Enhancement Enhancement to existing functionality Area: Calc Pertains to calculations labels Nov 21, 2024
@dopplershift dopplershift added this to the 1.7.0 milestone Nov 21, 2024
@dopplershift dopplershift merged commit 8251fed into Unidata:main Nov 21, 2024
39 checks passed
@dopplershift
Copy link
Member

Thanks @DWesl for the helpful comments that made the final version of this really nice.

@sgdecker sgdecker deleted the front_invest branch November 21, 2024 16:47
@@ -567,9 +567,18 @@ def frontogenesis(potential_temperature, u, v, dx=None, dy=None, x_dim=-1, y_dim

# Compute the angle (beta) between the wind field and the gradient of potential temperature
psi = 0.5 * np.arctan2(shrd, strd)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you're interested in eliminating all trig:

We have $\tan(2 \psi) =\frac{shrd}{strd}$, which implies $\cos(2\psi)=\frac{strd}{\sqrt{shrd^2+strd^2}}$. We use psi next line to find $\cos(\psi)$ and $\sin(\psi)$.

We show later that $\cos(2\psi) = 1 - 2 \sin^2(\psi) = 2 \cos^2(\psi) - 1$. Rearranging those identities gives
$\cos(\psi) = \sqrt{\frac{1 + \cos(2 \psi)}{2} }$
$\sin(\psi) = \sqrt{\frac{1 - \cos(2 \psi)}{2} } = \sqrt{1-\cos^2(\psi)}$
This lets you replace three trig functions and a multiplication with three or four multiplications, three additions, three square roots, and a division, and four trips through the data with eight. I have no idea whether the eliminating the trig calls makes up for the extra trips through the data.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I ran through that too and wasn't wild about all the square roots. I'm not passionate about trying it but if someone wants to see about performance/accuracy, I'd entertain it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Calc Pertains to calculations Type: Enhancement Enhancement to existing functionality
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Frontogenesis should not return nan
3 participants