Post-215 work: add support for retrieving liquid mass/saturation values directly from a model

[jeff-cohere]

These functions have been discarded as part of #215, since they're specific to the discretization. We can resurrect them in a more specific form if we can figure out how to articulate what kind of space these values live in (e.g. a finite element function space, or "on cells" in a finite volume method).

TDyGetLiquidMassValuesLocal
TDyGetSaturationValuesLocal

[jeff-cohere]

Hey @bishtgautam,

Do all of our "models" use liquid mass and saturation? If so, I think we can add these functions in the manner of PETSc's VecGetArray and VecRestoreArray methods. E.g.

• TDyGetLiquidMassArray: sets a pointer to the array containing the liquid mass values for the model, however it's represented
• TDyRestoreLiquidMassArray: unsets the pointer set by TDyGetLiquidMassArray

and the same for saturation. Does this make sense? If there are models that don't use liquid mass and saturation, we might have to use a more elaborate name.

[bishtgautam]

All models will have both of those variables.

[jeff-cohere]

Do you like this solution? If so, I can put together a PR. It will add one function to the set of virtual functions for our model interface.

[bishtgautam]

The PETSc-type approach looks good to me. For the application code, I think we should only expose TDyGetLiquidMassArrayRead (similar to VecGetArrayRead), so the application code does not modify the values.

[jeff-cohere]

Good point. We can add the const attribute to the pointer in the function to make it like VecGetArrayRead rather than VecGetArray. Do you think we should call it TDyGetLiquidMassArrayRead instead of just TDyGetLiquidMassArray, or do you think it's sufficient to say (in the comments describing it) that it's like VecGetArrayRead?

[jeff-cohere]

There's one issue with this approach: by allowing direct access to an array of liquid mass and saturation values, we require every implementation to store these things in their own arrays, which might not be necessary all the time.

For example, one calculates the liquid mass in each cell (for the MPFA-O method) by computing the product of the water density and the cell volume. This doesn't require the MPFA-O implementation(s) to store an extra array--we could actually compute the values on request if the user provides an array of the right length.

If we went this way, the interface would look like this:

• TDyGetLiquidMass(TDy tdy, PetscInt n, PetscReal liquid_mass[n]): computes the liquid mass, placing the results into the given array (of length n).
• TDyGetSaturation(TDy tdy, PetscInt n, PetscReal saturation[n]): computes the saturation, placing the results into the given array (of length n).

This might be simpler, since it lets different implementations decide how they want to provide the requested information.

[jedbrown]

The FE methods we're considering have some form of constant pressure so these are "per cell" similar to FV.

What does the caller need from these? A Get/Restore pair providing a Vec would have the TDy keep ownership, though it can allocate on demand. This has a potential advantage that the Vec can live on a GPU without changing the interface, and the caller can decide whether they need to access it on the CPU or GPU. Alternatively, the caller can provide a Vec to place the result in, in which case we check that it's the correct size and will write into the Vec on CPU or GPU memory depending on the preference of the Vec.

I imagine we'll want to plot these fields, in which case making it a field over an auxiliary DM would be desirable. (Then VecView(X, viewer) gives you many ways to visualize.) Cc: @knepley

[jeff-cohere]

As far as I know, these are diagnostic variables used to evaluate the quality of the solution in demo programs. So an interface using a Vec is probably fine, and probably at a better level of abstraction anyway. I'll work something up.

[jedbrown]

If it's really about a human monitoring a simulation, then it can be a TS (and SNES?) monitor that writes summary stats and/or fields to a file. There could be several such fields written by one monitor rather than a separate monitor per field (which, for many formats, will imply a separate file with separate mesh connectivity for each field).

[jeff-cohere]

Well, I suppose these variables might also be needed by (say) E3SM's coupler, so it's possible that monitors might not necessarily be the best fit. But @bishtgautam probably knows more.

[bishtgautam]

Certain fields (e.g. mass, saturation, temperature, etc) will be needed by ELM. I'm fine with either approach of

1. ELM/demo providing a Vec when calling TDy to get liquid mass, or
2. ELM/demo providing a pointer (like in get/resotre functions for Vec)

