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

Tue Dec 12 18:01:13 CST 2023

On Tue, Dec 12, 2023 at 2:36 PM MIGUEL MOLINOS PEREZ <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));
>
>
>
> [image: Screenshot 2023-12-12 at 19.42.13.png]
>
> On 4 Nov 2023, at 15:50, MIGUEL MOLINOS PEREZ <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> wrote:
>
> 
> On Sat, Nov 4, 2023 at 8:40 AM Mark Adams <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>
>> 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231212/a5d04bbe/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot 2023-12-12 at 19.42.13.png
Type: image/png
Size: 1252024 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231212/a5d04bbe/attachment-0001.png>