only average viscosity if the viscosity is computed in the material model

[jdannberg]

This pull request makes sure that the cell-wise averaging of material properties only averages the viscosity if the viscosity is actually computed in the material model (as indicted by in.requests_property(MaterialModel::MaterialProperties::viscosity)). Before this change, models with geometric averaging would crash in debug mode in case a material model did not compute the viscosity (for me, this happened with the steinberger model). This is because the viscosity is initialized with NaNs. But the average_property function only checks if values_out.size() == 0 (It has a comment that says: "if an output field has not been filled (because it was not requested), then simply do nothing -- no harm no foul", but that's not what happens in the other parts of the code, the size of in.viscosities is not actually zero if the viscosity is not requested).

I am not sure what I did is the best solution for this problem (I had to add a new argument to the average function, which needs to be handed over in many places in the code), so let me know if you have a better idea.

• [x] I have followed the instructions for indenting my code.

[tjhei]

It is a bit weird that the viscosity is now a special case here.

This is because the viscosity is initialized with NaNs

Did you check if it would be difficult to change this fact and have it size=0?

[jdannberg]

Did you check if it would be difficult to change this fact and have it size=0?

So there are two different ways that I see of doing this.

1. We change the constructor of the material model outputs. In this case, that would have to get the in.requests_property(MaterialModel::MaterialProperties::viscosity) variable and could then initialize the viscosities with size zero. I am not really sure that would make it better. The problem here is that requests_property is part of the material model inputs, but what we want to change are the material model outputs.

2. We do this it in each individual material model. So instead of just not touching the viscosity when in.requests_property(MaterialModel::MaterialProperties::viscosity) is false, the material model would have to resize it to zero. But that would only ever be a convention, and user-created models may not always follow it.

[tjhei]

okay, that is not ideal either.

How about adding an argument MaterialProperties::Property requested_properties to the averaging routine instead of the bool (you can pass in.requested_properties)?

[jdannberg]

I think that would not actually help, at least not at the moment.

In the constructor of the material model inputs, it's initialized as requested_properties(MaterialProperties::all_properties).

As far as I can tell, requested_properties is not changed anywhere in a way that the viscosity is removed from the list. At least, the implementation of requests_property at the moment is:

    template <int dim>
    bool
    MaterialModelInputs<dim>::requests_property(const MaterialProperties::Property &property) const
    {
      //TODO: Remove this once all callers set requested_properties correctly
      if ((property & MaterialProperties::Property::viscosity) != 0)
        return (strain_rate.size() != 0);

      ... some other stuff ...
    }

So it does not actually use what's stored in requested_properties.

Alternatively, should we just pass the material model inputs, and then the average() function can call in.requests_property(MaterialModel::MaterialProperties::viscosity)?

[tjhei]

I think that would not actually help, at least not at the moment.

That is annoying.

So it does not actually use what's stored in requested_properties.

We could change this. I don't think there are many places in user code that would need to be fixed. What do you think? I find it better than requiring a MaterialModelInputs object.

[jdannberg]

Yes, I think that would work, but would require going through all the places that resize in.strain_rate to 0. @anne-glerum made a pull request for this (#4248) that I totally missed and that @gassmoeller pointed me to, so we should get that merged first. I'll take a look at #4248 to see if I can figure out what might be the problem with the tests.

[tjhei]

Ah, nice! I would add the additional third argument to the material inputs and then go through all the compile errors to see if she missed a location...

[bangerth]

I don't particularly care for this patch because it treats the viscosity differently than other fields. What if we encounter that same issue with one of the other fields? Do we add more bool parameters?

An alternative would be to add code such as this to the average_property() function:

  if (std::any (viscosities.begin(), viscosities.end(), &std::is_nan))
    {
       // At least one of the elements of the array is a NaN. If so, we know that averaging
       // will result in NaNs for all elements; rather than compute those, just *set* them:
       std::fill (viscosities.begin(), viscosities.end(), numbers::signaling_nan<double());
       return;
    }

This implies a second read over all elements, every time we call the function. is_nan should be a relatively cheap function, and averaging probably doesn't cost very much compared to computing material properties in the first place, but it's a cost.

[bangerth]

What should we do here?

[jdannberg]

I am waiting for #4248. Then I can do what Timo suggested.

[gassmoeller]

#4248 has been merged, so when you have some time this can move forward.

[bangerth]

@jdannberg Want to pick this up again?