| /*========================================================================= |
| |
| Copyright (c) 2007, Los Alamos National Security, LLC |
| |
| All rights reserved. |
| |
| Copyright 2007. Los Alamos National Security, LLC. |
| This software was produced under U.S. Government contract DE-AC52-06NA25396 |
| for Los Alamos National Laboratory (LANL), which is operated by |
| Los Alamos National Security, LLC for the U.S. Department of Energy. |
| The U.S. Government has rights to use, reproduce, and distribute this software. |
| NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, |
| EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. |
| If software is modified to produce derivative works, such modified software |
| should be clearly marked, so as not to confuse it with the version available |
| from LANL. |
| |
| Additionally, redistribution and use in source and binary forms, with or |
| without modification, are permitted provided that the following conditions |
| are met: |
| - Redistributions of source code must retain the above copyright notice, |
| this list of conditions and the following disclaimer. |
| - Redistributions in binary form must reproduce the above copyright notice, |
| this list of conditions and the following disclaimer in the documentation |
| and/or other materials provided with the distribution. |
| - Neither the name of Los Alamos National Security, LLC, Los Alamos National |
| Laboratory, LANL, the U.S. Government, nor the names of its contributors |
| may be used to endorse or promote products derived from this software |
| without specific prior written permission. |
| |
| THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS |
| "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, |
| THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL SECURITY, LLC OR |
| CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; |
| OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, |
| WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR |
| OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF |
| ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| |
| =========================================================================*/ |
| |
| #include <iostream> |
| #include <fstream> |
| #include <sstream> |
| #include <iomanip> |
| |
| #include <sys/types.h> |
| #include <dirent.h> |
| |
| #include "Partition.h" |
| #include "InitialExchange.h" |
| |
| using namespace std; |
| |
| ///////////////////////////////////////////////////////////////////////// |
| // |
| // InitialExchange takes input from the particle initializer which originally |
| // contained alive particles on this processor, but after the particles move |
| // a little they might have moved to a neighbor processor. Input is in the |
| // form of arrays, but vectors are supplied which will be filled with the |
| // alive particles for this processor. As each particle is examined it is |
| // placed immediately in the vector if it is alive on this processor and the |
| // index is placed by neighbor if it is dead on this processor. After |
| // categorizing all particles, the dead are exchanged with neighbors where |
| // they become alive. |
| // |
| // After this initial exchange, ParticleExchange is called to place all |
| // dead particles (which will be the ones sent in this step being returned |
| // plus the others which were originally placed correctly on the neighbor). |
| // |
| ///////////////////////////////////////////////////////////////////////// |
| |
| InitialExchange::InitialExchange() |
| { |
| // Get the number of processors running this problem and rank |
| this->numProc = Partition::getNumProc(); |
| this->myProc = Partition::getMyProc(); |
| |
| // Get the number of processors in each dimension |
| Partition::getDecompSize(this->layoutSize); |
| |
| // Get my position within the Cartesian topology |
| Partition::getMyPosition(this->layoutPos); |
| |
| // Get neighbors of this processor including the wraparound |
| Partition::getNeighbors(this->neighbor); |
| |
| this->numberOfAliveParticles = 0; |
| } |
| |
| InitialExchange::~InitialExchange() |
| { |
| } |
| |
| ///////////////////////////////////////////////////////////////////////// |
| // |
| // Set parameters for particle distribution |
| // |
| ///////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::setParameters(POSVEL_T rL, POSVEL_T deadSz) |
| { |
| // Physical total space and amount of physical space to use for dead particles |
| this->boxSize = rL; |
| this->deadSize = deadSz; |
| |
| #ifdef DEBUG |
| if (this->myProc == MASTER) { |
| cout << endl << "------------------------------------" << endl; |
| cout << "boxSize: " << this->boxSize << endl; |
| cout << "deltaBox: " << this->deadSize << endl; |
| } |
| #endif |
| } |
| |
| ///////////////////////////////////////////////////////////////////////// |
| // |
| // Initialize the regions outside the alive area for this processor which |
| // contain the particles which must be sent away to neighbors |
| // |
| ///////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::initialize() |
| { |
| #ifdef DEBUG |
| if (this->myProc == MASTER) |
| cout << "Decomposition: [" << this->layoutSize[0] << ":" |
| << this->layoutSize[1] << ":" << this->layoutSize[2] << "]" << endl; |
| #endif |
| |
| // Set subextents on particle locations for this processor |
| POSVEL_T boxStep[DIMENSION]; |
| for (int dim = 0; dim < DIMENSION; dim++) { |
| boxStep[dim] = this->boxSize / this->layoutSize[dim]; |
| |
| // Particles in this region belong to this processor as alive |
| this->minAlive[dim] = this->layoutPos[dim] * boxStep[dim]; |
| this->maxAlive[dim] = this->minAlive[dim] + boxStep[dim]; |
| if (this->maxAlive[dim] > this->boxSize) |
| this->maxAlive[dim] = this->boxSize; |
| |
| // Particles in this region are dead on this processor but alive elsewhere |
| this->minDead[dim] = this->minAlive[dim] - this->deadSize; |
| this->maxDead[dim] = this->maxAlive[dim] + this->deadSize; |
| } |
| |
| // Set the ranges on the dead particles for each neighbor direction |
| calculateExchangeRegions(); |
| } |
| |
| ///////////////////////////////////////////////////////////////////////// |
| // |
| // Each of the 26 neighbors will be sent a rectangular region of particles |
| // on this processor where they are dead. The locations coming from the |
| // initializer's arrays are in the range [0:rL] and so particles to the |
| // left of the leftmost processor will have values at the far end of |
| // the box. When they are sent they won't need to be modified by offsets. |
| // |
| ///////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::calculateExchangeRegions() |
| { |
| // Initialize all neighbors to the entire available exchange range |
| for (int i = 0; i < NUM_OF_NEIGHBORS; i++) { |
| for (int dim = 0; dim < DIMENSION; dim++) { |
| this->minRange[i][dim] = this->minAlive[dim]; |
| this->maxRange[i][dim] = this->maxAlive[dim]; |
| } |
| } |
| |
| // Left face |
| this->minRange[X0][0] = this->minDead[0]; |
| this->maxRange[X0][0] = this->minAlive[0]; |
| |
| // Right face |
| this->minRange[X1][0] = this->maxAlive[0]; |
| this->maxRange[X1][0] = this->maxDead[0]; |
| |
| // Bottom face |
| this->minRange[Y0][1] = this->minDead[1]; |
| this->maxRange[Y0][1] = this->minAlive[1]; |
| |
| // Top face |
| this->minRange[Y1][1] = this->maxAlive[1]; |
| this->maxRange[Y1][1] = this->maxDead[1]; |
| |
| // Front face |
| this->minRange[Z0][2] = this->minDead[2]; |
| this->maxRange[Z0][2] = this->minAlive[2]; |
| |
| // Back face |
| this->minRange[Z1][2] = this->maxAlive[2]; |
| this->maxRange[Z1][2] = this->maxDead[2]; |
| |
| // Left bottom and top bars |
| this->minRange[X0_Y0][0] = this->minDead[0]; |
| this->maxRange[X0_Y0][0] = this->minAlive[0]; |
| this->minRange[X0_Y0][1] = this->minDead[1]; |
| this->maxRange[X0_Y0][1] = this->minAlive[1]; |
| |
| this->minRange[X0_Y1][0] = this->minDead[0]; |
| this->maxRange[X0_Y1][0] = this->minAlive[0]; |
| this->minRange[X0_Y1][1] = this->maxAlive[1]; |
| this->maxRange[X0_Y1][1] = this->maxDead[1]; |
| |
| // Right bottom and top bars |
| this->minRange[X1_Y0][0] = this->maxAlive[0]; |
| this->maxRange[X1_Y0][0] = this->maxDead[0]; |
| this->minRange[X1_Y0][1] = this->minDead[1]; |
| this->maxRange[X1_Y0][1] = this->minAlive[1]; |
| |
| this->minRange[X1_Y1][0] = this->maxAlive[0]; |
| this->maxRange[X1_Y1][0] = this->maxDead[0]; |
| this->minRange[X1_Y1][1] = this->maxAlive[1]; |
| this->maxRange[X1_Y1][1] = this->maxDead[1]; |
| |
| // Bottom front and back bars |
| this->minRange[Y0_Z0][1] = this->minDead[1]; |
| this->maxRange[Y0_Z0][1] = this->minAlive[1]; |
| this->minRange[Y0_Z0][2] = this->minDead[2]; |
| this->maxRange[Y0_Z0][2] = this->minAlive[2]; |
| |
| this->minRange[Y0_Z1][1] = this->minDead[1]; |
| this->maxRange[Y0_Z1][1] = this->minAlive[1]; |
| this->minRange[Y0_Z1][2] = this->maxAlive[2]; |
| this->maxRange[Y0_Z1][2] = this->maxDead[2]; |
| |
| // Top front and back bars |
| this->minRange[Y1_Z0][1] = this->maxAlive[1]; |
| this->maxRange[Y1_Z0][1] = this->maxDead[1]; |
| this->minRange[Y1_Z0][2] = this->minDead[2]; |
| this->maxRange[Y1_Z0][2] = this->minAlive[2]; |
| |
| this->minRange[Y1_Z1][1] = this->maxAlive[1]; |
| this->maxRange[Y1_Z1][1] = this->maxDead[1]; |
| this->minRange[Y1_Z1][2] = this->maxAlive[2]; |
| this->maxRange[Y1_Z1][2] = this->maxDead[2]; |
| |
| // Left front and back bars (vertical) |
| this->minRange[Z0_X0][0] = this->minDead[0]; |
| this->maxRange[Z0_X0][0] = this->minAlive[0]; |
| this->minRange[Z0_X0][2] = this->minDead[2]; |
| this->maxRange[Z0_X0][2] = this->minAlive[2]; |
| |
| this->minRange[Z1_X0][0] = this->minDead[0]; |
| this->maxRange[Z1_X0][0] = this->minAlive[0]; |
| this->minRange[Z1_X0][2] = this->maxAlive[2]; |
| this->maxRange[Z1_X0][2] = this->maxDead[2]; |
| |
| // Right front and back bars (vertical) |
| this->minRange[Z0_X1][0] = this->maxAlive[0]; |
| this->maxRange[Z0_X1][0] = this->maxDead[0]; |
| this->minRange[Z0_X1][2] = this->minDead[2]; |
| this->maxRange[Z0_X1][2] = this->minAlive[2]; |
| |
| this->minRange[Z1_X1][0] = this->maxAlive[0]; |
| this->maxRange[Z1_X1][0] = this->maxDead[0]; |
| this->minRange[Z1_X1][2] = this->maxAlive[2]; |
| this->maxRange[Z1_X1][2] = this->maxDead[2]; |
| |
| // Left bottom front corner |
| this->minRange[X0_Y0_Z0][0] = this->minDead[0]; |
| this->maxRange[X0_Y0_Z0][0] = this->minAlive[0]; |
| this->minRange[X0_Y0_Z0][1] = this->minDead[1]; |
| this->maxRange[X0_Y0_Z0][1] = this->minAlive[1]; |
| this->minRange[X0_Y0_Z0][2] = this->minDead[2]; |
| this->maxRange[X0_Y0_Z0][2] = this->minAlive[2]; |
| |
| // Left bottom back corner |
| this->minRange[X0_Y0_Z1][0] = this->minDead[0]; |
| this->maxRange[X0_Y0_Z1][0] = this->minAlive[0]; |
| this->minRange[X0_Y0_Z1][1] = this->minDead[1]; |
| this->maxRange[X0_Y0_Z1][1] = this->minAlive[1]; |
| this->minRange[X0_Y0_Z1][2] = this->maxAlive[2]; |
| this->maxRange[X0_Y0_Z1][2] = this->maxDead[2]; |
| |
| // Left top front corner |
| this->minRange[X0_Y1_Z0][0] = this->minDead[0]; |
| this->maxRange[X0_Y1_Z0][0] = this->minAlive[0]; |
| this->minRange[X0_Y1_Z0][1] = this->maxAlive[1]; |
| this->maxRange[X0_Y1_Z0][1] = this->maxDead[1]; |
| this->minRange[X0_Y1_Z0][2] = this->minDead[2]; |
| this->maxRange[X0_Y1_Z0][2] = this->minAlive[2]; |
| |
| // Left top back corner |
| this->minRange[X0_Y1_Z1][0] = this->minDead[0]; |
| this->maxRange[X0_Y1_Z1][0] = this->minAlive[0]; |
| this->minRange[X0_Y1_Z1][1] = this->maxAlive[1]; |
| this->maxRange[X0_Y1_Z1][1] = this->maxDead[1]; |
| this->minRange[X0_Y1_Z1][2] = this->maxAlive[2]; |
| this->maxRange[X0_Y1_Z1][2] = this->maxDead[2]; |
| |
| // Right bottom front corner |
| this->minRange[X1_Y0_Z0][0] = this->maxAlive[0]; |
| this->maxRange[X1_Y0_Z0][0] = this->maxDead[0]; |
| this->minRange[X1_Y0_Z0][1] = this->minDead[1]; |
| this->maxRange[X1_Y0_Z0][1] = this->minAlive[1]; |
| this->minRange[X1_Y0_Z0][2] = this->minDead[2]; |
| this->maxRange[X1_Y0_Z0][2] = this->minAlive[2]; |
| |
| // Right bottom back corner |
| this->minRange[X1_Y0_Z1][0] = this->maxAlive[0]; |
| this->maxRange[X1_Y0_Z1][0] = this->maxDead[0]; |
| this->minRange[X1_Y0_Z1][1] = this->minDead[1]; |
| this->maxRange[X1_Y0_Z1][1] = this->minAlive[1]; |
| this->minRange[X1_Y0_Z1][2] = this->maxAlive[2]; |
| this->maxRange[X1_Y0_Z1][2] = this->maxDead[2]; |
| |
| // Right top front corner |
| this->minRange[X1_Y1_Z0][0] = this->maxAlive[0]; |
| this->maxRange[X1_Y1_Z0][0] = this->maxDead[0]; |
| this->minRange[X1_Y1_Z0][1] = this->maxAlive[1]; |
| this->maxRange[X1_Y1_Z0][1] = this->maxDead[1]; |
| this->minRange[X1_Y1_Z0][2] = this->minDead[2]; |
| this->maxRange[X1_Y1_Z0][2] = this->minAlive[2]; |
| |
| // Right top back corner |
| this->minRange[X1_Y1_Z1][0] = this->maxAlive[0]; |
| this->maxRange[X1_Y1_Z1][0] = this->maxDead[0]; |
| this->minRange[X1_Y1_Z1][1] = this->maxAlive[1]; |
| this->maxRange[X1_Y1_Z1][1] = this->maxDead[1]; |
| this->minRange[X1_Y1_Z1][2] = this->maxAlive[2]; |
| this->maxRange[X1_Y1_Z1][2] = this->maxDead[2]; |
| |
| // Fix ranges for processors on a face in the decomposition |
| // Processor is on front edge in X dimension |
| if (this->layoutPos[0] == 0) { |
| this->minRange[X0][0] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0][0] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y1][0] = this->boxSize - this->deadSize; |
| this->minRange[Z0_X0][0] = this->boxSize - this->deadSize; |
| this->minRange[Z1_X0][0] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0_Z0][0] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0_Z1][0] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y1_Z0][0] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y1_Z1][0] = this->boxSize - this->deadSize; |
| |
| this->maxRange[X0][0] = this->boxSize; |
| this->maxRange[X0_Y0][0] = this->boxSize; |
| this->maxRange[X0_Y1][0] = this->boxSize; |
| this->maxRange[Z0_X0][0] = this->boxSize; |
| this->maxRange[Z1_X0][0] = this->boxSize; |
| this->maxRange[X0_Y0_Z0][0] = this->boxSize; |
| this->maxRange[X0_Y0_Z1][0] = this->boxSize; |
| this->maxRange[X0_Y1_Z0][0] = this->boxSize; |
| this->maxRange[X0_Y1_Z1][0] = this->boxSize; |
| } |
| |
| // Processor is on back edge in X dimension |
| if (this->layoutPos[0] == (this->layoutSize[0] - 1)) { |
| this->minRange[X1][0] = 0; |
| this->minRange[X1_Y1][0] = 0; |
| this->minRange[X1_Y0][0] = 0; |
| this->minRange[Z1_X1][0] = 0; |
| this->minRange[Z0_X1][0] = 0; |
| this->minRange[X1_Y1_Z1][0] = 0; |
| this->minRange[X1_Y1_Z0][0] = 0; |
| this->minRange[X1_Y0_Z1][0] = 0; |
| this->minRange[X1_Y0_Z0][0] = 0; |
| |
| this->maxRange[X1][0] = this->deadSize; |
| this->maxRange[X1_Y1][0] = this->deadSize; |
| this->maxRange[X1_Y0][0] = this->deadSize; |
| this->maxRange[Z1_X1][0] = this->deadSize; |
| this->maxRange[Z0_X1][0] = this->deadSize; |
| this->maxRange[X1_Y1_Z1][0] = this->deadSize; |
| this->maxRange[X1_Y1_Z0][0] = this->deadSize; |
| this->maxRange[X1_Y0_Z1][0] = this->deadSize; |
| this->maxRange[X1_Y0_Z0][0] = this->deadSize; |
| } |
| |
| // Processor is on front edge in Y dimension |
| if (this->layoutPos[1] == 0) { |
| this->minRange[Y0][1] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0][1] = this->boxSize - this->deadSize; |
| this->minRange[X1_Y0][1] = this->boxSize - this->deadSize; |
| this->minRange[Y0_Z0][1] = this->boxSize - this->deadSize; |
| this->minRange[Y0_Z1][1] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0_Z0][1] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0_Z1][1] = this->boxSize - this->deadSize; |
| this->minRange[X1_Y0_Z1][1] = this->boxSize - this->deadSize; |
| this->minRange[X1_Y0_Z0][1] = this->boxSize - this->deadSize; |
| |
| this->maxRange[Y0][1] = this->boxSize; |
| this->maxRange[X0_Y0][1] = this->boxSize; |
| this->maxRange[X1_Y0][1] = this->boxSize; |
| this->maxRange[Y0_Z0][1] = this->boxSize; |
| this->maxRange[Y0_Z1][1] = this->boxSize; |
| this->maxRange[X0_Y0_Z0][1] = this->boxSize; |
| this->maxRange[X0_Y0_Z1][1] = this->boxSize; |
| this->maxRange[X1_Y0_Z1][1] = this->boxSize; |
| this->maxRange[X1_Y0_Z0][1] = this->boxSize; |
| } |
| |
| // Processor is on back edge in Y dimension |
| if (this->layoutPos[1] == (this->layoutSize[1] - 1)) { |
| this->minRange[Y1][1] = 0; |
| this->minRange[X1_Y1][1] = 0; |
| this->minRange[X0_Y1][1] = 0; |
| this->minRange[Y1_Z1][1] = 0; |
| this->minRange[Y1_Z0][1] = 0; |
| this->minRange[X1_Y1_Z1][1] = 0; |
| this->minRange[X1_Y1_Z0][1] = 0; |
| this->minRange[X0_Y1_Z0][1] = 0; |
| this->minRange[X0_Y1_Z1][1] = 0; |
| |
| this->maxRange[Y1][1] = this->deadSize; |
| this->maxRange[X1_Y1][1] = this->deadSize; |
| this->maxRange[X0_Y1][1] = this->deadSize; |
| this->maxRange[Y1_Z1][1] = this->deadSize; |
| this->maxRange[Y1_Z0][1] = this->deadSize; |
| this->maxRange[X1_Y1_Z1][1] = this->deadSize; |
| this->maxRange[X1_Y1_Z0][1] = this->deadSize; |
| this->maxRange[X0_Y1_Z0][1] = this->deadSize; |
| this->maxRange[X0_Y1_Z1][1] = this->deadSize; |
| } |
| |
| // Processor is on front edge in Z dimension |
| if (this->layoutPos[2] == 0) { |
| this->minRange[Z0][2] = this->boxSize - this->deadSize; |
| this->minRange[Y0_Z0][2] = this->boxSize - this->deadSize; |
| this->minRange[Y1_Z0][2] = this->boxSize - this->deadSize; |
| this->minRange[Z0_X0][2] = this->boxSize - this->deadSize; |
| this->minRange[Z0_X1][2] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y0_Z0][2] = this->boxSize - this->deadSize; |
| this->minRange[X1_Y1_Z0][2] = this->boxSize - this->deadSize; |
| this->minRange[X0_Y1_Z0][2] = this->boxSize - this->deadSize; |
| this->minRange[X1_Y0_Z0][2] = this->boxSize - this->deadSize; |
| |
| this->maxRange[Z0][2] = this->boxSize; |
| this->maxRange[Y0_Z0][2] = this->boxSize; |
| this->maxRange[Y1_Z0][2] = this->boxSize; |
| this->maxRange[Z0_X0][2] = this->boxSize; |
| this->maxRange[Z0_X1][2] = this->boxSize; |
| this->maxRange[X0_Y0_Z0][2] = this->boxSize; |
| this->maxRange[X1_Y1_Z0][2] = this->boxSize; |
| this->maxRange[X0_Y1_Z0][2] = this->boxSize; |
| this->maxRange[X1_Y0_Z0][2] = this->boxSize; |
| } |
| |
| // Processor is on back edge in Z dimension |
| if (this->layoutPos[2] == (this->layoutSize[2] - 1)) { |
| this->minRange[Z1][2] = 0; |
| this->minRange[Y1_Z1][2] = 0; |
| this->minRange[Y0_Z1][2] = 0; |
| this->minRange[Z1_X1][2] = 0; |
| this->minRange[Z1_X0][2] = 0; |
| this->minRange[X1_Y1_Z1][2] = 0; |
| this->minRange[X0_Y0_Z1][2] = 0; |
| this->minRange[X1_Y0_Z1][2] = 0; |
| this->minRange[X0_Y1_Z1][2] = 0; |
| |
| this->maxRange[Z1][2] = this->deadSize; |
| this->maxRange[Y1_Z1][2] = this->deadSize; |
| this->maxRange[Y0_Z1][2] = this->deadSize; |
| this->maxRange[Z1_X1][2] = this->deadSize; |
| this->maxRange[Z1_X0][2] = this->deadSize; |
| this->maxRange[X1_Y1_Z1][2] = this->deadSize; |
| this->maxRange[X0_Y0_Z1][2] = this->deadSize; |
| this->maxRange[X1_Y0_Z1][2] = this->deadSize; |
| this->maxRange[X0_Y1_Z1][2] = this->deadSize; |
| } |
| } |
| |
| ///////////////////////////////////////////////////////////////////////// |
| // |
| // Set the particle vectors that have already been read and which |
| // contain only the alive particles for this processor |
| // |
| ///////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::setParticleArrays( |
| long count, |
| POSVEL_T* xLoc, |
| POSVEL_T* yLoc, |
| POSVEL_T* zLoc, |
| POSVEL_T* xVel, |
| POSVEL_T* yVel, |
| POSVEL_T* zVel, |
| POTENTIAL_T* potential, |
| ID_T* id, |
| MASK_T* maskData) |
| { |
| this->particleCount = count; |
| this->xxInit = xLoc; |
| this->yyInit = yLoc; |
| this->zzInit = zLoc; |
| this->vxInit = xVel; |
| this->vyInit = yVel; |
| this->vzInit = zVel; |
| this->potInit = potential; |
| this->tagInit = id; |
| this->maskInit = maskData; |
| } |
| |
| void InitialExchange::setParticleVectors( |
| vector<POSVEL_T>* xLoc, |
| vector<POSVEL_T>* yLoc, |
| vector<POSVEL_T>* zLoc, |
| vector<POSVEL_T>* xVel, |
| vector<POSVEL_T>* yVel, |
| vector<POSVEL_T>* zVel, |
| vector<POTENTIAL_T>* potential, |
| vector<ID_T>* id, |
| vector<MASK_T>* maskData, |
| vector<STATUS_T>* type) |
| { |
| this->xx = xLoc; |
| this->yy = yLoc; |
| this->zz = zLoc; |
| this->vx = xVel; |
| this->vy = yVel; |
| this->vz = zVel; |
| this->pot = potential; |
| this->tag = id; |
| this->mask = maskData; |
| this->status = type; |
| } |
| |
| ///////////////////////////////////////////////////////////////////////////// |
| // |
| // Identify the border particles which will be alive on other processors |
| // and send them and receive the particles which are dead on other processors |
| // but alive on this processor. Store all the newly acquired alive particles |
| // in the empty supplied vectors for the rest of the simulation. |
| // |
| ///////////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::exchangeParticles() |
| { |
| // Identify dead particles on this processor which must be sent |
| // because they are alive particles on neighbor processors |
| // x,y,z are still in physical units with wraparound included |
| identifyExchangeParticles(); |
| |
| // Exchange those particles with appropriate neighbors |
| // x,y,z are not in normalized units |
| exchangeNeighborParticles(); |
| |
| // Count the particles across processors |
| long totalAliveParticles = 0; |
| #ifndef USE_SERIAL_COSMO |
| MPI_Allreduce((void*) &this->numberOfAliveParticles, |
| (void*) &totalAliveParticles, |
| 1, MPI_LONG, MPI_SUM, MPI_COMM_WORLD); |
| #endif |
| |
| //cout << "InitialExchange Particles Rank " << setw(3) << this->myProc |
| // << " #alive = " << this->numberOfAliveParticles << endl; |
| |
| if (this->myProc == MASTER) { |
| cout << "InitialExchange TotalAliveParticles " |
| << totalAliveParticles << endl; |
| } |
| } |
| |
| ///////////////////////////////////////////////////////////////////////////// |
| // |
| // Iterate over all the particles on this processor and determine which are |
| // dead and must be sent away and add them to the vector for that neighbor |
| // |
| ///////////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::identifyExchangeParticles() |
| { |
| for (int i = 0; i < this->particleCount; i++) { |
| bool found = false; |
| |
| // Particle is alive on this processor so add to vectors |
| if ((this->xxInit[i] >= this->minAlive[0] && |
| this->xxInit[i] < this->maxAlive[0]) && |
| (this->yyInit[i] >= this->minAlive[1] && |
| this->yyInit[i] < this->maxAlive[1]) && |
| (this->zzInit[i] >= this->minAlive[2] && |
| this->zzInit[i] < this->maxAlive[2])) { |
| this->xx->push_back(this->xxInit[i]); |
| this->yy->push_back(this->yyInit[i]); |
| this->zz->push_back(this->zzInit[i]); |
| this->vx->push_back(this->vxInit[i]); |
| this->vy->push_back(this->vyInit[i]); |
| this->vz->push_back(this->vzInit[i]); |
| this->tag->push_back(this->tagInit[i]); |
| this->pot->push_back(this->potInit[i]); |
| this->mask->push_back(this->maskInit[i]); |
| this->numberOfAliveParticles++; |
| found = true; |
| |
| } else { |
| // Particle is dead here but which processor needs it as alive |
| for (int n = 0; n < NUM_OF_NEIGHBORS; n++) { |
| if ((this->xxInit[i] >= minRange[n][0]) && |
| (this->xxInit[i] < maxRange[n][0]) && |
| (this->yyInit[i] >= minRange[n][1]) && |
| (this->yyInit[i] < maxRange[n][1]) && |
| (this->zzInit[i] >= minRange[n][2]) && |
| (this->zzInit[i] < maxRange[n][2])) { |
| this->neighborParticles[n].push_back(i); |
| found = true; |
| } |
| } |
| } |
| if (found == false) { |
| cout << "Rank " << myProc << " Problem particle " << xxInit[i] |
| << " , " << yyInit[i] << " , " << zzInit[i] << endl; |
| } |
| } |
| } |
| |
| ///////////////////////////////////////////////////////////////////////////// |
| // |
| // Exchange the appropriate particles with neighbors |
| // Only the index of the particle to be exchanged is stored so fill out |
| // the message with location, velocity, tag. Status information doesn't |
| // have to be sent because all particles will be alive. |
| // Use the Cartesian communicator for neighbor exchange |
| // |
| ///////////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::exchangeNeighborParticles() |
| { |
| // Calculate the maximum number of particles to share for calculating buffer |
| int myShareSize = 0; |
| for (int n = 0; n < NUM_OF_NEIGHBORS; n++) |
| if (myShareSize < (int) this->neighborParticles[n].size()) |
| myShareSize = this->neighborParticles[n].size(); |
| |
| int maxDeadSize; |
| #ifndef USE_SERIAL_COSMO |
| MPI_Allreduce((void*) &myShareSize, |
| (void*) &maxDeadSize, |
| 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); |
| #endif |
| |
| // Allocate messages to send and receive MPI buffers |
| int bufferSize = (1 * sizeof(int)) + // number of particles |
| (maxDeadSize * |
| ((7 * sizeof(POSVEL_T)) + // location, velocity, potential |
| (1 * sizeof(ID_T)) + // id tag |
| (1 * sizeof(MASK_T)))); // mask |
| Message* sendMessage = new Message(bufferSize); |
| Message* recvMessage = new Message(bufferSize); |
| |
| //debug statement added by Adrian to see how much buffer space we're using |
| if(this->myProc == MASTER) { |
| printf("PXCH buffer = 2*%d = %f MB\n",bufferSize, |
| 2.0*bufferSize/1024.0/1024.0); |
| } |
| //MPI_Barrier(MPI_COMM_WORLD); |
| |
| // Exchange with each neighbor, with everyone sending in one direction and |
| // receiving from the other. Data corresponding to the particle index |
| // must be packed in the buffer. When the data is received it is unpacked |
| // into the location, velocity and tag vectors and the status is set |
| // to the neighbor who sent it |
| |
| for (int n = 0; n < NUM_OF_NEIGHBORS; n=n+2) { |
| // Neighbor pairs in Definition.h must match so that every processor |
| // sends and every processor receives on each exchange |
| exchange(n, n+1, sendMessage, recvMessage); |
| exchange(n+1, n, sendMessage, recvMessage); |
| } |
| |
| delete sendMessage; |
| delete recvMessage; |
| } |
| |
| ///////////////////////////////////////////////////////////////////////////// |
| // |
| // Pack particle data for the indicated neighbor into MPI message |
| // Send that message and receive from opposite neighbor |
| // Unpack the received particle data and add to particle buffers with |
| // an indication of dead and the neighbor on which particle is alive |
| // |
| ///////////////////////////////////////////////////////////////////////////// |
| |
| void InitialExchange::exchange( |
| int sendTo, |
| int recvFrom, |
| Message* sendMessage, |
| Message* recvMessage) |
| { |
| POSVEL_T posValue; |
| POTENTIAL_T potValue; |
| ID_T idValue; |
| MASK_T maskValue; |
| |
| #ifndef USE_SERIAL_COSMO |
| // Fill same message for each of the neighbors |
| sendMessage->reset(); |
| recvMessage->reset(); |
| |
| // Number of particles to share with neighbor |
| int sendParticleCount = this->neighborParticles[sendTo].size(); |
| |
| // Pack the number of particles being sent |
| sendMessage->putValue(&sendParticleCount); |
| |
| for (int i = 0; i < sendParticleCount; i++) { |
| int deadIndex = this->neighborParticles[sendTo][i]; |
| sendMessage->putValue(&this->xxInit[deadIndex]); |
| sendMessage->putValue(&this->yyInit[deadIndex]); |
| sendMessage->putValue(&this->zzInit[deadIndex]); |
| sendMessage->putValue(&this->vxInit[deadIndex]); |
| sendMessage->putValue(&this->vyInit[deadIndex]); |
| sendMessage->putValue(&this->vzInit[deadIndex]); |
| sendMessage->putValue(&this->potInit[deadIndex]); |
| sendMessage->putValue(&this->tagInit[deadIndex]); |
| sendMessage->putValue(&this->maskInit[deadIndex]); |
| this->particleCount--; |
| } |
| |
| // Send the message buffer |
| sendMessage->send(this->neighbor[sendTo]); |
| |
| // Receive the buffer from neighbor on other side |
| recvMessage->receive(this->neighbor[recvFrom]); |
| MPI_Barrier(Partition::getComm()); |
| |
| // Process the received buffer |
| int recvParticleCount; |
| recvMessage->getValue(&recvParticleCount); |
| |
| for (int i = 0; i < recvParticleCount; i++) { |
| recvMessage->getValue(&posValue); |
| this->xx->push_back(posValue); |
| recvMessage->getValue(&posValue); |
| this->yy->push_back(posValue); |
| recvMessage->getValue(&posValue); |
| this->zz->push_back(posValue); |
| recvMessage->getValue(&posValue); |
| this->vx->push_back(posValue); |
| recvMessage->getValue(&posValue); |
| this->vy->push_back(posValue); |
| recvMessage->getValue(&posValue); |
| this->vz->push_back(posValue); |
| recvMessage->getValue(&potValue); |
| this->pot->push_back(potValue); |
| recvMessage->getValue(&idValue); |
| this->tag->push_back(idValue); |
| recvMessage->getValue(&maskValue); |
| this->mask->push_back(maskValue); |
| |
| this->numberOfAliveParticles++; |
| this->particleCount++; |
| } |
| #endif |
| } |
