search

geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
218 stars 232 forks source link

remove checks for invalid strain_rate #4874

Open tjhei opened 2 years ago

tjhei commented 2 years ago

I am confused by the following: https://github.com/geodynamics/aspect/blob/0bb443da15a5eaf96baaa5a382100a5c0d959a75/source/material_model/rheology/visco_plastic.cc#L121-L122

I think only a SymmetricTensor filled with NaN would make this condition true. Why do we catch this situation here? Do we forget to fill it somewhere? If it is needed, we should use !std::isfinite() instead.

@naliboff Do you remember why this was added?

tjhei commented 1 year ago

I am looking back at this issue and seems to me that it is there so that at timestep 0 (where velocity and as such strainrate would be 0), we do not end up using 0 as a strainrate (or a minimum strainrate), but replace it with a reference strainrate. I guess this is ok?

bangerth commented 1 year ago

The second part of the condition is true if the strain rate is zero or very nearly so. (The condition is not a great choice because it is dimensionally not correct: the left hand side has units of strain, whereas the right hand side is an absolute number.)

If you had NaNs in the strain rate, then any comparison will return false, so the condition cannot be true.