Add ParticleData for MPM

[xuchenhan-tri]

Also extend a few matrix utility functions to support single precision floats.

---

This change is 

[xuchenhan-tri]

+@amcastro-tri for feature review please.

Reviewable status: LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

[amcastro-tri]

question on the single precision. This is to exploit SIMD instructions? is precision a problem?
Have you quantified any of these in your dev branch?

Reviewed 18 of 21 files at r1.
Reviewable status: 9 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

a discussion (no related file):
First pass. Working...

---

multibody/mpm/particle_data.h line 32 at r1 (raw file):

                 fem::internal::LinearConstitutiveModelData<T>>;

/* The collection of all physical attributes we care about for all particles.

nit, not sure if my correction is true, just a guess. Please clarify if this is meant to store all particles in an MPM sim, or just a particular subset.

Suggestion:

The collection of all physical attributes we care about for all particles in a full MPM model.

---

multibody/mpm/particle_data.h line 41 at r1 (raw file):

  DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(ParticleData);

  using Scalar = T;

nit, is this used?

---

multibody/mpm/particle_data.h line 131 at r1 (raw file):

   @param positions     The positions of the new particles in the world frame.
   @param total_volume  The per particle volume of the new particles is given by
                        total_volume divided by the number of new particles.

so this assumes that particles are more or less uniformly distributed in space, is that right?
If that's not true, what are the iplications?

---

multibody/mpm/particle_data.h line 135 at r1 (raw file):

                        of the new particles.
   @pre total_volume > 0. */
  void Sample(const std::vector<Vector3<double>>& positions,

could you say more about the semantics and usage of this function?
Why the name "sample", instead of something like say "Add"?
Right now there is only the default ctor for empty data and then the ideas is to add particles using this method? how many times do you think this will get called in your workflow? once? or once per deformable body (if you plan to have deformable bodies or the like)?

---

multibody/mpm/particle_data.h line 139 at r1 (raw file):

              const fem::DeformableBodyConfig<double>& config);

  /* Per particle state and data. All of the following fields have the same

what's the distinction between "state" and "data".
Maybe at least worth having them in separate clusters pre-fixed with a docs like (my guess):

/* Particle state. All other data is function of the state. */
std::vector<Vector3<T>> x_;  // positions
std::vector<Vector3<T>> v_;  // velocity
 std::vector<Matrix3<T>>
      F0_;                     // deformation gradient at the previous time step

/* Particle data, quantities function of the state. */
 std::vector<T> m_;           // mass  
  std::vector<Matrix3<T>> F_;  // deformation gradient 
  std::vector<Matrix3<T>> C_;  // affine velocity field
  std::vector<T> volume_;      // reference volume
  std::vector<ConstitutiveModelVariant<T>> constitutive_models_;
  std::vector<bool> in_constraint_;  // whether the particle is participating
                                     // in a constraint

Code quote:

state and data.

---

multibody/mpm/particle_data.h line 148 at r1 (raw file):

      F0_;                     // deformation gradient at the previous time step
  std::vector<Matrix3<T>> C_;  // affine velocity field
  std::vector<T> volume_;      // reference volume

nit, I'd suggest making this clear in the name. Otherwise easy to forget in a scope far away from this comment.

Suggestion:

td::vector<T> reference_volume_;      // reference volume

---

multibody/mpm/particle_data.cc line 139 at r1 (raw file):

    x_.emplace_back(positions[i].cast<T>());
    v_.emplace_back(Vector3<T>::Zero());
    F_.emplace_back(Matrix3<T>::Identity());

nit, I belive all of these can use push_back() and since they are r-values they'll get moved just fine.

---

multibody/fem/matrix_utilities.h line 34 at r1 (raw file):

 @tparam T The scalar type, can be a double, float, or AutoDiffXd. */
