Simplified Hopkins bonded contact model with orientation

[akturner]

Summary

Updates the Hopkins bonded contact model, so that rather than integrating the explicit bond positions in time, they are inferred from the element orientation and an initial bond angle. The concomitant DEMSI branch is https://github.com/akturner/DEMSI/tree/simplified_bonds

Related Issue(s)

None

Author(s)

Adrian K. Turner (ALNL)

[karapeterson]

Code changes look good both in LAMMPS and in the associated DEMSI branch.

[dsbolin]

@karapeterson @akturner @snikolov3 @jtclemm

There are some issues here pertaining to how fix_neigh_history operates. I think it's worth a discussion here before I try to suggest any code changes.

When copying history entries, the intended behavior is as follows:

• By default, history entries for pair i-j are the negative of history entries for pair j-i. This is because history entries in most granular pair styles are relative displacements between i and j. See
  https://github.com/E3SM-Project/DEMSI-LAMMPS/blob/1af703ae065b1960c7ca2037b6ddab5cc823f5ec/src/KOKKOS/fix_neigh_history_kokkos.cpp#L177
• If this is not the case (as with several entries in pair_gran_hopkins), the pair flag nondefault_history_transfer can be set by a pair style, and a method transfer_history can be implemented in the pair style to indicate correspondance between i-j and j-i values for each history entry. See two lines above previous link.
• Looks like for fix_neigh_history_kokkos, two issues are introduced:
  1.) The default behavior is no longer to flip the negative sign, see:
  https://github.com/E3SM-Project/DEMSI-LAMMPS/blob/1af703ae065b1960c7ca2037b6ddab5cc823f5ec/src/KOKKOS/fix_neigh_history_kokkos.cpp#L179
  This would impact pair_gran_hooke_thickness_kokkos, and may be a problem for PR Add kokkos versions of contact model and fix. #33.
  2.) The nondefault_history_transfer leads to a hard-coded transfer of histories, currently only applicable to my original implementation of pair_gran_hopkins. See:
  @https://github.com/E3SM-Project/DEMSI-LAMMPS/blob/1af703ae065b1960c7ca2037b6ddab5cc823f5ec/src/KOKKOS/fix_neigh_history_kokkos.cpp#L147
  This can be modified to match Adrian's new implementation that uses orientation, but would not be generic beyond that. That's fine if we only ever want a single non-default history transfer contact model.

I guess my question is: can fix_neigh_history_kokkos not make calls to pair->transfer_history(), as is the case for the non-Kokkos version? If not, we could hard-code a few options for history transfer in fix_neigh_history_kokkos, and indicate the option via some additional member of the pair style.

[karapeterson]

@dsbolin Thanks for looking into the fix_neighbor_history implementation! I have never been certain why the kokkos version did not have the negative sign (issue 1). That error is also in the master LAMMPS branch. I just tried testing the change on the vortex test case, which uses pair_gran_hooke_thickness, but that test case doesn't seem to use fix_neighbor_history at all. Does that make sense to you?

Regarding issue 2, it's not possible to call nonmember functions from a Kokkos for loop or lambda. That is why the non-default history transfer from pair_gran_hopkins was implemented within fix_neighbor_history_kokkos. There are probably better ways to do this and would be happy to hear your thoughts on it.

[dsbolin]

@karapeterson Thanks for the quick reply. Re: 1 - you're right, I see it in the master lammps master now too, not sure why that's the case, I'll look into this and ask Joel as well. And you're also right that pair_gran_hooke_thickness would not be affected, I forgot that that pair style does not use history-dependent friction.

Re: 2: Thanks, that makes sense now. What about if we use the nondefault_history_transfer to indicate different history transfer options, e.g.:

if (nondefault_history_transfer == 0){ //Default behavior
}
else if (nondefault_history_transfer == 1){ //pair_gran_hopkins using orientation
 //Some hard-coded transfer behavior
}
else if ( nondefault_history_transfer == 2){ //Some other option
}
...

[jtclemm]

I also do not understand what's going on in the master lammps version of fix_neigh_history_kokkos.cpp. But just to throw it out there, could we also generalize the sign of history transfers by storing a public array in the relevant Pair classes of length size_history that stores +1 or -1 (or just a boolean) which could be copied by the FixNeighHistoryKokkos class at setup? Just spitballing ideas to avoid having to hard code options into fix_neigh_history_kokkos.cpp.

[snikolov3]

