Postprocess

Postprocess#

Subsection: Postprocess#

Parameter name: List of postprocessors#

Default value:

Documentation: A comma separated list of postprocessor objects that should be run at the end of each time step. Some of these postprocessors will declare their own parameters which may, for example, include that they will actually do something only every so many time steps or years. Alternatively, the text ‘all’ indicates that all available postprocessors should be run after each time step.

The following postprocessors are available:

‘Stokes residual’: A postprocessor that outputs the Stokes residuals during the iterative solver algorithm into a file stokes_residuals.txt in the output directory.

‘basic statistics’: A postprocessor that outputs some simplified statistics like the Rayleigh number and other quantities that only make sense in certain model setups. The output is written after completing initial adaptive refinement steps. The postprocessor assumes a point at the surface at the adiabatic surface temperature and pressure is a reasonable reference condition for computing these properties. Furthermore, the Rayleigh number is computed using the model depth (i.e. not the radius of the Earth), as we need a definition that is geometry independent. Take care when comparing these values to published studies and make sure they use the same definitions.

‘boundary densities’: A postprocessor that computes the laterally averaged density at the top and bottom of the domain.

‘boundary pressures’: A postprocessor that computes the laterally averaged pressure at the top and bottom of the domain.

‘boundary strain rate residual statistics’: A postprocessor that computes some statistics about the surface strain rate residual along the top boundary. The residual is the difference between the second invariant of the model strain rate and the second strain rate invariant read from the input data file. Currently, the strain residual statistics, i.e., min, max and the rms magnitude, are computed at the top surface.

‘boundary velocity residual statistics’: A postprocessor that computes some statistics about the velocity residual along the top boundary. The velocity residual is the difference between the model solution velocities and the input velocities (GPlates model or ascii data). Currently, the velocity residual statistics, i.e., min, max and the rms magnitude, is computed at the top surface.

‘command’: A postprocessor that executes a command line process.

‘composition statistics’: A postprocessor that computes some statistics about the compositional fields, if present in this simulation. In particular, it computes maximal and minimal values of each field, as well as the total mass contained in this field as defined by the integral $m_i(t) = \int_\Omega c_i(\mathbf x,t) \; \text{d}x$.

