Skip to content

Commit

Permalink
Merge pull request #6065 from bangerth/doc-update
Browse files Browse the repository at this point in the history
Update the documentation of plugins (part 2).
  • Loading branch information
gassmoeller authored Oct 10, 2024
2 parents d769ee6 + da389bd commit b776b39
Show file tree
Hide file tree
Showing 7 changed files with 105 additions and 333 deletions.
107 changes: 22 additions & 85 deletions doc/sphinx/user/extending/plugin-types/material-models.md
Original file line number Diff line number Diff line change
Expand Up @@ -11,74 +11,20 @@ implementation of the new class should be in namespace
`aspect::MaterialModel`. An example of a material model implemented this way
is given in {ref}`sec:benchmarks:davies_et_al:case2.3`.

Specifically, your new class needs to implement the following interface:

```{code-block} c++
template <int dim>
class aspect::MaterialModel::Interface
{
public:
// Physical parameters used in the basic equations
virtual void evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const=0;
virtual bool is_compressible () const = 0;
// Reference quantities
virtual double reference_viscosity () const = 0;
// Functions used in dealing with run-time parameters
static void
declare_parameters (ParameterHandler &prm);
virtual void
parse_parameters (ParameterHandler &prm);
// Optional:
virtual void initialize ();
virtual void update ();
}
```

The main properties of the material are computed in the function evaluate()
that takes a struct of type MaterialModelInputs and is supposed to fill a
MaterialModelOutputs class. For performance reasons this function is
The main properties of the material are computed in the member function `evaluate()`
that takes a struct of type `MaterialModelInputs` and is supposed to fill a
`MaterialModelOutputs` object. For performance reasons this function is
handling lookups at an arbitrary number of positions, so for each variable
(for example viscosity), a std::vector is returned. Some of the MaterialModelOutputs
members that need to be filled are:
(for example viscosity), a `std::vector` is returned.

```{code-block} c++
class MaterialModelOutputs
{
public:
std::vector<double> viscosities;
std::vector<double> densities;
std::vector<double> thermal_expansion_coefficients;
std::vector<double> specific_heat;
// more here
}
```

The variables refer to the coefficients $\eta,C_p,k,\rho$ in equations
{math:numref}`eq:stokes-1`&ndash;{math:numref}`eq:temperature`, each as a function of
temperature, pressure, position, compositional fields and, in the case of the
viscosity, the strain rate (all handed in by MaterialModelInputs).
Implementations of evaluate() may of course choose to ignore dependencies on
any of these arguments. In writing a new material model, you should consider
Implementations of `evaluate()` may of course choose to ignore dependencies on
any of these material model inputs; for example, the `MaterialModels::Simpler` class
describes a constant viscosity and a density that only depends linearly on the temperature.
In writing a new material model, you should consider
coefficient self-consistency
({ref}`sec:coefficient_self_consistency`).

