Add a wrapper for Abaqus

[vigneshwaran-radhakrishnan]

I'm looking for a wrapper for calling ExaCMech from Abaqus/UMAT for two reasons.
First, the wrapper will make ExaCMech more accessible to many users of UMAT. Nowadays, UMAT can be called from several FE/FFT codes, including ExaConsist. Second, ExaConsist is still being improved to support periodic boundary conditions, controlled triaxiality, and Lode loading conditions. While these features are being implemented, a wrapper for ExaCMech will allow users to run jobs in Abaqus using ExaCMech's material models.

[rcarson3]

@vigneshwaran-radhakrishnan given I don't really use Abaqus and I last looked at the UMAT interface ~5 years ago when I added UMAT support into ExaConstit: https://github.com/LLNL/ExaConstit/blob/exaconstit-dev/src/mechanics_umat.cpp#L306-L578 . You will need to help me in creating a wrapper over ExaCMech given my lack of expertise around the Abaqus / UMAT interfaces. So, I would need to know what quantities specifically does Abaqus provide and what quantities does it expect back from the model. From there, we can try and map those to what ExaCMech expects.

Also, the initialization of models is another thing to be mindful of as best practice from the ExaCMech side of things would be to only create the object once and then re-use it for the lifetime of things rather than recreating the object every material point call. We would also need to instantiate the state / history variables during the very first time step / first call per quadrature point as well to make sure all of the lattice orientations and history variables are instantiated correctly. Hopefully, Abaqus has a way for that to be done for UMATs. If not additional logic will be required to guard against such things.

I will also mention that I'm not opposed to adding these wrappers to ExaCMech as we can have them only built when an optional CMake option is defined, but I probably wouldn't be the one to create them due to time constraints. However, I would be able to do of mentoring type thing to help guide someone into creating that initial prototype / implementation of a wrapper.

[vigneshwaran-radhakrishnan]

@rcarson3 Thanks. At each step, UMAT provides the following data for every integration point:

1. Initial data: All user-specified parameters
2. Previous step data: Cauchy stress, log strain, deformation gradient, and state variables (including all user-defined slip system orientations, hardening parameters, etc..)
3. Current step data: Deformation (strain) rate (D) and rotation rate

UMAT expects to return the updated Cauchy stress and the Jacobian matrix (del Stress/del Strain), and that's it.

I can put effort and time on to it. Guidance would be great. What is a good starting point? Start with calling pyecmech (or the jax version) fior every integration point from umat or something else?

[rcarson3]

You probably don't want to use the python bindings as that will be slower. Plus, the python stuff is more for prototyping new models or running material point simulation type things to get a better understanding of a material models behavior.

You can take a look at the miniapp we have for the library: https://github.com/LLNL/ExaCMech/blob/develop/miniapp/ and more specifically: https://github.com/LLNL/ExaCMech/blob/develop/miniapp/orientation_evolution.cxx and it's subsequent calls to see how an application code might call the library and what data it will provide the code. So, I'd say understanding the differences in the calls there to what UMAT expects would be a good first step.

I will admit the miniapp isn't the cleanest code right now, but it was meant to provide people with a model for how to use the library at a higher level.

Within the miniapp itself, you'll note there is an init_data call which setups all the state variables and this is only called once. In the orientation_evolution file, you'll notice we only create the material model once as we don't need to recreate our c++ classes constantly as all of their member variables stay the same through out the life of the program. For each time step, we have 3 steps we go through, a setup phase called setup_data found in the setup_kernels.cxx file, a batched material model invocation called from material_kernels.cxx and then a post-processing stage called retrieve_data in the retrieve_kernels.cxx file.

If I recall UMATs are per quadrature point calls right and are not batch quadrature point calls right? If so we'd probably just need to tell the underlying ECMech calls that it only needs to do 1 point for the calculation which is sub-optimal from a performance point of view but should be fine for now.

If Abaqus provides which newton iteration we're on for the nonlinear solve and we know the time step is 0 then that should be enough info to properly set-up the data only once ever. Since, you can put guards into the code to check for such things. I'm guessing the instantiation of the crystal orientation as it varies from element to element is an already understood solution there so I'm guessing that's something y'all already know how to handle. Outside of that, yeah I think what I mentioned up above might be the best 1st steps and we can go from there.

[rcarson3]

@vigneshwaran-radhakrishnan could you maybe drop this in Github GIST ( https://gist.github.com/ ) and then link that back here. It's a bit harder to read some of the above given the formatting issues that appear to be going on up above.

edit: Also, I think if you tell it that the file is a fortran / c++ file it should also do some auto-highlighting of the code which will likely make things easier to read.

[vigneshwaran-radhakrishnan]

@rcarson3
To answer your question in the last comment: Yes, umat is called for every quadrature point, not as a batch call.

[vigneshwaran-radhakrishnan]

@rcarson3

I have written a wrapper, which appears to be working. I successfully compiled the code with Abaqus, and it runs, though I still need to further evaluate the accuracy and convergence. As you said, I have uploaded the wrapper to this GitHub Gist: https://gist.github.com/vigneshwaran-radhakrishnan/97843d08a86fad66886f3fda1ac0a6ec

The Gist contains two files:

1. umatecmech.f90: The UMAT subroutine (written in Fortran 90) that Abaqus calls.
2. umat2exacmech.cxx: The C++ subroutine that interfaces with the ExaCMech constitutive model. (written based on orientation_evolution.cxx)

In the umatecmech.f90 file, I added an interface to call the C++ subroutine. For each increment, at every integration point, the previous stress, material properties, current state variables (including quaternions), and current velocity gradient are extracted and passed to umat2exacmech.cxx. The updated stress, state variables, and ddsdde are then returned from umat2exacmech.cxx to umatecmech.f90, which sends them back to Abaqus.

Can you please check if things are in correct place?

I have three questions for now:

1. The d_vgrad in the orientation_evolution.cxx code: Is it delta_vgrad or delta_vgrad*dt?
2. The material model updates quats based on the velocity gradient and stores it in state variables, correct?
3. I think orientation_evolution.cxx has only FCC, how do I add BCC material models to call from the wrapper?

[rcarson3]

@vigneshwaran-radhakrishnan thanks for the gist link that's much easier to read now. I have some thoughts on how things could probably be improved now that I have a better understanding of how everything works. Since, we can likely avoid creating a new material model every time step and a few other things.

As for your questions:

1. The d_vgrad in the orientation_evolution.cxx code: Is it delta_vgrad or delta_vgrad*dt?

Sorry about any confusion here that actually stands for device_velocity_gradient as it's used on the GPUs or CPUs depending on a users desire in the miniapp. So, you just need to pass in the current velocity gradient.

2. The material model updates quats based on the velocity gradient and stores it in state variables, correct?

Yes that is correct.

3. I think orientation_evolution.cxx has only FCC, how do I add BCC material models to call from the wrapper?

This is where we can do things a bit smarter as the miniapp / that orientation_evolution is the top level function and lives for the life of the program call. I think if we make use of C++ singleton patterns here and the makeMatModel() function: https://github.com/LLNL/ExaCMech/blob/develop/src/ecmech/ECMech_cases.h#L47 then supporting BCC materials should be "simple". It would also allow the material model to only be created once for the life of the program. I'll see if I can't sketch what that would look like over this week after work.

[vigneshwaran-radhakrishnan]

@rcarson3 Thanks for your input. I realized that in ExaCMech, v_grad represents delta_velocity/delta_time, whereas in Abaqus, it only refers to delta_velocity. I’ve updated the Fortran code accordingly. I also found another mistake—(dstran33 + drot) / dt is not the velocity gradient; it should be (dstran33 + omega) / dt, where omega is the spin derived from the incremental rotation tensor, drot.

https://gist.github.com/vigneshwaran-radhakrishnan/97843d08a86fad66886f3fda1ac0a6ec#file-umatecmech-f90

The code is converging and working well so far, but the accuracy still needs to be verified. Do you have any code snippets for adding BCC here?

[vigneshwaran-radhakrishnan]

@rcarson3 Please take a look at the gist again, I was able to add BCC (voce and MTS kinetics) into the wrapper. Let me run some test cases between Exaconstit and Abaqus, and I will update you.

[rcarson3]

@rcarson3 Thanks for your input. I realized that in ExaCMech, v_grad represents delta_velocity/delta_time, whereas in Abaqus, it only refers to delta_velocity. I’ve updated the Fortran code accordingly. I also found another mistake—(dstran33 + drot) / dt is not the velocity gradient; it should be (dstran33 + omega) / dt, where omega is the spin derived from the incremental rotation tensor, drot.

https://gist.github.com/vigneshwaran-radhakrishnan/97843d08a86fad66886f3fda1ac0a6ec#file-umatecmech-f90

The code is converging and working well so far, but the accuracy still needs to be verified. Do you have any code snippets for adding BCC here?

The velocity gradient term still seems like it might be off to me. Since $\nabla V = \dot{F} F^{-1}$ this shouldn't be equivalent to the (dstran33 + omega) / dt term you have above. Based on what I have in my ExaConstit UMAT wrapper and my notes here: https://github.com/LLNL/ExaConstit/blob/exaconstit-dev/src/mechanics_umat.cpp#L384-L387 I think a better calculation would be to do something like: (dfgrd1 - dfgrd0) / dt * inv(dfgrd1) as that essentially approximates the $\dot{F}$ term up above.

Next, I don't know if you've updated things in the ExaConstit post-processing / pre-processing function calls as it looks like you're borrowing the miniapp logic, but you'll likely want to be mindful of how Abaqus and standard Voigt notation differ which means some of the shear terms need to re-oriented. I had to do this in ExaConstit for my UMAT wrapper for both the stress and material tangent terms as seen in this function call which sets everything up and then returns things for ExaConstit: https://github.com/LLNL/ExaConstit/blob/exaconstit-dev/src/mechanics_umat.cpp#L306-L578

[rcarson3]

@vigneshwaran-radhakrishnan outside of my above comment I'll try and take a more in-depth look at the other portions of the code that you have tomorrow when I should have a bit more time in the afternoon.

[vigneshwaran-radhakrishnan]

@rcarson3 You are correct. I missed the voigt notation ordering part. I swapped the 4th (12 component), and 6th (23 compoents) columns and/or rows accordingly after I receive stress, ddsdde arrays back from exacmech into the fortran wrapper. I also had a mistake in the line do i=1,4; quats(i)=statev(9+i); enddo where I ealrier had do i=1,4; quats(i)=statev(10+i); enddo. I updated the code in the gist.

Regarding the velocity gradient: We know that \delta V = \dot F F^{-1}= D + \Omega. So, the velocity gradient can be obtained using either the deformation gradients (F_1, F_0) or the rate of deformations (D) and Spin (\Omega). I checked both approaches, and for the case I’m currently testing, they give the same result. However, as discussed by Nolan et. al 2022 (https://www.sciencedirect.com/science/article/pii/S1751616121005713), when using local orientations in Abaqus, the deformation gradient provided by Abaqus is not exactly what we expect. It is rotated to some intermediate frame. Because of this, I’m a bit hesitant towards using the deformation gradients than using the D and \Omega. What are your thoughts?

[rcarson3]

@vigneshwaran-radhakrishnan that's sounds about right that Abaqus would do things in an undocumented way which could cause less than ideal issues. I'd probably have to look at exactly what the dstran and drot terms that Abaqus is returning again as I do have a feeling that might be missing some terms as well. Since if I did my math right vgrad (L) can also be $L = \dot{V} V + V^{-1} + V \dot{R} R^{T} V^{-1}$ which due to this means the spin contribution has both the stretch and rotation contributions in it. However, it might be that it's okay to use your way up to a point as long as the time steps remain modest and there aren't any large increments in deformation in for a given step? It's been a while since I've worked through all the math related to these quantities though, so I could also be missing some other relationships here...

[vigneshwaran-radhakrishnan]

@rcarson3 Two questions:

1. How to pipe back when the material model did not converge in exacmech back to umat? I see it prints out "..failed to converge..after some iteration..", but I do not see any indicator/flag piping back.

2. What is the unit system? Between Exaconsist and Abaqus, one major difference I see is the saturation stress. In Exaconsit the (averaged) saturation stress is approx 0.35 Gpa for pulling along 100 direction, whereas in Abaqus it is 0.035Gpa. Am I missing something?

[rcarson3]

@vigneshwaran-radhakrishnan sorry for the slow reply I totally forgot to reply to this earlier in the week...

@rcarson3 Two questions:

1. How to pipe back when the material model did not converge in exacmech back to umat? I see it prints out "..failed to converge..after some iteration..", but I do not see any indicator/flag piping back.

So in the past for some of my auto-time stepping methods in ExaConstit (still need to add this to the public facing code...) I've used a try {...} catch {...} block:

   // We can now make the call to our material model set-up stage...
   // Everything else that we need should live on the class.
   // Within this function the model just needs to produce the Cauchy stress
   // and the material tangent matrix (d \sigma / d Vgrad_{sym})
   bool succeed_t = false;
   bool succeed = false;
   try{
      if (mech_type == MechType::UMAT) {
         model->ModelSetup(nqpts, nelems, space_dims, ndofs, el_jac, qpts_dshape, k);
      }
      else {
         // Takes in k vector and transforms into into our E-vector array
         P->Mult(k, px);
         elem_restrict_lex->Mult(px, el_x);
         model->ModelSetup(nqpts, nelems, space_dims, ndofs, el_jac, qpts_dshape, el_x);
      }
      succeed_t = true;
   }
   catch(const std::exception &exc) {
         // catch anything thrown within try block that derives from std::exception
         MFEM_WARNING(exc.what());
         succeed_t = false;
   }
   catch(...) {
      succeed_t = false;
   }

   MPI_Allreduce(&succeed_t, &succeed, 1, MPI_C_BOOL, MPI_LAND, MPI_COMM_WORLD);

   if (!succeed) {
      throw std::runtime_error(std::string("Material model setup portion of code failed for at least one integration point."));
   }

I'm not too sure the MPI_Reduce all would be necessary on your part of things as I'm guessing the UMAT just needs a flag passed up to signal that aspect of things.

2. What is the unit system? Between Exaconsist and Abaqus, one major difference I see is the saturation stress. In Exaconsit the (averaged) saturation stress is approx 0.35 Gpa for pulling along 100 direction, whereas in Abaqus it is 0.035Gpa. Am I missing something?

The units should be unitless for things. So, if we assume the meshes are the exact same dimensions and you're going out to the same strain and at the same strain rate you should be getting identical values. ExaConstit also reports the volume average value for this which maybe is causing issues? What value is Abaqus reporting / maybe a better question is from what is Abaqus reporting that value?