‘composition velocity statistics’: A postprocessor that computes the root mean square velocity over the area spanned by each compositional field (i.e. where the field values are larger or equal to 0.5.

‘core statistics’: A postprocessor that computes some statistics about the core evolution. (Working only with dynamic core boundary temperature plugin)

‘depth average’: A postprocessor that computes depth averaged quantities and writes them into a file <depth_average.ext> in the output directory, where the extension of the file is determined by the output format you select. In addition to the output format, a number of other parameters also influence this postprocessor, and they can be set in the section Postprocess/Depth average in the input file.

In the output files, the $x$-value of each data point corresponds to the depth, whereas the $y$-value corresponds to the simulation time. The time is provided in seconds or, if the global “Use years in output instead of seconds” parameter is set, in years.

‘domain volume statistics’: A postprocessor that computes the total area (in 2d) or volume (in 3d) of the computational domain.

‘dynamic topography’: A postprocessor that computes a measure of dynamic topography based on the stress at the surface and bottom. The data is written into text files named ‘dynamic_topography.NNNNN’ in the output directory, where NNNNN is the number of the time step.

The exact approach works as follows: At the centers of all cells that sit along the top surface, we evaluate the stress and evaluate the component of it in the direction in which gravity acts. In other words, we compute $\sigma_{rr}={\hat g}^T(2 \eta \varepsilon(\mathbf u)- \frac 13 (\textrm{div}\;\mathbf u)I)\hat g - p_d$ where $\hat g = \mathbf g/\|\mathbf g\|$ is the direction of the gravity vector $\mathbf g$ and $p_d=p-p_a$ is the dynamic pressure computed by subtracting the adiabatic pressure $p_a$ from the total pressure $p$ computed as part of the Stokes solve. From this, the dynamic topography is computed using the formula $h=\frac{\sigma_{rr}}{(\mathbf g \cdot \mathbf n) \rho}$ where $\rho$ is the density at the cell center. For the bottom surface we chose the convection that positive values are up (out) and negative values are in (down), analogous to the deformation of the upper surface. Note that this implementation takes the direction of gravity into account, which means that reversing the flow in backward advection calculations will not reverse the instantaneous topography because the reverse flow will be divided by the reverse surface gravity. The file format then consists of lines with Euclidean coordinates followed by the corresponding topography value.

(As a side note, the postprocessor chooses the cell center instead of the center of the cell face at the surface, where we really are interested in the quantity, since this often gives better accuracy. The results should in essence be the same, though.)

‘entropy viscosity statistics’: A postprocessor that computes the maximum and volume averagedentropy viscosity stabilization for the temperature field.

‘geoid’: A postprocessor that computes a representation of the geoid based on the density structure in the mantle, as well as the topography at the surface and core mantle boundary (CMB) if desired. The topography is based on the dynamic topography postprocessor in case of no free surface, and based on the real surface from the geometry model in case of a free surface. The geoid is computed from a spherical harmonic expansion, so the geometry of the domain must be a 3d spherical shell.

‘global statistics’: A postprocessor that outputs all the global statistics information, e.g. the time of the simulation, the timestep number, number of degrees of freedom and solver iterations for each timestep. The postprocessor can output different formats, the first printing one line in the statistics file per nonlinear solver iteration (if a nonlinear solver scheme is selected). The second prints one line per timestep, summing the information about all nonlinear iterations in this line. Note that this postprocessor is always active independent on whether or not it is selected in the parameter file.

‘gravity calculation’: A postprocessor that computes gravity, gravity anomalies, gravity potential and gravity gradients for a set of points (e.g. satellites) in or above the model surface for either a user-defined range of latitudes, longitudes and radius or a list of point coordinates.Spherical coordinates in the output file are radius, colatitude and colongitude. Gravity is here based on the density distribution from the material model (and non adiabatic). This means that the density may come directly from an ascii file. This postprocessor also computes theoretical gravity and its derivatives, which corresponds to the analytical solution of gravity in the same geometry but filled with a reference density. The reference density is also used to determine density anomalies for computing gravity anomalies. Thus one must carefully evaluate the meaning of the gravity anomaly output, because the solution may not reflect the actual gravity anomaly (due to differences in the assumed reference density). On way to guarantee correct gravity anomalies is to subtract gravity of a certain point from the average gravity on the map. Another way is to directly use density anomalies for this postprocessor.The average- minimum- and maximum gravity acceleration and potential are written into the statistics file.

‘heat flux densities’: A postprocessor that computes some statistics about the heat flux density for each boundary id. The heat flux density across each boundary is computed in outward direction, i.e., from the domain to the outside. The heat flux is computed as sum of advective heat flux and conductive heat flux through Neumann boundaries, both computed as integral over the boundary area, and conductive heat flux through Dirichlet boundaries, which is computed using the consistent boundary flux method as described in “Gresho, Lee, Sani, Maslanik, Eaton (1987). The consistent Galerkin FEM for computing derived boundary quantities in thermal and or fluids problems. International Journal for Numerical Methods in Fluids, 7(4), 371-394.”

Note that the “heat flux statistics” postprocessor computes the same quantity as the one here, but not divided by the area of the surface. In other words, it computes the total heat flux through each boundary.

‘heat flux map’: A postprocessor that computes the heat flux density across each boundary in outward direction, i.e., from the domain to the outside. The heat flux is computed as sum of advective heat flux and conductive heat flux through Neumann boundaries, both computed as integral over the boundary area, and conductive heat flux through Dirichlet boundaries, which is computed using the consistent boundary flux method as described in “Gresho, Lee, Sani, Maslanik, Eaton (1987). The consistent Galerkin FEM for computing derived boundary quantities in thermal and or fluids problems. International Journal for Numerical Methods in Fluids, 7(4), 371-394.”

‘heat flux statistics’: A postprocessor that computes some statistics about the heat flux density across each boundary in outward direction, i.e., from the domain to the outside. The heat flux is computed as sum of advective heat flux and conductive heat flux through Neumann boundaries, both computed as integral over the boundary area, and conductive heat flux through Dirichlet boundaries, which is computed using the consistent boundary flux method as described in “Gresho, Lee, Sani, Maslanik, Eaton (1987). The consistent Galerkin FEM for computing derived boundary quantities in thermal and or fluids problems. International Journal for Numerical Methods in Fluids, 7(4), 371-394.”The point-wise heat flux can be obtained from the heat flux map postprocessor, which outputs the heat flux to a file, or the heat flux map visualization postprocessor, which outputs the heat flux for visualization.

As stated, this postprocessor computes the outbound heat flux. If you are interested in the opposite direction, for example from the core into the mantle when the domain describes the mantle, then you need to multiply the result by -1.

\note{In geodynamics, the term “heat flux” is often understood to be the quantity $- k \nabla T$, which is really a heat flux density, i.e., a vector-valued field. In contrast to this, the current postprocessor only computes the integrated flux over each part of the boundary. Consequently, the units of the quantity computed here are $W=\frac{J}{s}$.}

The “heat flux densities” postprocessor computes the same quantity as the one here, but divided by the area of the surface.

‘heating statistics’: A postprocessor that computes some statistics about heating, averaged by volume.

‘load balance statistics’: A postprocessor that computes statistics about the distribution of cells, and if present particles across subdomains. In particular, it computes maximal, average and minimal number of cells across all ranks. If there are particles it also computes the maximal, average, and minimum number of particles across all ranks, and maximal, average, and minimal ratio between local number of particles and local number of cells across all processes. All of these numbers can be useful to assess the load balance between different MPI ranks, as the difference between the minimal and maximal load should be as small as possible.

‘mass flux statistics’: A postprocessor that computes some statistics about the mass flux across boundaries. For each boundary indicator (see your geometry description for which boundary indicators are used), the mass flux is computed in outward direction, i.e., from the domain to the outside, using the formula $\int_{\Gamma_i} \rho \mathbf v \cdot \mathbf n$ where $\Gamma_i$ is the part of the boundary with indicator $i$, $\rho$ is the density as reported by the material model, $\mathbf v$ is the velocity, and $\mathbf n$ is the outward normal.

As stated, this postprocessor computes the outbound mass flux. If you are interested in the opposite direction, for example from the core into the mantle when the domain describes the mantle, then you need to multiply the result by -1.

\note{In geodynamics, the term “mass flux” is often understood to be the quantity $\rho \mathbf v$, which is really a mass flux density, i.e., a vector-valued field. In contrast to this, the current postprocessor only computes the integrated flux over each part of the boundary. Consequently, the units of the quantity computed here are $\frac{kg}{s}$.}

‘material statistics’: A postprocessor that computes some statistics about the material properties. In particular, it computes the volume-averages of the density and viscosity, and the total mass in the model. Specifically, it implements the following formulas in each time step: $\left<\rho\right> = \frac{1}{|\Omega|} \int_\Omega \rho(\mathbf x) \, \text{d}x$, $\left<\eta\right> = \frac{1}{|\Omega|} \int_\Omega \eta(\mathbf x) \, \text{d}x$, $M = \int_\Omega \rho(\mathbf x) \, \text{d}x$, where $|\Omega|$ is the volume of the domain.

‘matrix statistics’: A postprocessor that computes some statistics about the matrices. In particular, it outputs total memory consumption, total non-zero elements, and non-zero elements per block, for system matrix and system preconditioner matrix.

‘maximum depth of field’: A postprocessor that for each compositional field outputs the largest depth at which a quadrature point is found where the field has a value of 0.5 or larger. For fields that do not represent materials, but for example track a certain quantity like strain, this value of 0.5 does not necessarily make sense.

‘melt statistics’: A postprocessor that computes some statistics about the melt (volume) fraction. If the material model does not implement a melt fraction function, the output is set to zero.

‘memory statistics’: A postprocessor that computes some statistics about the memory consumption. In particular, it computes the memory usage of the system matrix, triangulation, p4est, DoFHandler, current constraints, solution vector, and peak virtual memory usage, all in MB. It also outputs the memory usage of the system matrix to the screen.

‘mobility statistics’: A postprocessor that computes some statistics about mobility following Tackley (2000) and Lourenco et al. (2020).

‘particle count statistics’: A postprocessor that computes some statistics about the particle distribution, if present in this simulation. In particular, it computes minimal, average and maximal values of particles per cell in the global domain.

‘particles’: A Postprocessor that creates particles that follow the velocity field of the simulation. The particles can be generated and propagated in various ways and they can carry a number of constant or time-varying properties. The postprocessor can write output positions and properties of all particles at chosen intervals, although this is not mandatory. It also allows other parts of the code to query the particles for information.

‘point values’: A postprocessor that evaluates the solution (i.e., velocity, pressure, temperature, and compositional fields along with other fields that are treated as primary variables) at the end of every time step or after a user-specified time interval at a given set of points and then writes this data into the file <point_values.txt> in the output directory. The points at which the solution should be evaluated are specified in the section Postprocess/Point values in the input file.

In the output file, data is organized as (i) time, (ii) the 2 or 3 coordinates of the evaluation points, and (iii) followed by the values of the solution vector at this point. The time is provided in seconds or, if the global “Use years in output instead of seconds” parameter is set, in years. In the latter case, the velocity is also converted to meters/year, instead of meters/second.

\note{Evaluating the solution of a finite element field at arbitrarily chosen points is an expensive process. Using this postprocessor will only be efficient if the number of evaluation points or output times is relatively small. If you need a very large number of evaluation points, you should consider extracting this information from the visualization program you use to display the output of the ‘visualization’ postprocessor.}

‘pressure statistics’: A postprocessor that computes some statistics about the pressure field.

‘rotation statistics’: A postprocessor that computes some statistics about the rotational velocity of the model (i.e. integrated net rotation and angular momentum). In 2d we assume the model to be a cross-section through an infinite domain in z direction, with a zero z-velocity. Thus, the z-axis is the only possible rotation axis and both moment of inertia and angular momentum are scalar instead of tensor quantities.

‘spherical velocity statistics’: A postprocessor that computes radial, tangential and total RMS velocity.

‘temperature statistics’: A postprocessor that computes some statistics about the temperature field.

‘topography’: A postprocessor intended for use with a deforming top surface. After every step it loops over all the vertices on the top surface and determines the maximum and minimum topography relative to a reference datum (initial box height for a box geometry model or initial radius for a sphere/spherical shell geometry model). If ’Topography.Output to file’ is set to true, also outputs topography into text files named ‘topography.NNNNN’ in the output directory, where NNNNN is the number of the time step. The file format then consists of lines with Euclidean coordinates followed by the corresponding topography value.Topography is printed/written in meters.

‘velocity boundary statistics’: A postprocessor that computes some statistics about the velocity along the boundaries. For each boundary indicator (see your geometry description for which boundary indicators are used), the min and max velocity magnitude is computed.

‘velocity statistics’: A postprocessor that computes the root mean square and maximum velocity in the computational domain.

‘viscous dissipation statistics’: A postprocessor that outputs the viscous rate of dissipation of energy for each compositional field (where the field has a value of 0.5 or more) as well as over the whole domain. When all the fields represent lithologies and there is no background field, the sum of the individual field’s dissipation should equal that over the whole domain. The viscous dissipation is computed as: $\int_{V}\left(\sigma’ \dot{\epsilon}’ \right)$, where $\sigma’$ is the deviatoric stress and $\dot{\epsilon}’$ the deviatoric strain rate.Note then when shear heating is included in the temperature equation, it is better to use the ’heating statistics’ postprocessor.

‘visualization’: A postprocessor that takes the solution and writes it into files that can be read by a graphical visualization program. Additional run time parameters are read from the parameter subsection ’Visualization’.

‘volume of fluid statistics’: A postprocessor that computes some statistics about the volume-of-fluid fields.

Parameter name: Run postprocessors on nonlinear iterations#

Default value: false

Pattern: [Bool]

Documentation: Whether or not the postprocessors should be executed after each of the nonlinear iterations done within one time step. As this is mainly an option for the purposes of debugging, it is not supported when the ’Time between graphical output’ is larger than zero, or when the postprocessor is not intended to be run more than once per timestep.

Subsection: Postprocess / Boundary strain rate residual statistics#

Parameter name: Data directory#

Default value: $ASPECT_SOURCE_DIR/data/postprocess/boundary-strain-rate-residual/

Pattern: [DirectoryName]

Documentation: The name of a directory that contains the ascii data. This path may either be absolute (if starting with a ‘/’) or relative to the current directory. The path may also include the special text ‘$ASPECT_SOURCE_DIR’ which will be interpreted as the path in which the ASPECT source files were located when ASPECT was compiled. This interpretation allows, for example, to reference files located in the ‘data/’ subdirectory of ASPECT.

Parameter name: Data file name#

Default value: box_3d_boundary_strain_rate.txt

Pattern: [Anything]

Documentation: The file name of the input surface strain rate an ascii data. The file has one column in addition to the coordinate columns corresponding to the second invariant of strain rate.

Parameter name: Scale factor#

Default value: 1.

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Scalar factor, which is applied to the model data. You might want to use this to scale the input to a reference model.

Subsection: Postprocess / Boundary velocity residual statistics#

Parameter name: Data directory#

Default value: $ASPECT_SOURCE_DIR/data/boundary-velocity/gplates/

Pattern: [DirectoryName]

Documentation: The name of a directory that contains the GPlates model or the ascii data. This path may either be absolute (if starting with a ‘/’) or relative to the current directory. The path may also include the special text ‘$ASPECT_SOURCE_DIR’ which will be interpreted as the path in which the ASPECT source files were located when ASPECT was compiled. This interpretation allows, for example, to reference files located in the ‘data/’ subdirectory of ASPECT.

Parameter name: Data file name#

Default value: current_day.gpml

Pattern: [Anything]

Documentation: The file name of the input velocity as a GPlates model or an ascii data. For the GPlates model, provide file in the same format as described in the ’gplates’ boundary velocity plugin. For the ascii data, provide file in the same format as described in ’ascii data’ initial composition plugin.

Parameter name: Scale factor#

Default value: 1.

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Scalar factor, which is applied to the model data. You might want to use this to scale the input to a reference model. Another way to use this factor is to convert units of the input files. For instance, if you provide velocities in cm/year set this factor to 0.01.

Parameter name: Use ascii data#

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: The number of zones in depth direction within which we are to compute averages. By default, we subdivide the entire domain into 10 depth zones and compute temperature and other averages within each of these zones. However, if you have a very coarse mesh, it may not make much sense to subdivide the domain into so many zones and you may wish to choose less than this default. It may also make computations slightly faster. On the other hand, if you have an extremely highly resolved mesh, choosing more zones might also make sense.

Parameter name: Output format#

Default value: gnuplot, txt

Pattern: [MultipleSelection none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|tecplot_binary|vtk|vtu|hdf5|svg|deal.II intermediate|txt ]

Documentation: A list of formats in which the output shall be produced. The format in which the output is generated also determines the extension of the file into which data is written. The list of possible output formats that can be given here is documented in the appendix of the manual where the current parameter is described. By default the output is written as gnuplot file (for plotting), and as a simple text file.

Parameter name: Time between graphical output#

Default value: 1e8

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: The time interval between each generation of graphical output files. A value of zero indicates that output should be generated in each time step. Units: years if the ’Use years in output instead of seconds’ parameter is set; seconds otherwise.

Subsection: Postprocess / Dynamic core statistics#

Parameter name: Excess entropy only#

Default value: false

Pattern: [Bool]

Documentation: Output the excess entropy only instead the each entropy terms.

Subsection: Postprocess / Dynamic topography#

Parameter name: Density above#

Default value: 0.

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: Dynamic topography is calculated as the excess or lack of mass that is supported by mantle flow. This value depends on the density of material that is moved up or down, i.e. crustal rock, and the density of the material that is displaced (generally water or air). While the density of crustal rock is part of the material model, this parameter ‘Density above’ allows the user to specify the density value of material that is displaced above the solid surface. By default this material is assumed to be air, with a density of 0. Units: \si{\kilogram\per\meter\cubed}.

Parameter name: Density below#

Default value: 9900.

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: Dynamic topography is calculated as the excess or lack of mass that is supported by mantle flow. This value depends on the density of material that is moved up or down, i.e. mantle above CMB, and the density of the material that is displaced (generally outer core material). While the density of mantle rock is part of the material model, this parameter ‘Density below’ allows the user to specify the density value of material that is displaced below the solid surface. By default this material is assumed to be outer core material with a density of 9900. Units: \si{\kilogram\per\meter\cubed}.

Parameter name: Output bottom#

Default value: true

Pattern: [Bool]

Documentation: Whether to output a file containing the bottom (i.e., CMB) dynamic topography.

Parameter name: Output surface#

Default value: true

Pattern: [Bool]

Parameter name: Write statistics for each nonlinear iteration#

Default value: false

Pattern: [Bool]

Documentation: Whether to put every nonlinear iteration into a separate line in the statistics file (if true), or to output only one line per time step that contains the total number of iterations of the Stokes and advection linear system solver.

Subsection: Postprocess / Gravity calculation#

Parameter name: List of latitude#

Default value:

Pattern: [List of <[Double -90…90 (inclusive)]> of length 0…4294967295 (inclusive)]

Documentation: Parameter for the list of points sampling scheme: List of satellite latitude coordinates.

Parameter name: List of longitude#

Default value:

Pattern: [List of <[Double -180…180 (inclusive)]> of length 0…4294967295 (inclusive)]

Documentation: Parameter for the list of points sampling scheme: List of satellite longitude coordinates.

Parameter name: List of radius#

Default value:

Pattern: [List of <[Double 0…MAX_DOUBLE (inclusive)]> of length 0…4294967295 (inclusive)]

Documentation: Parameter for the list of points sampling scheme: List of satellite radius coordinates. Just specify one radius if all points values have the same radius. If not, make sure there are as many radius as longitude and latitude

Parameter name: Maximum latitude#

Default value: 90

Pattern: [Double -90…90 (inclusive)]

Documentation: Parameter for the uniform distribution sampling scheme: Gravity may be calculated for a sets of points along the latitude between a minimum and maximum latitude.

Parameter name: Maximum longitude#

Default value: 180.

Pattern: [Double -180…180 (inclusive)]

Documentation: Parameter for the uniform distribution sampling scheme: Gravity may be calculated for a sets of points along the longitude between a minimum and maximum longitude.

Parameter name: Maximum radius#

Default value: 0.

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: Parameter for the map sampling scheme: Maximum radius can be defined in or outside the model.

Parameter name: Minimum latitude#

Default value: -90.

Pattern: [Double -90…90 (inclusive)]

Documentation: Parameter for the uniform distribution sampling scheme: Gravity may be calculated for a sets of points along the latitude between a minimum and maximum latitude.

Parameter name: Minimum longitude#

Default value: -180.

Pattern: [Double -180…180 (inclusive)]

Documentation: Parameter for the uniform distribution sampling scheme: Gravity may be calculated for a sets of points along the longitude between a minimum and maximum longitude.

Parameter name: Minimum radius#

Default value: 0.

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: Parameter for the map sampling scheme: Minimum radius may be defined in or outside the model. Prescribe a minimum radius for a sampling coverage at a specific height.

Parameter name: Number points fibonacci spiral#

Default value: 200

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Parameter for the fibonacci spiral sampling scheme: This specifies the desired number of satellites per radius layer. The default value is 200. Note that sampling becomes more uniform with increasing number of satellites

Parameter name: Number points latitude#

Default value: 1

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Parameter for the map sampling scheme: This specifies the number of points along the latitude (e.g. gravity map) between a minimum and maximum latitude.

Parameter name: Number points longitude#

Default value: 1

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Parameter for the map sampling scheme: This specifies the number of points along the longitude (e.g. gravity map) between a minimum and maximum longitude.

Parameter name: Number points radius#

Default value: 1

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Parameter for the map sampling scheme: This specifies the number of points along the radius (e.g. depth profile) between a minimum and maximum radius.

Parameter name: Precision in gravity output#

Default value: 12

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: Set the precision of gravity acceleration, potential and gradients in the gravity output and statistics file.

Parameter name: Quadrature degree increase#

Default value: 0

Pattern: [Integer range -1…2147483647 (inclusive)]

Documentation: Quadrature degree increase over the velocity element degree may be required when gravity is calculated near the surface or inside the model. An increase in the quadrature element adds accuracy to the gravity solution from noise due to the model grid.

Parameter name: Reference density#

Default value: 3300.

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: Gravity anomalies may be computed using density anomalies relative to a reference density.

Parameter name: Sampling scheme#

Default value: map

Pattern: [Selection map|list|list of points|fibonacci spiral ]

Documentation: Choose the sampling scheme. By default, the map produces a grid of equally angled points between a minimum and maximum radius, longitude, and latitude. A list of points contains the specific coordinates of the satellites. The fibonacci spiral sampling scheme produces a uniformly distributed map on the surface of sphere defined by a minimum and/or maximum radius.

Parameter name: Time between gravity output#

Default value: 1e8

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: The time interval between each generation of gravity output files. A value of 0 indicates that output should be generated in each time step. Units: years if the ’Use years in output instead of seconds’ parameter is set; seconds otherwise.

Parameter name: Time steps between gravity output#

Default value: 2147483647

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: The maximum number of time steps between each generation of gravity output files.

Subsection: Postprocess / Memory statistics#

Parameter name: Output peak virtual memory (VmPeak)#

Default value: true

Pattern: [Bool]

Documentation: If set to ’true’, also output the peak virtual memory usage (computed as the maximum over all processors).

Subsection: Postprocess / Particles#

Parameter name: Allow cells without particles#

Default value: false

Pattern: [Bool]

Documentation: By default, every cell needs to contain particles to use this interpolator plugin. If this parameter is set to true, cells are allowed to have no particles. In case both the current cell and its neighbors are empty, the interpolator will return 0 for the current cell’s properties.

Parameter name: Data output format#

Default value: vtu

Pattern: [MultipleSelection none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|tecplot_binary|vtk|vtu|hdf5|svg|deal.II intermediate|ascii ]

Documentation: A comma separated list of file formats to be used for graphical output. The list of possible output formats that can be given here is documented in the appendix of the manual where the current parameter is described.

Parameter name: Exclude output properties#

Default value:

Pattern: [Anything]

Documentation: A comma separated list of particle properties that should not be output. If this list contains the entry ‘all’, only the id of particles will be provided in graphical output files.

Parameter name: Integration scheme#

Default value: rk2

Pattern: [Selection euler|rk2|rk4 ]

Documentation: This parameter is used to decide which method to use to solve the equation that describes the position of particles, i.e., $\frac{d}{dt}\mathbf x_k(t) = \mathbf u(\mathbf x_k(t),t)$, where $k$ is an index that runs over all particles, and $\mathbf u(\mathbf x,t)$ is the velocity field that results from the Stokes equations.

In practice, the exact velocity $\mathbf u(\mathbf x,t)$ is of course not available, but only a numerical approximation $\mathbf u_h(\mathbf x,t)$. Furthermore, this approximation is only available at discrete time steps, $\mathbf u^n(\mathbf x)=\mathbf u(\mathbf x,t^n)$, and these need to be interpolated between time steps if the integrator for the equation above requires an evaluation at time points between the discrete time steps. If we denote this interpolation in time by $\tilde{\mathbf u}_h(\mathbf x,t)$ where $\tilde{\mathbf u}_h(\mathbf x,t^n)=\mathbf u^n(\mathbf x)$, then the equation the differential equation solver really tries to solve is $\frac{d}{dt}\tilde{\mathbf x}_k(t) = \tilde{\mathbf u}_h(\mathbf x_k(t),t)$.

As a consequence of these considerations, if you try to assess convergence properties of an ODE integrator – for example to verify that the RK4 integrator converges with fourth order –, it is important to recall that the integrator may not solve the equation you think it solves. If, for example, we call the numerical solution of the ODE $\tilde{\mathbf x}_{k,h}(t)$, then the error will typically satisfy a relationship like [ | \tilde{\mathbf x}k(T) - \tilde{\mathbf x}{k,h}(T) | \le C(T) \Delta t^p] where $\Delta t$ is the time step and $p$ the convergence order of the method, and $C(T)$ is a (generally unknown) constant that depends on the end time $T$ at which one compares the solutions. On the other hand, an analytically computed trajectory would likely use the exact velocity, and one may be tempted to compute $\| \mathbf x_k(T) - \tilde{\mathbf x}_{k,h}(T) \|$, but this quantity will, in the best case, only satisfy an estimate of the form [ | \mathbf x_k(T) - \tilde{\mathbf x}_{k,h}(T) | \le C_1(T) \Delta t^p + C_2(T) | \mathbf u-\mathbf u_h | + C_3(T) | \mathbf u_h-\tilde{\mathbf u}_h |] with appropriately chosen norms for the second and third term. These second and third terms typically converge to zero at relatively low rates (compared to the order $p$ of the integrator, which can often be chosen relatively high) in the mesh size $h$ and the time step size $\\Delta t$, limiting the overall accuracy of the ODE integrator.

Select one of the following models:

‘euler’: Explicit Euler scheme integrator, where $y_{n+1} = y_n + \Delta t \, v(y_n)$. This requires only one integration substep per timestep.

‘rk2’: Second Order Runge Kutta integrator $y_{n+1} = y_n + \Delta t\, v(t_{n+1/2}, y_{n} + \frac{1}{2} k_1)$ where $k_1 = \Delta t\, v(t_{n}, y_{n})$

‘rk4’: Runge Kutta fourth order integrator, where $y_{n+1} = y_n + \frac{1}{6} k_1 + \frac{1}{3} k_2 + \frac{1}{3} k_3 + \frac{1}{6} k_4$ and $k_1$, $k_2$, $k_3$, $k_4$ are defined as usual.

Parameter name: Interpolation scheme#

Default value: cell average

Pattern: [Selection bilinear least squares|cell average|harmonic average|nearest neighbor|quadratic least squares ]

Documentation: Select one of the following models:

‘bilinear least squares’: Uses linear least squares to obtain the slopes and center of a 2d or 3d plane from the particle positions and a particular property value on those particles. Interpolate this property onto a vector of points. If the limiter is enabled then it will ensure the interpolated properties do not exceed the range of the minimum and maximum of the values of the property on the particles. Note that deal.II must be configured with BLAS and LAPACK to support this operation.

‘cell average’: Return the arithmetic average of all particle properties in the given cell, or in the neighboring cells if the given cell is empty. In case the neighboring cells are also empty, and ’Allow cells without particles’ is set to true, the interpolator returns 0. Otherwise, an exception is thrown.

‘harmonic average’: Return the harmonic average of all particle properties in the given cell. If the cell contains no particles, return the harmonic average of the properties in the neighboring cells. In case the neighboring cells are also empty, and ’Allow cells without particles’ is set to true, the interpolator returns 0. Otherwise, an exception is thrown.

‘nearest neighbor’: Return the properties of the nearest neighboring particle in the current cell, or nearest particle in nearest neighboring cell if current cell is empty. In case the neighboring cells are also empty, and ’Allow cells without particles’ is set to true, the interpolator returns 0. Otherwise, an exception is thrown.

‘quadratic least squares’: Interpolates particle properties onto a vector of points using a quadratic least squares method. Note that deal.II must be configured with BLAS/LAPACK.

Parameter name: List of particle properties#

Default value:

Documentation: A comma separated list of particle properties that should be tracked. By default none is selected, which means only position, velocity and id of the particles are output.

The following properties are available:

‘composition’: Implementation of a plugin in which the particle property is defined by the compositional fields in the model. This can be used to track solid compositionevolution over time.

‘cpo bingham average’: This is a particle property plugin which computes the Bingham average for the Crystal Preferred Orientation particle property plugin so that it can be visualized.

‘crystal preferred orientation’: The plugin manages and computes the evolution of Lattice/Crystal Preferred Orientations (LPO/CPO) on particles. Each ASPECT particle can be assigned many grains. Each grain is assigned a size and a orientation matrix. This allows for CPO evolution tracking with polycrystalline kinematic CrystalPreferredOrientation evolution models such as D-Rex (Kaminski and Ribe, 2001; Kaminski et al., 2004).

‘elastic stress’: A plugin in which the particle property tensor is defined as the total elastic stress a particle has accumulated. See the viscoelastic material model documentation for more detailed information.

‘function’: Implementation of a model in which the particle property is set by evaluating an explicit function at the initial position of each particle. The function is defined in the parameters in section “Particles|Function”. The format of these functions follows the syntax understood by the muparser library, see Section~\ref{sec:muparser-format}.

‘grain size’: A plugin in which the particle property is defined as the evolving grain size of a particle. See the grain_size material model documentation for more detailed information.

‘initial composition’: Implementation of a plugin in which the particle property is given as the initial composition at the particle’s initial position. The particle gets as many properties as there are compositional fields.

‘initial position’: Implementation of a plugin in which the particle property is given as the initial position of the particle. This property is vector-valued with as many components as there are space dimensions. In practice, it is often most useful to only visualize one of the components of this vector, or the magnitude of the vector. For example, in a spherical mantle simulation, the magnitude of this property equals the starting radius of a particle, and is thereby indicative of which part of the mantle a particle comes from.

‘integrated strain’: A plugin in which the particle property tensor is defined as the deformation gradient tensor $\mathbf F$ this particle has experienced. $\mathbf F$ can be polar-decomposed into the left stretching tensor $\mathbf L$ (the finite strain we are interested in), and the rotation tensor $\mathbf Q$. See the corresponding cookbook in the manual for more detailed information.

‘integrated strain invariant’: A plugin in which the particle property is defined as the finite strain invariant ($\varepsilon_{ii}$). This property is calculated with the timestep ($dt$) and the second invariant of the deviatoric strain rate tensor ($\dot{\varepsilon}_{ii}$), where the value at time step $n$ is $\varepsilon_{ii}^{n} = \varepsilon_{ii}^{n-1} + dt\dot{\varepsilon}_{ii}$.

‘melt particle’: Implementation of a plugin in which the particle property is defined as presence of melt above a threshold, which can be set as an input parameter. This property is set to 0 if melt is not present and set to 1 if melt is present.

‘pT path’: Implementation of a plugin in which the particle property is defined as the current pressure and temperature at this position. This can be used to generate pressure-temperature paths of material points over time.

‘position’: Implementation of a plugin in which the particle property is defined as the current position.

‘reference position’: Implementation of a plugin in which the particle property is defined as the current reference position.

‘strain rate’: Implementation of a plugin in which the time evolution of strain rate is saved and stored on the particles.

‘velocity’: Implementation of a plugin in which the particle property is defined as the recent velocity at this position.

‘viscoplastic strain invariants’: A plugin that calculates the finite strain invariant a particle has experienced and assigns it to either the plastic and/or viscous strain field based on whether the material is plastically yielding, or the total strain field used in the visco plastic material model. The implementation of this property is equivalent to the implementation for compositional fields that is located in the plugin in benchmarks/buiter\_et\_al\_2008\_jgr/plugin/,and is effectively the same as what the visco plastic material model uses for compositional fields.

Parameter name: Load balancing strategy#

Default value: repartition

Pattern: [MultipleSelection none|remove particles|add particles|remove and add particles|repartition ]

Documentation: Strategy that is used to balance the computational load across processors for adaptive meshes.

Parameter name: Maximum particles per cell#

Default value: 100

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Upper limit for particle number per cell. This limit is useful for adaptive meshes to prevent coarse cells from slowing down the whole model. It will be checked and enforced after mesh refinement, after MPI transfer of particles and after particle movement. If there are n\_number\_of\_particles $>$ max\_particles\_per\_cell particles in one cell then n\_number\_of\_particles - max\_particles\_per\_cell particles in this cell are randomly chosen and destroyed.

Parameter name: Minimum particles per cell#

Default value: 0

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Lower limit for particle number per cell. This limit is useful for adaptive meshes to prevent fine cells from being empty of particles. It will be checked and enforced after mesh refinement and after particle movement. If there are n\_number\_of\_particles $<$ min\_particles\_per\_cell particles in one cell then min\_particles\_per\_cell - n\_number\_of\_particles particles are generated and randomly placed in this cell. If the particles carry properties the individual property plugins control how the properties of the new particles are initialized.

Parameter name: Number of grouped files#

Default value: 16

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: VTU file output supports grouping files from several CPUs into a given number of files using MPI I/O when writing on a parallel filesystem. Select 0 for no grouping. This will disable parallel file output and instead write one file per processor. A value of 1 will generate one big file containing the whole solution, while a larger value will create that many files (at most as many as there are MPI ranks).

Parameter name: Number of particles#

Default value: 1000

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: Total number of particles to create (not per processor or per element). The number is parsed as a floating point number (so that one can specify, for example, ’1e4’ particles) but it is interpreted as an integer, of course.

Parameter name: Particle generator name#

Default value: random uniform

Documentation: Select one of the following models:

‘ascii file’: Generates a distribution of particles from coordinates specified in an Ascii data file. The file format is a simple text file, with as many columns as spatial dimensions and as many lines as particles to be generated. Initial comment lines starting with ‘#’ will be discarded. Note that this plugin always generates as many particles as there are coordinates in the data file, the “Postprocess/Particles/Number of particles” parameter has no effect on this plugin. All of the values that define this generator are read from a section “Postprocess/Particles/Generator/Ascii file” in the input file, see Section~\ref{parameters:Postprocess/Particles/Generator/Ascii_20file}.

‘probability density function’: Generate a random distribution of particles over the entire simulation domain. The probability density is prescribed in the form of a user-prescribed function. The format of this function follows the syntax understood by the muparser library, see Section~\ref{sec:muparser-format}. The return value of the function is always checked to be a non-negative probability density but it can be zero in parts of the domain.

‘quadrature points’: Generates particles at the quadrature points of each active cell of the triangulation. Here, Gauss quadrature of degree (velocity_degree + 1), is used similarly to the assembly of Stokes matrix.

‘random uniform’: Generates a random uniform distribution of particles over the entire simulation domain.

‘reference cell’: Generates a uniform distribution of particles per cell and spatial direction in the unit cell and transforms each of the particles back to real region in the model domain. Uniform here means the particles will be generated with an equal spacing in each spatial dimension.

‘uniform box’: Generate a uniform distribution of particles over a rectangular domain in 2d or 3d. Uniform here means the particles will be generated with an equal spacing in each spatial dimension. Note that in order to produce a regular distribution the number of generated particles might not exactly match the one specified in the input file.

‘uniform radial’: Generate a uniform distribution of particles over a spherical domain in 2d or 3d. Uniform here means the particles will be generated with an equal spacing in each spherical spatial dimension, i.e., the particles are created at positions that increase linearly with equal spacing in radius, colatitude and longitude around a certain center point. Note that in order to produce a regular distribution the number of generated particles might not exactly match the one specified in the input file.

Parameter name: Particle weight#

Default value: 10

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: Weight that is associated with the computational load of a single particle. The sum of particle weights will be added to the sum of cell weights to determine the partitioning of the mesh if the ‘repartition’ particle load balancing strategy is selected. The optimal weight depends on the used integrator and particle properties. In general for a more expensive integrator and more expensive properties a larger particle weight is recommended. Before adding the weights of particles, each cell already carries a weight of 1000 to account for the cost of field-based computations.

Parameter name: Temporary output location#

Default value:

Pattern: [Anything]

Documentation: On large clusters it can be advantageous to first write the output to a temporary file on a local file system and later move this file to a network file system. If this variable is set to a non-empty string it will be interpreted as a temporary storage location.

Parameter name: Time between data output#

Default value: 1e8

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: The time interval between each generation of output files. A value of zero indicates that output should be generated every time step.

Units: years if the ’Use years in output instead of seconds’ parameter is set; seconds otherwise.

Parameter name: Update ghost particles#

Default value: false

Pattern: [Bool]

Documentation: Some particle interpolation algorithms require knowledge about particles in neighboring cells. To allow this, particles in ghost cells need to be exchanged between the processes neighboring this cell. This parameter determines whether this transport is happening.

Parameter name: Write in background thread#

Default value: false

Pattern: [Bool]

Documentation: File operations can potentially take a long time, blocking the progress of the rest of the model run. Setting this variable to ‘true’ moves this process into a background thread, while the rest of the model continues.

Subsection: Postprocess / Particles / CPO Bingham Average#

Parameter name: Number of samples#

Default value: 0

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: This determines how many samples are taken when using the random draw volume averaging. Setting it to zero means that the number of samples is set to be equal to the number of grains.

Parameter name: Random number seed#

Default value: 1

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: The seed used to generate random numbers. This will make sure that results are reproducible as long as the problem is run with the same amount of MPI processes. It is implemented as final seed = Random number seed + MPI Rank.

Subsection: Postprocess / Particles / Crystal Preferred Orientation#

Parameter name: CPO derivatives algorithm#

Default value: Spin tensor

Pattern: [List of <[Anything]> of length 0…4294967295 (inclusive)]

Documentation: Options: Spin tensor

Parameter name: Number of grains per particle#

Default value: 50

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: The number of grains of each different mineral each particle contains.

Parameter name: Property advection max iterations#

Default value: 100

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: The Backward Euler property advection method involve internal iterations. This option allows for setting the maximum number of iterations. Note that when the iteration is ended by the max iteration amount an assert is thrown.

Parameter name: Property advection method#

Default value: Backward Euler

Pattern: [Anything]

Documentation: Options: Forward Euler, Backward Euler

Parameter name: Property advection tolerance#

Default value: 1e-10

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: The Backward Euler property advection method involve internal iterations. This option allows for setting a tolerance. When the norm of tensor new - tensor old is smaller than this tolerance, the iteration is stopped.

Parameter name: Random number seed#

Default value: 1

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: The seed used to generate random numbers. This will make sure that results are reproducible as long as the problem is run with the same number of MPI processes. It is implemented as final seed = user seed + MPI Rank.

Subsection: Postprocess / Particles / Crystal Preferred Orientation / D-Rex 2004#

Parameter name: Exponents p#

Default value: 1.5

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: This is exponent p as defined in equation 11 of Kaminski et al., 2004.

Parameter name: Mobility#

Default value: 50

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: The dimensionless intrinsic grain boundary mobility for both olivine and enstatite.

Parameter name: Nucleation efficiency#

Default value: 5

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: This is the dimensionless nucleation rate as defined in equation 8 of Kaminski et al., 2004.

Parameter name: Stress exponents#

Default value: 3.5

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: This is the power law exponent that characterizes the rheology of the slip systems. It is used in equation 11 of Kaminski et al., 2004.

Parameter name: Threshold GBS#

Default value: 0.3

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Documentation: The Dimensionless Grain Boundary Sliding (GBS) threshold. This is a grain size threshold below which grain deform by GBS and become strain-free grains.

Parameter name: Volume fractions minerals#

Default value: 0.5, 0.5

Pattern: [List of <[Double 0…MAX_DOUBLE (inclusive)]> of length 0…4294967295 (inclusive)]

Documentation: The volume fraction for the different minerals. There need to be the same amount of values as there are minerals

Subsection: Postprocess / Particles / Crystal Preferred Orientation / Initial grains#

Parameter name: Minerals#

Default value: Olivine: Karato 2008, Enstatite

Pattern: [List of <[Anything]> of length 0…4294967295 (inclusive)]

Documentation: This determines what minerals and fabrics or fabric selectors are used used for the LPO/CPO calculation. The options are Olivine: Passive, A-fabric, Olivine: B-fabric, Olivine: C-fabric, Olivine: D-fabric, Olivine: E-fabric, Olivine: Karato 2008 or Enstatite. Passive sets all RRSS entries to the maximum. The Karato 2008 selector selects a fabric based on stress and water content as defined in figure 4 of the Karato 2008 review paper (doi: 10.1146/annurev.earth.36.031207.124120).

Parameter name: Model name#

Default value: Uniform grains and random uniform rotations

Pattern: [Anything]

Documentation: The model used to initialize the CPO for all particles. Currently ’Uniform grains and random uniform rotations’ is the only valid option.

Parameter name: Volume fractions minerals#

Default value: 0.7, 0.3

Pattern: [List of <[Double 0…MAX_DOUBLE (inclusive)]> of length 0…4294967295 (inclusive)]

Documentation: The volume fractions for the different minerals. There need to be the same number of values as there are minerals.Note that the currently implemented scheme is incompressible and does not allow chemical interaction or the formation of new phases

Subsection: Postprocess / Particles / Function#

Parameter name: Function constants#

Default value:

Pattern: [Anything]

Documentation: Sometimes it is convenient to use symbolic constants in the expression that describes the function, rather than having to use its numeric value everywhere the constant appears. These values can be defined using this parameter, in the form ‘var1=value1, var2=value2, …’.

A typical example would be to set this runtime parameter to ‘pi=3.1415926536’ and then use ‘pi’ in the expression of the actual formula. (That said, for convenience this class actually defines both ‘pi’ and ‘Pi’ by default, but you get the idea.)

Parameter name: Function expression#

Default value: 0

Pattern: [Anything]

Documentation: The formula that denotes the function you want to evaluate for particular values of the independent variables. This expression may contain any of the usual operations such as addition or multiplication, as well as all of the common functions such as ‘sin’ or ‘cos’. In addition, it may contain expressions like ‘if(x>0, 1, -1)’ where the expression evaluates to the second argument if the first argument is true, and to the third argument otherwise. For a full overview of possible expressions accepted see the documentation of the muparser library at http://muparser.beltoforion.de/.

If the function you are describing represents a vector-valued function with multiple components, then separate the expressions for individual components by a semicolon.

Parameter name: Number of components#

Default value: 1

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: The number of function components where each component is described by a function expression delimited by a ’;’.

Parameter name: Variable names#

Default value: x,y,t

Pattern: [Anything]

Documentation: The names of the variables as they will be used in the function, separated by commas. By default, the names of variables at which the function will be evaluated are ‘x’ (in 1d), ‘x,y’ (in 2d) or ‘x,y,z’ (in 3d) for spatial coordinates and ‘t’ for time. You can then use these variable names in your function expression and they will be replaced by the values of these variables at which the function is currently evaluated. However, you can also choose a different set of names for the independent variables at which to evaluate your function expression. For example, if you work in spherical coordinates, you may wish to set this input parameter to ‘r,phi,theta,t’ and then use these variable names in your function expression.

Subsection: Postprocess / Particles / Generator#

Subsection: Postprocess / Particles / Generator / Ascii file#

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: The number of radial shells of particles that will be generated around the central point.

Subsection: Postprocess / Particles / Integrator#

Subsection: Postprocess / Particles / Integrator / RK2#

Parameter name: Higher order accurate in time#

Default value: true

Pattern: [Bool]

Documentation: Whether to correctly evaluate old and current velocity solution to reach higher-order accuracy in time. If set to ’false’ only the old velocity solution is evaluated to simulate a first order method in time. This is only recommended for benchmark purposes.

Subsection: Postprocess / Particles / Interpolator#

Subsection: Postprocess / Particles / Interpolator / Bilinear least squares#

Parameter name: Use boundary extrapolation#

Default value: false

Pattern: [List of <[Bool]> of length 0…4294967295 (inclusive)]

Documentation: Extends the range used by ’Use linear least squares limiter’ by linearly interpolating values at cell boundaries from neighboring cells. If more than one value is given, it will be treated as a list with one component per particle property. Enabling ’Use boundary extrapolation’ requires enabling ’Use linear least squares limiter’.

Parameter name: Use linear least squares limiter#

Default value: false

Pattern: [List of <[Bool]> of length 0…4294967295 (inclusive)]

Documentation: Limit the interpolation of particle properties onto the cell, so that the value of each property is no smaller than its minimum and no larger than its maximum on the particles of each cell, and the average of neighboring cells. If more than one value is given, it will be treated as a list with one component per particle property.

Subsection: Postprocess / Particles / Interpolator / Quadratic least squares#

Parameter name: Use boundary extrapolation#

Default value: false

Pattern: [List of <[Bool]> of length 0…4294967295 (inclusive)]

Documentation: Extends the range used by ’Use quadratic least squares limiter’ by linearly interpolating values at cell boundaries from neighboring cells. If more than one value is given, it will be treated as a list with one component per particle property. Enabling ’Use boundary extrapolation’ requires enabling ’Use quadratic least squares limiter’.

Default value: true

Pattern: [Bool]

Documentation: deal.II offers the possibility to linearly interpolate output fields of higher order elements to a finer resolution. This somewhat compensates the fact that most visualization software only offers linear interpolation between grid points and therefore the output file is a very coarse representation of the actual solution field. Activating this option increases the spatial resolution in each dimension by a factor equal to the polynomial degree used for the velocity finite element (usually 2). In other words, instead of showing one quadrilateral or hexahedron in the visualization per cell on which ASPECT computes, it shows multiple (for quadratic elements, it will describe each cell of the mesh on which we compute as $2\times 2$ or $2\times 2\times 2$ cells in 2d and 3d, respectively; correspondingly more subdivisions are used if you use cubic, quartic, or even higher order elements for the velocity).

The effect of using this option can be seen in the following picture showing a variation of the output produced with the input files from Section~\ref{sec:shell-simple-2d}:

\begin{center} \includegraphics[width=0.5\textwidth]{viz/parameters/build-patches}\end{center}Here, the left picture shows one visualization cell per computational cell (i.e., the option is switched off), and the right picture shows the same simulation with the option switched on (which is the default). The images show the same data, demonstrating that interpolating the solution onto bilinear shape functions as is commonly done in visualizing data loses information.

Of course, activating this option also greatly increases the amount of data ASPECT will write to disk: approximately by a factor of 4 in 2d, and a factor of 8 in 3d, when using quadratic elements for the velocity, and correspondingly more for even higher order elements.

Parameter name: List of output variables#

Default value:

Documentation: A comma separated list of visualization objects that should be run whenever writing graphical output. By default, the graphical output files will always contain the primary variables velocity, pressure, and temperature. However, one frequently wants to also visualize derived quantities, such as the thermodynamic phase that corresponds to a given temperature-pressure value, or the corresponding seismic wave speeds. The visualization objects do exactly this: they compute such derived quantities and place them into the output file. The current parameter is the place where you decide which of these additional output variables you want to have in your output file.

The following postprocessors are available:

‘ISA rotation timescale’: A visualization output object that generates output showing the timescale for the rotation of grains toward the infinite strain axis. Kaminski and Ribe (see \cite{Kaminski2002}) call this quantity $\tau_\text{ISA}$ and define it as $\tau_\text{ISA} \approx \frac{1}{\dot{\epsilon}}$ where $\dot{\epsilon}$ is the largest eigenvalue of the strain rate tensor. It can be used, along with the grain lag angle $\Theta$, to calculate the grain orientation lag parameter.

Physical units: \si{\second}.

‘Vp anomaly’: A visualization output object that generates output showing the percentage anomaly in the seismic compressional wave speed $V_p$ as a spatially variable function with one value per cell. This anomaly is either shown as a percentage anomaly relative to the reference profile given by adiabatic conditions (with the compositions given by the current composition, such that the reference could potentially change through time), or as a percentage change relative to the laterally averaged velocity at the depth of the cell. This velocity is calculated by linear interpolation between average values calculated within equally thick depth slices. The number of depth slices in the domain is user-defined. Typically, the best results will be obtained if the number of depth slices is balanced between being large enough to capture step changes in velocities, but small enough to maintain a reasonable number of evaluation points per slice. Bear in mind that lateral averaging subsamples the finite element mesh. Note that this plugin requires a material model that provides seismic velocities.

Physical units: None, the quantity being output is a fractional change provided as a percentage.

‘Vs anomaly’: A visualization output object that generates output showing the percentage anomaly in the seismic shear wave speed $V_s$ as a spatially variable function with one value per cell. This anomaly is either shown as a percentage anomaly relative to the reference profile given by adiabatic conditions (with the compositions given by the current composition, such that the reference could potentially change through time), or as a percentage change relative to the laterally averaged velocity at the depth of the cell. This velocity is calculated by linear interpolation between average values calculated within equally thick depth slices. The number of depth slices in the domain is user-defined. Typically, the best results will be obtained if the number of depth slices is balanced between being large enough to capture step changes in velocities, but small enough to maintain a reasonable number of evaluation points per slice. Bear in mind that lateral averaging subsamples the finite element mesh. Note that this plugin requires a material model that provides seismic velocities.

Physical units: None, the quantity being output is a fractional change provided as a percentage.

‘adiabat’: A visualization output object that generates adiabatic temperature, pressure, density, and density derivative (with regard to depth)as produced by the AdiabaticConditions class.

Physical units: \si{\kelvin}, \si{\pascal}, \si{\kilo\gram\per\meter\cubed\per\meter}, respectively, for the four components.

‘artificial viscosity’: A visualization output object that generates output showing the value of the artificial viscosity on each cell.

Physical units: \si{\watt\per\meter\per\kelvin}.

‘artificial viscosity composition’: A visualization output object that generates output showing the value of the artificial viscosity for a compositional field on each cell.

Physical units: \si{\meter\squared\per\second}.

‘boundary indicators’: A visualization output object that generates output about the used boundary indicators. In a loop over the active cells, if a cell lies at a domain boundary, the boundary indicator of the face along the boundary is requested. In case the cell does not lie along any domain boundary, the cell is assigned the value of the largest used boundary indicator plus one. When a cell is situated in one of the corners of the domain, multiple faces will have a boundary indicator. This postprocessor returns the value of the first face along a boundary that is encountered in a loop over all the faces.

Physical units: None.

‘boundary strain rate residual’: A visualization output object that generates output for the strain rate residual at the top surface. The residual is computed at each point at the surface as the difference between the strain rate invariant in the model and the input data, where the invariant is computed like in the ’strain rate’ postprocessor. The user chooses the input data as ascii data files with coordinate columns and column corresponding to the surface strain rate norm.

Physical units: $\frac{1}{\text{s}}$ or $\frac{1}{\text{year}}$, depending on settings in the input file.

‘boundary velocity residual’: A visualization output object that generates output for the velocity residual at the top surface. The residual is computed at each point at the surface as the difference between the modeled velocities and the input data velocities for each vector component. The user has an option to choose the input data as ascii data files (e.g. GPS velocities) with columns in the same format as described for the ’ascii data’ initial temperature plugin or a velocity field computed from the GPlates program as described in the gplates boundary velocity plugin.

Physical units: $\frac{\text{m}}{\text{s}}$ or $\frac{\text{m}}{\text{year}}$, depending on settings in the input file.

‘compositional vector’: A visualization output object that outputs vectors whose components are derived from compositional fields. Input parameters for this postprocessor are defined in section Postprocess/Visualization/Compositional fields as vectors.

Physical units: None.

‘depth’: A visualization output postprocessor that outputs the depth for all points inside the domain, as determined by the geometry model.

Physical units: \si{\meter}.

‘dynamic topography’: A visualization output object that generates output for the dynamic topography at the top and bottom of the model space. The approach to determine the dynamic topography requires us to compute the stress tensor and evaluate the component of it in the direction in which gravity acts. In other words, we compute $\sigma_{rr}={\hat g}^T(2 \eta \varepsilon(\mathbf u)-\frac 13 (\textrm{div}\;\mathbf u)I)\hat g - p_d$ where $\hat g = \mathbf g/\|\mathbf g\|$ is the direction of the gravity vector $\mathbf g$ and $p_d=p-p_a$ is the dynamic pressure computed by subtracting the adiabatic pressure $p_a$ from the total pressure $p$ computed as part of the Stokes solve. From this, the dynamic topography is computed using the formula $h=\frac{\sigma_{rr}}{(\mathbf g \cdot \mathbf n) \rho}$ where $\rho$ is the density at the cell center. For the bottom surface we chose the convection that positive values are up (out) and negative values are in (down), analogous to the deformation of the upper surface. Note that this implementation takes the direction of gravity into account, which means that reversing the flow in backward advection calculations will not reverse the instantaneous topography because the reverse flow will be divided by the reverse surface gravity.

Strictly speaking, the dynamic topography is of course a quantity that is only of interest at the surface. However, we compute it everywhere to make things fit into the framework within which we produce data for visualization. You probably only want to visualize whatever data this postprocessor generates at the surface of your domain and simply ignore the rest of the data generated.

Alternatively, consider using the “surface dynamic topography” visualization postprocessor to only output the dynamic topography at the boundary of the domain.

Physical units: \si{\meter}.

‘error indicator’: A visualization output object that generates output showing the estimated error or other mesh refinement indicator as a spatially variable function with one value per cell.

Physical units: None. (Strictly speaking, errors have physical units of course, but because error indicators can be computed from different solution components and other input, we consider error indicators unitless.)

‘geoid’: Visualization for the geoid solution. The geoid is given by the equivalent water column height due to a gravity perturbation.

Physical units: \si{\meter}.

‘grain lag angle’: A visualization output object that generates output showing the angle between the ~infinite strain axis and the flow velocity. Kaminski and Ribe (see \cite{Kaminski2002}) call this quantity $\Theta$ and define it as $\Theta = \cos^{-1}(\hat{u}\cdot\hat{e})$ where $\hat{u}=\vec{u}/|{u}|$, $\vec{u}$ is the local flow velocity, and $\hat{e}$ is the local infinite strain axis, which we calculate as the first eigenvector of the ’left stretch’ tensor. $\Theta$ can be used to calculate the grain orientation lag parameter.

Physical units: \si{\radian}.

‘gravity’: A visualization output object that outputs the gravity vector.

Physical units: \si {\meter\per\second\squared} .

‘heat flux map’: A visualization output object that generates output for the heat flux density across the top and bottom boundary in outward direction. The heat flux is computed as sum of advective heat flux and conductive heat flux through Neumann boundaries, both computed as integral over the boundary area, and conductive heat flux through Dirichlet boundaries, which is computed using the consistent boundary flux method as described in “Gresho, Lee, Sani, Maslanik, Eaton (1987). The consistent Galerkin FEM for computing derived boundary quantities in thermal and or fluids problems. International Journal for Numerical Methods in Fluids, 7(4), 371-394.” If only conductive heat flux through Dirichlet boundaries is of interest, the postprocessor can produce output of higher resolution by evaluating the CBF solution vector point-wise instead of computing cell-wise averaged values.

Physical units: \si{\watt\per\meter\squared}.

‘heating’: A visualization output object that generates output for all the heating terms used in the energy equation.

Physical units: \si{\watt\per\cubic\meter}.

‘material properties’: A visualization output object that generates output for the material properties given by the material model. The current postprocessor allows to output a (potentially large) subset of all of the information provided by material models at once, with just a single material model evaluation per output point. Although individual properties can still be listed in the “List of output variables”, this visualization plugin is called internally to avoid duplicated evaluations of the material model.

In almost all places inside ASPECT, the program can use “averaged” material properties, for example for the assembly of matrices and right hand side vectors. To accurately reflect the material parameters used internally, this visualization postprocessor averages in the same way as is used to do the assembly, and consequently the graphical output will reflect not pointwise properties, but averaged properties.

Physical units: Various.

‘maximum horizontal compressive stress’: A plugin that computes the direction and magnitude of the maximum horizontal component of the compressive stress as a vector field. The direction of this vector can often be used to visualize the principal mode of deformation (e.g., at normal faults or extensional margins) and can be correlated with seismic anisotropy. Recall that the compressive stress is simply the negative stress, $\sigma_c=-\sigma=-\left[ 2\eta (\varepsilon(\mathbf u) - \frac 13 (\nabla \cdot \mathbf u) I) + pI\right]$.

Following \cite{LundTownend07}, we define the maximum horizontal stress direction as that horizontal direction $\mathbf n$ that maximizes $\mathbf n^T \sigma_c \mathbf n$. We call a vector horizontal if it is perpendicular to the gravity vector $\mathbf g$.

In two space dimensions, $\mathbf n$ is simply a vector that is horizontal (we choose one of the two possible choices). This direction is then scaled by the size of the horizontal stress in this direction, i.e., the plugin outputs the vector $\mathbf w = (\mathbf n^T \sigma_c \mathbf n) \; \mathbf n$.

In three space dimensions, given two horizontal, perpendicular, unit length, but otherwise arbitrarily chosen vectors $\mathbf u,\mathbf v$, we can express $\mathbf n = (\cos \alpha)\mathbf u + (\sin\alpha)\mathbf v$ where $\alpha$ maximizes the expression \begin{align*} f(\alpha) = \mathbf n^T \sigma_c \mathbf n = (\mathbf u^T \sigma_c \mathbf u)(\cos\alpha)^2 +2(\mathbf u^T \sigma_c \mathbf v)(\cos\alpha)(\sin\alpha) +(\mathbf v^T \sigma_c \mathbf v)(\sin\alpha)^2.\end{align*}

The maximum of $f(\alpha)$ is attained where $f’(\alpha)=0$. Evaluating the derivative and using trigonometric identities, one finds that $\alpha$ has to satisfy the equation \begin{align*} \tan(2\alpha) = \frac{2.0\mathbf u^T \sigma_c \mathbf v} {\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v}.\end{align*}Since the transform $\alpha\mapsto\alpha+\pi$ flips the direction of $\mathbf n$, we only need to seek a solution to this equation in the interval $\alpha\in[0,\pi)$. These are given by $\alpha_1=\frac 12 \arctan \frac{\mathbf u^T \sigma_c \mathbf v}{\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v}$ and $\alpha_2=\alpha_1+\frac{\pi}{2}$, one of which will correspond to a minimum and the other to a maximum of $f(\alpha)$. One checks the sign of $f”(\alpha)=-2(\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v)\cos(2\alpha) - 2 (\mathbf u^T \sigma_c \mathbf v) \sin(2\alpha)$ for each of these to determine the $\alpha$ that maximizes $f(\alpha)$, and from this immediately arrives at the correct form for the maximum horizontal stress $\mathbf n$.

The description above computes a 3d direction vector $\mathbf n$. If one were to scale this vector the same way as done in 2d, i.e., with the magnitude of the stress in this direction, one will typically get vectors whose length is principally determined by the hydrostatic pressure at a given location simply because the hydrostatic pressure is the largest component of the overall stress. On the other hand, the hydrostatic pressure does not determine any principal direction because it is an isotropic, anti-compressive force. As a consequence, there are often points in simulations (e.g., at the center of convection rolls) where the stress has no dominant horizontal direction, and the algorithm above will then in essence choose a random direction because the stress is approximately equal in all horizontal directions. If one scaled the output by the magnitude of the stress in this direction (i.e., approximately equal to the hydrostatic pressure at this point), one would get randomly oriented vectors at these locations with significant lengths.

To avoid this problem, we scale the maximal horizontal compressive stress direction $\mathbf n$ by the difference between the stress in the maximal and minimal horizontal stress directions. In other words, let $\mathbf n_\perp=(\sin \alpha)\mathbf u - (\cos\alpha)\mathbf v$ be the horizontal direction perpendicular to $\mathbf n$, then this plugin outputs the vector quantity $\mathbf w = (\mathbf n^T \sigma_c \mathbf n -\mathbf n^T_\perp \sigma_c \mathbf n_\perp) \; \mathbf n$. In other words, the length of the vector produced indicates how dominant the direction of maximal horizontal compressive strength is.

Fig.~\ref{fig:max-horizontal-compressive-stress} shows a simple example for this kind of visualization in 3d.

\begin{figure} \includegraphics[width=0.3\textwidth] {viz/plugins/maximum_horizontal_compressive_stress/temperature.png} \hfill \includegraphics[width=0.3\textwidth] {viz/plugins/maximum_horizontal_compressive_stress/velocity.png} \hfill \includegraphics[width=0.3\textwidth] {viz/plugins/maximum_horizontal_compressive_stress/horizontal-stress.png} \caption{\it Illustration of the ‘maximum horizontal compressive stress’ visualization plugin. The left figure shows a ridge-like temperature anomaly. Together with no-slip boundary along all six boundaries, this results in two convection rolls (center). The maximal horizontal compressive strength at the bottom center of the domain is perpendicular to the ridge because the flow comes together there from the left and right, yielding a compressive force in left-right direction. At the top of the model, the flow separates outward, leading to a negative compressive stress in left-right direction; because there is no flow in front-back direction, the compressive strength in front-back direction is zero, making the along-ridge direction the dominant one. At the center of the convection rolls, both horizontal directions yield the same stress; the plugin therefore chooses an essentially arbitrary horizontal vector, but then uses a zero magnitude given that the difference between the maximal and minimal horizontal stress is zero at these points.} \label{fig:max-horizontal-compressive-stress}\end{figure}

Physical units: \si{\pascal}.

‘melt fraction’: A visualization output object that generates output for the melt fraction at the temperature and pressure of the current point. If the material model computes a melt fraction, this is the quantity that will be visualized. Otherwise, a specific parametrization for batch melting (as described in the following) will be used. It does not take into account latent heat. If there are no compositional fields, or no fields called ’pyroxenite’, this postprocessor will visualize the melt fraction of peridotite (calculated using the anhydrous model of Katz, 2003). If there is a compositional field called ’pyroxenite’, the postprocessor assumes that this compositional field is the content of pyroxenite, and will visualize the melt fraction for a mixture of peridotite and pyroxenite (using the melting model of Sobolev, 2011 for pyroxenite). All the parameters that were used in these calculations can be changed in the input file, the most relevant maybe being the mass fraction of Cpx in peridotite in the Katz melting model (Mass fraction cpx), which right now has a default of 15%. The corresponding $p$-$T$-diagrams can be generated by running the tests melt_postprocessor_peridotite and melt_postprocessor_pyroxenite.

Physical units: None.

‘melt material properties’: A visualization output object that generates output for melt related properties of the material model. Note that this postprocessor always outputs the compaction pressure, but can output a large range of additional properties, as selected in the “List of properties” parameter.

Physical units: Various, depending on what is being output.

‘named additional outputs’: Some material models can compute quantities other than those that typically appear in the equations that ASPECT solves (such as the viscosity, density, etc). Examples of quantities material models may be able to compute are seismic velocities, or other quantities that can be derived from the state variables and the material coefficients such as the stress or stress anisotropies. These quantities are generically referred to as ‘named outputs’ because they are given an explicit name different from the usual outputs of material models.

This visualization postprocessor outputs whatever quantities the material model can compute. What quantities these are is specific to the material model in use for a simulation, and for many models in fact does not contain any named outputs at all.

Physical units: Various, depending on what is being output.

‘nonadiabatic pressure’: A visualization output object that generates output for the non-adiabatic component of the pressure.

The variable that is outputted this way is computed by taking the pressure at each point and subtracting from it the adiabatic pressure computed at the beginning of the simulation. Because the adiabatic pressure is one way of defining a static pressure background field, what this visualization postprocessor therefore produces is one way to compute a dynamic pressure. There are, however, other ways as well, depending on the choice of the “background pressure”.

Physical units: \si{\pascal}.

‘nonadiabatic temperature’: A visualization output object that generates output for the non-adiabatic component of the temperature.

Physical units: \si{\kelvin}.

‘particle count’: A visualization output object that generates output about the number of particles per cell.

Physical units: None.

‘partition’: A visualization output object that generates output for the parallel partition that every cell of the mesh is associated with.

Physical units: None.

‘principal stress’: A visualization output object that outputs the principal stress values and directions, i.e., the eigenvalues and eigenvectors of the stress tensor. The postprocessor can either operate on the full stress tensor or only on the deviatoric stress tensor, depending on what run-time parameters are set.

Physical units: \si{\pascal}.

‘shear stress’: A visualization output object that generates output for the 3 (in 2d) or 6 (in 3d) components of the shear stress tensor, i.e., for the components of the tensor $-2\eta\varepsilon(\mathbf u)$ in the incompressible case and $-2\eta\left[\varepsilon(\mathbf u)-\tfrac 13(\textrm{tr}\;\varepsilon(\mathbf u))\mathbf I\right]$ in the compressible case. If elasticity is used, the elastic contribution is being accounted for. The shear stress differs from the full stress tensor by the absence of the pressure. Note that the convention of positive compressive stress is followed.

Physical units: \si{\pascal}.

‘spd factor’: A visualization output object that generates output for the spd factor. The spd factor is a factor which scales a part of the Jacobian used for the Newton solver to make sure that the Jacobian remains positive definite.

Physical units: None.

‘spherical velocity components’: A visualization output object that outputs the polar coordinates components $v_r$ and $v_\phi$ of the velocity field in 2d and the spherical coordinates components $v_r$, $v_{\phi}$ and $v_{\theta}$ of the velocity field in 3d.

Physical units: $\frac{\text{m}}{\text{s}}$ or $\frac{\text{m}}{\text{year}}$, depending on settings in the input file.

‘strain rate’: A visualization output object that generates output for the norm of the strain rate, i.e., for the quantity $\sqrt{\varepsilon(\mathbf u):\varepsilon(\mathbf u)}$ in the incompressible case and $\sqrt{[\varepsilon(\mathbf u)-\tfrac 13(\textrm{tr}\;\varepsilon(\mathbf u))\mathbf I]:[\varepsilon(\mathbf u)-\tfrac 13(\textrm{tr}\;\varepsilon(\mathbf u))\mathbf I]}$ in the compressible case.

Physical units: \si{\per\second}.

‘strain rate tensor’: A visualization output object that generates output for the 4 (in 2d) or 9 (in 3d) components of the strain rate tensor, i.e., for the components of the tensor $\varepsilon(\mathbf u)$ in the incompressible case and $\varepsilon(\mathbf u)-\tfrac 13(\textrm{tr}\;\varepsilon(\mathbf u))\mathbf I$ in the compressible case.

Physical units: \si{\per\second}.

‘stress’: A visualization output object that generates output for the 3 (in 2d) or 6 (in 3d) components of the stress tensor, i.e., for the components of the tensor $-2\eta\varepsilon(\mathbf u)+pI$ in the incompressible case and $-2\eta\left[\varepsilon(\mathbf u)-\tfrac 13(\textrm{tr}\;\varepsilon(\mathbf u))\mathbf I\right]+pI$ in the compressible case. If elasticity is used, the elastic contribution is being accounted for. Note that the convention of positive compressive stress is followed.

Physical units: \si{\pascal}.

‘stress second invariant’: A visualization output object that outputs the second moment invariant of the deviatoric stress tensor.

Physical units: \si{\pascal}.

‘surface dynamic topography’: A visualization output object that generates output for the dynamic topography at the top and bottom of the model space. The approach to determine the dynamic topography requires us to compute the stress tensor and evaluate the component of it in the direction in which gravity acts. In other words, we compute $\sigma_{rr}={\hat g}^T(2 \eta \varepsilon(\mathbf u)-\frac 13 (\textrm{div}\;\mathbf u)I)\hat g - p_d$ where $\hat g = \mathbf g/\|\mathbf g\|$ is the direction of the gravity vector $\mathbf g$ and $p_d=p-p_a$ is the dynamic pressure computed by subtracting the adiabatic pressure $p_a$ from the total pressure $p$ computed as part of the Stokes solve. From this, the dynamic topography is computed using the formula $h=\frac{\sigma_{rr}}{(\mathbf g \cdot \mathbf n) \rho}$ where $\rho$ is the density at the cell center. For the bottom surface we chose the convection that positive values are up (out) and negative values are in (down), analogous to the deformation of the upper surface. Note that this implementation takes the direction of gravity into account, which means that reversing the flow in backward advection calculations will not reverse the instantaneous topography because the reverse flow will be divided by the reverse surface gravity.

In contrast to the ‘dynamic topography’ visualization postprocessor, this plugin really only evaluates the dynamic topography at faces of cells that are adjacent to ‘bottom’ and ‘top’ boundaries, and only outputs information on the surface of the domain, rather than padding the information with zeros in the interior of the domain.

Physical units: \si{\meter}.

‘surface stress’: A visualization output object that generates output on the surface of the domain for the 3 (in 2d) or 6 (in 3d) components of the stress tensor, i.e., for the components of the tensor $-2\eta\varepsilon(\mathbf u)+pI$ in the incompressible case and $-2\eta\left[\varepsilon(\mathbf u)-\tfrac 13(\textrm{tr}\;\varepsilon(\mathbf u))\mathbf I\right]+pI$ in the compressible case. If elasticity is included, its contribution is accounted for. Note that the convention of positive compressive stress is followed.

Physical units: \si{\pascal}.

‘temperature anomaly’: A visualization output postprocessor that outputs the temperature minus the depth-average of the temperature.The average temperature is calculated using the lateral averaging function from the “depth average” postprocessor and interpolated linearly between the layers specified through “Number of depth slices”.

Physical units: \si{\kelvin}.

‘vertical heat flux’: A visualization output object that generates output for the heat flux in the vertical direction, which is the sum of the advective and the conductive heat flux, with the sign convention of positive flux upwards.

Physical units: \si{\watt\per\square\meter}.

‘volume of fluid values’: A visualization output object that outputs the volume fraction and optionally a level set field and the interface normal vectors of volume of fluid fields.

Physical units: None.

‘volumetric strain rate’: A visualization output object that generates output for the volumetric strain rate, i.e., for the quantity $\nabla\cdot\mathbf u = \textrm{div}\; \mathbf u = \textrm{trace}\; \varepsilon(\mathbf u)$. This should be zero (in some average sense) in incompressible convection models, but can be non-zero in compressible models and models with melt transport.

Physical units: \si{\per\second}.

Default value:

Pattern: [Anything]

Parameter name: Time between graphical output#

Default value: 1e8

Pattern: [Double 0…MAX_DOUBLE (inclusive)]

Parameter name: Time steps between graphical output#

Default value: 2147483647

Pattern: [Integer range 0…2147483647 (inclusive)]

Documentation: The maximum number of time steps between each generation of graphical output files.

Parameter name: Write higher order output#

Default value: false

Pattern: [Bool]

Documentation: deal.II offers the possibility to write vtu files with higher order representations of the output data. This means each cell will correctly show the higher order representation of the output data instead of the linear interpolation between vertices that ParaView and VisIt usually show. Note that activating this option is safe and recommended, but requires that (i) “Output format” is set to “vtu”, (ii) “Interpolate output” is set to true, (iii) you use a sufficiently new version of Paraview or VisIt to read the files (Paraview version 5.5 or newer, and VisIt version to be determined), and (iv) you use deal.II version 9.1.0 or newer. The effect of using this option can be seen in the following picture:

\begin{center} \includegraphics[width=0.5\textwidth]{viz/parameters/higher-order-output}\end{center}The top figure shows the plain output without interpolation or higher order output. The middle figure shows output that was interpolated as discussed for the “Interpolate output” option. The bottom panel shows higher order output that achieves better accuracy than the interpolated output at a lower memory cost.

Parameter name: Write in background thread#

Documentation: Names of vectors as they will appear in the output.

Subsection: Postprocess / Visualization / Heat flux map#

Parameter name: Output point wise heat flux#

Default value: false

Pattern: [Bool]

Documentation: A boolean flag that controls whether to output the heat flux map as a point wise value, or as a cell-wise averaged value. The point wise output is more accurate, but it currently omits prescribed heat flux values at boundaries and advective heat flux that is caused by velocities non-tangential to boundaries. If you do not use these two features it is recommended to switch this setting on to benefit from the increased output resolution.

Subsection: Postprocess / Visualization / Material properties#

Parameter name: List of material properties#

Default value: density,thermal expansivity,specific heat,viscosity

Documentation: A comma separated list of material properties that should be written whenever writing graphical output. By default, the material properties will always contain the density, thermal expansivity, specific heat and viscosity. The following material properties are available:

Subsection: Postprocess / Visualization / Melt fraction#

Parameter name: A1#

Default value: 1085.7

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Constant parameter in the quadratic function that approximates the solidus of peridotite. Units: \si{\degreeCelsius}.

Parameter name: A2#

Default value: 1.329e-7

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the linear pressure term in the quadratic function that approximates the solidus of peridotite. \si{\degreeCelsius\per\pascal}.

Parameter name: A3#

Default value: -5.1e-18

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the quadratic pressure term in the quadratic function that approximates the solidus of peridotite. \si{\degreeCelsius\per\pascal\squared}.

Parameter name: B1#

Default value: 1475.0

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Constant parameter in the quadratic function that approximates the lherzolite liquidus used for calculating the fraction of peridotite-derived melt. Units: \si{\degreeCelsius}.

Parameter name: B2#

Default value: 8.0e-8

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the linear pressure term in the quadratic function that approximates the lherzolite liquidus used for calculating the fraction of peridotite-derived melt. \si{\degreeCelsius\per\pascal}.

Parameter name: B3#

Default value: -3.2e-18

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the quadratic pressure term in the quadratic function that approximates the lherzolite liquidus used for calculating the fraction of peridotite-derived melt. \si{\degreeCelsius\per\pascal\squared}.

Parameter name: C1#

Default value: 1780.0

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Constant parameter in the quadratic function that approximates the liquidus of peridotite. Units: \si{\degreeCelsius}.

Parameter name: C2#

Default value: 4.50e-8

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the linear pressure term in the quadratic function that approximates the liquidus of peridotite. \si{\degreeCelsius\per\pascal}.

Parameter name: C3#

Default value: -2.0e-18

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the quadratic pressure term in the quadratic function that approximates the liquidus of peridotite. \si{\degreeCelsius\per\pascal\squared}.

Parameter name: D1#

Default value: 976.0

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Constant parameter in the quadratic function that approximates the solidus of pyroxenite. Units: \si{\degreeCelsius}.

Parameter name: D2#

Default value: 1.329e-7

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the linear pressure term in the quadratic function that approximates the solidus of pyroxenite. Note that this factor is different from the value given in Sobolev, 2011, because they use the potential temperature whereas we use the absolute temperature. \si{\degreeCelsius\per\pascal}.

Parameter name: D3#

Default value: -5.1e-18

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the quadratic pressure term in the quadratic function that approximates the solidus of pyroxenite. \si{\degreeCelsius\per\pascal\squared}.

Parameter name: E1#

Default value: 663.8

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the linear depletion term in the quadratic function that approximates the melt fraction of pyroxenite. \si{\degreeCelsius\per\pascal}.

Parameter name: E2#

Default value: -611.4

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the quadratic depletion term in the quadratic function that approximates the melt fraction of pyroxenite. \si{\degreeCelsius\per\pascal\squared}.

Parameter name: Mass fraction cpx#

Default value: 0.15

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Mass fraction of clinopyroxene in the peridotite to be molten. Units: non-dimensional.

Parameter name: beta#

Default value: 1.5

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Exponent of the melting temperature in the melt fraction calculation. Units: non-dimensional.

Parameter name: r1#

Default value: 0.5

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Constant in the linear function that approximates the clinopyroxene reaction coefficient. Units: non-dimensional.

Parameter name: r2#

Default value: 8e-11

Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]

