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

Tue Dec 12 18:43:29 CST 2023

On Tue, Dec 12, 2023 at 7:28 PM MIGUEL MOLINOS PEREZ <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> wrote:
>
> 
> 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));
>>
>>
>>
>> <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/>
>
>

-- 
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/52990b36/attachment.html>