On Tue, Dec 12, 2023 at 7:28 PM MIGUEL MOLINOS PEREZ <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 <knep...@gmail.com> wrote: > > > On Tue, Dec 12, 2023 at 2:36 PM MIGUEL MOLINOS PEREZ <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, &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 <mmoli...@us.es> wrote: >> >> Thank you Mark! I will have a look to it. >> >> Best, >> Miguel >> >> >> On 4 Nov 2023, at 13:54, Matthew Knepley <knep...@gmail.com> wrote: >> >> >> On Sat, Nov 4, 2023 at 8:40 AM Mark Adams <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 <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/ >> <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/>