Artificial viscosity for composition if u=0

[tjhei]

compute_viscosity contains the following code:

    // If the velocity is 0 we have to assume a sensible velocity to calculate
    // an artificial diffusion. We choose similar to nondimensional
    // formulations: v ~ thermal_diffusivity / length_scale, which cancels
    // the density and specific heat from the entropy formulation. It seems
    // surprising at first that only the conductivity remains, but remember
    // that this actually *is* an additional artificial diffusion.
    if (std::abs(global_u_infty) < 1e-50)
      return parameters.stabilization_beta *
             max_conductivity / geometry_model->length_scale() *
             cell_diameter;

but max_conductivity is 0 if this is a compositional field meaning we return 0. This can't be right.

[jdannberg]

Oh, you are right, I didn't even think about this. (And at some point I also wanted to make a test case that verified that we compute the artificial viscosity correctly for the case of compositions, but I had no idea how to do that, so if you have any ideas we could add one. For now we only have a test where the results look about right.)

But should we use the same value as for the temperature here? Even if that might not be exactly what we want, I don't really see any other parameter we could use...

[bangerth]

I think I stumbled across this before as well. If the velocity is zero, the evolution of the advected fields reverts to a pointwise ODE. I recall that this leads to instabilities if the velocity becomes zero only in some parts of the domain, or at individual points. We can "fix" this by smoothing the artificial viscosity as implemented by @rrgrove6 .

But the code block above is only relevant if the velocity is essentially zero everywhere. In that case, we have an ODE at every node of the domain, and I think we should be able to solve this in a stable way. This may be why we've never seen any issue with this.

(I will note that I'm not fond of hard coded numbers such as the 1e-50 above. My preference would be if we had something like

  abs(global_u_infty) < 1e-12*GridTools::diameter(triangulation)/(end_time-start_time)

or similar.)

[bangerth]

This is really a duplicate of #293.