Mdsio section

[timothyas]

What changes does this PR introduce?

(Bug fix, feature, docs update, ...)

The current mdsio read/write routines for YZ or XZ slices don't take useSingleCPUIO into account, and also mds_write_sec_[xz/yz] do not write meta files. This PR makes some additions so that both of these are addressed: the mdsio read/write routines follow the same single cpu io logic as mds_write_field, and prints meta files accordingly.

Perhaps the first issue, whether mds_[read/write]_sec uses singlecpuio logic or not, does not matter. I found it to be annoying that this was turned on, yet all the xx_obcs files always print out as tiled files (since pkg/autodiff always uses non global files). I dug into addressing this (which turned out to be more tedious than expected...) and figured I would open up the PR to see what people think. I think that at least writing meta files would be a useful addition, although relatively minor.

What is the current behaviour?

(You can also link to an open issue here)

It does not matter if singlecpuio is turned on for sliced fields. The time where I noticed this is with the autodiff+obcs+ctrl packages in use.

What is the new behaviour

(if this is a feature change)?

Now if singlecpuio is turned on, xz or yz slices are read/written with the master process via some new scatter/gather routines.

Does this PR introduce a breaking change?

(What changes might users need to make in their application due to this PR?)

Honestly, it might. I did not wade through the exch2 logic dutifully since it is very confusing and I'm not familiar with it. At least something that needs to change is the size of the x_buffer and y_buffer fields in mdsio_[read/write]_section, and the indexing that happens in the [scatter/gather]_[xz/yz] routines could be incorrect. However, I wanted to see if anyone thinks this is worth keeping before digging into that...

Another catch is with the gather_[xz/yz] routines. In this routine it's not clear which processor has the sliced information. For instance, if the slice is on one of the boundaries, then it would live on all the processors (and tiles) that intersect with the boundary. This would be difficult and perhaps overly intrusive to determine or pass as an input to the gather routine, so I solved this by simply only passing values to the global buffer if they are nonzero. Presumably, the sliced data are zero on all processors and tiles except where that slice is "active", and nonzero only where it is active, and the gather routine would then successfully grab the nonzero information. This was just a first pass so could be faulty (and all of this not worth it at all), but I'd be curious what the gurus think...

Other information:

Suggested addition to tag-index

