Writing to the Dependency Graph and PK tree

ATS builds on the multi-physics tool Arcos, which set out the design of a multi-physics framework. Since the initial paper in Environmental Modeling & Software, advances by the development team of Amanzi and ATS have evolved the concepts. This lists out best practices and design patterns for using Arcos in writing ATS Evaluators and PKs.

If you are unsure of what an Evaluator or PK is, please first read the original paper.

Names vs Types

Individual variables in ATS are referred to and accessed by their name. Names are strings, and are always be user-defined. On the other hand, types are strings that are defined in the code, and are used to map from the user spec to specific classes in the code. They are developer-defined, and are effectively keys into dictionaries whose value are classes that can be instantiated.

For instance, one might define an evaluator named my_saturation that is of type WRM evaluator.

Naming conventions

See https://amanzi.github.io/ats/stable/input_spec/introduction.html#naming

Don't hard-code names

Don't hard-code names in your evaluators or PKs. For instance, there is a temptation to write something like (pseudo-code):

class WRM : public Evaluator
{
  WRM() {
    dependencies_.push_back( "pressure" );
  }
};

But this has many bad implications:

1. Your evaluator may only be used on the subsurface domain.
2. Your evaluator cannot be generic -- what if a user is doing a dual-permeability media problem and wants to write your WRM as a function of some partial pressure, or macropore pressure?
3. The user cannot customize this name -- they should be able to call it whatever they want. (Ok, only maybe bad.)

Instead, use Keys::readKey() from src/utils/Key.hh to get the variable name. Note that this includes a domain prefix, provided by the user, and a default value, where the developer can provide the "standard" name (e.g. "pressure" in our example).

For instance, most classes that inherit from EvaluatorSecondaryMonotypeCV have, in their constructor, something like:

WRM::WRM(Teuchos::ParameterList& plist)
{
  std::string domain_name = Keys::getDomain(my_keys_.front().first);
  pres_key_ = Keys::readKey(plist_, domain_name, "pressure", "pressure");
}

This will allow the user to provide either:

<Parameter name="pressure key" type="string" value="my_pressure_name"/> OR <Parameter name="pressure key suffix" type="string" value="my_pressure_name"/>

which will result in a user-provided name or suffix being used to compute MyEvaluator.

Tags

Tags are conceptually a way of distinguishing between two copies of the same variable at different times in the unit time increment. The most commonly used tags are Tags::CURRENT and Tags::NEXT. These refer to variables at the beginning and end of the global timestep.

All data is defined at a tag. All output is written at the Tags::NEXT tag (unless this tag doesn't exist, then we write Tags::CURRENT). All evaluators are defined at a tag, and may be different for the same variable at different tags.

Conceptually, the unique identifier for a piece of data is a KeyTag, which (to extend the above naming convention) are sometimes separated by the @ sign, e.g. DOMAIN-ROOT_NAME@TAG, and is simply a std::pair<Key,Tag>.

Evaluators can span tags

By default, most evaluators work at a single tag, and take their dependencies as the same tag as they compute. But they don't have to.

See, for instance, an EvaluatorTemporalInterpolation.

Don't hard-code tags.

All evaluators compute data at a single tag. All PKs advance their physics from a tag_current to a tag_next. If you hard-code either of these, you limit the reusability of your Evaluator/PK.

Prefer to write evaluators that infer their tag from the computed quantity's tag:

WRM::WRM(Teuchos::ParameterList& plist)
{
  std::string domain_name = Keys::getDomain(my_keys_.front().first);
  pres_key_ = Keys::readKey(plist_, domain_name, "pressure", "pressure");
  Tag tag = my_keys_.front().second;
  dependencies_.insert(KeyTag{ pres_key_, tag });
}

Prefer to write your PK's variables with respect to PK::tag_current_ rather than Tags::CURRENT and PK::tag_next_ rather than Tags::NEXT. This allows, for instance, automated subcycling of your PK via MPCSubcycled with no extra work. If you hard-code things like:

S_->Require<CompositeVector, CompositeVectorSpace>(pres_key_, Tags::NEXT);

you cannot be subcycled, instead prefer:

S_->Require<CompositeVector, CompositeVectorSpace>(pres_key_, tag_next_);

Evaluators

Evaluators are relatively easy to write, and are a very common extension point.

Start from the right base class.

The hardest part is starting from the right base class. Most new evaluators are "secondary" evaluators, meaning they depend upon other data. The most common base classes to start from are:

1. EvaluatorSecondaryMonotypeCV These are for evaluators that compute a field (CompositeVector) as a function of only other fields.

The most common types of these include:

• Algebraic evaluators, which contain all of the same set of components (e.g. if the computed quantity is on cells, all dependencies are also defined on cells) and are defined on the same mesh. Think of these as constitutive relationships, e.g. density as a function of temperature and pressure.
• Functions of fields on different meshes (e.g. column integrators, which compute a quantity on the surface as a function of data on the subsurface below)
• Functions of fields with different components (e.g. computing cell values as a function of the neighboring faces)

2. EvaluatorSecondary These are for evaluators that may be functions of different data types (e.g. EvaluatorTemporalInterpolation, which computes a field as a function of other fields and also of time, which is a double).

EnsureCompatibility

Evaluator::EnsureCompatibility() is one of the most confusing functions to write for an evaluator. What this does is check as many things as possible about this data and its dependencies to make sure they are consistent.

All the possible things one might want to do in this method are documented in EvaluatorSecondaryMonotypeCV::EnsureCompatibility(). It is very common to override one or many of the sub-methods called in that method. For algebraic evaluators that compute only one thing, no overrides are needed. Common overrides include:

• override EnsureCompatibility_Structure_() -- for evaluators that compute multiple fields, typically the must have the same mesh and same components. In that case just:

void EnsureCompatibility_Structure_() { EnsureCompatibility_StructureSame_(); }

If the fields have different components, then ensure those are correct in how they are different.

• override EnsureCompatibility_ToDeps_() to require dependencies with different meshes, components, etc.

PKs

PKs are much harder to write. They must include all fields they use, time integration, the ability to recover from a failed timestep, etc. To write a new PK, it is recommended that the developer first read and understand another PK that is very similar to what they want to do, e.g. implicit vs explicit time integration, etc.

Choose the right base class, and use it!

There are a lot of base classes that provide some basic functionality:

• PK is the fundamental base class. It provides a verbose object for logging, a name, a type, and a TreeVector solution_. All PKs will (directly or indirectly) derive from this base class, and all other base classes also derive from this base class. However, due to the "dreaded diamond" of multiple inheritance in C++, most base classes derive from this virtually, meaning all PKs should explicitly call this constructor. See any other PK as an example of how to do this.

• PK_Physical is a base class for individual, uncoupled PDEs. They are "physical" in the sense that they have a single primary variable, and are defined on a mesh. They provide a domain name, mesh, primary variable key_, and (in the case of PK_Physical_Default) provide methods to commit a step, fail a step, read the initial conditions, etc.

• PK_BDF is a base class for