Yea I think storing a public array that fix_neigh_history_kokkos.cpp inherits should be straightforward. For example in the lines below (DEMSI-LAMMPS/src/KOKKOS/fix_neigh_history_kokkos.cpp:147-176:

      if (nondefault_history_transfer) {
        if (d_firstvalue(i,dnum*jj+8) < d_firstvalue(i,dnum*jj+9)) {
          d_valuepartner(j,dnum*m) = d_firstvalue(i,dnum*jj+4);
          d_valuepartner(j,dnum*m+1) = d_firstvalue(i,dnum*jj+5);
          d_valuepartner(j,dnum*m+2) = d_firstvalue(i,dnum*jj+6);
          d_valuepartner(j,dnum*m+3) = d_firstvalue(i,dnum*jj+7);

          d_valuepartner(j,dnum*m+4) = d_firstvalue(i,dnum*jj);
          d_valuepartner(j,dnum*m+5) = d_firstvalue(i,dnum*jj+1);
          d_valuepartner(j,dnum*m+6) = d_firstvalue(i,dnum*jj+2);
          d_valuepartner(j,dnum*m+7) = d_firstvalue(i,dnum*jj+3);

          d_valuepartner(j,dnum*m+8)  = d_firstvalue(i,dnum*jj+8);
          d_valuepartner(j,dnum*m+9)  = d_firstvalue(i,dnum*jj+9);
          d_valuepartner(j,dnum*m+10) = d_firstvalue(i,dnum*jj+10);
          d_valuepartner(j,dnum*m+11) = d_firstvalue(i,dnum*jj+11);
        } else {
          d_valuepartner(j,dnum*m)   = -d_firstvalue(i,dnum*jj);
          d_valuepartner(j,dnum*m+1) = -d_firstvalue(i,dnum*jj+1);
          d_valuepartner(j,dnum*m+2) = -d_firstvalue(i,dnum*jj+2);
          d_valuepartner(j,dnum*m+3) = -d_firstvalue(i,dnum*jj+3);

          d_valuepartner(j,dnum*m+4) = d_firstvalue(i,dnum*jj+4);
          d_valuepartner(j,dnum*m+5) = d_firstvalue(i,dnum*jj+5);
          d_valuepartner(j,dnum*m+8) = d_firstvalue(i,dnum*jj+8);
          d_valuepartner(j,dnum*m+9) = d_firstvalue(i,dnum*jj+9);
          d_valuepartner(j,dnum*m+10) = d_firstvalue(i,dnum*jj+10);
          d_valuepartner(j,dnum*m+11) = d_firstvalue(i,dnum*jj+11);
        }
      }

We would store a sign array (of length size_history) in fix_neigh_history.cpp which has different -1/1 values depending on the nondefault_history_transfer type (0,1,2... etc.). Is that what you had in mind? The only issue I see is that changes in the order of assignment in fix_neigh_history_kokkos.cpp could break things if these are not reflected in the indices of the sign array. Also I was wondering, in the else{} statement above, why are the signs between d_valuepartner(j,dnumm+3) and d_valuepartner(j,dnumm+4) different? Aren't these the x-y components of the s1 and s2 displacement vectors farthest away from the current element? The else statement is for the case when the bond is broken, because I think indices dnumjj+8 and dnumjj+9 correspond to chi1 and chi2?

[karapeterson]

Agree @jtclemm, @snikolov3 - having an array associated with a pair style that indicates the sign and index for the transfer may be a good solution.

[dsbolin]

in the else{} statement above, why are the signs between d_valuepartner(j,dnumm+3) and d_valuepartner(j,dnumm+4) different?

The else block here corresponds to non-bonded interactions. Entries 0, 1 are normal force components, entries 2 and 3 are tangential displacement components, so all of those get a sign flip if the order of the elements changes. Entry 4 is contact thickness, 5 is reference overlap, etc., so those don't get sign flips.
See comments here:
https://github.com/E3SM-Project/DEMSI-LAMMPS/blob/1af703ae065b1960c7ca2037b6ddab5cc823f5ec/src/USER-DEMSI/pair_gran_hopkins.cpp#L112-L126

[dsbolin]

Agree @jtclemm, @snikolov3 - having an array associated with a pair style that indicates the sign and index for the transfer may be a good solution.

In the original implementation, the if/else construction above was part of the pair style's transfer_history function, so that the indices/ordering are different depending on bonded/non-bonded interactions (which are also stored in the history entries, indices 8 and 9 in this case). So I don't think an array of indices would work in general.

Does kokkos allow function pointers? If so, could we have a function pointer that gets assigned to the pair's member function outside of the parallel for loop, then called within the for loop? If the history transfer functionality can't be part of the pair style, I'm not sure we can avoid hard-coding a different option for each 'nondefault_history_transfer' pair style. I don't think we'd have very many of these, so it wouldn't be too bad.

[karapeterson]

In the original implementation, the if/else construction above was part of the pair style's transfer_history function, so that the indices/ordering are different depending on bonded/non-bonded interactions (which are also stored in the history entries, indices 8 and 9 in this case). So I don't think an array of indices would work in general.

Good point @dsbolin.

I don't think we can use a function pointer - but will try to find out for sure.

[jtclemm]

In the original implementation, the if/else construction above was part of the pair style's transfer_history function, so that the indices/ordering are different depending on bonded/non-bonded interactions (which are also stored in the history entries, indices 8 and 9 in this case). So I don't think an array of indices would work in general.

Ah, I failed to understand that.

If it ends up that we can't use a function pointer, could we also just expand the history array by 4 and move non-bonded data to the end? Then the first 12 values are always transferred with no sign flip and the last 4 are always inverted. Not the most efficient solution though.

[dsbolin]

In the original implementation, the if/else construction above was part of the pair style's transfer_history function, so that the indices/ordering are different depending on bonded/non-bonded interactions (which are also stored in the history entries, indices 8 and 9 in this case). So I don't think an array of indices would work in general.

Ah, I failed to understand that.

If it ends up that we can't use a function pointer, could we also just expand the history array by 4 and move non-bonded data to the end? Then the first 12 values are always transferred with no sign flip and the last 4 are always inverted. Not the most efficient solution though.

That could work, and we could require that only two types of history transfers exist, corresponding to sign flip/no sign flip. My original pair_gran_hopkins implementation had order switches as well (history[0] = history[4]), but that can probably be avoided on the pair style side (e.g. Adrian's new pair_gran_hopkins implementation in this PR checks for atom global id order in the pair style, so that the order wouldn't need to be switched in fix_neigh_history). But, like you said, not the most efficient or general solution - though I'm not sure how much cost is incurred for having superfluous history entries.