/* * vector_dist_performance_common.hpp * * Created on: Dec 25, 2016 * Author: i-bird */ #ifndef SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_ #define SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_ #include "Vector/vector_dist.hpp" /*! \brief Calculate and put particles' forces * * \param NN Cell list * \param vd Distributed vector * \param r_cut Cut-off radius */ template void calc_forces(T & NN, V & vd, float r_cut) { auto it_v = vd.getDomainIterator(); float sum[dim]; for (size_t i = 0; i < dim; i++) sum[i] = 0; while (it_v.isNext()) { //key vect_dist_key_dx key = it_v.get(); // Get the position of the particles Point p = vd.getPos(key); for (size_t i = 0; i < dim; i++) sum[i] = 0; // Get the neighborhood of the particle auto cell_it = NN.template getNNIterator(NN.getCell(p)); while(cell_it.isNext()) { auto nnp = cell_it.get(); // p != q if (nnp == key.getKey()) { ++cell_it; continue; } Point q = vd.getPos(nnp); if (p.distance2(q) < r_cut*r_cut) { //Calculate the forces float num[dim]; for (size_t i = 0; i < dim; i++) num[i] = vd.getPos(key)[i] - vd.getPos(nnp)[i]; float denom = 0; for (size_t i = 0; i < dim; i++) denom += num[i] * num[i]; float res[dim]; for (size_t i = 0; i < dim; i++) res[i] = num[i] / denom; for (size_t i = 0; i < dim; i++) sum[i] += res[i]; } //Next particle in a cell ++cell_it; } //Put the forces for (size_t i = 0; i < dim; i++) vd.template getProp(key)[i] += sum[i]; //Next particle in cell list ++it_v; } } /*! \brief For each particle of vd calculate the accumulation of the distances of the neighborhood * particles inside vd2 * * * \param NN Cell list vd * \param NN2 Cell list vd2 * \param vd Distributed vector * \param vd2 Distributed vector 2 * \param r_cut Cut-off radius * */ template void cross_calc(T & NN, T & NN2, V & vd, V & vd2) { auto it_v = vd.getDomainIterator(); while (it_v.isNext()) { //key vect_dist_key_dx key = it_v.get(); // Get the position of the particles Point p = vd.getPos(key); // Get the neighborhood of the particle auto cell_it = NN2.template getNNIterator(NN2.getCell(p)); double sum = 0.0; while(cell_it.isNext()) { auto nnp = cell_it.get(); Point q = vd2.getPos(nnp); sum += norm(p - q); //Next particle in a cell ++cell_it; } vd.template getProp(key) = sum; //Next particle in cell list ++it_v; } } /*! \brief Initialize a distributed vector * * \param vd Distributed vector */ template void vd_initialize(v_dist & vd, Vcluster & v_cl) { // The random generator engine std::default_random_engine eg(v_cl.getProcessUnitID()*4313); std::uniform_real_distribution ud(0.0f, 1.0f); //! [Create a vector of random elements on each processor 2D] auto it = vd.getIterator(); while (it.isNext()) { auto key = it.get(); for (size_t i = 0; i < dim; i++) {vd.getPos(key)[i] = ud(eg);} ++it; } vd.map(); } /*! \brief Initialize a distributed vector * * \param vd Distributed vector * \param box where to initialize * \param npart number of particles to add * */ template void vd_initialize_box_nomap(v_dist & vd, const Box & box, Vcluster & v_cl, int npart) { // The random generator engine std::default_random_engine eg(v_cl.getProcessUnitID()*4313); std::uniform_real_distribution ud(0.0f, 1.0f); size_t n = vd.size_local(); size_t n_stop = n + npart; auto it = vd.getIterator(n,n_stop); while (it.isNext()) { auto key = it.get(); for (size_t i = 0; i < dim; i++) {vd.getPos(key)[i] = (box.getHigh(i) - box.getLow(i))*ud(eg) + box.getLow(i);} ++it; } } /*! \brief Initialize 2 distributed vectors with equally positioned particles * * \param vd, vd2 Distributed vectors * \param v_cl Global vcluster * \param k_int Number of particles */ template void vd_initialize_double(v_dist & vd,v_dist & vd2, Vcluster<> & v_cl, size_t k_int) { // The random generator engine std::default_random_engine eg(v_cl.getProcessUnitID()*4313); std::uniform_real_distribution ud(0.0f, 1.0f); //! [Create a vector of random elements on each processor 2D] auto it = vd.getIterator(); while (it.isNext()) { auto key = it.get(); for (size_t i = 0; i < dim; i++) { vd.getPos(key)[i] = ud(eg); vd2.getPos(key)[i] = vd.getPos(key)[i]; } ++it; } vd.map(); vd2.map(); } /*! \brief Calculate and put particles' forces * * \param NN Cell list hilbert * \param vd Distributed vector * \param r_cut Cut-off radius */ template void calc_forces_hilb(T & NN, V & vd, float r_cut) { auto it_cl = NN.getIterator(); float sum[dim]; for (size_t i = 0; i < dim; i++) sum[i] = 0; while (it_cl.isNext()) { //key auto key = it_cl.get(); // Get the position of the particles Point p = vd.getPos(key); for (size_t i = 0; i < dim; i++) sum[i] = 0; // Get the neighborhood of the particle auto cell_it = NN.template getNNIterator(NN.getCell(p)); while(cell_it.isNext()) { auto nnp = cell_it.get(); // p != q if (nnp == key) { ++cell_it; continue; } Point q = vd.getPos(nnp); if (p.distance2(q) < r_cut*r_cut) { //Calculate the forces float num[dim]; for (size_t i = 0; i < dim; i++) num[i] = vd.getPos(key)[i] - vd.getPos(nnp)[i]; float denom = 0; for (size_t i = 0; i < dim; i++) denom += num[i] * num[i]; float res[dim]; for (size_t i = 0; i < dim; i++) res[i] = num[i] / denom; for (size_t i = 0; i < dim; i++) sum[i] += res[i]; } //Next particle in a cell ++cell_it; } //Put the forces for (size_t i = 0; i < dim; i++) vd.template getProp(key)[i] += sum[i]; //Next particle in cell list ++it_cl; } } #endif /* SRC_VECTOR_PERFORMANCE_VECTOR_DIST_PERFORMANCE_COMMON_HPP_ */