Mack and I are testing the quadratic particle interpolation on problems in solid mechanics / plate tectonics that John Naliboff has added to the benchmarks directory in the ASPECT repository. In particular, we are computing the following benchmark from
The analytical solution is given in a comment at the top of the parameter file
The analytical solution for viscoelastic stress build-up in an incompressible medium with a constant
strain-rate is: simga_ij = 2 * edot_ii * eta * (1 - e^(-mu*t/eta)), where sigma_ij is the elastic deviatoric stress tensor, edot_ij is the deviatoric strain-rate, eta is the background viscosity, mu is the shear modulus and t is time.
My thought is to essentially copy the structure of the computational version of the exact solution for SolCx, which is actually in SolCx.h not SolCx.cc, a portion of which I quote here (sorry for the length):
* @name Reference quantities
virtual double reference_viscosity () const
* Returns the viscosity value on the right half of the domain,
* typically 1 or 1e6
double get_eta_B() const
* Returns the background density of this model. See the
* corresponding member variable of this class for more information.
double get_background_density() const
* Viscosity value on the right half of the domain, typically 1 or
* A constant background density over which the density variations
* are overlaid. This constant density has no effect on the dynamic
* pressure and consequently on the flow field, but it contributes
* to the total pressure via the adiabatic pressure. We use this
* field to support our claim in the first ASPECT paper that the
* accuracy of the solutions is guaranteed even if we don’t subtract
* the adiabatic pressure in our computations.
My question is this: "Are all of the quantities in the exact (analytical) solution: simga_ij, edot_ii, eta, and mu stored somewhere where we can easily access them, or will we have to compute some them ourselves in our solution.cc file?