The remaining functions are used in postprocessing as well as handling
run-time parameters. The exact meaning of these member functions is documented
in the [aspect::MaterialModel::Interface](https://aspect.geodynamics.org/doc/doxygen/namespaceaspect_1_1MaterialModel.html)
class documentation. Note that some of the
functions listed above have a default implementation, as discussed on the
documentation page just mentioned.

The function `is_compressible` returns whether we should consider the material
The function `is_compressible()` returns whether we should consider the material
as compressible or not, see {ref}`sec:methods:approximate-equations:ba` on the
Boussinesq model. As discussed there, incompressibility as described by this
function does not necessarily imply that the density is constant; rather, it
Expand All @@ -87,31 +33,22 @@ compressibility simply means whether we should solve the continuity equation
as $\nabla \cdot (\rho \mathbf u)=0$ (compressible Stokes) or as
$\nabla \cdot \mathbf{u}=0$ (incompressible Stokes).

The purpose of the parameter handling functions has been discussed in the
general overview of plugins above.

The functions initialize() and update() can be implemented if desired (the
default implementation does nothing) and are useful if the material model has
internal state. The function initialize() is called once during the
initialization of ASPECT and can be used to
allocate memory, initialize state, or read information from an external file.
The function update() is called at the beginning of every time step.
The other functions derived classes can overload are used in postprocessing as well as handling
run-time parameters. The exact meaning of these member functions is documented
in the [aspect::MaterialModel::Interface](https://aspect.geodynamics.org/doc/doxygen/namespaceaspect_1_1MaterialModel.html)
class documentation. Note that some of these
functions have a default implementation, as discussed on the
documentation page just mentioned.

Additionally, every material model has a member variable
"modeldependence," declared in the Interface class, which can be
accessed from the plugin as `this$\rightarrow$modeldependence`.
This structure describes the nonlinear dependence of the various coefficients
on pressure, temperature, composition or strain rate. This information will be
used in future versions of ASPECT to implement
a fully nonlinear solution scheme based on, for example, a Newton iteration.
"model_dependence," declared in the `Interface` class, which can be
accessed from the plugin as `this->model_dependence`.
This object describes the nonlinear dependence of the various coefficients
on pressure, temperature, composition, or strain rate. This information is
used to implement
fully nonlinear solution schemes based on, for example, a Newton iteration.
The initialization of this variable is optional, but only plugins that declare
correct dependencies can benefit from these solver types. All packaged
material models declare their dependencies in the parseparameters() function
material models declare their dependencies in the `parse_parameters()` function
and can be used as a starting point for implementations of new material
models.

Older versions of ASPECT used to have
individual functions like `viscosity()` instead of the `evaluate()` function
discussed above. This old interface is no longer supported, restructure your
plugin to implement `evaluate()` instead (even if this function only calls the
old functions).
34 changes: 4 additions & 30 deletions doc/sphinx/user/extending/plugin-types/mesh-refinement.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ refinement criteria (see the parameter "Strategy" in
{ref}`parameters:Mesh_20refinement/Strategy`). Each such plugin is responsible for producing a
vector of values (one per active cell on the current processor, though only
those values for cells that the current processor owns are used) with an
indicator of how badly this cell needs to be refined: large values mean that
indicator of how badly this cell needs to be refined: Large values mean that
the cell should be refined, small values that the cell may be coarsened away.

To implement a new mesh refinement criterion, you need to overload the
Expand All @@ -33,35 +33,9 @@ To implement a new mesh refinement criterion, you need to overload the
The implementation of the new class should be in namespace
`aspect::MeshRefinement`.

Specifically, your new class needs to implement the following basic interface:

```{code-block} c++
template <int dim>
class aspect::MeshRefinement::Interface
{
public:
virtual
void
execute (Vector<float> &error_indicators) const = 0;
static
void
declare_parameters (ParameterHandler &prm);
virtual
void
parse_parameters (ParameterHandler &prm);
};
```

The first of these functions computes the set of refinement criteria (one per
cell) and returns it in the given argument. Typical examples can be found in
The primary function of this interface, `execute()`, computes the set of refinement criteria (one per
cell) and returns through its function argument. Typical examples can be found in
the existing implementations in the `source/mesh_refinement` directory. As
usual, your termination criterion implementation will likely need to be
derived from the `SimulatorAccess` to get access to the current state of the
derived from the `SimulatorAccess` class to get access to the current state of the
simulation.

The remaining functions are obvious, and are also discussed in the
documentation of this interface class at
`aspect::MeshRefinement::Interface`. The purpose of the last two functions
has been discussed in the general overview of plugins above.
97 changes: 31 additions & 66 deletions doc/sphinx/user/extending/plugin-types/postprocessors.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,11 @@

Postprocessors are arguably the most complex and powerful of the plugins
available in ASPECT since they do not only
passively provide any information but can actually compute quantities derived
from the solution. They are executed once at the end of each time step and,
unlike all the other plugins discussed above, there can be an arbitrary number
of active postprocessors in the same program (for the plugins discussed in
previous sections it was clear that there is always exactly one material
model, geometry model, etc.).
passively provide information that enters the simulation as parameters or boundary
conditions, but can actually compute quantities derived
from the solution. They are executed once at the end of each time step (or at the
end of nonlinear iterations); there can be an arbitrary number
of active postprocessors in the same program.

## Motivation.

Expand All @@ -24,15 +23,17 @@ solution. Examples for already existing postprocessors are:
- Computing statistics about the velocity field (e.g., computing minimal,
maximal, and average velocities), temperature field (minimal, maximal, and
average temperatures), or about the heat fluxes across boundaries of the
domain. This is provided by the
domain. These are provided by the
`aspect::Postprocess::VelocityStatistics`,
`aspect::Postprocess::TemperatureStatistics`,
`aspect::Postprocess::HeatFluxStatistics` classes, respectively.

Since writing this text, there may have been other additions as well.
In addition to these examples, there are many other postprocessors that are part of ASPECT.

However, postprocessors can be more powerful than this. For example, while the
ones listed above are by and large stateless, i.e., they do not carry
Postprocessors can be more powerful than the examples: They do not only need
to passively take the solution and transform it into something that ends up
in output files or on the screen. For example, while the
ones listed above are by and large "stateless", i.e., they do not carry
information from one invocation at one timestep to the next invocation,[^footnote1]
there is nothing that prohibits postprocessors from doing so. For example, the
following ideas would fit nicely into the postprocessor framework:
Expand All @@ -47,8 +48,8 @@ following ideas would fit nicely into the postprocessor framework:
carry any information that feeds into material properties; in other words,
they are *passive*), their location could well be stored in a
postprocessor object and then be output in periodic intervals for
visualization. In fact, such a passive particle postprocessor is already
available.
visualization. In fact, this is how ASPECT's particle infrastructure
was originally implemented.

- *Surface or crustal processes:* Another possibility would be to keep track
of surface or crustal processes induced by mantle flow. An example would
Expand All @@ -60,12 +61,19 @@ following ideas would fit nicely into the postprocessor framework:
eventually forming a trench; if the divergence is negative, a mountain
belt eventually forms.

In all of these cases, the essential limitation is that postprocessors are
*passive*, i.e., that they do not affect the simulation but only observe it.
While these examples suggest that postprocessors are
*passive*, i.e., do not affect the simulation but only observe it, this is not
an essential limitation. In fact, some of the existing postprocessors perform
actions, store the results, and these results can then be queried by other
parts of ASPECT or other postprocessors. For example, one could imagine letting
a postprocessor compute a depth average at the end of each time step, and
then letting a material model query this information the next time it is
asked to provide material properties for a specific point.


## The statistics file.

Postprocessors fall into two categories: ones that produce lots of output
Postprocessors tend to fall into two categories: ones that produce lots of output
every time they run (e.g., the visualization postprocessor), and ones that
only produce one, two, or in any case a small and fixed number of often
numerical results (e.g., the postprocessors computing velocity, temperature,
Expand All @@ -75,7 +83,7 @@ allows the latter class of postprocessors to store their data into a central
file that is updated at the end of each time step, after all postprocessors
are run.

To this end, the function that executes each of the postprocessors is given a
To this end, the `execute()` function that postprocessors implement and that performs their "action" is given a
reference to a `dealii::TableHandler` object that allows to store data in
named columns, with one row for each time step. This table is then stored in
the `statistics` file in the directory designated for output in the input
Expand All @@ -88,6 +96,7 @@ type, though it often is. An example of text-based entries in this table is
the visualization class that stores the name of the graphical output file
written in a particular time step.


## Implementing a postprocessor.

Ultimately, implementing a new postprocessor is no different than any of the
Expand Down Expand Up @@ -123,56 +132,12 @@ Primarily, this difficulty results from two facts:
certainly good examples to start from in trying to understand how to do
this.

Given these comments, the interface a postprocessor class has to implement is
rather basic:

```{code-block} c++
template <int dim>
class aspect::Postprocess::Interface
{
public:
virtual
std::pair<std::string,std::string>
execute (TableHandler &statistics) = 0;
virtual
void
save (std::map<std::string, std::string> &status_strings) const;
virtual
void
load (const std::map<std::string, std::string> &status_strings);
static
void
declare_parameters (ParameterHandler &prm);
virtual
void
parse_parameters (ParameterHandler &prm);
};
```

The purpose of these functions is described in detail in the documentation of
the `aspect::Postprocess::Interface` class. While the first one is
responsible for evaluating the solution at the end of a time step, the
`save/load` functions are used in checkpointing the program and restarting it
at a previously saved point during the simulation. The first of these
functions therefore, needs to store the status of the object as a string under
a unique key in the database described by the argument, while the latter
function restores the same state as before by looking up the status string
under the same key. The default implementation of these functions is to do
nothing; postprocessors that do have non-static member variables that contain
a state need to overload these functions.

There are numerous postprocessors already implemented. If you want to
implement a new one, it would be helpful to look at the existing ones to see
how they implement their functionality.

## Postprocessors and checkpoint/restart.

Postprocessors have `save()` and `load()` functions that are used to write the
data a postprocessor has into a checkpoint file, and to load it again upon
Postprocessors have `save()` and `load()` functions (inherited from
`aspect::Postprocess::Interface`) that are used to write the
data a postprocessor stores between successive invocations into a checkpoint file, and to load it again upon
restart. This is important since many postprocessors store some state &ndash;
say, a temporal average over all the time steps seen so far, or the number of
the last graphical output file generated so that we know how the next one
Expand All @@ -191,13 +156,13 @@ is the postprocessor that supports advecting passive particles along the
velocity field: on every processor, it handles only those particles that lie
inside the part of the domain that is owned by this MPI rank. The
serialization approach outlined above can not work in this case, for obvious
reasons. In cases like this, one needs to implement the `save()` and `load()`
reasons. In cases like this, one needs to implement `save()` and `load()`
differently than usual: one needs to put all variables that are common across
processors into the maps of string as usual, but one then also needs to save
all state that is different across processors, from all processors. There are
two ways: If the amount of data is small, you can use MPI communications to
two ways: If the amount of data is small, you can use MPI communication to
send the state of all processors to processor zero, and have processor zero
store it in the result so that it gets written into the checkpoint file; in
return it so that it gets written into the checkpoint file; in
the `load()` function, you will then have to identify which part of the text
written by processor 0 is relevant to the current processor. Or, if your
postprocessor stores a large amount of data, you may want to open a restart
Expand Down
41 changes: 6 additions & 35 deletions doc/sphinx/user/extending/plugin-types/temp-bc.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,41 +12,12 @@ To implement a new boundary conditions model, you need to overload the
The implementation of the new class should be in namespace
`aspect::BoundaryTemperature`.

Specifically, your new class needs to implement the following basic interface:

```{code-block} c++
template <int dim>
class aspect::BoundaryTemperature::Interface
{
public:
virtual
double
temperature (const GeometryModel::Interface<dim> &geometry_model,
const unsigned int boundary_indicator,
const Point<dim> &location) const = 0;
virtual
double minimal_temperature () const = 0;
virtual
double maximal_temperature () const = 0;
static
void
declare_parameters (ParameterHandler &prm);
virtual
void
parse_parameters (ParameterHandler &prm);
};
```

The first of these functions needs to provide the fixed temperature at the
given point. The geometry model and the boundary indicator of the particular
The principle function you need to overload needs to provide the fixed temperature at a
given point. The boundary indicator of the particular
piece of boundary on which the point is located is also given as a hint in
determining where this point may be located; this may, for example, be used to
determine if a point is on the inner or outer boundary of a spherical shell.
The remaining functions are obvious, and are also discussed in the
documentation of this interface class at
`aspect::BoundaryTemperature::Interface`. The purpose of the last two
functions has been discussed in the general overview of plugins above.
Models may also need to query the geometry model (via the `SimulatorAccess` class)
to understand what a specific boundary indicator is supposed to mean &ndash; e.g.,
if the given boundary indicator corresponds to an inner or outer surface of a shell
geometry.
Loading

0 comments on commit b776b39

Please sign in to comment.