New differential operators: Laplacian, grad, div, and curl #293
Conversation
|
I'm not sure why the test action is failing for Python 3.7, it works for me at Python 3.8. Any ideas? |
As you likely saw it was one failure due to a partitioning reshape issue, which seems funky. I also do not know immediately what might be up. I have re-run the Actions jobs to see if we get the same result again, in case there is some test flakiness involved (also not good and would need fixing, of course)... |
|
After re-running, the funky test failure has not emerged, but of the two jobs that ran before cancelling the others, both Python 3.9, the Ubuntu passed but the Mac OS failed on some some genuine failures relating to this PR: I am not immediately sure why these are apparently failing for only Mac OS - something you might want to ponder whist I review the PR? |
|
Thanks, Sadie - I'll investigate the new error. I'll probably commit some print statements, just for the purpose of diagnosis. |
Codecov Report
@@ Coverage Diff @@
## master #293 +/- ##
==========================================
- Coverage 73.80% 71.59% -2.20%
==========================================
Files 85 88 +3
Lines 19179 19778 +599
==========================================
+ Hits 14153 14159 +6
- Misses 5026 5619 +593
Flags with carried forward coverage won't be shown. Click here to find out more.
Continue to review full report at Codecov.
|
|
Hi Sadie - I give up. I can't fathom the automated tests - they seem either pass or time out, now. So I'm putting this up for review, as it seems OK, as best I can make it. Thanks. |
|
Thanks David, I have noticed the automated test jobs hanging on other PRs so I doubt it is an issue stemming from this PR. Sorry for the trouble, which I will investigate shortly, though I will review this first... |
There was a problem hiding this comment.
I'm mid-review, now reviewing the details of the method implementations, but have checked the tests and the new documentation (with much consultation of Wolfram MathWorld and working around their usage of a different angle symbol convention for spherical polar coordinates!), both of which look great and more than suitably comprehensive.
Just noting some initial comments and suggestions to consider, in the meantime, including the following general comment...
The unit tests look very good, and are certainly sufficient for testing, though I was wondering if we could test on some simple vector calculus identities to act as a basic integration test and further validate some, if not all of, the new operators, as that would be simple and quick to add.
For example whilst playing around with the new methods I was checking, based on two of the identities listed here, whether the divergence of the curl of some (of any) field ended up being close enough (considering floating point calculations) to zero, and likewise with the curl of the gradient, like so:
In [64]: import cf
...: a = cf.example_field(0)
...:
...: a_x, a_y = a.grad_xy(radius="earth", one_sided_at_boundary=True)
...: c = cf.curl_xy(a_x, a_y, radius="earth", one_sided_at_boundary=True)
...: d = cf.div_xy(a_x, a_y, radius="earth", one_sided_at_boundary=True)
...:
...: dc = cf.div_xy(c, c, radius="earth", one_sided_at_boundary=True)
...: print(dc.identity()) # check name reflects context
...: print(dc.array) # should be zero within calc approximation
...: zeros = dc.copy()
...: zeros[...] = 0
...: print(cf.equals(dc, zeros, verbose=3))
...:
...: cg = cf.curl_xy(a_x, a_y, radius="earth", one_sided_at_boundary=True)
...: print(cg.identity()) # check name reflects context
...: print(cg.array) # should be zero within calc approximation
...: zeros = cg.copy()
...: zeros[...] = 0
...: print(cf.equals(cg, zeros, verbose=3))
long_name=Horizontal divergence of (long_name=Horizontal curl of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity), long_name=Horizontal curl of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity))
[[ 3.00738521e-37 -7.51846302e-38 0.00000000e+00 4.94399257e-37
8.17934420e-37 -1.25764388e-36 -5.17195900e-37 6.88059993e-37]
[-6.48640473e-38 -2.71237017e-38 5.93601504e-38 2.20155877e-38
-5.81820676e-38 -6.76206288e-38 8.29496038e-38 5.34651036e-38]
[-5.35128390e-38 1.96619143e-38 -8.10800592e-39 -2.10808154e-38
1.62160118e-38 1.13512083e-38 -1.94592142e-38 -1.29728095e-38]
[-4.81471642e-38 6.09340161e-38 -2.73335270e-38 5.18912379e-38
1.42671699e-38 -1.53289750e-38 3.10583659e-38 -3.14495159e-38]
[ 1.42683698e-37 6.13375424e-37 -1.42683698e-37 -4.22686361e-38
1.21925124e-38 8.05513260e-37 6.37863089e-38 -8.71473313e-37]]
True
long_name=Horizontal curl of (long_name=X gradient of specific_humidity, long_name=Y gradient of specific_humidity)
[[ 0.00000000e+00 -5.01626996e-31 0.00000000e+00 -5.01626996e-31
5.01626996e-31 2.00650799e-30 -1.00325399e-30 -5.01626996e-31]
[ 3.67216448e-31 -2.75412336e-31 0.00000000e+00 1.83608224e-31
0.00000000e+00 -9.18041120e-32 1.83608224e-31 1.83608224e-31]
[ 0.00000000e+00 1.62288275e-32 -3.24576551e-32 0.00000000e+00
1.62288275e-32 -3.24576551e-32 0.00000000e+00 3.24576551e-32]
[-3.67216448e-31 4.87709345e-32 -9.18041120e-32 -1.83608224e-31
2.75412336e-31 9.18041120e-32 -1.83608224e-31 0.00000000e+00]
[-1.25406749e-30 5.64330371e-31 0.00000000e+00 0.00000000e+00
-4.98977600e-31 1.00325399e-30 1.00325399e-30 0.00000000e+00]]
TrueSo my suggestion is a small test method called test_diff_operators or similar which tests some basic combinations, like the above (not essential, but could be useful).
cf/field.py
Outdated
| """ | ||
| from numpy import pi | ||
|
|
||
| f = _inplace_enabled_define_and_cleanup(self) |
There was a problem hiding this comment.
In-place Laplacian calculation doesn't seem to work properly, though I can't immediately see why since the standard in-place logic and decorator are all there - any ideas what is going on? I see:
In [15]: import cf
...: a = cf.example_field(0)
...: print("A IS", a, a.array)
...:
...: b = a.laplacian_xy(radius="earth")
...: print("B IS", b, b.array)
...:
...: c = a.laplacian_xy(radius="earth", inplace=True)
...: print("IN-PLACE A IS", a, a.array, "C IS", c)
...: cf.equals(a, b, verbose=3)
A IS Field: specific_humidity (ncvar%q)
----------------------------------
Data : specific_humidity(latitude(5), longitude(8)) 1
Cell methods : area: mean
Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
: longitude(8) = [22.5, ..., 337.5] degrees_east
: time(1) = [2019-01-01 00:00:00] [[0.007 0.034 0.003 0.014 0.018 0.037 0.024 0.029]
[0.023 0.036 0.045 0.062 0.046 0.073 0.006 0.066]
[0.11 0.131 0.124 0.146 0.087 0.103 0.057 0.011]
[0.029 0.059 0.039 0.07 0.058 0.072 0.009 0.017]
[0.006 0.036 0.019 0.035 0.018 0.037 0.034 0.013]]
B IS Field: long_name=Horizontal Laplacian of specific_humidity (ncvar%q)
--------------------------------------------------------------------
Data : long_name=Horizontal Laplacian of specific_humidity(latitude(5), longitude(8)) m-2.rad-2
Cell methods : area: mean
Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
: longitude(8) = [22.5, ..., 337.5] degrees_east
: time(1) = [2019-01-01 00:00:00] [[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-2.143064530824368e-15 -2.6749442771355633e-15 -2.4238409002370553e-15
-2.6377524366954036e-15 -1.0992306748059246e-15 -1.6075152632499628e-15
3.5426406417251047e-16 2.2860907541031372e-15]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]]
IN-PLACE A IS Field: specific_humidity (ncvar%q)
----------------------------------
Data : specific_humidity(latitude(5), longitude(8)) 1
Cell methods : area: mean
Dimension coords: latitude(5) = [-1.3089969389957472, ..., 1.3089969389957472] radians
: longitude(8) = [0.39269908169872414, ..., 5.8904862254808625] radians
: time(1) = [2019-01-01 00:00:00] [[0.007 0.034 0.003 0.014 0.018 0.037 0.024 0.029]
[0.023 0.036 0.045 0.062 0.046 0.073 0.006 0.066]
[0.11 0.131 0.124 0.146 0.087 0.103 0.057 0.011]
[0.029 0.059 0.039 0.07 0.058 0.072 0.009 0.017]
[0.006 0.036 0.019 0.035 0.018 0.037 0.034 0.013]] C IS None
Field: Different Units: <Units: 1> != <Units: m-2.rad-2>
Out[15]: FalseThere was a problem hiding this comment.
Ah - got it. The culprit is that f is reassigned:
f = _inplace_enabled_define_and_cleanup(self)
<snip>
f = d2f_dx2 + d2f_dy2
return fI'll fix that and add the test that caught it.
There was a problem hiding this comment.
Hmm. Weird things started happening with NaNs when it was properly in-place. I think the answer is to remove inplace from laplacian_xy for now, and see if anyone ever notices (there would be no real performance benefit from it, in this case, anyway)
There was a problem hiding this comment.
(Also, none of the other new methods/function have inplace, so it sort of consistent)
There was a problem hiding this comment.
I think the answer is to remove inplace from laplacian_xy for now, and see if anyone ever notices (there would be no real performance benefit from it, in this case, anyway)
(Also, none of the other new methods/function have inplace, so it sort of consistent)
Sounds good to me! If anyone asks it is for consistency...
Co-authored-by: Sadie L. Bartholomew <[email protected]>
…l/cf-python into differential-operators
|
New test |
|
Great - all done, then. Thanks for your review that has improved things a lot. |
Fixes #292
Fixes #178