Documentation: Prefactor of the linear pressure term in the linear function that approximates the clinopyroxene reaction coefficient. Units: \si{\per\pascal}.

Subsection: Postprocess / Visualization / Melt material properties#

Parameter name: List of properties#

Default value: compaction viscosity,permeability

Documentation: A comma separated list of melt properties that should be written whenever writing graphical output. The following material properties are available:

Subsection: Postprocess / Visualization / Principal stress#

Parameter name: Use deviatoric stress#

Default value: false

Pattern: [Bool]

Documentation: Whether to use the deviatoric stress tensor instead of the full stress tensor to compute principal stress directions and values.

Subsection: Postprocess / Visualization / Temperature anomaly#

Parameter name: Number of depth slices#

Default value: 20

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: Number of depth slices used to define average temperature.

Parameter name: Use maximal temperature for bottom#

Default value: true

Pattern: [Bool]

Documentation: If true, use the specified boundary temperatures as average temperatures at the surface. If false, extrapolate the temperature gradient between the first and second cells to the surface. This option will only work for models with a fixed surface temperature.

Parameter name: Use minimal temperature for surface#

Default value: true

Pattern: [Bool]

Documentation: Whether to use the minimal specified boundary temperature as the bottom boundary temperature. This option will only work for models with a fixed bottom boundary temperature.

