6.3 Materials, geometries, gravitation and other plugin types

6.3.1 Material models

The material model is responsible for describing the various coefficients in the equations that ASPECT solves. To implement a new material model, you need to overload the aspect::MaterialModel::Interface class and use the ASPECT_REGISTER_MATERIAL_MODEL macro to register your new class. The implementation of the new class should be in namespace aspect::MaterialModel. An example of a material model implemented this way is given in Section 5.4.12.

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

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 structure. 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. The following members of MaterialModelOutputs need to be filled:

The variables refer to the coefficients $\eta ,C_p,k,\rho$ in equations (1)–(3), 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.

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 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 as compressible or not, see Section 2.10.3 on the Boussinesq model. As discussed there, incompressibility as described by this function does not necessarily imply that the density is constant; rather, it may still depend on temperature or pressure. In the current context, compressibility simply means whether we should solve the continuity equation as $\nabla \cdot \left(\rho u\right)=0$ (compressible Stokes) or as $\nabla \cdot 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.

Additionally, every material model has a member variable “model_dependence”, declared in the Interface class, which can be accessed from the plugin as “this$\to$model_dependence”. 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. 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 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).

6.3.2 Heating models

The heating model is responsible for describing the various terms in the energy equation (3), using the coefficients provided by the material model. These can be source terms such as radiogenic heat production or shear heating, they can be terms on the left-hand side of the equation, such as part of the latent heating terms, or they can be heating processes related to reactions. Each of these terms is described by a “heating model”, and a simulation can have none, one, or many heating models that are active throughout a simulation, with each heating model usually only implementing the terms for one specific heating process. One can then decide in the input file which heating processes should be included in the computation by providing a list of heating models in the input file.

When the equations are assembled and solved, the heating terms from all heating models used in the computation are added up.

To implement a new heating model, you need to overload the aspect::HeatingModel::Interface class and use the ASPECT_REGISTER_HEATING_MODEL macro to register your new class. The implementation of the new class should be in namespace aspect::HeatingModel.

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

The main properties of the material are computed in the function evaluate() that takes references to MaterialModelInputs and MaterialModelOutputs objects and is supposed to fill the HeatingModelOutputs structure. As in the material model, this function is handling lookups at an arbitrary number of positions, so for each heating term (for example the heating source terms), a std::vector is returned. The following members of HeatingModelOutputs need to be filled:

Heating source terms are terms on the right-hand side of the equations, such as the adiabatic heating $\alpha T\left(u\cdot \nabla p\right)$ in equation (3). An example for a left-hand side heating term is the temperature-derivative term $\rho T\Delta S\frac{\partial X}{\partial T}$ that is part of latent heat production (see equation (5)).⁴⁰ Rates of temperature change⁴¹ are used when the heating term is related to a reaction process, happening on a faster time scale than the temperature advection. All of these terms can depend on any of the material model inputs or outputs. Implementations of evaluate() may of course choose to ignore dependencies on any of these arguments.

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::HeatingModel::Interface class documentation. Note that some of the functions listed above have a default implementation, as discussed on the documentation page just mentioned.

Just like for material models, the functions initialize() and update() can be implemented if desired (the default implementation does nothing) and are useful if the heating model has an internal state. The function initialize() is called once during the initialization of ASPECT and can be used to allocate memory for the heating model, initialize state, or read information from an external file. The function update() is called at the beginning of every time step.

6.3.3 Geometry models

The geometry model is responsible for describing the domain in which we want to solve the equations. A domain is described in deal.II by a coarse mesh and, if necessary, an object that characterizes the boundary. Together, these two suffice to reconstruct any domain by adaptively refining the coarse mesh and placing new nodes generated by refining cells onto the surface described by the boundary object. The geometry model is also responsible for marking different parts of the boundary with different boundary indicators for which one can then, in the input file, select whether these boundaries should be Dirichlet-type (fixed temperature) or Neumann-type (no heat flux) boundaries for the temperature, and what kind of velocity conditions should hold there. In deal.II, a boundary indicator is a number of type types::boundary_id, but since boundaries are hard to remember and get right in input files, geometry models also have a function that provide a map from symbolic names that can be used to describe pieces of the boundary to the corresponding boundary indicators. For example, the simple box geometry model in 2d provides the map {"left"$\to$0, "right"$\to$1, "bottom"$\to$2,"top"$\to$3}, and we have consistently used these symbolic names in the input files used in this manual.

