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

Added option for curvilinear grid, including rotation angles of grid #22

Merged
merged 8 commits into from
Jan 19, 2017

Conversation

sambarluc
Copy link
Contributor

I added an option to separately deal with Cartesian curvilinear grids. The only difference with standard Cartesian is that sine/cosine defining the grid rotation at each point are now being read.

@codecov-io
Copy link

codecov-io commented Dec 14, 2016

Current coverage is 91.01% (diff: 71.42%)

Merging #22 into master will decrease coverage by 0.29%

@@             master        #22   diff @@
==========================================
  Files             4          4          
  Lines           529        534     +5   
  Methods           0          0          
  Messages          0          0          
  Branches        118        120     +2   
==========================================
+ Hits            483        486     +3   
- Misses           31         32     +1   
- Partials         15         16     +1   

Powered by Codecov. Last update 21be463...6d4c115

Copy link
Member

@rabernat rabernat left a comment

Choose a reason for hiding this comment

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

This is great! Thanks @sambarluc so much for your contribution!

There are two points that would be good to address:

  • Other geometries (e.g. LLC) also have these angle parameters. Should we enable reading of those as well?
  • We will need a test case for curvilinear geometry. This takes the form of a tarred gzipped directory of files in the xmitgcm/tests directory and some modifications to test_mds_store.py. Ideally this will be as small as possible (a 2D case would be great).

filename='AngleCS'),
SN = dict(dims=["j", "i"], attrs=dict(standard_name="Sin of grid orientation angle",
long_name="AngleSN", units=" ", coordinate="YC XC"),
filename='AngleSN')
Copy link
Member

Choose a reason for hiding this comment

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

Is there a reason to call these variables SN and CS rather than just AngleSN and AngleCS?

In the only cases where I change the grid names from their filename, I did it to match what the MITgcm mnc output does.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

SN and CS are the netcdf names, yes, I tried to follow what you had done for the other variables.

Concerning the test, I noticed the complaint in the pull request, but there are no verifications available using curvilinear grid but not LLC (happy to be proven wrong). I will try to make a small test case based on my configuration but I am afraid that will take a little time.

Copy link
Member

Choose a reason for hiding this comment

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

The test example for the llc geometry is described here:
https://github.com/xgcm/xmitgcm/blob/master/xmitgcm/test/test_mds_store.py#L90
(The files are stored in xmitgcm/test/global_oce_llc90.tar.gz)

I call the llc geometry just llc, but perhaps curvilinear would be a better description. Maybe llc means something else, which is related to exch2. In that case, we already have a test case we can use.

Copy link
Member

Choose a reason for hiding this comment

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

SN and CS are the netcdf names, yes, I tried to follow what you had done for the other variables.

perfect

Copy link
Member

Choose a reason for hiding this comment

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

Following up on our earlier discussion, it seems like we should use this pull request to add reading of CS and SN for the LLC geometry as well.

@rabernat
Copy link
Member

Ping! Any progress on this @sambarluc?

@sambarluc
Copy link
Contributor Author

I'm afraid, not at all, and won't make any before next year, I'm very sorry. A very busy end of year between flu, proposals and much needed holidays.

@rabernat
Copy link
Member

I think we should probably replace geometry='llc' with geometry='curvilinear', since the llc geometry is in fact curvilinear. Since llc is not actually a geometry, perhaps we need a separate option for llc, like simply llc=True (which would automatically imply `geometry='curvilinear').

@sambarluc
Copy link
Contributor Author

Hi there, I am finally working on a test for the curvilinear grid. I am struggling with a last failure in pytest. In test_values_and_endianness[curvilinear_leman] I get the assertion error:

E           AssertionError: 
E           Not equal to tolerance rtol=1e-07, atol=0
E           
E           (mismatch 100.0%)
E            x: array(501919.21875, dtype=float32)
E            y: array(501919.0)

test_mds_store.py:376: AssertionError

The value I believe is the correct one is 501919.0. How should I interpret x and y arrays? This should be the check in the endianness, is that 0.22 some sort of rounding error?
(sorry, I know I am very ignorant on these things, please do laugh if I'm saying bullshit)

Thanks, Andrea

@rabernat
Copy link
Member

This is a strange error. Could you just push your current branch so that we can see your full code? It doesn't have to be finished.

You might have to rebase because there have been lots of changes to master recently (e.g. #23, #32).

@sambarluc
Copy link
Contributor Author

This should be it. At least on my machine, one failure as described above.

@codecov-io
Copy link

codecov-io commented Jan 18, 2017

Current coverage is 92.11% (diff: 100%)

Merging #22 into master will increase coverage by 0.80%

@@             master        #22   diff @@
==========================================
  Files             4          4          
  Lines           529        596    +67   
  Methods           0          0          
  Messages          0          0          
  Branches        118        128    +10   
==========================================
+ Hits            483        549    +66   
  Misses           31         31          
- Partials         15         16     +1   

Powered by Codecov. Last update 21be463...72fda6e

'shape': (35, 64, 340),
'test_iternum': 6,
'dtype': np.dtype('f4'),
'expected_values': {'XC': ((0,0), 501919.0)},
Copy link
Member

Choose a reason for hiding this comment

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

This is where the value is set.

@rabernat
Copy link
Member

Ok @sambarluc thanks for this.

Yes, I can see that Travis CI is getting the same error. The error arises due to a mismatch between the value that xmitgcm reads from the file 'XC.data' and the expected value, which you specified to be 501919.0 in the test metadata.

To check this out, I downloaded curvilinear_lmean.tar.gz from your branch and manually examined the file:

>>> data = np.fromfile('XC.data', dtype='>f4')
>>> print('%10f' % data[0])
501919.218750

So it looks like your expected value is actually incorrect.

Curiously, there seems to be some weird behavior involving print. If I just call print with no formatting,

>>> print(data[0])
501919.0

I get the number with no decimals. So this might be the reason you put in that expected value. That is really strange and feels like a bug to me. But there is no doubt that the correct value is 501919.21875.

It seems to print correctly with a slice

>>> print(data[:1])
[ 501919.21875]

@sambarluc
Copy link
Contributor Author

Ok, I fixed this and now it completes the tests successfully. I was tricked by the issue you mention. I really don't get it but at least it's fixed.

@rabernat
Copy link
Member

I wanted to get your opinion on the relationship between curvilinear and llc. If you look at the llc data files, you see that usingCurvilinearGrid=.TRUE.. However, it is a curvilinear grid on the sphere, while yours is cartesian.

Do you think it makes sense to just allow geometry={'cartesian', 'sphericalpolar', 'curvilinear'} and make llc into a separate parameter (related to the presence of multiple facets)? And if so, how do we distinguish from cartesian vs spherical curvilinear?

@sambarluc
Copy link
Contributor Author

I am not familiar with "llc" so the opinion is not based on direct experience. Personally, I see no added value in reproducing the MITgcm grid hierarchy, while having two separate option "llc" and "curvilinear" is in my view easier for the user.

@rabernat
Copy link
Member

Ok, let's leave it as is for now.

One final request before merging: could you edit the documentation to describe the new option for grid geometry?
https://github.com/xgcm/xmitgcm/blob/master/doc/index.rst#grid-geometries

@rabernat
Copy link
Member

Great! Thanks so much for this contribution

@rabernat rabernat merged commit 0a13f62 into MITgcm:master Jan 19, 2017
@sambarluc
Copy link
Contributor Author

Thank you, I'm finding xmitgcm very useful in my work!

@sambarluc sambarluc deleted the add_grid_rotation branch January 20, 2017 07:31
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.

3 participants