Couple C++ SHOC to AD

[PeterCaldwell]

Port shoc_interface.F90 to C++ and couple it to our atmosphere driver and the C++ version of SHOC we've created.

Depends on #831

[tcclevenger]

shoc_main is ported (See https://github.com/E3SM-Project/scream/issues/831). Beginning to work here now. Basing on p3-ad branch: https://github.com/E3SM-Project/scream/tree/aarondonahue/p3_interface_ad_cxx_redo

A few comments:

1. Unlike p3, the Fortran impl of SHOC in the AD seems incomplete. For instance, there is no call to shoc_main() in shoc_main_f90() in scream_shoc_interface.F90. This means that in the shoc stand alone test and shoc-p3 coupled test never actually call shoc_main(). Fortran is getting replaced, but the tests need to be completed and input files (I believe) need to be created and parsed correctly.
2. When calling SHOCMacrophysics::set_grids(), I (think) I will need a way to access num_tracers for the input/output fields. Not sure if that should be stored in the GridManager. (like with columns and levels) or hardcoded or elsewhere..
3. shoc_init(), which I believe only calculates npbl, needs to be implemented.

[AaronDonahue]

"When calling SHOCMacrophysics::set_grids(), I (think) I will need a way to access num_tracers for the input/output fields. Not sure if that should be stored in the GridManager. (like with columns and levels) or hardcoded or elsewhere.."

@tcclevenger , in a departure from EAM and the F90 implementation, SCREAM will define each tracer specifically instead of keep the Q array. So presumably we won't have to worry about something like num_tracers. That being said, if for some reason SHOC needs to see the tracers as a monolithic array then it's possible that information could be gathered specifically for SHOC.

[tcclevenger]

@tcclevenger , in a departure from EAM and the F90 implementation, SCREAM will define each tracer specifically instead of keep the Q array. So presumably we won't have to worry about something like num_tracers. That being said, if for some reason SHOC needs to see the tracers as a monolithic array then it's possible that information could be gathered specifically for SHOC.

@AaronDonahue Currently shoc_main() (C++ and Fortran) takes in a 3d view of size (shcol, nlev, num_tracers). Will this change? Are the total number of tracers supposed to be constant?

[AaronDonahue]

I can't speak to SHOC in particular since I haven't looked at that code directly. I guess the main question is, could the interface be changed to take those tracers that SHOC needs? This might be a good question for @bogensch.

To answer your question, in general in EAM there is a maximum tracer size that is defined at compile-time, but on a smaller scale physics keeps track of the total number of tracers that have been registered at runtime. In short, any process can add a new tracer. In SCREAM right now that happens in p3 and shoc, and this is done using cnst_add. Everytime cnst_add is issued the total number of tracers is bumped up by 1.

In SCREAM we want to change this approach. Instead of having a 3D variable state%q we want each tracer to be it's own individual 2D tracer. So processes can explicitly request the tracers they care about instead of having to take the full q variable and only use a subset of them.

[bogensch]

I would totally vote for SHOC to keep the tracer array as a monolithic array, as they are all transported the same and would require substantial modifications to the code that I see as unnecessary (unless I'm missing something here).

The tracers that need to be transport should be gathered in the qtracer array in the interface. This is already done in SCREAMv0 as the num_qtracers does not necessarily equal "pcnst" and only the tracers we want SHOC to be transported are fed into this array ("edsclr_in" in the inteface).

[AaronDonahue]

@bogensch , thanks for responding. So it looks like SHOC does have a treatment for a monolithic q array. I wasn't sure if q was passed as input and then broken into it's individual parts within the main body of code. I think this might be a good thing to get comments on from others in the group. My initial thought would be to do something like what @jeff-cohere is working on for remapping to the dynamics, which would be to gather all the "tracers" into an array locally. Instead of maintaining a 3D q array in the field manager.

@tcclevenger , what that would look like in implementation would be a local 3D view in run_impl that is populated with all of the 2D views of each tracer. It would be a bit more complicated than that, but that's the gist.

Alternatively, @jeff-cohere and @bartgol is there any merit to having the field manager have its own get_tracer_view command that would grab all of the 2d tracer fields and assign them to a 3d view? Follow-up, would any changes to values in the 3d view automatically be reflected in the 2d fields themselves, in other words, could we have the 3d view point to the memory addresses of the 2d views? Maybe as an array of 2d views, or an array of pointers? Just thinking outside of the box here.

[bartgol]

This complicates things. I suggest we discuss this at the v1 telecon. It is of course possible to have a monolithic Q in the field manager, and have SHOC use that, but there are a few problems:

• the dimension of Q is not known until all tracers have been declared. One could do what homme does, and use a compile-time upper bound.
• we'd need a way to pass to all processes an information regarding the idx of each tracer in the Q array.

• from what @tcclevenger said, it seems that SHOC expects tracers shaped as (ncols,nlevs,ntracers), which is not what we usually do for physics layouts (we usually put nlevs as the fastest striding dim).

  In general, this peculiar way of treating tracers (i.e., I use only a few, but I want the whole big array) is going to require some changes to the design. I would vote against keeping both approaches (single Q array vs separate tracers), since that greatly complicates things. On the other hand, if we decide to go down the route of a single Q array, then the phys-dyn tracers remapping task @jeff-cohere is working on might no longer be needed.

From the design point of view, I'm not a fan of a monolithic Q array, for a couple of reasons:

• it hides the info of who updates what. With a single Q field, every process will likely update Q, but who updates Q1,...Qn would be harder to track).
• it is a big bucket to pass around (somewhat like pbuf), which hides actual context.

• it creates horizontal dependency between atm process. By this I mean that, as @tcclevenger pointed out in his original comment, one needs to know how many tracers there are in order to create Q; but that number depends on what other atm process need.

  Long story short: anything can be done, so if the consensus is that a big Q is better, we'll work on it. But I want to make sure that there is a reason for that other than "that's what it is right now". Scream is doing things from scratch, so it can be a chance to change things that have been kept the way they are for lack of resources/will to change/upgrade them.

[jeff-cohere]

One additional wrinkle here is that the "3D" tracer array is laid out differently in HOMME than in SHOC, for performance reasons. So the representation of tracers is likely to be different between those two packages anyway, unless we can figure out a single representation for both. Maybe I'm missing something.

I'm with @bartgol that treating the tracers as a big monolith adds a lot of hidden assumptions and requires more expertise to manipulate. I agree we should discuss this at the telecon.

[bartgol]

I believe this is completed, so I'm closing. @tcclevenger reopen if you think this is not yet done.