Stress output multiplies with viscosity, not averaged viscosity

[bangerth]

When outputting the stress tensor, we multiply the strain rate tensor by the viscosity. This leads to trouble when the velocity field is computed with an averaged viscosity field, because in that case we compute the stress not with the viscosity that was used, but something entirely unrelated. The resulting plots are often wildly wrong with over- and undershoots.

Would y'all agree that this is wrong?

[dsarahstamps]

[bangerth]

First is the strain rate, roughly piecewise constant as expected. The second is the stress as computed, which doesn't work because it multiplies the piecewise constant strain rate by a function that has a jump somewhere within the lowest row of red cells in the first picture. In other words, the greenish line that winds its way through the second picture appears completely artificial to me.

[gassmoeller]

I agree that this is strange. What do you propose to do? Should the stress output use the averaged material properties? This would make it slightly inconsistent with other postprocessors (e.g. the material properties), but I guess the averaging matters most for viscosity and stress anyway, so I would be ok with that.

But additionally for the call to average() you need a cell iterator, which is not available in the visualization postprocessors. That hit me and @jdannberg already at several occasions, but we usually worked around it somehow, e.g. by providing fallbacks in the material model. This will not be possible for you though, since you need the cell reference at least for the projection averaging I guess.

[bangerth]

That's indeed a good point, @gassmoeller . No cell no averaging. I didn't think of this. We do have an open issue in deal.II (see https://github.com/dealii/dealii/issues/391) that would allow us to pass the cell, but this isn't the case right now. This issue may simply have to wait till we figure this out in deal.II.

In general, I don't think all postprocessors should use the averaged fields. This one seems to be special in my mind, and the pictures clearly show why the current approach isn't working.

For the moment we (try to) work around this by using the "averaging" material model that does the averaging as part of the evaluate() function. Consequently, the postprocessors get to also see the averaged material properties.

[jdannberg]

I think it would be great to have the cell in the visualization postprocessors once this option is available in deal.ii. I have several material and heating models where I need the cell reference (to get the solution from an old time step, or to get some input that is not part of the material model inputs), and right now there is no way to visualize these properties.

For averaging in the postprocessors: There are a few others where it would probably make sense to have averaged material properties. One example is the shear heating, also because the stress is included in the computation of the heating rate. But otherwise I agree that not all of the postprocessors have to use the averaged properties.

[gassmoeller]

Since we now have the cell reference in the postprocessors, a call to MaterialModel::MaterialAveraging::average should do the trick. In case someone is interested this should be a good starter project to learn about the material averaging.

[jaustermann]

@gassmoeller and @bangerth - I'd be happy to give this a try. I looked at visualization/stress.cc and it doesn't seem addressed yet but wanted to check just in case given this is an old issue. If this issue is superseded by other work let me know, otherwise, I'd be happy to work on changing the viscosity from being extracted at each quadrature point to using the cell average.

[gassmoeller]

Yes that would be a nice project. Take a look at https://github.com/geodynamics/aspect/pull/4076 where Wolfgang implemented this for the material properties postprocessor, there is a slight catch for project_to_Q1 averaging, but everything else should be pretty straightforward.

[jaustermann]

I wanted to set up an easy test case to recreate this in order to check that adjusting the code actually fixes it. Below is a rising sphere in a box with viscosity contrasts and the shear stress does have these spikes at the viscosity boundary. Just double-checking - is that where the issue shows? I don't see this in the total stress magnitude though.

[bangerth]

In your visualization, both the strain and stress are shown as piecewise constants. Are you setting Point-wise stress and strain to true in your input file? What happens if you set it to false?

Separately, it might be easier to understand the output if you set Interpolate output to false.

[jaustermann]

Hi Wolfgang. The Point-wise stress and strain was already set to false in the above runs. The polynomial degree in the velocity is the default (2). The Interpolate output was set to true but if that is changed to false, the result doesn't change too much - the strain rate looks pretty much the same, the shear stress looks like this:

I'm not sure why the stress and strain are shown as piecewise constants. I'm attaching the parameter file here: depth_dependent_box_none_simple.txt And the output parameter file with all variables is here: parameters.txt

[jaustermann]

Ping!