API for Coupling with GCM

[citibeth]

We would like to two-way couple PISM to a GCM step by step, compiling the two into one executable. This introduces the following requirements, for which we propose a draft API (almost certainly subject to change as this turns into a working system):

1. It must be possible to timestep PISM explicitly. Rather than calling IceModel::run(), we need a way to step forward to a particular time. We propose:

    /* Steps forward to the time t */
    IceModel::run_to(double t)

1. We need to give our own atmosphere climate to PISM, rather than reading it from any files. We propose ConstSGivenClimate to replace PSGivenClimate. ConstSGivenClimate is essentially a buffer that stores ONE set of SMB and energy fluxes, and then reports that as the flux whenever asked (I believe the asking happens in PISM through the methods ice_surface_mass_flux() and ice_surface_temperature()). No interpolation, averaging or anything else fancy should be needed. The GCM will change the fluxes provided by ConstSGivenClimate at the beginning of each timestep. This will be done via set_surface_mass_flux() and set_surface_temperature() methods.
2. There needs to be a way to initialize PETSc and PISM outside of a main() program. For now, PISM is the only user of PETSc. We believe that no additional functionality is needed here. We intend to construct argv[] and argc variables programmatically, and then send them to PetscInitialize(). We will use information in a related recent issue (#189) to get the communicators set up correctly. There does not seem to be an obvious way to initialize PISM without "pretending" to set up command line options. It feels a bit kludgy. But no matter, the existing ways should work.
3. There needs to be a way to construct the PISM classes DIFFERENTLY from how they are constructed through the normal command-line initialization. In particular, we would like to instantiate an instance of ConstSGivenClimate, and use that in place of PSGivenClimate. We don't know exactly how this is best accomplished, or what ConstSGivenClimate should subclass. We suggest a method IceModel::attach_sclimate(), but that is very much just a guess at how it might work.
4. We will need to replace POGivenClimate as well. But first things first...

Below is a sample of a test program we would envision to write with this API. The program should do approximately the same thing as pismr.cc --- except it would read its forcings out of a file of my own specification. (This file is just for testing, it would be replaced by no file when integrated with the GCM).

class ConstSGivenClimate : public PSModifier
{
public:
    ConstSGivenClimate(IceGrid &g, const NCConfigVariable &conf);

    /** Does nothing */
    void attach_atmosphere_model(PISMAtmosphereModel *input);

    /** Does nothing */
    PetscErrorCode init(PISMVars &vars);

    /** Does nothing */
    PetscErrorCode update(PetscReal my_t, PetscReal my_dt);

    /** ModelE will call this with the SMB it has computed over a
    coupling timestep.  The call will be made from ONLY ONE MPI node,
    the array must be distributed properly according to PISM's domain
    decomposition. */
    void set_surface_mass_flux(blitz::Array<double,2> const &surface);
    void set_surface_temperature(blitz::Array<double,2> const &surface);

    /** Used within PISM to retrieve values stored by
    set_surface_mass_flux(), etc. */
    PetscErrorCode ice_surface_mass_flux(IceModelVec2S &result);
    PetscErrorCode ice_surface_temperature(IceModelVec2S &result);

    // ------------------------------------

private:
    /** Stores the values set by set_xyz() methods above */
    IceModelVec2S surface_mass_flux, surface_temperature;
};

int main(int argc, char *argv[]) {
...
     IceGrid g(com, rank, size, config);
     IceModel m(g, config, overrides);

printf("start = %f\n", g.time->start());
printf("end = %f\n", g.time->end());

     ierr = m.setExecName("pismr"); CHKERRQ(ierr);

     ierr = m.init(); CHKERRQ(ierr);

    std::unique_ptr<ConstSGivenClimate> climate(new ConstSGivenClimate);
    ConstSGivenClimate *pclimate = climate.get();
    m.attach_sclimate(std::move(climate));

    auto nc(NcFile("bob_pism_forcings.nc"));

    blitz:array<double,1> times;    // Will read these from forcing file
    blitz::Array<double,2> mass_flux, temperature;

    m.set_time(times(0));
    for (int i=0; i<times.size(0)-1; ++i) {
        read_next_record(nc, i, mass_flux, temperature);
        pclimate->set_surface_mass_flux(mass_flux);
        pclimate->set_surface_temperature(temperature);
        ierr = m.run_to(times(i+1)); CHKERRQ(ierr);
    }

[ckhroulev]

General comments

First of all: I agree that developing an API will take several iterations, so please go ahead and put together a first approximation of the coupling code. I'll try to help.

One point I cannot overemphasize: if you want me to be able to help, always use the dev branch to base your code on and merge often. (Do not start with stable0.5.) If merging does not seem to work, ask for help. This will almost certainly require less effort than working independently and merging the code later.

I recommend creating a "topical" branch and pushing new code to that branch.

I'm sure I don't need to stress this, but here goes anyway: write test scripts as you go.

Below are my comments regarding individual questions raised in the issue. All these have to do with passing data from the GCM to the ISM.

We need to think about the other direction, too.

"Manual" time-stepping of PISM

I agree that the GCM will have to manage time-stepping and specify the starting time and the duration of the "step". PISM will have to decide whether to take one or multiple time-steps to arrive at the final time, though.

I will add IceModel::run_to(double time_in_seconds).

Providing surface mass balance and surface temperature

Implementing something like the proposed ConstSGivenClimate class should work. The ideal "parent" class for it would be PSConstant, but it was removed from PISM to reduce code duplication. (In other words: we don't need it any more.)

One should be able to copy/paste the implementation of PSConstantPIK to get all the trivial bits.

I will add the IceModel::attach_surface_model(PISMSurfaceModel *) method. (IceModel used to have such a method in the past; it's perfectly OK to put it back in.)

As mentioned in the issue, ConstSGivenClimate will have to scatter SMB and temperature fields according to PISM's grid distribution. See PBLingleClark::transfer_from_proc0(Vec source, IceModelVec2S *result) for the sequence of PETSc calls that performs a scatter like this one.

Initializing PISM outside of main()

I am working of this. See #175, #176, #150, #144.

So far creating argv and argc "programmatically" sounds like a good idea. In the future PISM will take all options and parameters from a provided JSON string. PETSc will always use command-line options, but all PETSc options in PISM will have appropriate prefixes.

Regarding setting up communicators and #189: sets of processors used by PISM and by the GCM can overlap in the setup described in a comment to #189. (I made an assumption that they do not for simplicity.)

Attaching a custom surface model to IceModel

See above.

Attaching a different ocean model

No different from a surface model.

[citibeth]

Constantine,

I believe that last night, I checked in this run_to() functionality on the dev branch. So it would need to be reviewed, but not written from scratch.

https://github.com/pism/pism/commit/e0ef0f33111e51e6dfb5178c9f1217a701e3a745

I'll look at other aspects of your email a little later.

Thank you, -- Bob Fischer

"Manual" time-stepping of PISM I agree that the GCM will have to manage time-stepping and specify the starting time and the duration of the "step". PISM will have to decide whether to take one or multiple time-steps to arrive at the final time, though. I will add IceModel::run_to(double time_in_seconds).

-- Bob

---

From: Constantine Khroulev [notifications@github.com] Sent: Monday, November 04, 2013 10:30 PM To: pism/pism Cc: Fischer, Robert P (GISS-6110)[COLUMBIA UNIVERSITY] Subject: Re: [pism] API for Coupling with GCM (#196)

General comments

First of all: I agree that developing an API will take several iterations, so please go ahead and put together a first approximation of the coupling code. I'll try to help.

One point I cannot overemphasize: if you want me to be able to help, always use the dev branch to base your code on and merge often. (Do not start with stable0.5.) If merging does not seem to work, ask for help. This will almost certainly require less effort than working independently and merging the code later.

I recommend creating a "topical" branch and pushing new code to that branch.

I'm sure I don't need to stress this, but here goes anyway: write test scripts as you go.

Below are my comments regarding individual questions raised in the issue. All these have to do with passing data from the GCM to the ISM.

We need to think about the other direction, too.

"Manual" time-stepping of PISM

I agree that the GCM will have to manage time-stepping and specify the starting time and the duration of the "step". PISM will have to decide whether to take one or multiple time-steps to arrive at the final time, though.

I will add IceModel::run_to(double time_in_seconds).

Providing surface mass balance and surface temperature

Implementing something like the proposed ConstSGivenClimate class should work. The ideal "parent" class for it would be PSConstant, but it was removed from PISM to reduce code duplication. (In other words: we don't need it any more.)

One should be able to copy/paste the implementation of PSConstantPIK to get all the trivial bits.

I will add the IceModel::attach_surface_model(PISMSurfaceModel *) method. (IceModel used to have such a method in the past; it's perfectly OK to put it back in.)

As mentioned in the issue, ConstSGivenClimate will have to scatter SMB and temperature fields according to PISM's grid distribution. See PBLingleClark::transfer_from_proc0(Vec source, IceModelVec2S *result) for the sequence of PETSc calls that performs a scatter like this one.

Initializing PISM outside of main()

I am working of this. See #175https://github.com/pism/pism/issues/175, #176https://github.com/pism/pism/issues/176, #150https://github.com/pism/pism/issues/150, #144https://github.com/pism/pism/issues/144.

So far creating argv and argc "programmatically" sounds like a good idea. In the future PISM will take all options and parameters from a provided JSONhttp://www.json.org string. PETSc will always use command-line options, but all PETSc options in PISM will have appropriate prefixes.

Regarding setting up communicators and #189https://github.com/pism/pism/issues/189: sets of processors used by PISM and by the GCM can overlap in the setup described in a comment to #189https://github.com/pism/pism/issues/189. (I made an assumption that they do not for simplicity.)

Attaching a custom surface model to IceModel

See above.

Attaching a different ocean model

No different from a surface model.

— Reply to this email directly or view it on GitHubhttps://github.com/pism/pism/issues/196#issuecomment-27745172.

[ckhroulev]

Bob,

Yes, I saw that (I get notified when someone pushes to this repository). It looks good to me; I also added a couple of lines of doxygen comments.

Note that (as mentioned above) I also added IceModel::attach_{surface,ocean}_model() to IceModel.

Let me know if this works for you.

[citibeth]

What are the units that PISM expects for SMB and top surface T? Looking at PSConstantPIK, I surmise they are:

ice_surface_temp: K climatic_mass_balance: m s-1

For climatic_mass_balance: is this meters/second of water equivalent? That would be equivalent to 1000 kg/(s m^2)?

[ckhroulev]

ice_surface_temp: K

Yes; all temperatures are in Kelvin.

climatic_mass_balance: m s-1 For climatic_mass_balance: is this meters/second of water equivalent? That would be equivalent to 1000 kg/(s m^2)?

PISM expects PISMSurfaceModel::ice_surface_mass_flux(result) to set result to the SMB in m/second ice equivalent.

A derived class of PISMSurfaceModel will have to do the conversion if necessary.

[citibeth]

It looks like the density of ice is obtained as the "ice_density" configuration parameter:

 PetscReal ice_density = config.get("ice_density")

Thanks, -- Bob

https://github.com/pism/pism/search?q=rho&ref=cmdform https://github.com/pism/pism/search?q=ice_density&ref=cmdform

---

From: Constantine Khroulev [notifications@github.com] Sent: Thursday, November 21, 2013 6:17 PM To: pism/pism Cc: Fischer, Robert P (GISS-6110)[COLUMBIA UNIVERSITY] Subject: Re: [pism] API for Coupling with GCM (#196)

ice_surface_temp: K

Yes; all temperatures are in Kelvin.

climatic_mass_balance: m s-1 For climatic_mass_balance: is this meters/second of water equivalent? That would be equivalent to 1000 kg/(s m^2)?

PISM expects PISMSurfaceModel::ice_surface_mass_flux(result) to set result to the SMB in m/second ice equivalent.

An derived class of PISMSurfaceModel will have to do the conversion if necessary.

— Reply to this email directly or view it on GitHubhttps://github.com/pism/pism/issues/196#issuecomment-29034541.

[ckhroulev]

It looks like the density of ice is obtained as the "ice_density" configuration parameter: PetscReal ice_density = config.get("ice_density")

Yes, this is the way to go in PISM. All physical constants are obtained using the config.get("...") mechanism to ensure consistency across all of PISM.

This does not ensure consistency between PISM and the code it is coupled to. The coupling code might have to use config.set("...") while initializing PISM to take care of this.

Default values of configuration flags and parameters can be found here.

[citibeth]

I need to obtain the following four fields out of PISM. I will root around the code, looking for how to obtain them...

geothermal heat flux internal deformational heat basal runoff elevation changes

[bueler]

Right. Constantine will probably add to this answer, because at least it is likely to be incomplete. But,

• geothermal heat flux: We call it "geothermal flux" because "geothermal heat flux" is redundant. ;-) IceModelVec2S IceModel::vGhf; "upward geothermal flux at bedrock surface"; definition at line 240 of file src/base/iceModel.hh.
• internal deformational heat: Call PISMStressBalance::get_volumetric_strain_heating(), definition at line 97 of file src/base/stressbalance/PISMStressBalance.hh
• basal runoff: Not sure of your definition of the this term, but closest match is probably IceModelVec2S IceModel::basal_melt_rate; "ice basal melt rate from energy conservation and subshelf melt, in ice thickness per time"; definition at line 240 of file src/base/iceModel.hh. There is also class IceModel_grounded_basal_flux_2D_cumulative defined at line 435 of src/base/iceModel_diagnostics.hh
• elevation changes: The surface elevation itself is computed, and you can compute its changes. It is IceModelVec2S IceModel::vh defined at line 240 of src/base/IceModel.hh. [NOTE: Now the variable is called ice_surface_elevation]

Note IceModel::vh has this allocation/set_attrs() call at line 235 of src/base/IceModel.cc:

ierr = vh.create(grid, "usurf", true, WIDE_STENCIL); CHKERRQ(ierr);
ierr = vh.set_attrs("diagnostic", "ice upper surface elevation",
                    "m", "surface_altitude"); CHKERRQ(ierr);
ierr = variables.add(vh); CHKERRQ(ierr);

which illustrates that we have too many names for things:

• Glimmer-matching short names ("usurf"),
• CF-convention standard names for a subset of variables ("surface_altitude"),
• a long name ("ice upper surface elevation"),
• and a historical internal PISM name ("vh"). We know that is a pain! At least we must keep some internal variable name, a standard name, and a long name, so just complying with CF means some of such duplication.

Line numbers above refer to current dev branch.

By the way, I did

$ cd PISM-BUILD-DIR
$ make browser_base
$ chromium-browser doc/browser/base/html/index.html

and looked at the list of IceModelVec data members of IceModel in the browser file:///home/bueler/pism/builddev/doc/browser/base/html/classIceModel.html

Note that one (preferred) way to access IceModel members from PISMComponents is through a PISMVars object; do

$ make browser_util
$ chromium-browser doc/browser/util/html/index.html

and look at file:///home/bueler/pism/builddev/doc/browser/util/html/classPISMVars.html

[bueler]

Regarding "basal runoff": The physics itself is, IMHO, subtle.

There is a many-term energy balance at the base which determines the local rate of conversion of ice to-and-from liquid water, plus there is (lateral) transport of liquid water (in general). For example, there is

• liquid water generated at the ice-rock interface by heating of the ice, which might occur with or without the movement of the liquid water there (e.g. geothermal and friction sources),
• liquid water generated at the ice-subglacialwater interface generated by turbulent dissipation of heat in the moving water there already,
• and water which appears at the subglacier through drainage through the ice from above (including both supraglacial and englacial sources).

Which do you mean? The latter two of the above are certainly "hydrology modeling". Roughly-speaking, the first of the above list is the only contributor to IceModel::basal_melt_rate at this time in PISM. How much of this physics do you want to model?

See subsection 3.4 of A. Aschwanden, E. Bueler, C. Khroulev, H. Blatter (2012) An enthalpy formulation for glaciers and ice sheets. Journal of Glaciology 58 (209) pp. 441–457

[citibeth]

Which do you mean? How much of this physics do you want to model?

I would like to start off simple, I'm sure there are different options that will involve more or less hydrology physics? At the end of the day, I need to know about any water that LEAVES the ice sheet as runoff, so I can give that water to the ocean. Water that's running around above, inside or below the ice, but still under the ice sheet "footprint," I'm happy to let the ice model deal with.

See subsection of FIXME

Not sure what document. I will grep around for it in the documentation source.

---

From: Ed Bueler [notifications@github.com] Sent: Monday, January 06, 2014 7:15 PM To: pism/pism Cc: Fischer, Robert P (GISS-6110)[COLUMBIA UNIVERSITY] Subject: Re: [pism] API for Coupling with GCM (#196)

Regarding "basal runoff": The physics itself is, IMHO, subtle.

There is a many-term energy balance at the base which determines the local rate of conversion of ice to-and-from liquid water, plus there is (lateral) transport of liquid water (in general). For example, there is

• liquid water generated at the ice-rock interface by heating of the ice,
• liquid water generated at the ice-subglacialwater interface generated by turbulent dissipation of heat in the water there already,
• and water which appears at the subglacier through drainage through the ice from above (including both supraglacial and englacial sources).

Which do you mean? How much of this physics do you want to model?

See subsection of FIXME

— Reply to this email directly or view it on GitHubhttps://github.com/pism/pism/issues/196#issuecomment-31701740.

[ckhroulev]

This is ancient and almost certainly no longer relevant. (Coupling is an important topic but PISM's code changed over the last 8 years).