[mpm] Add SparseGrid

[xuchenhan-tri]

In addition:

• Fix a couple of stale documentation in SpGrid
• Add a sugar for GridState.
• Add MassAndMomentum struct.

---

This change is 

[xuchenhan-tri]

+@amcastro-tri for feature review (or delegation) please. This PR is mostly introducing a wrapper around SpGrid; it's not crucial for understanding the core algorithm, so it's fine to delegate to other reviewers if needed.

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

[amcastro-tri]

checkpoint before lunch. A bunch of API design/questions.

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

---

multibody/mpm/sparse_grid.h line 17 at r1 (raw file):

namespace internal {

/* SparseGrid is a 3D grid that is sparsely populated with GridData (see

minor, consider SparseMpmGrid, since having both SpGrid and SparseGrid in the same namespace is a bit confusing.

I suggested that name based on my understanding which is:

1. SpGrid is the true wrapper around the SPGrid:: namespace. This is generic code, not MPM specific.
2. The new class in this file is now specific to MPM, storing (MPM) GridData.

The docs should probably state this distinction to introduce those not familiar with the design.

---

multibody/mpm/sparse_grid.h line 47 at r1 (raw file):

  /* Constructs a SparseGrid with grid spacing `dx` in meters.
   @pre dx > 0. */
  explicit SparseGrid(double dx);

any reason not to pass T dx? and avoid the multiple static_cast<T>(dx_) in the source?

---

multibody/mpm/sparse_grid.h line 64 at r1 (raw file):

                         SparseGrid's SPGrid and dx and the particle positions.
  */
  void Allocate(const ParticleSorter& particles);

I have a question in the order of operations here.
I imagine you want to pass an already "sorted" sorter here. That is, I already called Sort() with a valid SpGrid. However, this class is the owner of the underlying SpGrid, and pressumably the idea was to hide it inside this class?

Also, I realize I have new questions that did not come up in the review of ParticleSorter. Well, one question. Are methods data_indices(),base_node_offsets(), colored_ranges() and GetActiveBlockOffsets() conditioned to be called "after" the call to Sort()? or is their return useful if I didn't call Sort() yet?

---

multibody/mpm/sparse_grid.h line 73 at r1 (raw file):

  Pad<Vector3<T>> GetPadNodes(const Vector3<T>& q_WP) const;

  /* Given the SpGrid offset of a grid node, returns the grid data in the pad

question on the design here. As before, is the idea to "hide" the SpGrid under the hood of this class?
That is, should I expect mentions of SpGrid in the docs or elements such as offsets in the API? In this case, how would a user obtain the offset?

---

multibody/mpm/sparse_grid.h line 80 at r1 (raw file):

  }

  /* Given the offset (see multibody::mpm::internal::SpGrid) of a grid node,

that still doesn't tell me the recommended way to iterate through nodes in this class.
Maybe an example would help?

Code quote:

 (see multibody::mpm::internal::SpGrid)

---

multibody/mpm/sparse_grid.h line 96 at r1 (raw file):

  /* Returns grid data from all non-inactive grid nodes (see
   GridData::is_inactive) as a pair of world space coordinate and grid data of

is this concept of "active" somewhat related to the concept of "supported" introduced at the top of this file? or are those orthogonal?

Code quote:

GridData::is_inactive)

---

multibody/mpm/sparse_grid.h line 103 at r1 (raw file):

  /* Computes the mass, linear momentum, and angular momentum (about world
   origin) on the grid.
   @note Testing only. */

minor,

I could see this useful for logging or other users. I'd personally remoe the @note. Or maybe @note Useuful for testing to leave the door open to other uses.

Code quote:

 @note Testing only.

---

multibody/mpm/sparse_grid.h line 107 at r1 (raw file):

