Sporadic occurrences of stresses equal to NaN

[thiagordonho]

Describe the bug The computed stresses at some particles will sometimes be nan when using the Mohr-Coulomb material. Further inspection of the issue lead to the conclusion that this happens when (1-r*r) is negative and std::abs(1-r*r) > tolerance in material_utility.tcc where the invariants and Lode parameters are calculated. This particular condition results in dtheta_dr = nan because a square root of a negative number is used in its calculation:

    // Update derivative of theta in terms of R, check for sqrt of zero
    const double factor =
      (std::abs(1 - r * r) < tolerance) ? tolerance : (1 - r * r);
    dtheta_dr = -1.0 / (3.0 * sqrt(factor));

To Reproduce Steps to reproduce the behavior:

1. Compile 'cmake -DCMAKE_CXX_COMPILER=g++ .. && make -j8'
2. Run on one node with './mpm/build/mpm -f \path\to\input\file -i input-file.json'
3. On condition "locate_particles": true Using the Mohr-Coulomb
4. See error: particle outside of mesh domain

Expected behavior A clear and concise description of what you expected to happen.

Screenshots As aforementioned, if locate_particles is set to true, the simulation will stop due to a particle outside the mesh domain. Otherwise, the material point with stresses equal to nan, along with all the other material points in its cell and neighboring cells, will be removed as exemplified below:

This happens because the nan values are used to compute the internal force, nodal velocities and acceleration which are later mapped to all material points, leaving their parameters "corrupted" as well.

Runtime environment (please complete the following information):

• Any OS
• Branch develop

Additional context

• Debugging the problem, it was evidenced that the value of r that leads to this problem is usually very close to -1 (e.g. r = -1.000000000000029).
• Changing the computation of factor to use the absolute value of 1-r*r fixes the problem numerically.
• Since this issue happens in material_utility.tcc, any other constitutive model that uses the functions in this file are also susceptible to the same issue.

[bodhinandach]

@thiagordonho Thanks for figuring this out. @ezrayst we might have this issue as well.

[ezrayst]

@thiagordonho thanks so much for this. Can you review #653 and approve it? We can merge it quickly. J3 should be positive and J2 is always positive, though it is not the case numerically. Thus, I think your suggestion is good in making sure it's positive by std::abs()