[petsc-users] Domain decomposition in PETSc for Molecular Dynamics

Wed Dec 13 04:18:50 CST 2023

Thank you for the feedback. Seems like the second option is the right one (according to MD literature).

Best,
Miguel

On 13 Dec 2023, at 01:43, Matthew Knepley <knepley at gmail.com> wrote:

On Tue, Dec 12, 2023 at 7:28 PM MIGUEL MOLINOS PEREZ <mmolinos at us.es<mailto:mmolinos at us.es>> wrote:
I meant the list of atoms which lies inside of a sphere of radius R_cutoff centered at the mean position of a given atom.

Okay, this is possible in parallel, but would require hard work and I have not done it yet, although I think all the tools are coded.

In serial, there are at least two ways to do it. First, you can use a k-d tree implementation, since they usually have the radius query. I have not put one of these in, because I did not like any implementation, but Underworld and Firedrake have and it works fine. Second, you can choose a grid size of R_cutoff for the background grid, and then check neighbors. This is probably how I would start.

  Thanks,

     Matt

Best,
Miguel

On 13 Dec 2023, at 01:14, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Tue, Dec 12, 2023 at 2:36 PM MIGUEL MOLINOS PEREZ <mmolinos at us.es<mailto:mmolinos at us.es>> wrote:
Dear Matthew and Mark,

Thank you for four useful guidance.  I have taken as a starting point the example in "dm/tutorials/swarm_ex3.c" to build a first approximation for domain decomposition in my molecular dynamics code (diffusive molecular dynamic to be more precise :-) ). And I must say that I am very happy with the result. However, in my journey integrating domain decomposition into my code, I am facing some new (and expected) problems.  The first is in the implementation of the nearest neighbor algorithm (list of atoms closest to each atom).

Can you help me understand this? For a given atom, there should be a single "closest" atom (barring
degeneracies in distance). What do you mean by the list of closest atoms?

  Thanks,

     Matt

My current approach to the problem is a brute force algorithm (double loop through the list of atoms and calculate the distance). However, it seems that if I call the "neighbours" function after the "DMSwarmMigrate" function the search algorithm does not work correctly. My thoughts / hints are:

  *   The two nested for loops should be done on the global indexing of the atoms instead of the local one (I don't know how to get this number).
  *   If I print the mean position of atom #0 (for example) each range prints a different value of the average position. One of them is the correct position corresponding to site #0, the others are different (but identically labeled) atomic sites. Which means that the site_i index is not bijective.

I believe that solving this problem will increase my understanding of the domain decomposition approach and may allow me to fix the remaining parts of my code.

Any additional comments are greatly appreciated. For instance, I will be happy to be pointed to any piece of code (petsc examples for example) with solves a similar problem in order to self-learn learn by example.

Many thanks in advance.

Best,
Miguel

This is the piece of code (simplified) which computes the list of neighbours for each atomic site. DMD is a structure which contains the atomistic information (DMSWARM), and the background mesh and bounding cell (DMDA and DMShell)

int neighbours(DMD* Simulation) {

PetscFunctionBegin;
PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

PetscCall(DMSwarmGetLocalSize(Simulation->atomistic_data, &n_atoms_local));

//! Get array with the mean position of the atoms
DMSwarmGetField(Simulation->atomistic_data, DMSwarmPICField_coor, &blocksize, NULL,
(void**)&mean_q_ptr);
Eigen::Map<MatrixType> mean_q(mean_q_ptr, n_atoms_local, dim);

int* neigh = Simulation->neigh;
int* numneigh = Simulation->numneigh;

for (unsigned int site_i = 0; site_i < n_atoms_local; site_i++) {

//! Get mean position of site i
Eigen::Vector3d mean_q_i = mean_q.block<1, 3>(site_i, 0);

//! Search neighbourhs in the main cell (+ periodic cells)
for (unsigned site_j = 0; site_j < n_atoms_local; site_j++) {
if (site_i != site_j) {
//! Get mean position of site j in the periodic box
Eigen::Vector3d mean_q_j = mean_q.block<1, 3>(site_j, 0);

//! Check is site j is the neibourhood of the site i
double norm_r_ij = (mean_q_i - mean_q_j).norm();
if ((norm_r_ij <= r_cutoff_ADP) && (numneigh[site_i] < maxneigh)) {
neigh[site_i * maxneigh + numneigh[site_i]] = site_j;
numneigh[site_i] += 1;
}
}
}

} // MPI for loop (site_i)

DMSwarmRestoreField(Simulation->atomistic_data, DMSwarmPICField_coor, &blocksize,
NULL, (void**)&mean_q_ptr);

return EXIT_SUCCESS;
}

This is the piece of code that I use to read the atomic positions (mean_q) from a file:
//! @brief mean_q: Mean value of each atomic position
double* mean_q;
PetscCall(DMSwarmGetField(atomistic_data, DMSwarmPICField_coor, &blocksize,
NULL, (void**)&mean_q));

cnt = 0;
for (PetscInt site_i = 0; site_i < n_atoms_local; site_i++) {
if (cnt < n_atoms) {
mean_q[blocksize * cnt + 0] = Simulation_file.mean_q[cnt * dim + 0];
mean_q[blocksize * cnt + 1] = Simulation_file.mean_q[cnt * dim + 1];
mean_q[blocksize * cnt + 2] = Simulation_file.mean_q[cnt * dim + 2];

cnt++;
}
}
PetscCall(DMSwarmRestoreField(atomistic_data, DMSwarmPICField_coor,
&blocksize, NULL, (void**)&mean_q));

<Screenshot 2023-12-12 at 19.42.13.png>

On 4 Nov 2023, at 15:50, MIGUEL MOLINOS PEREZ <mmolinos at us.es<mailto:mmolinos at us.es>> wrote:

Thank you Mark! I will have a look to it.

Best,
Miguel

On 4 Nov 2023, at 13:54, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Sat, Nov 4, 2023 at 8:40 AM Mark Adams <mfadams at lbl.gov<mailto:mfadams at lbl.gov>> wrote:
Hi MIGUEL,

This might be a good place to start: https://petsc.org/main/manual/vec/
Feel free to ask more specific questions, but the docs are a good place to start.

Thanks,
Mark

On Fri, Nov 3, 2023 at 5:19 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es<mailto:mmolinos at us.es>> wrote:
Dear all,

I am currently working on the development of a in-house molecular dynamics code using PETSc and C++. So far the code works great, however it is a little bit slow since I am not exploiting MPI for PETSc vectors. I was wondering if there is a way to perform the domain decomposition efficiently using some PETSc functionality. Any feedback is highly appreciated.

It sounds like you mean "is there a way to specify a communication construct that can send my particle
information automatically". We use PetscSF for that. You can see how this works with the DMSwarm class, which represents a particle discretization. You can either use that, or if it does not work for you, do the same things with your class.

  Thanks,

     Matt

Best regards,
Miguel

--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>

--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>

--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231213/b49a4bef/attachment-0001.html>