Subsection: Postprocess / Visualization / Volume of Fluid#

Parameter name: Output interface normals#

Default value: false

Pattern: [Bool]

Documentation: Include the internal data for the interface normal on the unit cells.

Parameter name: Output interface reconstruction contour#

Default value: false

Pattern: [Bool]

Documentation: Include fields defined such that the 0 contour is the fluid interface.

Subsection: Postprocess / Visualization / Vp anomaly#

Parameter name: Average velocity scheme#

Default value: reference profile

Pattern: [Selection reference profile|lateral average ]

Documentation: Scheme to compute the average velocity-depth profile. The reference profile option evaluates the conditions along the reference adiabat according to the material model. The lateral average option instead calculates a lateral average from subdivision of the mesh. The lateral average option may produce spurious results where there are sharp velocity changes.

Parameter name: Number of depth slices#

Default value: 50

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: Number of depth slices used to define average seismic compressional wave velocities from which anomalies are calculated. Units: non-dimensional.

Subsection: Postprocess / Visualization / Vs anomaly#

Parameter name: Average velocity scheme#

Default value: reference profile

Pattern: [Selection reference profile|lateral average ]

Parameter name: Number of depth slices#

Default value: 50

Pattern: [Integer range 1…2147483647 (inclusive)]

Documentation: Number of depth slices used to define average seismic shear wave velocities from which anomalies are calculated. Units: non-dimensional.

