<div dir="ltr"><div dir="ltr">On Tue, Dec 12, 2023 at 7:28 PM MIGUEL MOLINOS PEREZ <<a href="mailto:mmolinos@us.es">mmolinos@us.es</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">



<div dir="auto">
<div dir="ltr">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.</div></div></blockquote><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto">
<div dir="ltr">Best,</div>
<div dir="ltr">Miguel</div>
<div dir="ltr"><br>
<blockquote type="cite">On 13 Dec 2023, at 01:14, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
<br>
</blockquote>
</div>
<blockquote type="cite">
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">On Tue, Dec 12, 2023 at 2:36 PM MIGUEL MOLINOS PEREZ <<a href="mailto:mmolinos@us.es" target="_blank">mmolinos@us.es</a>> wrote:<br>
</div>
<div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>Dear Matthew and Mark,
<div><br>
</div>
<div>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). </div>
</div>
</blockquote>
<div><br>
</div>
<div>Can you help me understand this? For a given atom, there should be a single "closest" atom (barring</div>
<div>degeneracies in distance). What do you mean by the list of closest atoms?</div>
<div><br>
</div>
<div>  Thanks,</div>
<div><br>
</div>
<div>     Matt</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<div>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:</div>
<div>
<ul>
<li>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).  </li><li>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.</li></ul>
</div>
<div><br>
</div>
<div>
<div>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.</div>
</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>Many thanks in advance.</div>
<div><br>
</div>
<div>
<div style="color:rgb(0,0,0)">Best,</div>
<div style="color:rgb(0,0,0)">Miguel</div>
</div>
<div><font color="#000000"><span style="white-space:pre-wrap"><br>
</span></font></div>
<div><font color="#000000"><span style="white-space:pre-wrap">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</span>
<span style="white-space:pre-wrap">(DMSWARM), and the background mesh and bounding cell (DMDA and DMShell)</span></font></div>
<div><br>
</div>
<div>
<div style="color:rgb(54,54,54);background-color:rgb(255,255,255);font-family:Menlo,Monaco,"Courier New",monospace;line-height:18px;white-space:pre-wrap">
<div><span style="color:rgb(63,151,223)">int</span> <span style="color:rgb(99,99,36)">
neighbours</span>(<span style="color:rgb(70,224,192)">DMD</span><span style="color:rgb(63,151,223)">*</span>
<span style="color:rgb(9,89,132)">Simulation</span>) {</div>
<br>
<div><span style="color:rgb(63,151,223)">PetscFunctionBegin</span>;</div>
<div><span style="color:rgb(63,151,223)">PetscCallMPI</span>(<span style="color:rgb(99,99,36)">MPI_Comm_rank</span>(PETSC_COMM_WORLD, &<span style="color:rgb(9,89,132)">rank</span>));</div>
<br>
<div><span style="color:rgb(63,151,223)">PetscCall</span>(<span style="color:rgb(99,99,36)">DMSwarmGetLocalSize</span>(<span style="color:rgb(9,89,132)">Simulation</span>-><span style="color:rgb(9,89,132)">atomistic_data</span>, &<span style="color:rgb(9,89,132)">n_atoms_local</span>));</div>
<br>
<div><span style="color:rgb(146,205,120)">//! Get array with the mean position of the atoms</span></div>
<div><span style="color:rgb(99,99,36)">DMSwarmGetField</span>(<span style="color:rgb(9,89,132)">Simulation</span>-><span style="color:rgb(9,89,132)">atomistic_data</span>,
<span style="color:rgb(15,157,235)">DMSwarmPICField_coor</span>, &<span style="color:rgb(9,89,132)">blocksize</span>,
<span style="color:rgb(63,151,223)">NULL</span>,</div>
<div>(<span style="color:rgb(63,151,223)">void</span>**)&<span style="color:rgb(9,89,132)">mean_q_ptr</span>);</div>
<div><span style="color:rgb(70,224,192)">Eigen</span>::Map<MatrixType> <span style="color:rgb(99,99,36)">
mean_q</span>(mean_q_ptr, n_atoms_local, dim);</div>
<br>
<div><span style="color:rgb(63,151,223)">int</span>* <span style="color:rgb(9,89,132)">
neigh</span> = <span style="color:rgb(9,89,132)">Simulation</span>-><span style="color:rgb(9,89,132)">neigh</span>;</div>
<div><span style="color:rgb(63,151,223)">int</span>* <span style="color:rgb(9,89,132)">
numneigh</span> = <span style="color:rgb(9,89,132)">Simulation</span>-><span style="color:rgb(9,89,132)">numneigh</span>;</div>
<div><br>
</div>
<div><span style="color:rgb(157,78,150)">for</span> (<span style="color:rgb(63,151,223)">unsigned</span>
<span style="color:rgb(63,151,223)">int</span> <span style="color:rgb(9,89,132)">
site_i</span> = <span style="color:rgb(73,104,57)">0</span>; <span style="color:rgb(9,89,132)">
site_i</span> < <span style="color:rgb(9,89,132)">n_atoms_local</span>; <span style="color:rgb(9,89,132)">
site_i</span>++) {</div>
<br>
<div><span style="color:rgb(146,205,120)">//! Get mean position of site i</span></div>
<div><span style="color:rgb(70,224,192)">Eigen</span>::Vector3d <span style="color:rgb(9,89,132)">
mean_q_i</span> = <span style="color:rgb(9,89,132)">mean_q</span>.<span style="color:rgb(9,89,132)">block</span><<span style="color:rgb(73,104,57)">1</span>,
<span style="color:rgb(73,104,57)">3</span>>(site_i, <span style="color:rgb(73,104,57)">
0</span>);</div>
<br>
<div><span style="color:rgb(146,205,120)">//! Search neighbourhs in the main cell (+ periodic cells)</span></div>
<div><span style="color:rgb(157,78,150)">for</span> (<span style="color:rgb(63,151,223)">unsigned</span>
<span style="color:rgb(9,89,132)">site_j</span> = <span style="color:rgb(73,104,57)">
0</span>; <span style="color:rgb(9,89,132)">site_j</span> < <span style="color:rgb(9,89,132)">
n_atoms_local</span>; <span style="color:rgb(9,89,132)">site_j</span>++) {</div>
<div><span style="white-space:pre-wrap"></span><span style="color:rgb(157,78,150)">if</span> (<span style="color:rgb(9,89,132)">site_i</span> !=
<span style="color:rgb(9,89,132)">site_j</span>) {</div>
<div><span style="color:rgb(146,205,120)">//! Get mean position of site j in the periodic box</span></div>
<div><span style="white-space:pre-wrap"></span><span style="color:rgb(70,224,192)">Eigen</span>::Vector3d mean_q_j =
<span style="color:rgb(9,89,132)">mean_q</span>.<span style="color:rgb(9,89,132)">block</span><<span style="color:rgb(73,104,57)">1</span>,
<span style="color:rgb(73,104,57)">3</span>>(site_j, <span style="color:rgb(73,104,57)">
0</span>);</div>
<br>
<div><span style="color:rgb(146,205,120)"><span style="white-space:pre-wrap"></span>//! Check is site j is the neibourhood of the site i</span></div>
<div><span style="white-space:pre-wrap"></span><span style="color:rgb(63,151,223)">double</span>
<span style="color:rgb(9,89,132)">norm_r_ij</span> = (mean_q_i - mean_q_j).<span style="color:rgb(99,99,36)">norm</span>();</div>
<div><span style="white-space:pre-wrap"></span><span style="color:rgb(157,78,150)">if</span> ((<span style="color:rgb(9,89,132)">norm_r_ij</span> <=
<span style="color:rgb(63,151,223)">r_cutoff_ADP</span>) && (<span style="color:rgb(9,89,132)">numneigh</span>[<span style="color:rgb(9,89,132)">site_i</span>] <
<span style="color:rgb(63,151,223)">maxneigh</span>)) {</div>
<div><span style="white-space:pre-wrap"></span><span style="color:rgb(9,89,132)">neigh</span>[<span style="color:rgb(9,89,132)">site_i</span> *
<span style="color:rgb(63,151,223)">maxneigh</span> + <span style="color:rgb(9,89,132)">
numneigh</span>[<span style="color:rgb(9,89,132)">site_i</span>]] = <span style="color:rgb(9,89,132)">
site_j</span>;</div>
<div><span style="white-space:pre-wrap"></span><span style="color:rgb(9,89,132)">numneigh</span>[<span style="color:rgb(9,89,132)">site_i</span>] +=
<span style="color:rgb(73,104,57)">1</span>;</div>
<div><span style="white-space:pre-wrap"></span>} </div>
<div>}</div>
<div>}</div>
<br>
<div>}<span style="color:rgb(146,205,120)"> // MPI for loop (site_i)</span></div>
<br>
<div><span style="color:rgb(99,99,36)">DMSwarmRestoreField</span>(<span style="color:rgb(9,89,132)">Simulation</span>-><span style="color:rgb(9,89,132)">atomistic_data</span>,
<span style="color:rgb(15,157,235)">DMSwarmPICField_coor</span>, &<span style="color:rgb(9,89,132)">blocksize</span>,</div>
<div><span style="color:rgb(63,151,223)">NULL</span>, (<span style="color:rgb(63,151,223)">void</span>**)&<span style="color:rgb(9,89,132)">mean_q_ptr</span>);</div>
<br>
<div><span style="color:rgb(157,78,150)">return</span> <span style="color:rgb(63,151,223)">
EXIT_SUCCESS</span>;</div>
<div>}</div>
</div>
</div>
<div><br>
</div>
<div><br>
</div>
<div>This is the piece of code that I use to read the atomic positions (mean_q) from a file:</div>
<div>
<div style="color:rgb(54,54,54);background-color:rgb(255,255,255);font-family:Menlo,Monaco,"Courier New",monospace;line-height:18px;white-space:pre-wrap">
<div><span style="color:rgb(146,205,120)">//! </span><span style="color:rgb(63,151,223)">@brief</span><span style="color:rgb(146,205,120)"> mean_q: Mean value of each atomic position</span></div>
<div><span style="color:rgb(63,151,223)">double</span>* <span style="color:rgb(9,89,132)">
mean_q</span>;</div>
<div><span style="color:rgb(63,151,223)">PetscCall</span>(<span style="color:rgb(99,99,36)">DMSwarmGetField</span>(<span style="color:rgb(9,89,132)">atomistic_data</span>,
<span style="color:rgb(15,157,235)">DMSwarmPICField_coor</span>, &<span style="color:rgb(9,89,132)">blocksize</span>,</div>
<div><span style="color:rgb(63,151,223)">NULL</span>, (<span style="color:rgb(63,151,223)">void</span>**)&<span style="color:rgb(9,89,132)">mean_q</span>));</div>
<br>
<div><span style="color:rgb(9,89,132)">cnt</span> = <span style="color:rgb(73,104,57)">
0</span>;</div>
<div><span style="color:rgb(157,78,150)">for</span> (<span style="color:rgb(70,224,192)">PetscInt</span>
<span style="color:rgb(9,89,132)">site_i</span> = <span style="color:rgb(73,104,57)">
0</span>; <span style="color:rgb(9,89,132)">site_i</span> < <span style="color:rgb(9,89,132)">
n_atoms_local</span>; <span style="color:rgb(9,89,132)">site_i</span>++) {</div>
<div><span style="color:rgb(157,78,150)">if</span> (<span style="color:rgb(9,89,132)">cnt</span> <
<span style="color:rgb(9,89,132)">n_atoms</span>) {</div>
<div><span style="color:rgb(9,89,132)">mean_q</span>[<span style="color:rgb(9,89,132)">blocksize</span> *
<span style="color:rgb(9,89,132)">cnt</span> + <span style="color:rgb(73,104,57)">
0</span>] = <span style="color:rgb(9,89,132)">Simulation_file</span>.<span style="color:rgb(9,89,132)">mean_q</span>[<span style="color:rgb(9,89,132)">cnt</span> *
<span style="color:rgb(9,89,132)">dim</span> + <span style="color:rgb(73,104,57)">
0</span>];</div>
<div><span style="color:rgb(9,89,132)">mean_q</span>[<span style="color:rgb(9,89,132)">blocksize</span> *
<span style="color:rgb(9,89,132)">cnt</span> + <span style="color:rgb(73,104,57)">
1</span>] = <span style="color:rgb(9,89,132)">Simulation_file</span>.<span style="color:rgb(9,89,132)">mean_q</span>[<span style="color:rgb(9,89,132)">cnt</span> *
<span style="color:rgb(9,89,132)">dim</span> + <span style="color:rgb(73,104,57)">
1</span>];</div>
<div><span style="color:rgb(9,89,132)">mean_q</span>[<span style="color:rgb(9,89,132)">blocksize</span> *
<span style="color:rgb(9,89,132)">cnt</span> + <span style="color:rgb(73,104,57)">
2</span>] = <span style="color:rgb(9,89,132)">Simulation_file</span>.<span style="color:rgb(9,89,132)">mean_q</span>[<span style="color:rgb(9,89,132)">cnt</span> *
<span style="color:rgb(9,89,132)">dim</span> + <span style="color:rgb(73,104,57)">
2</span>];</div>
<br>
<div><span style="color:rgb(9,89,132)">cnt</span>++;</div>
<div>}</div>
<div>}</div>
<div><span style="color:rgb(63,151,223)">PetscCall</span>(<span style="color:rgb(99,99,36)">DMSwarmRestoreField</span>(<span style="color:rgb(9,89,132)">atomistic_data</span>,
<span style="color:rgb(15,157,235)">DMSwarmPICField_coor</span>,</div>
<div>&<span style="color:rgb(9,89,132)">blocksize</span>, <span style="color:rgb(63,151,223)">
NULL</span>, (<span style="color:rgb(63,151,223)">void</span>**)&<span style="color:rgb(9,89,132)">mean_q</span>));</div>
</div>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br id="m_4736401739173453802m_7237655108684911485lineBreakAtBeginningOfMessage">
<div>
<div><Screenshot 2023-12-12 at 19.42.13.png></div>
<br>
<blockquote type="cite">
<div>On 4 Nov 2023, at 15:50, MIGUEL MOLINOS PEREZ <<a href="mailto:mmolinos@us.es" target="_blank">mmolinos@us.es</a>> wrote:</div>
<br>
<div>
<div dir="auto">
<div dir="ltr">Thank you Mark! I will have a look to it. <br>
<div dir="ltr">
<div dir="ltr">
<div><br>
</div>
<div>Best,</div>
<div>Miguel</div>
<div><br id="m_4736401739173453802m_7237655108684911485lineBreakAtBeginningOfSignature">
<div dir="ltr">
<blockquote type="cite"><br>
On 4 Nov 2023, at 13:54, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
<br>
</blockquote>
</div>
<blockquote type="cite">
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">On Sat, Nov 4, 2023 at 8:40 AM Mark Adams <<a href="mailto:mfadams@lbl.gov" target="_blank">mfadams@lbl.gov</a>> wrote:<br>
</div>
<div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr">Hi MIGUEL,
<div><br>
</div>
<div>This might be a good place to start: <a href="https://petsc.org/main/manual/vec/" target="_blank">https://petsc.org/main/manual/vec/</a><br>
</div>
<div>Feel free to ask more specific questions, but the docs are a good place to start.</div>
<div><br>
</div>
<div>Thanks,</div>
<div>Mark</div>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Fri, Nov 3, 2023 at 5:19 AM MIGUEL MOLINOS PEREZ <<a href="mailto:mmolinos@us.es" target="_blank">mmolinos@us.es</a>> wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Dear all, <br>
<br>
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.<br>
</blockquote>
</div>
</blockquote>
<div><br>
</div>
<div>It sounds like you mean "is there a way to specify a communication construct that can send my particle</div>
<div>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.</div>
<div><br>
</div>
<div>  Thanks,</div>
<div><br>
</div>
<div>     Matt</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Best regards,<br>
Miguel </blockquote>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
<span class="gmail_signature_prefix">-- </span><br>
<div dir="ltr" class="gmail_signature">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener</div>
<div><br>
</div>
<div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
<span class="gmail_signature_prefix">-- </span><br>
<div dir="ltr" class="gmail_signature">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener</div>
<div><br>
</div>
<div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>

</blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>