template <typename T>
void PolarDecompose(const Matrix3<T>& F, EigenPtr<Matrix3<T>> R,

any thoughts, comments and/or even warning we should place in these docs regarding precision?
I wonder how sensitive (or not) the JacobiSVD might be in the polar decomposition when using "double"

[xuchenhan-tri]

I wrote some answers to this comment and inline comments below. I'll push the changes once you are done reviewing in case they mess up the diffs.

Reviewable status: 9 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/fem/matrix_utilities.h line 34 at r1 (raw file):
This uses Eigen::JacobiSVD, which is stable (if that's the concern). The relative tolerance (for off-diagonal entry being zero in Σ for A = UΣV^T) is 2.0*numeric_limits::epsilon(), so it's tight. Obviously, we won't get the same number of digits as we do when using double, but I think that's not something particularly worth documenting about.

question on the single precision. This is to exploit SIMD instructions? is precision a problem?

The bigger picture is that performance of MPM is memory-bound, so double vs. float has a significant impact on the RTR when the MPM part of the sim is dominating the runtime. Supporting floats give users a knob to make a trade-off between accuracy and speed.

Have you quantified any of these in your dev branch?
For performance, I'll later push a benchmark for the transfer kernel. For precision, no; but we can discuss what's a good benchmark problem for that is, and I'd be happy to add it.

---

multibody/mpm/particle_data.h line 32 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit, not sure if my correction is true, just a guess. Please clarify if this is meant to store all particles in an MPM sim, or just a particular subset.

Working

Yes this is the full model.

---

multibody/mpm/particle_data.h line 41 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit, is this used?

Not in this PR, but in future classes that depend on ParticleData.

---

multibody/mpm/particle_data.h line 131 at r1 (raw file):
Working

so this assumes that particles are more or less uniformly distributed in space, is that right?

Yes, they are sampled using poisson disc.

---

multibody/mpm/particle_data.h line 135 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

could you say more about the semantics and usage of this function?
Why the name "sample", instead of something like say "Add"?
Right now there is only the default ctor for empty data and then the ideas is to add particles using this method? how many times do you think this will get called in your workflow? once? or once per deformable body (if you plan to have deformable bodies or the like)?

working

---

multibody/mpm/particle_data.h line 139 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

what's the distinction between "state" and "data".
Maybe at least worth having them in separate clusters pre-fixed with a docs like (my guess):

/* Particle state. All other data is function of the state. */
std::vector<Vector3<T>> x_;  // positions
std::vector<Vector3<T>> v_;  // velocity
 std::vector<Matrix3<T>>
      F0_;                     // deformation gradient at the previous time step

/* Particle data, quantities function of the state. */
 std::vector<T> m_;           // mass  
  std::vector<Matrix3<T>> F_;  // deformation gradient 
  std::vector<Matrix3<T>> C_;  // affine velocity field
  std::vector<T> volume_;      // reference volume
  std::vector<ConstitutiveModelVariant<T>> constitutive_models_;
  std::vector<bool> in_constraint_;  // whether the particle is participating
                                     // in a constraint

working

---

multibody/mpm/particle_data.h line 148 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit, I'd suggest making this clear in the name. Otherwise easy to forget in a scope far away from this comment.

working

---

multibody/mpm/particle_data.cc line 139 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit, I belive all of these can use push_back() and since they are r-values they'll get moved just fine.

working

The return type for Matrix3::Identity() is not a Matrix3, but an expression instead. So they are not necessarily equivalent, depending on whether Matrix3 has a constructor that takes that expression (which I suspect exists).

Although looking at this again, it's better to do F_.insert(F_.end(), num_new_particles, Matrix3<T>::Identity()) in one go.

[amcastro-tri]

Reviewed 3 of 21 files at r1, all commit messages.
Reviewable status: 16 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

multibody/fem/matrix_utilities.h line 34 at r1 (raw file):

Previously, xuchenhan-tri wrote…

This uses Eigen::JacobiSVD, which is stable (if that's the concern). The relative tolerance (for off-diagonal entry being zero in Σ for A = UΣV^T) is 2.0*numeric_limits::epsilon(), so it's tight. Obviously, we won't get the same number of digits as we do when using double, but I think that's not something particularly worth documenting about.

question on the single precision. This is to exploit SIMD instructions? is precision a problem?

The bigger picture is that performance of MPM is memory-bound, so double vs. float has a significant impact on the RTR when the MPM part of the sim is dominating the runtime. Supporting floats give users a knob to make a trade-off between accuracy and speed.

Have you quantified any of these in your dev branch?
For performance, I'll later push a benchmark for the transfer kernel. For precision, no; but we can discuss what's a good benchmark problem for that is, and I'd be happy to add it.

Sounds good thanks.

---

multibody/mpm/particle_data.h line 41 at r1 (raw file):

Previously, xuchenhan-tri wrote…

Not in this PR, but in future classes that depend on ParticleData.

okay, non blocking. But usually a good practice to push in that future PR instead.

---

multibody/mpm/test/particle_data_test.cc line 65 at r1 (raw file):

    EXPECT_EQ(this->dut_.v()[i], Vector3<TypeParam>::Zero());
    EXPECT_EQ(this->dut_.C()[i], Matrix3<TypeParam>::Zero());
    EXPECT_EQ(this->dut_.F()[i], Matrix3<TypeParam>::Identity());

missing test for tau_volume(), constitutive_models(), in_constraint(),, and their mutable versions.

---

multibody/mpm/test/particle_data_test.cc line 92 at r1 (raw file):

  const TypeParam nonzero_energy =
      this->dut_.ComputeTotalEnergy(nontrivial_F, F0);
  EXPECT_GT(nonzero_energy, 0);

couldn't we compare against an exact value?
ditto below, we do know the volume per particle.

---

multibody/mpm/test/particle_data_test.cc line 222 at r1 (raw file):

      Eigen::Map<const Matrix3d>(energy.derivatives().data(), 3, 3);
  const Matrix3d tau_volume_expected = P * particle_F[0].transpose();
  EXPECT_TRUE(CompareMatrices(tau_volume, tau_volume_expected, 1e-14,

btw, cool!

---

multibody/mpm/particle_data.h line 82 at r1 (raw file):

   @note The deformation gradients (and their previous time step values) stored
   in `this` ParticleData are not necessarily used for the computation; the
   passed in F and F0 are used instead. */

without more context, I find this semantics a bit strange.
Could you help me understand why this design choice?

ditto below.

Code quote:

   @note The deformation gradients (and their previous time step values) stored
   in `this` ParticleData are not necessarily used for the computation; the
   passed in F and F0 are used instead. */

---

multibody/mpm/particle_data.h line 91 at r1 (raw file):

   @param[in] F            The deformation gradients of the particles.
   @param[in] F0           The deformation gradients at the previous time step.
   @param[out] tau_volume  The Kirchhoff stress of each particle scaled by the

ditto here, why not to modify the tau_volume stored by this data?

---

multibody/mpm/particle_data.h line 92 at r1 (raw file):

   @param[in] F0           The deformation gradients at the previous time step.
   @param[out] tau_volume  The Kirchhoff stress of each particle scaled by the
                           volume of the particle.

nit,

Suggestion:

   @param[out] tau_volume  The Kirchhoff stress of each particle scaled by the
                           reference volume of the particle.

---

multibody/mpm/particle_data.h line 100 at r1 (raw file):

   Kirchhoff stress τ. In other words, `tau_volume` is given by

    volume * τ = volume * ∂Ψ(F, F0)/∂F * Fpᵀ

this is even sranger than in the previous method. Maybe there is another design pattern we could use?
Do you have example code where you use this?

Code quote:

   @warn F and F0 are used to compute the strain, but the deformation gradients
   stored in `this` ParticleData are used to convert First-Piola stress into
   Kirchhoff stress τ. In other words, `tau_volume` is given by

    volume * τ = volume * ∂Ψ(F, F0)/∂F * Fpᵀ

---

multibody/mpm/particle_data.cc line 29 at r1 (raw file):

template <typename T>
DeformationGradientDataVariant<T> MakeStrainData(

minor, probably a good idea to keep sayiing "gradient datt" here and in other places in this PR, instead of changing it to "strain data".

---

multibody/mpm/particle_data.cc line 56 at r1 (raw file):

          const Matrix3<T>& F_p = F[p];
          const Matrix3<T>& F0_p = F0[p];
          using StrainDataType = typename std::decay_t<decltype(model)>::Data;

I think std::remove_reference_t here does the same thing. Should we use that instead to be more precise on what exactly the intent of this line is?

[xuchenhan-tri]

Reviewable status: 16 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

a discussion (no related file):

Previously, amcastro-tri (Alejandro Castro) wrote…

First pass. Working...

It looks like (from reviewable icons) you've done a full pass. Let me know if that's the case, then I'll push changes I have based on the comments.

---

multibody/mpm/particle_data.h line 41 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

okay, non blocking. But usually a good practice to push in that future PR instead.

working. I'll cover this in the unit test.

---

multibody/mpm/particle_data.h line 82 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

without more context, I find this semantics a bit strange.
Could you help me understand why this design choice?

ditto below.

There are a few threads with essentially the same question. I'll answer below in a central location.

---

multibody/mpm/particle_data.h line 91 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

ditto here, why not to modify the tau_volume stored by this data?

There are a few threads with essentially the same question. I'll answer below in a central location.

---

multibody/mpm/particle_data.h line 92 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

nit,

working

---

multibody/mpm/particle_data.h line 100 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

this is even sranger than in the previous method. Maybe there is another design pattern we could use?
Do you have example code where you use this?

The reason I have the existing signature is because I need to compute energy/stress/stress derivatives at different "states".

In addition to the use case you have in mind (i.e. given the current F, compute the current stress), I also need to compute the stress, for example, at the "next F" (i.e. the deformation gradient in Newton iterations that's a function of the unknown variables being solved).

The obvious solution is to just have two copies of ParticleData; but the states for the newton solve and the states for particles aren't exactly the same because the newton solve is really happening on grid dofs.

Looking at this a again, a less confusing way to do this is probably a free function like

Computes the stresses given the stress strain relationship and the strain.
@param[in] constitutive_models  The stress strain relationships.
@param[in] F                    Deformation gradients
@param[out] tau_volume          Volume scaled stresses.
  void ComputeKirchhoffStress(const std::vector<ConstitutiveModelVariant<T>>& constitutive_models,
                              const std::vector<Matrix3<T>>& F,
                              const std::vector<Matrix3<T>>& F0,
                              std::vector<Matrix3<T>>* tau_volume,
                              Parallelism parallelism = false) const;

instead of a class member function.

What do you think?

(FWIW, here is the callsite I'm referring to in the newton solve https://github.com/xuchenhan-tri/drake/blob/mpm-restructure/multibody/mpm/mpm_state.cc#L324-L327).

---

multibody/mpm/particle_data.cc line 29 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

minor, probably a good idea to keep sayiing "gradient datt" here and in other places in this PR, instead of changing it to "strain data".

Working

---

multibody/mpm/particle_data.cc line 56 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

I think std::remove_reference_t here does the same thing. Should we use that instead to be more precise on what exactly the intent of this line is?

Working

---

multibody/mpm/test/particle_data_test.cc line 65 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

missing test for tau_volume(), constitutive_models(), in_constraint(),, and their mutable versions.

working

---

multibody/mpm/test/particle_data_test.cc line 92 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

couldn't we compare against an exact value?
ditto below, we do know the volume per particle.

I assume you are suggesting to compute the analytic solution with some easy-to-work-with numbers. We already go through that exercise in the constitutive model unit tests. Not sure if it's worth it here.

[amcastro-tri]

Reviewable status: 16 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers

---

a discussion (no related file):

Previously, xuchenhan-tri wrote…

It looks like (from reviewable icons) you've done a full pass. Let me know if that's the case, then I'll push changes I have based on the comments.

oh yes, sorry, feature review pass is done. Thanks!

---

multibody/mpm/particle_data.h line 100 at r1 (raw file):

but the states for the newton solve and the states for particles aren't exactly the same because the newton solve is really happening on grid dofs.

Humm... I think probably easier in a f2f, but if that's the case, why is it that the Newton solver needs to know about particles at all? shouldn't the implicit Newton work on a FEM model on the support grid?

Looking at this a again, a less confusing way to do this is probably a free function like

Without understanding yet your comment on the Newton solver, this does at least match my mental picture a bit better. If we think ala system framework, we'd usually have the nice separation model/state. It'd seem this ParticleData class is more related to the concept of state rather than model. Obviously I am open to other more common C++-ish patterns where the distinction is not as clear. But from what you tell me here, where someone else (the Newton solver) needs only to update quantities based off particle data, then probably free functions (or even as part of anMpmModel class) would seem to lead to an easier to follow data flow.

[xuchenhan-tri]

Dismissed @amcastro-tri from a discussion.
Reviewable status: 10 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @xuchenhan-tri)

---

a discussion (no related file):

Previously, amcastro-tri (Alejandro Castro) wrote…

oh yes, sorry, feature review pass is done. Thanks!

Done.

---

multibody/mpm/particle_data.h line 32 at r1 (raw file):

Previously, xuchenhan-tri wrote…

Working

Yes this is the full model.

Done

---

multibody/mpm/particle_data.h line 41 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working. I'll cover this in the unit test.

Done

---

multibody/mpm/particle_data.h line 92 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working

Done

---

multibody/mpm/particle_data.h line 100 at r1 (raw file):

..., why is it that the Newton solver needs to know about particles at all? shouldn't the implicit Newton work on a FEM model on the support grid?

Unfortunately, even though the states in the newton solve lives on the grid (i.e grid velocities), there are state-dependent data that lives on particles. In particular, the equation being solved by Newton (i.e. f(x) = 0) involves the internal forces on the grid nodes, which is a function of the stress on particles, which is a function of the deformation gradient (on particle), which is a function of grid velocities (on grid). So, we need machinery to propagate a change in grid velocity through this chain to get a change in f(x).

The ComputeKirchhoffStress function is a part of that chain (input deformation gradient and outputs stress).

Hopefully this helps. Reading the code I linked above on MpmState should also help somewhat.

We should discuss f2f if this is still unclear.

---

multibody/mpm/particle_data.h line 131 at r1 (raw file):

Previously, xuchenhan-tri wrote…

Working

so this assumes that particles are more or less uniformly distributed in space, is that right?

Yes, they are sampled using poisson disc.

Done. Added a note.

---

multibody/mpm/particle_data.h line 135 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working

Done. Added a note this function can be called in multiple passes. (So, it supports the workflow of adding particles once per MPM geometry).

---

multibody/mpm/particle_data.h line 139 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working

Done.

Note that F and C are actually states in MPM.

---

multibody/mpm/particle_data.h line 148 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working

I tried splattering the world "reference" in front of all instances of "volume" and I don't like the result. It makes variables like dPdF_reference_volume too long and thus less readable in places where easily matching code with mathematical equations are useful.

---

multibody/mpm/particle_data.cc line 29 at r1 (raw file):

Previously, xuchenhan-tri wrote…

Working

Done

---

multibody/mpm/particle_data.cc line 56 at r1 (raw file):

Previously, xuchenhan-tri wrote…

Working

Done

---

multibody/mpm/particle_data.cc line 139 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working

The return type for Matrix3::Identity() is not a Matrix3, but an expression instead. So they are not necessarily equivalent, depending on whether Matrix3 has a constructor that takes that expression (which I suspect exists).

Although looking at this again, it's better to do F_.insert(F_.end(), num_new_particles, Matrix3<T>::Identity()) in one go.

Done

---

multibody/mpm/test/particle_data_test.cc line 65 at r1 (raw file):

Previously, xuchenhan-tri wrote…

working

Done.

I removed a few mutable getters that don't make sense (because the data are constant).