diff --git a/doc/sphinx/user/extending/plugin-types/material-models.md b/doc/sphinx/user/extending/plugin-types/material-models.md index 592491c1ca2..f447067221b 100644 --- a/doc/sphinx/user/extending/plugin-types/material-models.md +++ b/doc/sphinx/user/extending/plugin-types/material-models.md @@ -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 - 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 viscosities; - std::vector densities; - std::vector thermal_expansion_coefficients; - std::vector specific_heat; - // more here -} -``` - -The variables refer to the coefficients $\eta,C_p,k,\rho$ in equations -{math:numref}`eq:stokes-1`–{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 @@ -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). diff --git a/doc/sphinx/user/extending/plugin-types/mesh-refinement.md b/doc/sphinx/user/extending/plugin-types/mesh-refinement.md index 15b65a1fcfc..558814836f1 100644 --- a/doc/sphinx/user/extending/plugin-types/mesh-refinement.md +++ b/doc/sphinx/user/extending/plugin-types/mesh-refinement.md @@ -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 @@ -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 - class aspect::MeshRefinement::Interface - { - public: - virtual - void - execute (Vector &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. diff --git a/doc/sphinx/user/extending/plugin-types/postprocessors.md b/doc/sphinx/user/extending/plugin-types/postprocessors.md index 9eed908e215..05023e87164 100644 --- a/doc/sphinx/user/extending/plugin-types/postprocessors.md +++ b/doc/sphinx/user/extending/plugin-types/postprocessors.md @@ -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. @@ -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: @@ -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 @@ -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, @@ -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 @@ -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 @@ -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 - class aspect::Postprocess::Interface - { - public: - virtual - std::pair - execute (TableHandler &statistics) = 0; - - virtual - void - save (std::map &status_strings) const; - - virtual - void - load (const std::map &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 – 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 @@ -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 diff --git a/doc/sphinx/user/extending/plugin-types/temp-bc.md b/doc/sphinx/user/extending/plugin-types/temp-bc.md index f04d740f676..a56f4798492 100644 --- a/doc/sphinx/user/extending/plugin-types/temp-bc.md +++ b/doc/sphinx/user/extending/plugin-types/temp-bc.md @@ -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 - class aspect::BoundaryTemperature::Interface - { - public: - virtual - double - temperature (const GeometryModel::Interface &geometry_model, - const unsigned int boundary_indicator, - const Point &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 – e.g., +if the given boundary indicator corresponds to an inner or outer surface of a shell +geometry. diff --git a/doc/sphinx/user/extending/plugin-types/terminate-criteria.md b/doc/sphinx/user/extending/plugin-types/terminate-criteria.md index 7b23e628b95..d20ca273ac8 100644 --- a/doc/sphinx/user/extending/plugin-types/terminate-criteria.md +++ b/doc/sphinx/user/extending/plugin-types/terminate-criteria.md @@ -15,34 +15,8 @@ To implement a termination criterion, you need to overload the implementation of the new class should be in namespace `aspect::TerminationCriteria`. -Specifically, your new class needs to implement the following basic interface: - -```{code-block} c++ -template - class aspect::TerminationCriteria::Interface - { - public: - virtual - bool - execute () const = 0; - - static - void - declare_parameters (ParameterHandler &prm); - - virtual - void - parse_parameters (ParameterHandler &prm); - }; -``` - -The first of these functions returns a value that indicates whether the +The principal function that needs to be overloaded 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 +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. diff --git a/doc/sphinx/user/extending/plugin-types/velocity-bc.md b/doc/sphinx/user/extending/plugin-types/velocity-bc.md index f1b2b2ab7e5..4b434f68895 100644 --- a/doc/sphinx/user/extending/plugin-types/velocity-bc.md +++ b/doc/sphinx/user/extending/plugin-types/velocity-bc.md @@ -1,5 +1,7 @@ (sec:extending:plugin-types:velocity-bc)= -# Prescribed velocity boundary conditions +# Velocity boundary conditions + +## 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 @@ -16,50 +18,38 @@ 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](https://aspect.geodynamics.org/doc/doxygen/classaspect_1_1BoundaryVelocity_1_1Interface.html) class and use the -`ASPECT_REGISTER_VELOCITY_BOUNDARY_CONDITIONS` macro to register your new +`ASPECT_REGISTER_BOUNDARY_VELOCITY_MODEL` 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: - -```{code-block} c++ -template - class aspect::VelocityBoundaryConditions::Interface - { - public: - virtual - Tensor<1,dim> - boundary_velocity (const Point &position) const = 0; - - virtual - void - initialize (const GeometryModel::Interface &geometry_model); - - virtual - void - update (); - - static - void - declare_parameters (ParameterHandler &prm); - - virtual - void - parse_parameters (ParameterHandler &prm); - }; -``` +`aspect::BoundaryVelocity`. -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 +In essence, the main function you have to implement for this plugin system is one +that, given a point returns the prescribed velocity at that point. +There are also member functions in the base class you can overload that are called +at the beginning of each time step; this is useful if one needs to perform some +operation once for 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 +The remaining functions are discussed in the documentation of this interface class at [aspect::VelocityBoundaryConditions::Interface](https://aspect.geodynamics.org/doc/doxygen/classaspect_1_1BoundaryVelocity_1_1Interface.html). -The purpose of the last two -functions has been discussed in the general overview of plugins above. + + +## Prescribed traction boundary conditions + +Alternatively, at a boundary one can prescribe the traction (i.e., a force density) +that *drives* the velocity, rather than the velocity itself. As for prescribed +velocities, a plugin system allows to describe this information, with the key +function being one that for a given point returns the prescribed traction. + +To implement a new boundary traction model, you need to overload the +[aspect::BoundaryTraction::Interface](https://aspect.geodynamics.org/doc/doxygen/classaspect_1_1BoundaryTraction_1_1Interface.html) +class and use the +`ASPECT_REGISTER_BOUNDARY_TRACTION_MODEL` macro to register your new +class. The implementation of the new class should be in namespace +`aspect::BoundaryTraction`. +The member functions that can be overloaded are discussed in the +documentation of this interface class at +[aspect::BoundaryTraction::Interface](https://aspect.geodynamics.org/doc/doxygen/classaspect_1_1BoundaryTraction_1_1Interface.html). diff --git a/doc/sphinx/user/extending/plugin-types/vis-postprocessors.md b/doc/sphinx/user/extending/plugin-types/vis-postprocessors.md index 871c995df58..9b59070cc2e 100644 --- a/doc/sphinx/user/extending/plugin-types/vis-postprocessors.md +++ b/doc/sphinx/user/extending/plugin-types/vis-postprocessors.md @@ -7,17 +7,17 @@ already implemented in ASPECT is the outputs it as a collection of files that can then be visualized graphically, see {ref}`sec:run-aspect:visualizing-results`. 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. +velocity, pressure, temperature, and any compositional fields there may be; on the other hand, we are frequently +interested in *derived*, spatially and temporally variable quantities such as +the viscosity, temperature and strain rate at a given +location, the difference between the actual and the adiabatic pressure, 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, +furthermore not scale very well as the list of things that can be output 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 +computed and output is implemented in a number of plugins, and what *is* computed is selected in the input parameter file (see {ref}`parameters:Postprocess/Visualization`). Defining visualization postprocessors works in much the same way as for the @@ -48,7 +48,7 @@ Visualization plugins can come in two flavors: class (and possibly the [SimulatorAccess][aspect::SimulatorAccess class] 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 + classes can be thought of as transformers: 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 @@ -70,53 +70,14 @@ Visualization plugins can come in two flavors: `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: - - ```{code-block} c++ - class Interface - { - public: - static - void - declare_parameters (ParameterHandler &prm); - - virtual - void - parse_parameters (ParameterHandler &prm); - - virtual - void save (std::map &status_strings) const; - - virtual - void load (const std::map &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: - - ```{code-block} c++ - class dealii::DataPostprocessorScalar - { - public: - virtual - void - compute_derived_quantities_vector - (const std::vector > &uh, - const std::vector > > &duh, - const std::vector > > &dduh, - const std::vector > &normals, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const; - }; - ``` - - What this function does is described in detail in the deal.II + `dealii::DataPostprocessorVector` classes. You will want to read up + on the specifics of what function you have to implement in the deal.II documentation. In addition, one has to write a suitable constructor to - call `dealii::DataPostprocessorScalar::DataPostprocessorScalar`. + call `dealii::DataPostprocessorScalar::DataPostprocessorScalar` or + the constructor other deal.II classes. - *Plugins that compute things from the solution in a cell-wise way:* The second possibility is for a class to not derive from