(To avoid unnecessary merge conflicts, please don't update tag-index. One of the admins will do that when merging your pull request.)

o mdsio_[read/write]section rewritten to include useSingleCPUIO logic, mirroring mdsio[read/write]field
o [scatter/gather][xz/yz]_[r4/r8] routines added

[timothyas]

I have to say I don't really understand why the tutorial_plume_on_slope check is failing. I ran the testreport on the cluster at our institute and without turning on mpi it passed, with output matching to at least 16 digits of accuracy. I believe the verification results are without mpi, but running with mpi I get the same result as travis, where it simply says that the run failed... But even running with mpi on our cluster, diffing the STDOUT.0000 and results/output.txt files shows no differences in the monitor output. So I'm not really sure where I've gone wrong.

[mjlosch]

@timothyas First of all, I like your effort. I too am a little annoyed by the multiple files and the missing meta information.

I am not too much of an MDSIO expert but I can reproduce the failing experiment on my laptop. The numbers are OK (I guess), but it fails at the end in MDS_WRITE_SEC_YZ, because after calling MDS_PASS_YZ_R8toRL, the variable nNz is suddenly zero and then the length_of_rec = 0 and corresponding open statement returns an error ("Fortran runtime error: RECL parameter is non-positive in OPEN statement"). It is totally unclear to me, why nNz is zero, but it happens in or after line 297: https://github.com/timothyas/MITgcm/blob/f3530ccf6a512f61194fd16fcdac79ad8b19937b/pkg/mdsio/mdsio_pass_r8torl.F#L297
(by primitive "hello debugging"). I have no idea why ... but I am guessing it has to do with inconsistent field sizes? E.g. r8seg is defined with size (sNy,nSx,nSy), but the corresponding field in S/R MDS_PASS_YZ_R8toRL is buffer(1-oL:sNy+oL,nNz,nSx,nSy) and nNz > 1 on entry?

[mjlosch]

@timothyas I think I now have a fix. Would you like me to check it in for you to see? Not sure that it does the right thing (whether the pickup file is correct, contains only zeros), but at least the test does not break.

[timothyas]

Hey @mjlosch, thanks so much for the feedback and helpful debugging! Yes I'd be interested to see your fix. I'm re-wrapping my head around all these buffers, segs, and pass routines... Now that you mention it, this does look suspect. Here I was mirroring (or, trying to mirror) the style of mdsio_write_field, which calls mds_pass_r8torl with sharedLocBuf(1:sNx,1:sNy,nSx,nSy) in the place of buffer(1-oLi:sNx+oLi,1-oLj:sNy+oLj,nNz,nSx,nSy) (using oLi=oLj=0, just like I was setting oL=0 for the YZ/XZ routines). It seemed like passing 1 to nNz would make all of this kosher but perhaps not. All that to say, I'm looking forward to seeing your fix :)

[mjlosch]

Hi Tim, just added my fix, but be aware that this only fixes the runtime error. The pickup_orlanskiE files have 26 records of 60x1 in them now, but with useSingleCpuIO = .TRUE. they are non-zero, with .FALSE. (default of this experiment), they are all zero. I am not sure what to expect here.

[timothyas]

Thanks @mjlosch this is super helpful. I bet this has to do with the corresponding read routines, I'll check it out

[mjlosch]

@timothyas you are right. With useSingleCpuIO=.FALSE., there are 4 orlanski-pickup files, and the fourth (pickup_orlanskiE.ckptA.004.001.data) contains the correct data. I guess "rdmds" fails because it does not know how to concatenate 4 fields of the same shape, when there should only be one. Not sure how to fix that, but that may not be necessary anyway ...

[timothyas]

Hey @mjlosch and @jm-c I think this is now in very good shape. I have made some fixes related to exch2 and tested all possible options (SingleCPUIO = T/F, globalFiles =T, with/without MPI, with or without exch2) and in every case I'm getting tests passing to 16 digits, same as master.

Additionally, @antnguyen13 gave me a lightweight ASTE setup to test with. This has been a good test because it uses a variety of options: obcs, exch2, and obcs loaded via the exf package (i.e. exf_set_obcs). Also seaice is in there, it's just a big ol' party. I hooked this up with testreport and get the comparison shown below, where the only difference is this PR, which actually uses SingleCPUIO (and thus MPI scatter) to read in the OBCS fields at each time step. In the current master branch, these are read in as global files.

Some seaice related details on the ASTE setup that I think @mjlosch will find interesting (big thanks to @antnguyen13 for some coaching here):

• with SEAICEscaleSurfStress=.TRUE., something weird happens after 16 iterations and differences in the ~10-11 decimal place of uice/vice turn into differences in the 3rd/4th place
• with this parameter set to .FALSE. I can safely make it through 20 time steps with the output shown below

However, almost no matter what the options were, I could always get through 10 iterations with output similar to what is shown below (i.e. passing with pretty good agreement). It seems like these stability issues are a matter of seaice being a tricky beast, rather than MDSIO problems, so I feel pretty confident about the code.

I'd be curious for any feedback and I'm ready to do what I can to help get this merged.

"miniASTE" testreport output after 20 iterations (note it reads in single precision OBCS fields):

There were 10 decimal places of similarity for "Press. Solver (cg2d)"
There were 16 decimal places of similarity for "Theta minimum"
There were 16 decimal places of similarity for "Theta maximum"
There were 16 decimal places of similarity for "Theta mean"
There were 16 decimal places of similarity for "Theta Std.Dev"
There were 16 decimal places of similarity for "Salt minimum"
There were 14 decimal places of similarity for "Salt maximum"
There were 16 decimal places of similarity for "Salt mean"
There were 14 decimal places of similarity for "Salt Std.Dev"
There were 14 decimal places of similarity for "U minimum"
There were 13 decimal places of similarity for "U maximum"
There were 10 decimal places of similarity for "U mean"
There were 11 decimal places of similarity for "U Std.Dev"
There were 13 decimal places of similarity for "V minimum"
There were 13 decimal places of similarity for "V maximum"
There were 10 decimal places of similarity for "V mean"
There were 13 decimal places of similarity for "V Std.Dev"
There were 12 decimal places of similarity for "Eta minimum"
There were 14 decimal places of similarity for "Eta maximum"
There were 13 decimal places of similarity for "Eta mean"
There were 13 decimal places of similarity for "Eta Std.Dev"
There were 22 decimal places of similarity for "SIce Area min"
There were 16 decimal places of similarity for "SIce Area max"
There were 13 decimal places of similarity for "SIce Area mean"
There were 14 decimal places of similarity for "SIce Area StDv"
There were 22 decimal places of similarity for "SIce Heff min"
There were 13 decimal places of similarity for "SIce Heff max"
There were 16 decimal places of similarity for "SIce Heff mean"
There were 13 decimal places of similarity for "SIce Heff StDv"
There were 13 decimal places of similarity for "SIce Uice min"
There were 13 decimal places of similarity for "SIce Uice max"
There were 10 decimal places of similarity for "SIce Uice mean"
There were 11 decimal places of similarity for "SIce Uice StDv"
There were 13 decimal places of similarity for "SIce Vice min"
There were 13 decimal places of similarity for "SIce Vice max"
There were 9 decimal places of similarity for "SIce Vice mean"
There were 11 decimal places of similarity for "SIce Vice StDv"

[jm-c]

@timothyas I did not check all what I wanted to but have few comments:

1. looking at mdsio_read_section.F, there is likely a problem with W2_useE2ioLayOut and multiple facets (like, e.g., full 6 facet cube), because either exch2_xStack_Nx > Nx or exch2_yStack_Ny > Ny so that an out-of-bounds index problem is expected either in MDS_READ_SEC_XZ, lines 152-155 or in MDS_READ_SEC_YZ, lines 712-715. The problem here is not
   so much on your side but the fact that we don't have any experiment that test OBCS with pkg/exch2 + W2_useE2ioLayOut + multiple facets. @jahn do you still have a set-up like this to test this ?
   I think it's an important feature and we want to keep it working.
2. trying to read/write local array (i.e., not in common block) sections in multi-threaded run is likely not going to work since some buffers are not in common block but should be (like in 2d & 3d I/O routines). Note that for x-z and y-z sections this is
   not currently supported (and only works in multi-threaded for variables in common block). may be we could leave this multi-threads issue for later since it's not the most important right now and see how to fix it a bit later.
3. a suggestion: for "slice", we could consider using a single character string in N,S,E,W instead of a 5-character one.
   From my point of view, this would have 2 advantages:

1. simplify some pieces of code, e.g.:
   a) pkg/autodiff/active_file_gen.F, lines 88-89
   b) exf/exf_set_obcs.F, lines 71-72
2. less prone to coding mistake, e.g., using 'east' or 'west' (4 characters only) instead of 'east ' (with a blank) and 'west ' (with a blank, 5 characters in total).

I will continue to review the other bits but it's taking me some time ...

[timothyas]

Thanks so much for the feedback @jm-c. I agree with your points and I will work on your list in reverse order, starting with 3. Hopefully this will simplify things a bit.

[jm-c]

@timothyas I wrote a new issue #753 that proposes an alternative pkg/rw interface. Please take a look and let us know what you think.

Note that it we go this route (which would imply also postponing to a future PR many of the changes to S/R calls) the list of files with potential merge-conflicts would go down to a single one (mdsio_read_section.F).