  /* Returns the SpGrid underlying this SparseGrid. */
  const SpGrid<GridData<T>>& spgrid() const { return spgrid_; }

more about the workflow. Is this something I want? or this is meant only for testing?

[xuchenhan-tri]

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

---

multibody/mpm/sparse_grid.h line 17 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

minor, consider SparseMpmGrid, since having both SpGrid and SparseGrid in the same namespace is a bit confusing.

I suggested that name based on my understanding which is:

1. SpGrid is the true wrapper around the SPGrid:: namespace. This is generic code, not MPM specific.
2. The new class in this file is now specific to MPM, storing (MPM) GridData.

The docs should probably state this distinction to introduce those not familiar with the design.

I updated the doc to state the distinction, but I'm not convinced yet about the name change. The spelling of SpGrid and SparseGrid are different enough that no one will mistake them for the same class.

Shoving the word MPM into a class inside the mpm namespace doesn't add much in explaining that this the class used in the MPM transfer instead of the class that's a container.

---

multibody/mpm/sparse_grid.h line 47 at r1 (raw file):

any reason not to pass T dx?

There's no good reason to use a float for dx to me. The reason I allow float for GridData is to reduce the memory footprint which has a significant impact on performance in transfers.

and avoid the multiple static_cast<T>(dx_) in the source?

I don't understand the question. There's no static_cast<T>(dx_) in the source.

---

multibody/mpm/sparse_grid.h line 64 at r1 (raw file):

I imagine you want to pass an already "sorted" sorter here. That is, I already called Sort() with a valid SpGrid. However, this class is the owner of the underlying SpGrid, and pressumably the idea was to hide it inside this class?

I don't think I fully understand the question here. But I'll attempt an answer. Yes, ParticleSorter::Sort() needs to be called before calling Allocate() here, as documented. However, the quantifier you added ("with a valid SpGrid") is ambiguous. To be clear, an SpGrid is always "valid" for the purpose of this discussion. In particular, you don't have to call Allocate on SparseGrid to obtain a "valid" SpGrid to be passed to ParticleSorter::Sort() (because, obviously, that creates an impossible cycle).

I'll work on making that clear. Perhaps there's a piece of SpGrid that I can split off (similar to SpGridFlags) so that I don't need to pass the entire SpGrid to ParticleSorter for the sort.

Are methods data_indices(),base_node_offsets(), colored_ranges() and GetActiveBlockOffsets() conditioned to be called "after" the call to Sort()? or is their return useful if I didn't call Sort() yet?

Yes, you need to call Sort() for those accessors to return meaningful result, according to the doc for ParticleSorter.

After `Sort()` is called, the ParticleSorter provides information about the particles and its background grid based on the result of the sort.

---

multibody/mpm/sparse_grid.h line 73 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

question on the design here. As before, is the idea to "hide" the SpGrid under the hood of this class?
That is, should I expect mentions of SpGrid in the docs or elements such as offsets in the API? In this case, how would a user obtain the offset?

Yes, I want to hide SpGrid as much as possible, but we need to expose the concept of offsets for performance reasons.

All of this is internal code and the end user would obviously never interact with the concept of an offset. For a developer-as-a-user of this code, they can obtain offsets through SpGrid and ParticleSorter as I explain in the thread immediately below this one.

I will, however, add sugar in future PRs to help loop through the grid nodes people care about in the most efficient way so that the exposure of "offsets" to algorithm developers are minimal too.

---

multibody/mpm/sparse_grid.h line 80 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

that still doesn't tell me the recommended way to iterate through nodes in this class.
Maybe an example would help?

I don't think this is the place to teach developers how to iterate through the grid nodes. This class really only does three things stated in the class documentation.

1. allocate memory for the grid data (if necessary).
2. access the grid data (in chunk through GetPadData).
3. modify the grid data (in chunk, through SetPadData).

If someone wants to iterate the grid, they'd have to obtain the list of node offsets somewhere else. One way to do that is through SpGrid APIs (e.g. SpGrid::IterateGridWithOffset) or use the list of base node offsets curated by ParticleSorter::base_node_offsets. The latter pattern is used heavily in the transfers, so I'll later add a sugar for that.

---

multibody/mpm/sparse_grid.h line 96 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

is this concept of "active" somewhat related to the concept of "supported" introduced at the top of this file? or are those orthogonal?

Being "supported" is a necessary but not sufficient condition for a grid node to be "active".

Particle data only gets transferred to supported grid nodes, but before the transfer happens, the grid nodes are not active.

After grid allocation, some allocated grid nodes are supported, the others are not, but all are inactive. After particle to grid transfer (P2G) happens, the supported grid nodes are no longer inactive, but the non-supported grid nodes remain inactive.

---

multibody/mpm/sparse_grid.h line 103 at r1 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

minor,

I could see this useful for logging or other users. I'd personally remoe the @note. Or maybe @note Useuful for testing to leave the door open to other uses.

Done

---

multibody/mpm/sparse_grid.h line 107 at r1 (raw file):

or this is meant only for testing?

This is the entry point into SpGrid if you want to get your hand dirty with the grid data structure. I could label it as advanced if you think that's appropriate. I don't particularly see the appeal since this is all internal.

Another use of this function is in ParticleSorter::Sort which needs some information about how the data is laid out. You already pointed this out in an earlier thread.

Is this something I want?

This is really a question about your personal preference, so I won't be able to tell :)

[amcastro-tri]

Reviewed 1 of 1 files at r2.
Reviewable status: 7 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)