To implement a new geometry model, you need to overload the aspect::GeometryModel::Interface class and use the ASPECT_REGISTER_GEOMETRY_MODEL macro to register your new class. The implementation of the new class should be in namespace aspect::GeometryModel.

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

The kind of information these functions need to provide is extensively discussed in the documentation of this interface class at aspect::GeometryModel::Interface. The purpose of the last two functions has been discussed in the general overview of plugins above.

The create_coarse_mesh function does not only create the actual mesh (i.e., the locations of the vertices of the coarse mesh and how they connect to cells) but it must also set the boundary indicators for all parts of the boundary of the mesh. The deal.II glossary describes the purpose of boundary indicators as follows:

Two comments are in order here. First, if a coarse triangulation’s faces already accurately represent where you want to pose which boundary condition (for example to set temperature values or determine which are no-flow and which are tangential flow boundary conditions), then it is sufficient to set these boundary indicators only once at the beginning of the program since they will be inherited upon mesh refinement to the child faces. Here, at the beginning of the program is equivalent to inside the create_coarse_mesh()) function of the geometry module shown above that generates the coarse mesh.

Secondly, however, if you can only accurately determine which boundary indicator should hold where on a refined mesh – for example because the coarse mesh is the cube ${\left[0,L\right]}^3$ and you want to have a fixed velocity boundary describing an extending slab only for those faces for which $z>L-L_{\text{slab}}$ – then you need a way to set the boundary indicator for all boundary faces either to the value representing the slab or the fluid underneath after every mesh refinement step. By doing so, child faces can obtain boundary indicators different from that of their parents. deal.II triangulations support this kind of operations using a so-called post-refinement signal. In essence, what this means is that you can provide a function that will be called by the triangulation immediately after every mesh refinement step.

The way to do this is by writing a function that sets boundary indicators and that will be called by the Triangulation class. The triangulation does not provide a pointer to itself to the function being called, nor any other information, so the trick is to get this information into the function. C++ provides a nice mechanism for this that is best explained using an example:

What the call to std_cxx1x::bind does is to produce an object that can be called like a function with no arguments. It does so by taking the address of a function that does, in fact, take an argument but permanently fix this one argument to a reference to the coarse grid triangulation. After each refinement step, the triangulation will then call the object so created which will in turn call set_boundary_indicators<dim> with the reference to the coarse grid as argument.

This approach can be generalized. In the example above, we have used a global function that will be called. However, sometimes it is necessary that this function is in fact a member function of the class that generates the mesh, for example because it needs to access run-time parameters. This can be achieved as follows: assuming the set_boundary_indicators() function has been declared as a (non-static, but possibly private) member function of the MyGeometry class, then the following will work:

Here, like any other member function, set_boundary_indicators implicitly takes a pointer or reference to the object it belongs to as first argument. std::bind again creates an object that can be called like a global function with no arguments, and this object in turn calls set_boundary_indicators with a pointer to the current object and a reference to the triangulation to work on. Note that because the create_coarse_mesh function is declared as const, it is necessary that the set_boundary_indicators function is also declared const.

6.3.4 Gravity models

The gravity model is responsible for describing the magnitude and direction of the gravity vector at each point inside the domain. To implement a new gravity model, you need to overload the aspect::GravityModel::Interface class and use the ASPECT_REGISTER_GRAVITY_MODEL macro to register your new class. The implementation of the new class should be in namespace aspect::GravityModel.

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

The kind of information these functions need to provide is discussed in the documentation of this interface class at aspect::GravityModel::Interface. The first needs to return a gravity vector at a given position, whereas the second is called at the beginning of each time step, for example to allow a model to update itself based on the current time or the solution of the previous time step. The purpose of the last two functions has been discussed in the general overview of plugins above.

6.3.5 Initial conditions

The initial conditions model is responsible for describing the initial temperature distribution throughout the domain. It essentially has to provide a function that for each point can return the initial temperature. Note that the model (1)–(3) does not require initial values for the pressure or velocity. However, if coefficients are nonlinear, one can significantly reduce the number of initial nonlinear iterations if a good guess for them is available; consequently, ASPECT initializes the pressure with the adiabatically computed hydrostatic pressure, and a zero velocity. Neither of these two has to be provided by the objects considered in this section.

To implement a new initial conditions model, you need to overload the aspect::InitialConditions::Interface class and use the ASPECT_REGISTER_INITIAL_CONDITIONS macro to register your new class. The implementation of the new class should be in namespace aspect::InitialConditions.

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

The meaning of the first class should be clear. The purpose of the last two functions has been discussed in the general overview of plugins above.