• If we go with the second approach, we should add the suffix Read as it will be consistent with the PETSc terminology.
• If we go with the first approach, we don't need to add the suffix Read because it will be a TDy-specific terminology.

[jedbrown]

Although I don't expect this to be a bottleneck, populating all of these fields in one shot is likely to be more efficient than one at a time. One could make an auxiliary DM and add the desired fields, then TDy would examine and populate the desired fields.

[jeff-cohere]

Thanks for the discussion, guys. I'll come back with a couple of interface ideas, as seen from the perspective of the ELM: a Vec-based interface and a DM-based one, so we can see how they look.

[jeff-cohere]

@jedbrown , a question about your "auxiliary DM" idea: What's the easiest way to create one of these from an existing DMPlex with the minimum needed information?

For example, I see that the structured grid implementations of DM can be "compatible" if they represent the same domain and have the same parallel layout, and there's a way to construct a "compatible" structured DM from another one. But we're using DMPlex, which doesn't seem to have this concept of "compatibility", at least not yet. Is there a way to set up an "auxiliary DM" that relates to the dycore's DM(Plex) in such a way that it can properly manage fields but doesn't result in all the unstructured mesh data being duplicated? That is: if all we care about is storing a set of fields that are defined properly on our domain, is there a "minimal" auxiliary DM that we can create using our DMPlex as a model?

[jedbrown]

DMClone can be used to make a new DM with the same topology, but new fields. For example, see dmAux usage: https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex12.c#L760-765

[jeff-cohere]

Okay, how about this:

• PetscErrorCode TDyCreateDiagnostics(TDy tdy, DM *diagnostics): creates an auxiliary DM (using DMClone under the hood) with diagnostic fields defined, based on the selected model and discretization.
• PetscErrorCode TDyComputeDiagnostics(TDy tdy, DM diagnostics): given an auxiliary DM created by TDyCreateDiagnostics, this function populates the diagnostic fields present with up-to-date values based on the dycore's state.

Then the demo in question (or ELM) can access the diagnostic fields with DMGetNumFields and DMGetField. I admit I'm not very familiar with PETSc's discretization objects, but this seems like a nice abstraction in the sense that it doesn't leave much room for error.

If we wanted to use a Vec interface instead, we could offer an interface to query the number and names of diagnostic fields, and an interface to populate a set of Vecs for those fields. But the DM approach guarantees that the bookkeeping is done properly, so if we like the above interface, let's give it a try.

[jedbrown]

TDyComputeDiagnostics(TDy, DM, Vec) because DM is like a "function space" (plus optional specification of a problem) so we still need a Vec to store the data (an element of the function space).

[jeff-cohere]

Okay. And the single Vec stores all the diagnostics using an indexing scheme that is defined by the DM, right? Is there a nice one-liner for creating this Vec from the DM?

[knepley]

DMCreateGlobalVector()

[jeff-cohere]

Thanks, Jed and Matt. @bishtgautam , I'll put together a PR that illustrates how this might work in the demo with the old stuff commented out.

[bishtgautam]

So, the Vec created by DMCreateGlobalVector() would be a vector with a stride corresponding to the following variables?

https://github.com/TDycores-Project/TDycore/blob/master/include/private/tdympfaoimpl.h#L55-L57 https://github.com/TDycores-Project/TDycore/blob/master/include/private/tdympfaoimpl.h#L72-L76

[jedbrown]

Yeah, it's one value per field per cell and can either use VecStride or DMPlex depending on how you want to access subsets of the fields. VecView can be used to write a file for visualization that will include all the named fields.

[jeff-cohere]

Hey @jedbrown . Suppose I wanted to extract a single field from a diagnostics vector created from DMCreateGlobalVector(diag_dm, &diag_vec) and then extract one of the diagnostic fields into a separate vector (S, say for the saturation field), to write it out to an HDF5 file using our own I/O interface. Is there a way to create the vector S so that it has the proper block size for storing just one field so that the data can be extracted with a call to VecStrideGather(diag_vec, DIAG_SATURATION, S, INSERT_VALUES)? Or is there a better way to do this altogether?

[jedbrown]

The general way is to DMCreateSubDM for the subset of fields that you want. It's flexible and used in places like PCFieldSplit.