Prevent frontogenesis from returning nans #3696
Conversation
|
Ideas on how to make my test pass would be greatly appreciated! |
sgdecker
force-pushed
the
front_invest
branch
2 times, most recently
from
a2f713c
to
98d632d
Compare
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.
sgdecker
force-pushed
the
front_invest
branch
from
98d632d
to
12f4285
Compare
| @@ -567,7 +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) | |||
| beta = np.arcsin((-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta) | |||
| arg = (-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta | |||
There was a problem hiding this comment.
One option would be something like this:
| arg = (-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta | |
| arg = (-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) | |
| nonzero_denominator = mag_theta != 0 | |
| arg[nonzero_denominator] /= mag_theta[nonzero_denominator] |
Given how mag_theta is calculated, arg should be already zero where mag_theta is zero.
While searching for the where argument to np.divide, I found a StackOverflow answer suggesting
| arg = (-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta | |
| arg = np.divide(-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)), mag_theta, out=np.zeros_like(mag_theta), where=mag_theta!=0) |
which is a single expresion and seems somewhat more explicit about what it is doing and why.
| beta = np.arcsin(arg) | ||
|
|
||
| return 0.5 * mag_theta * (tdef * np.cos(2 * beta) - div) |
There was a problem hiding this comment.
If I remember the cosine double-angle identities correctly, this should be equivalent:
| beta = np.arcsin(arg) | |
| return 0.5 * mag_theta * (tdef * np.cos(2 * beta) - div) | |
| # cos(2*theta) = cos(theta)**2 - sin(theta)**2 = 1 - 2 * sin(theta)**2 | |
| # second equality uses cos(theta)**2 + sin(theta)**2 = 1 | |
| # sin(arcsin(arg)) == arg for all -1 <= arg <= 1 | |
| return 0.5 * mag_theta * (tdef * (1 - 2 * arg**2) - div) |
The np.clip provides an opportunity to correct roundoff error, but we might be able to skip it and still get decent answers. (Fortunately the errors are near one rather than zero, so going a bit past doesn't introduce sign errors, so the relevant thing is whether this roundoff error is smaller than the discretization error in mag_theta and div).
If you want to declare the trigonometric identities out-of-scope for this PR, that is entirely understandable.
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