ASPECT 2.5.0

Postprocess

Contents

Postprocess#

Subsection: Postprocess#

Parameter name: List of postprocessors#

Parameter name: Run postprocessors on nonlinear iterations#

Subsection: Postprocess / Boundary strain rate residual statistics#

Parameter name: Data directory#

Parameter name: Data file name#

Parameter name: Scale factor#

Subsection: Postprocess / Boundary velocity residual statistics#

Parameter name: Data directory#

Parameter name: Data file name#

Parameter name: Scale factor#

Parameter name: Use ascii data#

Parameter name: Use spherical unit vectors#

Subsection: Postprocess / Command#

Parameter name: Command#

Parameter name: Run on all processes#

Parameter name: Terminate on failure#

Subsection: Postprocess / Composition velocity statistics#

Parameter name: Names of selected compositional fields#

Subsection: Postprocess / Depth average#

Parameter name: Depth boundaries of zones#

Parameter name: List of output variables#

Parameter name: Number of zones#

Parameter name: Output format#

Parameter name: Time between graphical output#

Subsection: Postprocess / Dynamic core statistics#

Parameter name: Excess entropy only#

Subsection: Postprocess / Dynamic topography#

Parameter name: Density above#

Parameter name: Density below#

Parameter name: Output bottom#

Parameter name: Output surface#

