I recently encountered a bug in the MFIX version 23.1 when I ran a DEM simulation involving superquadratic particles and a gravity vector set to zero in all directions.
The bug occurred in the subroutines Sq_angle_with_gravity in the file des/sq_properties.f.
Using zero gravity led to a division by 0.

I solved the problem by modifying the line:
“gravity=gravity/gravity_norm”
to:
“IF(gravity_norm == 0.0) THEN
gravity=(/0,1,0/)
ELSE
gravity=gravity/gravity_norm
ENDIF”

I don’t know if it’s the most appropriate solution for everyone, but I wanted to inform you of this little bug.

Thanks for the suggestion. Your change should be safe - when the magnitute of the g-vector is zero, the angle is undefined, so it shouldn’t matter what value the function Sq_angle_with_gravity returns in that case.

I would probably write:

IF(gravity_norm == 0.0) THEN
Sq_angle_with_gravity = 0
ELSE

i.e. skip the rest of the calculation and just return 0 if |g| is 0.

But, I took a closer look at the code to see where Sq_angle_with_gravity is used, and it seems that this value is not used at all. The call sites look like this:

There is a commented-out line storing cos_angle into des_usr_var but other than this, the variable cos_angle is not used. So I think the right fix is to remove these unneeded calls completely. I will do this for the next MFiX release.