I have a question about interpreting the maximum_horizontal_compressive_stress output in ASPECT. I’m trying to compute the direction of the maximum horizontal compressive stress (SHmax) from two different outputs, but I’m getting inconsistent results.
Method 1: Using maximum_horizontal_compressive_stress:0 and :1
I assumed that:
maximum_horizontal_compressive_stress:0 → East component
maximum_horizontal_compressive_stress:1 → North component
Then I computed the direction as: \theta =arctan2(max\_comp\_stress:0, max\_comp\_stress:1)=35.09°
But if I mistakenly reverse the order: \theta=arctan2(1.7242e+08,1.2112e+08)=54.9°
The following is a brief explanation of section 2.1 of the article
1. Input
A 3×3 symmetric stress tensor \mathbf{S} in geographic coordinates (North, East, Down).
2. Principal Stress Decomposition
The stress tensor is decomposed via eigenvalue analysis:
\mathbf{S} \mathbf{v}_i = \lambda_i \mathbf{v}_i
Eigenvalues (\lambda_1 \geq \lambda_2 \geq \lambda_3) represent the principal stresses.
Corresponding eigenvectors \mathbf{v}_i give the principal stress directions.
3. Objective: Maximize Normal Stress on Vertical Planes
The algorithm assumes that SHmax corresponds to the orientation of a vertical fault plane that experiences the maximum normal stress.
A trial normal vector on the horizontal plane is defined as:
\mathbf{n}(\alpha) = (\cos\alpha,\ \sin\alpha,\ 0)
where \alpha is the azimuth angle (clockwise from North).
4. Normal Stress Expression
The normal stress S_n on a plane with normal \mathbf{n} is computed using:
S_n = \sum_{i=1}^3 \lambda_i (\mathbf{v}_i \cdot \mathbf{n})^2
5. Optimization via Calculus
The first derivative \frac{dS_n}{d\alpha} is symbolically computed and set to zero to find stationary points.
The second derivative \frac{d^2S_n}{d\alpha^2} is evaluated to classify each point:
If negative → local maximum → candidate for SHmax.
If positive → local minimum → corresponds to Shmin; SHmax is orthogonal (\alpha + 90^\circ).
6. Output
All valid solutions are converted to degrees, normalized to [0^\circ, 180^\circ), deduplicated, and returned as a list of possible SHmax orientations.
And my result of the angle is 65.245°
question
The two methods give very different results
Is maximum_horizontal_compressive_stress:0,1 meant to represent the vector components of SHmax?If so, why does its direction not match the analytical principal stress direction?
Does ASPECT directly output the SHmax direction, or do I need to compute it from the stress tensor?
Any help would be greatly appreciated.
Best,
Haolin
The maximum horizontal compressive stress post processor follows the methods of Lund and Townsend 2007, and the exact output is described in the manual postprocessor visualization section (search for maximum horizontal compressive stress within it). Can you review that documentation and let us know if you have any further questions?
I looked into the code and have a few observations and recommendations.
pressure sign
The code builds compressive_stress as compressive_stress = in.pressure[q] * I; compressive_stress -= 2 * eta * deviatoric_strain_rate; so hydrostatic pressure is added as +p I. That sign only changes magnitudes, not principal directions — so it does not affect the computed SHmax angle (it only affects vector lengths). You can ignore this for angle comparisons.
quadrature vs node and averaging
ASPECT computes stresses at quadrature (Gauss) points. Many postprocessors (including the stress output) then transform/average those quad-point values to node/output points (DataOut / average_quantities).
The stress field you compare from the stress postprocessor is effectively the quad→node averaged field, but maximum_horizontal_compressive_stress computes the direction per quad point and then the code writes angle/values into computed_quantities. The code then may apply a simple arithmetic average to those components before final output.
why you cannot match Python node-based results
The crucial problem: the plugin computes angles α at quad points and then those angles (or raw α components) are averaged arithmetically. Averaging angles is not a valid way to combine directions — it produces intermediate values and is especially wrong when some quad points choose α and some choose α+90° (the plugin’s stationary-point logic can yield orthogonal solutions). Therefore even if you enable averaging, the averaged “angle” from the plugin will generally NOT match an angle computed from the stress tensor evaluated at nodes (your Python calculation).
In short: (a) quad-point → node arithmetic averaging of angles distorts directions, and (b) different quad points may pick orthogonal stationary solutions, so a simple average can’t recover the nodal principal direction.
Recommendations
For correct, comparable results use the stress tensor output and compute SHmax in Python at the nodes (or perform a proper projection of quad-point stresses to nodes first). That is the simplest and most robust approach.
If you want the plugin to be useful directly in ASPECT: implement a nodal projection (local vertex interpolation or global L2 projection) of the stress tensor, then compute angles from these nodal stresses before writing output.
I’m also a beginner with ASPECT and deal.II, so I can be mistaken on low-level API details, but from the code view above this is why plugin-angle output and Python-node output differ. My suggestion: postprocess the stress tensor in Python (node-based) to compute SHmax for reliable comparison.
@tiannh7 These are all excellent points that should probably be fixed in the code base. You could open github issues for each thing that could be improved/fixed, and we would of course be excited if you also wrote patches given that you probably know the code better than anyone at this point Some feedback about some of the points:
You’re entirely right that computing angles and then averaging is wrong. The averaging should happen first.
Adding the pressure doesn’t affect the angle, or the selection of which direction is dominant, but at the end of the day I believe you also want to output the magnitude of the max horizontal stress, and for that you do need the pressure.
Thank you very much for your confirmation and valuable feedback! I’m glad the observations make sense from the development perspective.
You’re absolutely right — while the pressure sign doesn’t affect the direction, including it is essential for correctly computing the magnitude of the maximum horizontal compressive stress, which is indeed important for quantitative analysis.
Following up on your suggestion, I’ve opened an issue to track this problem. In the issue, I summarized the core problem — that averaging angles instead of stresses leads to incorrect SHmax directions.
I’ve also started working on a PR Improvements to maximum_horizontal_compressive_stress Plugin to implement a more robust approach, which will aim to compute SHmax based on properly interpolated nodal stresses rather than averaged angles from quadrature points.
@XHL-CUG By the way, I’ve performed a Python-based postprocessing script to compute SHmax from nodal-averaged stress tensors, and the results show good agreement with the expected orientations. When comparing these node-based results to the raw stress output from ASPECT, the angles are much more consistent and physically reasonable.
I’m still learning the deeper aspects of deal.II and ASPECT’s postprocessing pipeline, so any guidance or code review would be greatly appreciated as I refine the implementation.
Some feedback about some of the points: