Porting a constitutive modeling library from MOOSE to MFEM

[hugary1995]

Hi MFEM developers,

MOOSE has integrated with a constitutive modeling library named NEML2. I would like to explore the possibility of porting it over to MFEM because MFEM has better support for programming models on different platforms, mainly CUDA and HIP.

In a nutshell, NEML2 takes a vector of strain and maps it to a vector stress, along with a vector of dstress/dstrain.

My questions are

1. How do I get this "vector of strain" in MFEM? I see that MFEM does assembly using various "integrators" which seem to be element-level routines. It seems to me (without much further detailed reading) that these element-level routines maps from element degrees of freedoms from the solution vector all the way to element residual/Jacobian. I wonder if it is possible to break up the integrator "in the middle", by which I'm imagining one integrator to get me the vector of strain, and another integrator to map from vector of stress to the linear/bilinear form.
2. I guess depends on the answer to 1, where should I put my customization? I'd like to contribute back to MFEM, so I'd like to start thinking about this.

Thanks,
Gary

[jandrej]

Hi @hugary1995. We're in the process of introducing DifferentiableOperator in MFEM. This interface will allow any user to ad-hoc describe FE operators by describing inputs, outputs and physical models on the quadrature point.

An example of a function on a quadrature point for elasticity looks like

   MFEM_HOST_DEVICE inline
   auto operator()(
      const tensor<real_t, dim, dim> &dudxi,
      const tensor<real_t, dim, dim> &J,
      const double &w) const
   {
      constexpr auto I = mfem::internal::IsotropicIdentity<dim>();
      constexpr real_t nu = 0.4;
      constexpr real_t mu = 0.5 * 1e6;
      constexpr real_t lambda = 2.0 * mu * nu / (1.0 - 2.0 * nu);

      auto invJ = inv(J);
      auto dudx = dudxi * invJ;
      auto F = I + dudx;
      // St. Venant-Kirchhoff model
      auto C = transpose(F) * F;
      auto E = 0.5 * (C - I);
      // second Piola–Kirchhoff stress tensor
      auto PK2 = lambda * tr(E) * I + 2.0 * mu * E;
      auto JxW = det(J) * w * transpose(invJ);
      return mfem::tuple{F * PK2 * JxW};
   }

where DifferentiableOperator calls this function individually for each quadrature point. The user is also welcome to break this computation up into several different operations if a batch computation of material responses is desired. This functionality is planned to be upstreamed in the next months.

Simply defining an integrator won't give you ability to run on accelerator devices.

[hugary1995]

That sounds great!

If I break up the operator, what would be the data type of the batched strain tensor? Would it be tensor<real_t, dim, dim, dim, dim> where the additional dimensions are respectively for elements and quadrature points? Ideally, I would like to directly map the memory (between MFEM and NEML2 data structures) without copy. I think that's possible as long as the batched strain tensor is contiguous in memory.

Anyways, I was just thinking ahead about the potential implementation. If you don't mind, would you ping me when this capability is in the upstream? I'll check back next FY in case I don't hear back.

[jandrej]

Keep in mind that all quadrature points are already batched/threaded internally for performance reasons. So you may not need any more fine grained control as we already try to saturate the full accelerator.

You can certainly output contiguous memory associated with the strain tensor from such a kernel, transform it from a material response and then integrate it in a FE sense.

[rcarson3]

@hugary1995, So it's been a while but back when both @reverendbedford and those of use at LLNL were creating our various own material model libraries, we (@nrbarton more than me) pointed @reverendbedford to some of our libraries ExaCMech (https://github.com/LLNL/ExaCMech) and ExaConstit (https://github.com/LLNL/ExaConstit) for ways we we're getting things up and running on the GPUs. ExaConstit is our MFEM based nonlinear solid mechanics code. It's an updated lagrangian code, so the primal variable is velocity rather than displacement. You might find these headers useful to see how you could do things in MFEM: https://github.com/LLNL/ExaConstit/blob/exaconstit-dev/src/mechanics_ecmech.cpp (sets up all of our velocity-based models) and https://github.com/LLNL/ExaConstit/blob/exaconstit-dev/src/mechanics_operator.cpp#L283-L416 (sets up all the quantities needed to run a material model).

Feel free to reach out to me if you'd like to talk more about the above as I'm happy to help given I worked through all of this years ago.