---

multibody/mpm/sparse_grid.h line 64 at r1 (raw file):

Previously, xuchenhan-tri wrote…

I imagine you want to pass an already "sorted" sorter here. That is, I already called Sort() with a valid SpGrid. However, this class is the owner of the underlying SpGrid, and pressumably the idea was to hide it inside this class?

I don't think I fully understand the question here. But I'll attempt an answer. Yes, ParticleSorter::Sort() needs to be called before calling Allocate() here, as documented. However, the quantifier you added ("with a valid SpGrid") is ambiguous. To be clear, an SpGrid is always "valid" for the purpose of this discussion. In particular, you don't have to call Allocate on SparseGrid to obtain a "valid" SpGrid to be passed to ParticleSorter::Sort() (because, obviously, that creates an impossible cycle).

I'll work on making that clear. Perhaps there's a piece of SpGrid that I can split off (similar to SpGridFlags) so that I don't need to pass the entire SpGrid to ParticleSorter for the sort.

Are methods data_indices(),base_node_offsets(), colored_ranges() and GetActiveBlockOffsets() conditioned to be called "after" the call to Sort()? or is their return useful if I didn't call Sort() yet?

Yes, you need to call Sort() for those accessors to return meaningful result, according to the doc for ParticleSorter.

After `Sort()` is called, the ParticleSorter provides information about the particles and its background grid based on the result of the sort.

I am trying to undertand the purpose of yet another separate class. For instance, why didn't we just extend the SpGrid to implement the suggar here to work with pads?

That being said, what would this look in an actual program? my try:

// .... Load/compute particle locations ...
std::vector<Vector3d> particle_positions;

// Setup grid
SparseGrid grid(1e-2);
const SpGrid& spgrid = grid.spgrid();
ParticleSorter sorter;
sorter.Sort(spgrid, grid.dx(), particle_positions);
grid.Allocate(sorter);

// From now on we are ready to perform grid operations ...

Per code above it'd seem that the state of the class depends on the input sorter at Allocate() time. Would it make sense for the particle sorter to be a member of this class? So that the code looks like:

// .... Load/compute particle locations ...
std::vector<Vector3d> particle_positions;

// Setup grid
SparseGrid grid(1e-2);
grid.SetPositions(particle_positions);  // Class Sort() and Allocate() under the hood.

Say now you want to write your P2G, do you work with SparseGrid or with SpGrid? I am confused about when to use one or the other.

---

multibody/mpm/sparse_grid.h line 96 at r1 (raw file):

Previously, xuchenhan-tri wrote…

Being "supported" is a necessary but not sufficient condition for a grid node to be "active".

