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

$ASPECT_DIR/benchmarks/viscoelastic_stress_build-up

The analytical solution is given in a comment at the top of the parameter file

viscoelastic_stress_build-up.prm:

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

{

return 1;

}

/

* @}

/*

/

* Returns the viscosity value on the right half of the domain,

* typically 1 or 1e6

/*

double get_eta_B() const

{

return eta_B;

}

/

* Returns the background density of this model. See the

* corresponding member variable of this class for more information.

/*

double get_background_density() const

{

return background_density;

}

private:

/

* Viscosity value on the right half of the domain, typically 1 or

* 1e6

/*

double eta_B;

/

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

*/

double background_density;

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?