6.3.6 Prescribed velocity boundary conditions

Most of the time, one chooses relatively simple boundary values for the velocity: either a zero boundary velocity, a tangential flow model in which the tangential velocity is unspecified but the normal velocity is zero at the boundary, or one in which all components of the velocity are unspecified (i.e., for example, an outflow or inflow condition where the total stress in the fluid is assumed to be zero). However, sometimes we want to choose a velocity model in which the velocity on the boundary equals some prescribed value. A typical example is one in which plate velocities are known, for example their current values or historical reconstructions. In that case, one needs a model in which one needs to be able to evaluate the velocity at individual points at the boundary. This can be implemented via plugins.

To implement a new boundary velocity model, you need to overload the aspect::VelocityBoundaryConditions::Interface class and use the ASPECT_REGISTER_VELOCITY_BOUNDARY_CONDITIONS macro to register your new class. The implementation of the new class should be in namespace aspect::VelocityBoundaryConditions.

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

The first of these functions needs to provide the velocity at the given point. The next two are other member functions that can (but need not) be overloaded if a model wants to do initialization steps at the beginning of the program or at the beginning of each time step. Examples are models that need to call an external program to obtain plate velocities for the current time, or from historical records, in which case it is far cheaper to do so only once at the beginning of the time step than for every boundary point separately. See, for example, the aspect::VelocityBoundaryConditions::GPlates class.

The remaining functions are obvious, and are also discussed in the documentation of this interface class at aspect::VelocityBoundaryConditions::Interface. The purpose of the last two functions has been discussed in the general overview of plugins above.

6.3.7 Temperature boundary conditions

The boundary conditions are responsible for describing the temperature values at those parts of the boundary at which the temperature is fixed (see Section 6.3.3 for how it is determined which parts of the boundary this applies to).

To implement a new boundary conditions model, you need to overload the aspect::BoundaryTemperature::Interface class and use the ASPECT_REGISTER_BOUNDARY_TEMPERATURE_MODEL macro to register your new class. The implementation of the new class should be in namespace aspect::BoundaryTemperature.

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

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 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.

6.3.8 Postprocessors: Evaluating the solution after each time step

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.).

Motivation. The original motivation for postprocessors is that the goal of a simulation is of course not the simulation itself, but that we want to do something with the solution. Examples for already existing postprocessors are:

• Generating output in file formats that are understood by visualization programs. This is facilitated by the aspect::Postprocess::Visualization class and a separate class of visualization postprocessors, see Section 6.3.9.
• 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 aspect::Postprocess::VelocityStatistics, aspect::Postprocess::TemperatureStatistics, aspect::Postprocess::HeatFluxStatistics classes, respectively.

Since writing this text, there may have been other additions as well.

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 information from one invocation at one timestep to the next invocation,⁴² there is nothing that prohibits postprocessors from doing so. For example, the following ideas would fit nicely into the postprocessor framework:

• Passive particles: If one would like to follow the trajectory of material as it is advected along with the flow field, one technique is to use particles. To implement this, one would start with an initial population of particles distributed in a certain way, for example close to the core-mantle boundary. At the end of each time step, one would then need to move them forward with the flow field by one time increment. As long as these particles do not affect the flow field (i.e., they do not 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.
• Surface or crustal processes: Another possibility would be to keep track of surface or crustal processes induced by mantle flow. An example would be to keep track of the thermal history of a piece of crust by updating it every time step with the heat flux from the mantle below. One could also imagine integrating changes in the surface topography by considering the surface divergence of the surface velocity computed in the previous time step: if the surface divergence is positive, the topography is lowered, 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.

The statistics file. Postprocessors 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, or heat flux statistics). While the former are on their own in implementing how they want to store their data to disk, there is a mechanism in place that 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 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 parameter file. It allows for easy visualization of trends over all time steps. To see how to put data into this statistics object, take a look at the existing postprocessor objects.

Note that the data deposited into the statistics object need not be numeric in 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 other plugins. Specifically, you’ll have to write a class that overloads the aspect::Postprocess::Interface base class and use the ASPECT_REGISTER_POSTPROCESSOR macro to register your new class. The implementation of the new class should be in namespace aspect::Postprocess.

In reality, however, implementing new postprocessors is often more difficult. Primarily, this difficulty results from two facts:

• Postprocessors are not self-contained (only providing information) but in fact need to access the solution of the model at each time step. That is, of course, the purpose of postprocessors, but it requires that the writer of a plugin has a certain amount of knowledge of how the solution is computed by the main Simulator class, and how it is represented in data structures. To alleviate this somewhat, and to insulate the two worlds from each other, postprocessors do not directly access the data structures of the simulator class. Rather, in addition to deriving from the aspect::Postprocess::Interface base class, postprocessors also derive from the SimulatorAccess class that has a number of member functions postprocessors can call to obtain read-only access to some of the information stored in the main class of ASPECT. See the documentation of this class to see what kind of information is available to postprocessors. See also Section 6.1 for more information about the SimulatorAccess class.
• Writing a new postprocessor typically requires a fair amount of knowledge how to leverage the deal.II library to extract information from the solution. The existing postprocessors are 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:

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 restart. This is important since many postprocessors store some state – 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 needs to be numbered.

The typical case is that this state is the same across all processors of a parallel computation. Consequently, what ASPECT writes into the checkpoint file is only the state obtained from the postprocessors on processor 0 of a parallel computation. On restart, all processors read from the same file and the postprocessors on all processors will be initialized by what the same postprocessor on processor 0 wrote.

There are situations where postprocessors do in fact store complementary information on different processors. At the time of writing this, one example 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() 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 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 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 file specifically for this postprocessor, use MPI I/O or other ways to write into it, and do the reverse operation in load().

Note that this approach requires that ASPECT actually calls the save() function on all processors. This in fact happens – though ASPECT also discards the result on all but processor zero.

6.3.9 Visualization postprocessors

As mentioned in the previous section, one of the postprocessors that are already implemented in ASPECT is the aspect::Postprocess::Visualization class that takes the solution and outputs it as a collection of files that can then be visualized graphically, see Section 4.4. The question is which variables to output: the solution of the basic equations we solve here is characterized by the velocity, pressure and temperature; on the other hand, we are frequently interested in derived, spatially and temporally variable quantities such as the viscosity for the actual pressure, temperature and strain rate at a given location, or seismic wave speeds.

ASPECT already implements a good number of such derived quantities that one may want to visualize. On the other hand, always outputting all of them would yield very large output files, and would furthermore not scale very well as the list continues to grow. Consequently, as with the postprocessors described in the previous section, what can be computed is implemented in a number of plugins and what is computed is selected in the input parameter file (see Section ??).

Defining visualization postprocessors works in much the same way as for the other plugins discussed in this section. Specifically, an implementation of such a plugin needs to be a class that derives from interface classes, should by convention be in namespace aspect::Postprocess::VisualizationPostprocessors, and is registered using a macro, here called ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR. Like the postprocessor plugins, visualization postprocessors can derive from class aspect::Postprocess::SimulatorAccess if they need to know specifics of the simulation such as access to the material models and to get access to the introspection facility outlined in Section 6.1. A typical example is the plugin that produces the viscosity as a spatially variable field by evaluating the viscosity function of the material model using the pressure, temperature and location of each visualization point (implemented in the aspect::Postprocess::VisualizationPostprocessors::Viscosity class). On the other hand, a hypothetical plugin that simply outputs the norm of the strain rate $\sqrt{\epsilon \left(u\right):\epsilon \left(u\right)}$ would not need access to anything but the solution vector (which the plugin’s main function is given as an argument) and consequently is not derived from the aspect::Postprocess::SimulatorAccess class.⁴³

Visualization plugins can come in two flavors:

• Plugins that compute things from the solution in a point-wise way: The classes in this group are derived not only from the respective interface class (and possibly the SimulatorAccess class) but also from the deal.II class DataPostprocessor or any of the classes like DataPostprocessorScalar or DataPostprocessorVector. These classes can be thought of as filters: DataOut will call a function in them for every cell and this function will transform the values or gradients of the solution and other information such as the location of quadrature points into the desired quantity to output. A typical case would be if the quantity $g\left(x\right)$ you want to output can be written as a function $g\left(x\right)=G\left(u\left(x\right),\nabla u\left(x\right),x,...\right)$ in a point-wise sense where $u\left(x\right)$ is the value of the solution vector (i.e., the velocities, pressure, temperature, etc) at an evaluation point. In the context of this program an example would be to output the density of the medium as a spatially variable function since this is a quantity that for realistic media depends point-wise on the values of the solution.

  To sum this, slightly confusing multiple inheritance up, visualization postprocessors do the following:

  • If necessary, they derive from aspect::Postprocess::SimulatorAccess.
  • They derive from aspect::Postprocess::VisualizationPostprocessors::Interface. The functions of this interface class are all already implemented as doing nothing in the base class but can be overridden in a plugin. Specifically, the following functions exist:
        class Interface 
        { 
          public: 
            static 
            void 
            declare_parameters (ParameterHandler &prm); 
     
            virtual 
            void 
            parse_parameters (ParameterHandler &prm); 
     
            virtual 
            void save (std::map<std::string, std::string> &status_strings) const; 
     
            virtual 
            void load (const std::map<std::string, std::string> &status_strings); 
        }; 
      
  • They derive from either the dealii::DataPostprocessor class, or the simpler to use dealii::DataPostprocessorScalar or dealii::DataPostprocessorVector classes. For example, to derive from the second of these classes, the following interface functions has to be implemented:
        class dealii::DataPostprocessorScalar 
        { 
          public: 
            virtual 
            void 
            compute_derived_quantities_vector 
              (const std::vector<Vector<double> >              &uh, 
               const std::vector<std::vector<Tensor<1,dim> > > &duh, 
               const std::vector<std::vector<Tensor<2,dim> > > &dduh, 
               const std::vector<Point<dim> >                  &normals, 
               const std::vector<Point<dim> >                  &evaluation_points, 
               std::vector<Vector<double> >                    &computed_quantities) const; 
        }; 
      

    What this function does is described in detail in the deal.II documentation. In addition, one has to write a suitable constructor to call dealii::DataPostprocessorScalar::DataPostprocessorScalar.

• Plugins that compute things from the solution in a cell-wise way: The second possibility is for a class to not derive from dealii::DataPostprocessor but instead from the aspect::Postprocess::VisualizationPostprocessors::CellDataVectorCreator class. In this case, a visualization postprocessor would generate and return a vector that consists of one element per cell. The intent of this option is to output quantities that are not point-wise functions of the solution but instead can only be computed as integrals or other functionals on a per-cell basis. A typical case would be error estimators that do depend on the solution but not in a point-wise sense; rather, they yield one value per cell of the mesh. See the documentation of the CellDataVectorCreator class for more information.

If all of this sounds confusing, we recommend consulting the implementation of the various visualization plugins that already exist in the ASPECT sources, and using them as a template.

6.3.10 Mesh refinement criteria

Despite research since the mid-1980s, it isn’t completely clear how to refine meshes for complex situations like the ones modeled by ASPECT. The basic problem is that mesh refinement criteria either can refine based on some variable such as the temperature, the pressure, the velocity, or a compositional field, but that oftentimes this by itself is not quite what one wants. For example, we know that Earth has discontinuities, e.g., at 440km and 610km depth. In these places, densities and other material properties suddenly change. Their resolution in computation models is important as we know that they affect convection patterns. At the same time, there is only a small effect on the primary variables in a computation – maybe a jump in the second or third derivative, for example, but not a discontinuity that would be clear to see. As a consequence, automatic refinement criteria do not always refine these interfaces as well as necessary.

To alleviate this, ASPECT has plugins for mesh refinement. Through the parameters in Section ??, one can select when to refine but also which refinement criteria should be used and how they should be combined if multiple refinement criteria are selected. Furthermore, through the usual plugin mechanism, one can extend the list of available mesh refinement criteria (see the parameter “Strategy” in Section ??). 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 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 aspect::MeshRefinement::Interface class and use the ASPECT_REGISTER_MESH_REFINEMENT_CRITERION macro to register your new class. The implementation of the new class should be in namespace aspect::MeshRefinement.

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

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 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 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.

6.3.11 Criteria for terminating a simulation

ASPECT allows for different ways of terminating a simulation. For example, the simulation may have reached a final time specified in the input file. However, it also allows for ways to terminate a simulation when it has reached a steady state (or, rather, some criterion determines that it is close enough to steady state), or by an external action such as placing a specially named file in the output directory. The criteria determining termination of a simulation are all implemented in plugins. The parameters describing these criteria are listed in Section ??.

To implement a termination criterion, you need to overload the aspect::TerminationCriteria::Interface class and use the ASPECT_REGISTER_TERMINATION_CRITERION macro to register your new class. The implementation of the new class should be in namespace aspect::TerminationCriteria.

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

The first of these functions returns a value that indicates whether the simulation should be terminated. Typical examples can be found in the existing implementations in the source/termination_criteria 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 simulation.

The remaining functions are obvious, and are also discussed in the documentation of this interface class at aspect::TerminationCriteria::Interface. The purpose of the last two functions has been discussed in the general overview of plugins above.