Particle data only gets transferred to supported grid nodes, but before the transfer happens, the grid nodes are not active.

After grid allocation, some allocated grid nodes are supported, the others are not, but all are inactive. After particle to grid transfer (P2G) happens, the supported grid nodes are no longer inactive, but the non-supported grid nodes remain inactive.

great, thanks for the explanation. Could you place this definition of "active" somewhere? I don't find it in GridData either.

---

multibody/mpm/sparse_grid.cc line 30 at r2 (raw file):

      for (int k = 0; k < 3; ++k) {
        const Vector3<int> offset(i - 1, j - 1, k - 1);
        result[i][j][k] = static_cast<T>(dx_) * (base_node + offset).cast<T>();

example of static_cast mentioned above

[xuchenhan-tri]

@amcastro-tri I'll move the discussion from slack to here for more visibility.

You mentioned

 IMO, even though I see you are trying to split particles from grid,  there's really no obvious reason to split between SparseGrid, Particles and ParticlesSorter. To me all those are tightly intertwined. They all depend on the particles positions, and therefore having sparate classes (in my mind thus far) seems to require the careful coordination of these classes to be in sync.I think (again, very high level view), that SparseGrid could sotre the Particles (i.e postions, data and sorter) so that the coordination is done with a single entry point. Then, the ParticleSorter::Iterate() could be exposed through a SparseGrid::Iterate().

I gave this careful thought but decided against taking the suggested pattern over the current implementation. First, I think shoving particles into the grid class is a hard no. Particles and grid are two related but absolutely distinct concepts, and anyone with working knowledge of MPM would find it confusing to have to reach for particles inside a grid class. As it is right now, the SparseGrid class is essentially SpGrid instantiated on actual MPM GridData, and it's a very nice abstraction. You can query data on grid nodes, talk about grid spacing, and meaningfully discuss concepts like "grid momentum" (through grid.ComputeTotalMassAndMomentum()). Imagine the confusions sowed when there's a competing grid.particles().ComputeTotalMassAndMomentum().

The next step down is hiding ParticleSorter inside SparseGrid. Ostensibly it brings two benefits:

1. it removes the coordination problem between sorting particles and allocating for grid memory.
2. it hides some internal details about SpGrid; we can hide the 1d indexing of the grid.

Re 1: It seems to me that the coordination problem between particle sorting and grid allocation can simply be done with doing them in a lock step. Something like

void SparseGrid::Allocate(const ParticleData<T>& particle_data, ParticleSorter* sorter) {
  sorter->Sort(spgrid_, dx_, particle_data.x());
  spgrid_.Allocate(sorter->GetActiveBlockOffsets());
}

Re 2: For this to happen, the turn any 1d index APIs (e.g. SparseGrid::SetPadData() and SparseGrid::GetPadData()) into implementation details and implement the sugar for iteration over particles in SparseGrid like you suggested. However, that would prevent reusing the iteration helpers. For example, to facilitate autodiff (e.g. for testing), I have a different grid data structure that can accommodate AutoDiffXd as a scalar (SpGrid doesn't), and I'd need exactly the same type of iterations. In that sense, separating the iteration implementation from the data structure itself seems like a good idea.

I'm starting to think that trying to eliminate the 1d indexing from interface of SparseGrid should be a non-goal for this class. It's not hard to wrap one's head around that there's linear indexing for a grid structure, and with the iteration sugar, the transfer algorithm doesn't need to bother with the indexing anyway (see https://github.com/xuchenhan-tri/drake/blob/a43facf6162e12da2dc8592cb12194c3ac24c29f/multibody/mpm/transfer.cc#L36-L97). OTOH, I'm fine with cutting SparseGrid::spgrid() and SparseGrid::mutable_sparse_grid() loose.

[amcastro-tri]

I'd imagine you are in paper crunch mode, but I just wanted to double check you are not waiting for me on this one?

Reviewable status: 8 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)

[xuchenhan-tri]

Sorry about the delay, but no, this PR is not waiting on you.