Confusion about deviators

[bangerth]

In the documentation, we explain how when we deal with compressible strains, we always compute

  eps - 1/3 trace(eps)*I

even in for 2d models where one might think that we should be using a factor of 1/2 instead.

But we have plenty of places in the code where we do things such as the following (here in stress_second_invariant.cc):

                const SymmetricTensor<2, dim> deviatoric_strain_rate = (this->get_material_model().is_compressible()
                                                                        ? strain_rate - 1. / 3 * trace(strain_rate) * unit_symmetric_tensor<dim>()   // *** Note: 1/3
                                                                        : strain_rate);

                const double eta = out.viscosities[q];

                stress += -2. * eta * deviatoric_strain_rate;
              }

            const SymmetricTensor<2, dim> deviatoric_stress = deviator(stress);                       // *** Note: 1/2
            const double stress_invariant = std::sqrt(second_invariant(deviatoric_stress));

The issue here is that deviator() is defined as

SymmetricTensor<2, dim, Number>
deviator(const SymmetricTensor<2, dim, Number> &t)
{
  SymmetricTensor<2, dim, Number> tmp = t;

  // subtract scaled trace from the diagonal
  const Number tr = trace(t) / dim;                                          // ********* Note: 1/dim
  for (unsigned int i = 0; i < dim; ++i)
    tmp.data[i] -= tr;

  return tmp;
}

and for 2d models, this means that we're using a different deviator definition.

This raises the question of whether we are doing something wrong when we're in 2d by calling deviator()? There are, as far as I can see, 12 places in source/ where we call that function.

[anne-glerum]

That does seem inconsistent. I also found places in the material models where we do:

const SymmetricTensor<2,dim> strain_rate_deviator = deviator(in.strain_rate[i]);

So not just for deviatoric stresses/postprocessors, but also for the strain rate we use 1/2 in 2D where it should be 1/3.

[gassmoeller]

I agree that this looks suspicious, but I also vaguely remember that we had a discussion about it at some point already. I think we need to take some time and talk about the correct solution at the hackathon and then document that solution, or make sure we cannot accidentally use the wrong option again.