Subsection: Postprocess / Geoid#

Parameter name: Also output the gravity anomaly#

Parameter name: Also output the spherical harmonic coefficients of CMB dynamic topography contribution#

Parameter name: Also output the spherical harmonic coefficients of density anomaly contribution#

Parameter name: Also output the spherical harmonic coefficients of geoid anomaly#

Parameter name: Also output the spherical harmonic coefficients of surface dynamic topography contribution#

Parameter name: Density above#

Parameter name: Density below#

Parameter name: Include CMB topography contribution#

Parameter name: Include surface topography contribution#

Parameter name: Include the contributon from dynamic topography#

Parameter name: Maximum degree#

Parameter name: Minimum degree#

Parameter name: Output CMB topography contribution coefficients#

Parameter name: Output data in geographical coordinates#

Parameter name: Output density anomaly contribution coefficients#

Parameter name: Output geoid anomaly coefficients#

Parameter name: Output gravity anomaly#

Parameter name: Output surface topography contribution coefficients#

Subsection: Postprocess / Global statistics#

Parameter name: Write statistics for each nonlinear iteration#

Subsection: Postprocess / Gravity calculation#

Parameter name: List of latitude#

Parameter name: List of longitude#

Parameter name: List of radius#

Parameter name: Maximum latitude#

Parameter name: Maximum longitude#

Parameter name: Maximum radius#

Parameter name: Minimum latitude#

Parameter name: Minimum longitude#

Parameter name: Minimum radius#

Parameter name: Number points fibonacci spiral#

Parameter name: Number points latitude#

Parameter name: Number points longitude#

Parameter name: Number points radius#

Parameter name: Precision in gravity output#

Parameter name: Quadrature degree increase#

Parameter name: Reference density#

Parameter name: Sampling scheme#

Parameter name: Time between gravity output#

Parameter name: Time steps between gravity output#

Subsection: Postprocess / Memory statistics#

Parameter name: Output peak virtual memory (VmPeak)#

Subsection: Postprocess / Particles#