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

[MIGUEL MOLINOS PEREZ]

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  wrote:

On Tue, Dec 12, 2023 at 7:28 PM MIGUEL MOLINOS PEREZ 
mailto:mmoli...@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 
mailto:knep...@gmail.com>> wrote:


On Tue, Dec 12, 2023 at 2:36 PM MIGUEL MOLINOS PEREZ 
mailto:mmoli...@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, ));

PetscCall(DMSwarmGetLocalSize(Simulation->atomistic_data, _atoms_local));

//! Get array with the mean position of the atoms
DMSwarmGetField(Simulation->atomistic_data, DMSwarmPICField_coor, , 
NULL,
(void**)_q_ptr);
Eigen::Map 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, 
,
NULL, (void**)_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, ,
NULL, (void**)_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] = 

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

[Matthew Knepley]

On Tue, Dec 12, 2023 at 7:28 PM MIGUEL MOLINOS PEREZ  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  wrote:
>
> 
> On Tue, Dec 12, 2023 at 2:36 PM MIGUEL MOLINOS PEREZ 
> 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, ));
>>
>> PetscCall(DMSwarmGetLocalSize(Simulation->atomistic_data, _atoms_local
>> ));
>>
>> //! Get array with the mean position of the atoms
>> DMSwarmGetField(Simulation->atomistic_data, DMSwarmPICField_coor, &
>> blocksize, NULL,
>> (void**)_q_ptr);
>> Eigen::Map 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**)_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**)_q));
>>
>> cnt = 0;
>> for (PetscInt site_i = 0; site_i < n_atoms_local; site_i++) {
>> if (cnt < n_atoms) {
>> 

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

[MIGUEL MOLINOS PEREZ]

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

Best,
Miguel


On 4 Nov 2023, at 13:54, Matthew Knepley  wrote:


On Sat, Nov 4, 2023 at 8:40 AM Mark Adams 
mailto:mfad...@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 
mailto:mmoli...@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/

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

[Matthew Knepley]

On Sat, Nov 4, 2023 at 8:40 AM Mark Adams  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 
> 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/ 

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

[Mark Adams]

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  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.
>
> Best regards,
> Miguel

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

[MIGUEL MOLINOS PEREZ]

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.

Best regards,
Miguel