OOPSE/libmdtools/DumpWriter.cpp

#define _FILE_OFFSET_BITS 64

#include <string.h>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <utility>

#ifdef IS_MPI
#include <mpi.h>
#include "mpiSimulation.hpp"

namespace dWrite{
  void DieDieDie( void );
}

using namespace dWrite;
#endif //is_mpi

#include "ReadWrite.hpp"
#include "simError.h"

DumpWriter::DumpWriter( SimInfo* the_entry_plug ){

  entry_plug = the_entry_plug;

#ifdef IS_MPI
  if(worldRank == 0 ){
#endif // is_mpi

    dumpFile.open(entry_plug->sampleName, ios::out | ios::trunc );

    if( !dumpFile ){

      sprintf( painCave.errMsg,
               "Could not open \"%s\" for dump output.\n",
               entry_plug->sampleName);
      painCave.isFatal = 1;
      simError();
    }

#ifdef IS_MPI
  }

  //sort the local atoms by global index
  sortByGlobalIndex();
  
  sprintf( checkPointMsg,
           "Sucessfully opened output file for dumping.\n");
  MPIcheckPoint();
#endif // is_mpi
}

DumpWriter::~DumpWriter( ){

#ifdef IS_MPI
  if(worldRank == 0 ){
#endif // is_mpi

    dumpFile.close();

#ifdef IS_MPI
  }
#endif // is_mpi
}

#ifdef IS_MPI

/**
 * A hook function to load balancing
 */

void DumpWriter::update(){
  sortByGlobalIndex();          
}
  
/**
 * Auxiliary sorting function
 */
 
bool indexSortingCriterion(const pair<int, int>& p1, const pair<int, int>& p2){
  return p1.second < p2.second;
}

/**
 * Sorting the local index by global index
 */
 
void DumpWriter::sortByGlobalIndex(){
  Atom** atoms = entry_plug->atoms;
  
  indexArray.clear();
  
  for(int i = 0; i < mpiSim->getMyNlocal();i++) 
    indexArray.push_back(make_pair(i, atoms[i]->getGlobalIndex()));
  
  sort(indexArray.begin(), indexArray.end(), indexSortingCriterion);    

  //for (int i = 0; i < mpiSim->getMyNlocal(); i++) {
  //  printf("node %d has global %d at local %d\n", worldRank, indexArray[i].second, indexArray[i].first);
  //}
    
}

#endif

void DumpWriter::writeDump(double currentTime){

  ofstream finalOut;
  vector<ofstream*> fileStreams;

#ifdef IS_MPI
  if(worldRank == 0 ){
#endif    
    finalOut.open( entry_plug->finalName, ios::out | ios::trunc );
    if( !finalOut ){
      sprintf( painCave.errMsg,
               "Could not open \"%s\" for final dump output.\n",
               entry_plug->finalName );
      painCave.isFatal = 1;
      simError();
    }
#ifdef IS_MPI
  }
#endif // is_mpi

  fileStreams.push_back(&finalOut); 
  fileStreams.push_back(&dumpFile);

  writeFrame(fileStreams, currentTime);

#ifdef IS_MPI
  finalOut.close();
#endif
        
}

void DumpWriter::writeFinal(double currentTime){

  ofstream finalOut;
  vector<ofstream*> fileStreams;

#ifdef IS_MPI
  if(worldRank == 0 ){
#endif // is_mpi

    finalOut.open( entry_plug->finalName, ios::out | ios::trunc );

    if( !finalOut ){
      sprintf( painCave.errMsg,
               "Could not open \"%s\" for final dump output.\n",
               entry_plug->finalName );
      painCave.isFatal = 1;
      simError();
    }

#ifdef IS_MPI
  }
#endif // is_mpi
  
  fileStreams.push_back(&finalOut);  
  writeFrame(fileStreams, currentTime);

#ifdef IS_MPI
  finalOut.close();
#endif
  
}

void DumpWriter::writeFrame( vector<ofstream*>& outFile, double currentTime ){

  const int BUFFERSIZE = 2000;
  const int MINIBUFFERSIZE = 100;

  char tempBuffer[BUFFERSIZE];  
  char writeLine[BUFFERSIZE];

  int i, k;

#ifdef IS_MPI
  
  /*********************************************************************
   * Documentation?  You want DOCUMENTATION?
   * 
   * Why all the potatoes below?  
   *
   * To make a long story short, the original version of DumpWriter
   * worked in the most inefficient way possible.  Node 0 would 
   * poke each of the node for an individual atom's formatted data 
   * as node 0 worked its way down the global index. This was particularly 
   * inefficient since the method blocked all processors at every atom 
   * (and did it twice!).
   *
   * An intermediate version of DumpWriter could be described from Node
   * zero's perspective as follows:
   * 
   *  1) Have 100 of your friends stand in a circle.
   *  2) When you say go, have all of them start tossing potatoes at
   *     you (one at a time).
   *  3) Catch the potatoes.
   *
   * It was an improvement, but MPI has buffers and caches that could 
   * best be described in this analogy as "potato nets", so there's no 
   * need to block the processors atom-by-atom.
   * 
   * This new and improved DumpWriter works in an even more efficient 
   * way:
   * 
   *  1) Have 100 of your friend stand in a circle.
   *  2) When you say go, have them start tossing 5-pound bags of 
   *     potatoes at you.
   *  3) Once you've caught a friend's bag of potatoes,
   *     toss them a spud to let them know they can toss another bag.
   *
   * How's THAT for documentation?
   *
   *********************************************************************/

  int *potatoes;
  int myPotato;

  int nProc;
  int j, which_node, done, which_atom, local_index, currentIndex;
  double atomData6[6];
  double atomData13[13];
  int isDirectional;
  char* atomTypeString;
  char MPIatomTypeString[MINIBUFFERSIZE];

#else //is_mpi
  int nAtoms = entry_plug->n_atoms;
#endif //is_mpi

  double q[4];
  DirectionalAtom* dAtom;
  Atom** atoms = entry_plug->atoms;
  double pos[3], vel[3];

#ifndef IS_MPI
  
  for(k = 0; k < outFile.size(); k++){
    *outFile[k] << nAtoms << "\n";

    *outFile[k] << currentTime << ";\t"
               << entry_plug->Hmat[0][0] << "\t"
                     << entry_plug->Hmat[1][0] << "\t"
                     << entry_plug->Hmat[2][0] << ";\t"
               
               << entry_plug->Hmat[0][1] << "\t"
                     << entry_plug->Hmat[1][1] << "\t"
                     << entry_plug->Hmat[2][1] << ";\t"

                     << entry_plug->Hmat[0][2] << "\t"
                     << entry_plug->Hmat[1][2] << "\t"
                     << entry_plug->Hmat[2][2] << ";";

    //write out additional parameters, such as chi and eta
    *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
  }
  
  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getPos(pos);
    atoms[i]->getVel(vel);

    sprintf( tempBuffer,
             "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
             atoms[i]->getType(),
             pos[0],
             pos[1],
             pos[2],
             vel[0],
             vel[1],
             vel[2]);
    strcpy( writeLine, tempBuffer );

    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];
      dAtom->getQ( q );

      sprintf( tempBuffer,
               "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
               q[0],
               q[1],
               q[2],
               q[3],
               dAtom->getJx(),
               dAtom->getJy(),
               dAtom->getJz());
      strcat( writeLine, tempBuffer );
    }
    else
      strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );

    for(k = 0; k < outFile.size(); k++)
      *outFile[k] << writeLine;
  }

#else // is_mpi

  /* code to find maximum tag value */
  
  int *tagub, flag, MAXTAG;
  MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
  if (flag) {
    MAXTAG = *tagub;
  } else {
    MAXTAG = 32767;
  }  

  int haveError;

  MPI_Status istatus;
  int *AtomToProcMap = mpiSim->getAtomToProcMap();

  // write out header and node 0's coordinates

  if( worldRank == 0 ){

    // Node 0 needs a list of the magic potatoes for each processor;

    nProc = mpiSim->getNumberProcessors();
    potatoes = new int[nProc];

    //write out the comment lines
    for (i = 0; i < nProc; i++) 
      potatoes[i] = 0;
    
      for(k = 0; k < outFile.size(); k++){
        *outFile[k] << mpiSim->getTotAtoms() << "\n";

        *outFile[k] << currentTime << ";\t"
                         << entry_plug->Hmat[0][0] << "\t"
                         << entry_plug->Hmat[1][0] << "\t"
                         << entry_plug->Hmat[2][0] << ";\t"

                         << entry_plug->Hmat[0][1] << "\t"
                         << entry_plug->Hmat[1][1] << "\t"
                         << entry_plug->Hmat[2][1] << ";\t"

                         << entry_plug->Hmat[0][2] << "\t"
                         << entry_plug->Hmat[1][2] << "\t"
                         << entry_plug->Hmat[2][2] << ";";
  
        *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
    }

    currentIndex = 0;

    for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
      
      // Get the Node number which has this atom;
      
      which_node = AtomToProcMap[i];
      
      if (which_node != 0) {

        if (potatoes[which_node] + 3 >= MAXTAG) {
          // The potato was going to exceed the maximum value, 
          // so wrap this processor potato back to 0:         

          potatoes[which_node] = 0;          
          MPI_Send(0, 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
          
        }

        myPotato = potatoes[which_node];        
        
        MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
                 myPotato, MPI_COMM_WORLD, &istatus);
        
        atomTypeString = MPIatomTypeString;
        
        myPotato++;

        MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
                 myPotato, MPI_COMM_WORLD, &istatus);
               
        myPotato++;

        if (isDirectional) {          
          MPI_Recv(atomData13, 13, MPI_DOUBLE, which_node,
                   myPotato, MPI_COMM_WORLD, &istatus);
        } else {
          MPI_Recv(atomData6, 6, MPI_DOUBLE, which_node,
                   myPotato, MPI_COMM_WORLD, &istatus);          
        }
        
        myPotato++;
        potatoes[which_node] = myPotato;

      } else {
        
        haveError = 0;
        which_atom = i;
        
        //local_index = -1;

        //for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
        //  if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
        //}
        
        //if (local_index != -1) {
          
          local_index = indexArray[currentIndex].first;        
          
          if (which_atom == indexArray[currentIndex].second) {
            
            atomTypeString = atoms[local_index]->getType();
            
          atoms[local_index]->getPos(pos);
          atoms[local_index]->getVel(vel);          
          
          atomData6[0] = pos[0];
          atomData6[1] = pos[1];
          atomData6[2] = pos[2];

          atomData6[3] = vel[0];
          atomData6[4] = vel[1];
          atomData6[5] = vel[2];
          
          isDirectional = 0;

          if( atoms[local_index]->isDirectional() ){

            isDirectional = 1;
            
            dAtom = (DirectionalAtom *)atoms[local_index];
            dAtom->getQ( q );

            for (int j = 0; j < 6 ; j++)
              atomData13[j] = atomData6[j];            
            
            atomData13[6] = q[0];
            atomData13[7] = q[1];
            atomData13[8] = q[2];
            atomData13[9] = q[3];
            
            atomData13[10] = dAtom->getJx();
            atomData13[11] = dAtom->getJy();
            atomData13[12] = dAtom->getJz();
          }
          
        } else {
          sprintf(painCave.errMsg,
                  "Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
                  which_atom, worldRank, currentIndex, local_index );
          haveError= 1;
          simError();
        }
        
        if(haveError) DieDieDie();
        
        currentIndex++;
      }
      // If we've survived to here, format the line:
      
      if (!isDirectional) {
        
        sprintf( writeLine,
                 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
                 atomTypeString,
                 atomData6[0],
                 atomData6[1],
                 atomData6[2],
                 atomData6[3],
                 atomData6[4],
                 atomData6[5]);
        
        strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
        
      } else {
        
        sprintf( writeLine,
                 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
                 atomTypeString,
                 atomData13[0],
                 atomData13[1],
                 atomData13[2],
                 atomData13[3],
                 atomData13[4],
                 atomData13[5],
                 atomData13[6],
                 atomData13[7],
                 atomData13[8],
                 atomData13[9],
                 atomData13[10],
                 atomData13[11],
                 atomData13[12]);
        
      }
      
      for(k = 0; k < outFile.size(); k++)
        *outFile[k] << writeLine;
    }
    
    for(k = 0; k < outFile.size(); k++)
      outFile[k]->flush();
    
    sprintf( checkPointMsg,
             "Sucessfully took a dump.\n");
    
    MPIcheckPoint();        
    
    delete[] potatoes;
    
  } else {

    // worldRank != 0, so I'm a remote node.  

    // Set my magic potato to 0:

    myPotato = 0;
    currentIndex = 0;
    
    for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
      
      // Am I the node which has this atom?
      
      if (AtomToProcMap[i] == worldRank) {

        if (myPotato + 3 >= MAXTAG) {
          
          // The potato was going to exceed the maximum value, 
          // so wrap this processor potato back to 0 (and block until
          // node 0 says we can go:
          
          MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
          
        }
        which_atom = i; 

        //local_index = -1;

        //for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
        // if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
        //}
        
        //if (local_index != -1) {

        local_index = indexArray[currentIndex].first;        
                
        if (which_atom == indexArray[currentIndex].second) {
         
          atomTypeString = atoms[local_index]->getType();
          
          atoms[local_index]->getPos(pos);
          atoms[local_index]->getVel(vel);
          
          atomData6[0] = pos[0];
          atomData6[1] = pos[1];
          atomData6[2] = pos[2];

          atomData6[3] = vel[0];
          atomData6[4] = vel[1];
          atomData6[5] = vel[2];
          
          isDirectional = 0;

          if( atoms[local_index]->isDirectional() ){

            isDirectional = 1;
            
            dAtom = (DirectionalAtom *)atoms[local_index];
            dAtom->getQ( q );
            
            for (int j = 0; j < 6 ; j++)
              atomData13[j] = atomData6[j];
            
            atomData13[6] = q[0];
            atomData13[7] = q[1];
            atomData13[8] = q[2];
            atomData13[9] = q[3];
  
            atomData13[10] = dAtom->getJx();
            atomData13[11] = dAtom->getJy();
            atomData13[12] = dAtom->getJz();
          }

        } else {
          sprintf(painCave.errMsg,
                  "Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
                  which_atom, worldRank, currentIndex, local_index );
          haveError= 1;
          simError();
        }
        
        strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);

        // null terminate the string before sending (just in case):
        MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';

        MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
                             myPotato, MPI_COMM_WORLD);
        
        myPotato++;

        MPI_Send(&isDirectional, 1, MPI_INT, 0,
                             myPotato, MPI_COMM_WORLD);
        
        myPotato++;
        
        if (isDirectional) {

          MPI_Send(atomData13, 13, MPI_DOUBLE, 0,
                   myPotato, MPI_COMM_WORLD);
          
        } else {

          MPI_Send(atomData6, 6, MPI_DOUBLE, 0,
                   myPotato, MPI_COMM_WORLD);
        }

        myPotato++;  
        currentIndex++;    
      }
    }

    sprintf( checkPointMsg,
             "Sucessfully took a dump.\n");
    MPIcheckPoint();        
    
  } 
  
#endif // is_mpi
}

#ifdef IS_MPI

// a couple of functions to let us escape the write loop

void dWrite::DieDieDie( void ){

  MPI_Finalize();
  exit (0);
}

#endif //is_mpi
Revision:	952
Committed:	Fri Jan 16 21:55:39 2004 UTC (20 years, 5 months ago) by tim
File size:	15432 byte(s)
Log Message:	fix a bug in creating eor file
#	User	Rev	Content
1	mmeineke	723	#define _FILE_OFFSET_BITS 64
2
3	gezelter	829	#include <string.h>
4	mmeineke	377	#include <iostream>
5			#include <fstream>
6	tim	929	#include <algorithm>
7			#include <utility>
8	mmeineke	377
9			#ifdef IS_MPI
10			#include <mpi.h>
11			#include "mpiSimulation.hpp"
12	mmeineke	440
13			namespace dWrite{
14	gezelter	907	void DieDieDie( void );
15	mmeineke	440	}
16
17			using namespace dWrite;
18	mmeineke	377	#endif //is_mpi
19
20			#include "ReadWrite.hpp"
21			#include "simError.h"
22
23			DumpWriter::DumpWriter( SimInfo* the_entry_plug ){
24
25			entry_plug = the_entry_plug;
26
27			#ifdef IS_MPI
28			if(worldRank == 0 ){
29			#endif // is_mpi
30	tim	837
31	tim	929	dumpFile.open(entry_plug->sampleName, ios::out \| ios::trunc );
32	tim	837
33	tim	929	if( !dumpFile ){
34	tim	837
35	mmeineke	377	sprintf( painCave.errMsg,
36			"Could not open \"%s\" for dump output.\n",
37	tim	929	entry_plug->sampleName);
38	mmeineke	377	painCave.isFatal = 1;
39			simError();
40			}
41	mmeineke	469
42	mmeineke	377	#ifdef IS_MPI
43			}
44
45	tim	929	//sort the local atoms by global index
46			sortByGlobalIndex();
47
48	mmeineke	377	sprintf( checkPointMsg,
49			"Sucessfully opened output file for dumping.\n");
50			MPIcheckPoint();
51			#endif // is_mpi
52			}
53
54			DumpWriter::~DumpWriter( ){
55
56			#ifdef IS_MPI
57			if(worldRank == 0 ){
58			#endif // is_mpi
59
60	tim	929	dumpFile.close();
61	mmeineke	377
62			#ifdef IS_MPI
63			}
64			#endif // is_mpi
65			}
66
67	tim	929	#ifdef IS_MPI
68	tim	837
69	tim	929	/**
70			* A hook function to load balancing
71			*/
72
73			void DumpWriter::update(){
74			sortByGlobalIndex();
75			}
76
77			/**
78			* Auxiliary sorting function
79			*/
80
81			bool indexSortingCriterion(const pair<int, int>& p1, const pair<int, int>& p2){
82			return p1.second < p2.second;
83			}
84
85			/**
86			* Sorting the local index by global index
87			*/
88
89			void DumpWriter::sortByGlobalIndex(){
90			Atom** atoms = entry_plug->atoms;
91
92			indexArray.clear();
93
94	chuckv	949	for(int i = 0; i < mpiSim->getMyNlocal();i++)
95	tim	929	indexArray.push_back(make_pair(i, atoms[i]->getGlobalIndex()));
96
97			sort(indexArray.begin(), indexArray.end(), indexSortingCriterion);
98	chuckv	949
99			//for (int i = 0; i < mpiSim->getMyNlocal(); i++) {
100			// printf("node %d has global %d at local %d\n", worldRank, indexArray[i].second, indexArray[i].first);
101			//}
102
103	tim	929	}
104	chuckv	949
105	tim	929	#endif
106
107			void DumpWriter::writeDump(double currentTime){
108
109	tim	936	ofstream finalOut;
110	tim	934	vector<ofstream*> fileStreams;
111
112			#ifdef IS_MPI
113			if(worldRank == 0 ){
114	tim	952	#endif
115	tim	936	finalOut.open( entry_plug->finalName, ios::out \| ios::trunc );
116			if( !finalOut ){
117			sprintf( painCave.errMsg,
118			"Could not open \"%s\" for final dump output.\n",
119			entry_plug->finalName );
120			painCave.isFatal = 1;
121			simError();
122			}
123	tim	952	#ifdef IS_MPI
124	tim	934	}
125			#endif // is_mpi
126
127			fileStreams.push_back(&finalOut);
128			fileStreams.push_back(&dumpFile);
129
130			writeFrame(fileStreams, currentTime);
131	tim	936
132			#ifdef IS_MPI
133			finalOut.close();
134			#endif
135	tim	929
136			}
137
138			void DumpWriter::writeFinal(double currentTime){
139
140	tim	936	ofstream finalOut;
141	tim	934	vector<ofstream*> fileStreams;
142
143	tim	929	#ifdef IS_MPI
144			if(worldRank == 0 ){
145	mmeineke	951	#endif // is_mpi
146	tim	936
147			finalOut.open( entry_plug->finalName, ios::out \| ios::trunc );
148
149			if( !finalOut ){
150			sprintf( painCave.errMsg,
151			"Could not open \"%s\" for final dump output.\n",
152			entry_plug->finalName );
153			painCave.isFatal = 1;
154			simError();
155			}
156
157	mmeineke	951	#ifdef IS_MPI
158	tim	934	}
159	tim	929	#endif // is_mpi
160
161	tim	934	fileStreams.push_back(&finalOut);
162			writeFrame(fileStreams, currentTime);
163	tim	936
164			#ifdef IS_MPI
165			finalOut.close();
166			#endif
167	tim	929
168			}
169
170	tim	934	void DumpWriter::writeFrame( vector<ofstream*>& outFile, double currentTime ){
171	tim	929
172	mmeineke	377	const int BUFFERSIZE = 2000;
173	gezelter	912	const int MINIBUFFERSIZE = 100;
174	gezelter	907
175	tim	936	char tempBuffer[BUFFERSIZE];
176	mmeineke	377	char writeLine[BUFFERSIZE];
177
178	tim	934	int i, k;
179	gezelter	916
180	mmeineke	787	#ifdef IS_MPI
181	gezelter	916
182	gezelter	947	/*********************************************************************
183			* Documentation? You want DOCUMENTATION?
184			*
185			* Why all the potatoes below?
186			*
187			* To make a long story short, the original version of DumpWriter
188			* worked in the most inefficient way possible. Node 0 would
189			* poke each of the node for an individual atom's formatted data
190			* as node 0 worked its way down the global index. This was particularly
191			* inefficient since the method blocked all processors at every atom
192			* (and did it twice!).
193			*
194			* An intermediate version of DumpWriter could be described from Node
195			* zero's perspective as follows:
196			*
197			* 1) Have 100 of your friends stand in a circle.
198			* 2) When you say go, have all of them start tossing potatoes at
199			* you (one at a time).
200			* 3) Catch the potatoes.
201			*
202			* It was an improvement, but MPI has buffers and caches that could
203			* best be described in this analogy as "potato nets", so there's no
204			* need to block the processors atom-by-atom.
205			*
206			* This new and improved DumpWriter works in an even more efficient
207			* way:
208			*
209			* 1) Have 100 of your friend stand in a circle.
210			* 2) When you say go, have them start tossing 5-pound bags of
211			* potatoes at you.
212			* 3) Once you've caught a friend's bag of potatoes,
213			* toss them a spud to let them know they can toss another bag.
214			*
215			* How's THAT for documentation?
216			*
217			*********************************************************************/
218
219	gezelter	916	int *potatoes;
220			int myPotato;
221
222			int nProc;
223	tim	929	int j, which_node, done, which_atom, local_index, currentIndex;
224	gezelter	916	double atomData6[6];
225			double atomData13[13];
226	gezelter	907	int isDirectional;
227			char* atomTypeString;
228	gezelter	910	char MPIatomTypeString[MINIBUFFERSIZE];
229	gezelter	916
230	mmeineke	787	#else //is_mpi
231			int nAtoms = entry_plug->n_atoms;
232			#endif //is_mpi
233
234	mmeineke	377	double q[4];
235			DirectionalAtom* dAtom;
236			Atom** atoms = entry_plug->atoms;
237	mmeineke	670	double pos[3], vel[3];
238	tim	837
239	mmeineke	377	#ifndef IS_MPI
240	tim	934
241			for(k = 0; k < outFile.size(); k++){
242			*outFile[k] << nAtoms << "\n";
243	tim	837
244	tim	934	*outFile[k] << currentTime << ";\t"
245			<< entry_plug->Hmat[0][0] << "\t"
246			<< entry_plug->Hmat[1][0] << "\t"
247			<< entry_plug->Hmat[2][0] << ";\t"
248
249			<< entry_plug->Hmat[0][1] << "\t"
250			<< entry_plug->Hmat[1][1] << "\t"
251			<< entry_plug->Hmat[2][1] << ";\t"
252	tim	837
253	tim	934	<< entry_plug->Hmat[0][2] << "\t"
254			<< entry_plug->Hmat[1][2] << "\t"
255			<< entry_plug->Hmat[2][2] << ";";
256	mmeineke	572
257	tim	934	//write out additional parameters, such as chi and eta
258			*outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
259			}
260
261	mmeineke	377	for( i=0; i<nAtoms; i++ ){
262	tim	837
263	mmeineke	670	atoms[i]->getPos(pos);
264			atoms[i]->getVel(vel);
265	mmeineke	377
266			sprintf( tempBuffer,
267			"%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
268			atoms[i]->getType(),
269	mmeineke	670	pos[0],
270			pos[1],
271			pos[2],
272			vel[0],
273			vel[1],
274			vel[2]);
275	mmeineke	377	strcpy( writeLine, tempBuffer );
276
277			if( atoms[i]->isDirectional() ){
278	tim	837
279	mmeineke	377	dAtom = (DirectionalAtom *)atoms[i];
280			dAtom->getQ( q );
281	tim	837
282	mmeineke	377	sprintf( tempBuffer,
283			"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
284			q[0],
285			q[1],
286			q[2],
287			q[3],
288			dAtom->getJx(),
289			dAtom->getJy(),
290			dAtom->getJz());
291			strcat( writeLine, tempBuffer );
292			}
293			else
294			strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
295	tim	837
296	tim	934	for(k = 0; k < outFile.size(); k++)
297			*outFile[k] << writeLine;
298	mmeineke	377	}
299
300			#else // is_mpi
301	gezelter	416
302	chuckv	913	/* code to find maximum tag value */
303	tim	919
304	tim	920	int *tagub, flag, MAXTAG;
305	chuckv	913	MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
306			if (flag) {
307	tim	920	MAXTAG = *tagub;
308	chuckv	913	} else {
309			MAXTAG = 32767;
310	gezelter	916	}
311	mmeineke	440
312			int haveError;
313
314	mmeineke	447	MPI_Status istatus;
315	gezelter	416	int *AtomToProcMap = mpiSim->getAtomToProcMap();
316	tim	837
317	mmeineke	377	// write out header and node 0's coordinates
318	tim	837
319	mmeineke	377	if( worldRank == 0 ){
320	gezelter	916
321			// Node 0 needs a list of the magic potatoes for each processor;
322
323			nProc = mpiSim->getNumberProcessors();
324			potatoes = new int[nProc];
325
326	tim	934	//write out the comment lines
327	gezelter	916	for (i = 0; i < nProc; i++)
328			potatoes[i] = 0;
329
330	tim	934	for(k = 0; k < outFile.size(); k++){
331			*outFile[k] << mpiSim->getTotAtoms() << "\n";
332	tim	837
333	tim	934	*outFile[k] << currentTime << ";\t"
334			<< entry_plug->Hmat[0][0] << "\t"
335			<< entry_plug->Hmat[1][0] << "\t"
336			<< entry_plug->Hmat[2][0] << ";\t"
337	tim	837
338	tim	934	<< entry_plug->Hmat[0][1] << "\t"
339			<< entry_plug->Hmat[1][1] << "\t"
340			<< entry_plug->Hmat[2][1] << ";\t"
341	tim	837
342	tim	934	<< entry_plug->Hmat[0][2] << "\t"
343			<< entry_plug->Hmat[1][2] << "\t"
344			<< entry_plug->Hmat[2][2] << ";";
345
346			*outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
347			}
348	tim	837
349	tim	929	currentIndex = 0;
350	tim	934
351	gezelter	416	for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
352	chuckv	913
353	gezelter	417	// Get the Node number which has this atom;
354	chuckv	913
355	tim	837	which_node = AtomToProcMap[i];
356	chuckv	913
357	gezelter	907	if (which_node != 0) {
358	gezelter	916
359			if (potatoes[which_node] + 3 >= MAXTAG) {
360			// The potato was going to exceed the maximum value,
361			// so wrap this processor potato back to 0:
362
363			potatoes[which_node] = 0;
364			MPI_Send(0, 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
365
366			}
367
368			myPotato = potatoes[which_node];
369	gezelter	907
370	gezelter	910	MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
371	gezelter	916	myPotato, MPI_COMM_WORLD, &istatus);
372	gezelter	907
373	tim	920	atomTypeString = MPIatomTypeString;
374
375	gezelter	916	myPotato++;
376
377	gezelter	907	MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
378	gezelter	916	myPotato, MPI_COMM_WORLD, &istatus);
379
380			myPotato++;
381	gezelter	907
382	gezelter	916	if (isDirectional) {
383			MPI_Recv(atomData13, 13, MPI_DOUBLE, which_node,
384			myPotato, MPI_COMM_WORLD, &istatus);
385			} else {
386			MPI_Recv(atomData6, 6, MPI_DOUBLE, which_node,
387			myPotato, MPI_COMM_WORLD, &istatus);
388	gezelter	907	}
389	gezelter	916
390			myPotato++;
391			potatoes[which_node] = myPotato;
392	gezelter	907
393			} else {
394
395	tim	934	haveError = 0;
396	chuckv	949	which_atom = i;
397	gezelter	916
398	chuckv	949	//local_index = -1;
399
400			//for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
401			// if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
402			//}
403
404			//if (local_index != -1) {
405
406			local_index = indexArray[currentIndex].first;
407	gezelter	916
408	chuckv	949	if (which_atom == indexArray[currentIndex].second) {
409
410			atomTypeString = atoms[local_index]->getType();
411
412			atoms[local_index]->getPos(pos);
413			atoms[local_index]->getVel(vel);
414
415	gezelter	916	atomData6[0] = pos[0];
416			atomData6[1] = pos[1];
417			atomData6[2] = pos[2];
418	tim	837
419	gezelter	916	atomData6[3] = vel[0];
420			atomData6[4] = vel[1];
421			atomData6[5] = vel[2];
422	gezelter	907
423			isDirectional = 0;
424
425	chuckv	436	if( atoms[local_index]->isDirectional() ){
426	tim	837
427	gezelter	907	isDirectional = 1;
428
429	chuckv	436	dAtom = (DirectionalAtom *)atoms[local_index];
430			dAtom->getQ( q );
431	gezelter	916
432			for (int j = 0; j < 6 ; j++)
433			atomData13[j] = atomData6[j];
434	gezelter	907
435	gezelter	916	atomData13[6] = q[0];
436			atomData13[7] = q[1];
437			atomData13[8] = q[2];
438			atomData13[9] = q[3];
439
440			atomData13[10] = dAtom->getJx();
441			atomData13[11] = dAtom->getJy();
442			atomData13[12] = dAtom->getJz();
443	gezelter	907	}
444	gezelter	916
445	gezelter	907	} else {
446	chuckv	949	sprintf(painCave.errMsg,
447			"Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
448			which_atom, worldRank, currentIndex, local_index );
449			haveError= 1;
450			simError();
451			}
452	gezelter	916
453	tim	934	if(haveError) DieDieDie();
454	gezelter	916
455	chuckv	949	currentIndex++;
456	tim	926	}
457			// If we've survived to here, format the line:
458
459			if (!isDirectional) {
460
461	tim	934	sprintf( writeLine,
462	chuckv	949	"%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
463			atomTypeString,
464			atomData6[0],
465			atomData6[1],
466			atomData6[2],
467			atomData6[3],
468			atomData6[4],
469			atomData6[5]);
470	tim	926
471	chuckv	949	strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
472
473	tim	926	} else {
474
475	chuckv	949	sprintf( writeLine,
476			"%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
477			atomTypeString,
478			atomData13[0],
479			atomData13[1],
480			atomData13[2],
481			atomData13[3],
482			atomData13[4],
483			atomData13[5],
484			atomData13[6],
485			atomData13[7],
486			atomData13[8],
487			atomData13[9],
488			atomData13[10],
489			atomData13[11],
490			atomData13[12]);
491	gezelter	907
492			}
493	tim	926
494	tim	934	for(k = 0; k < outFile.size(); k++)
495			*outFile[k] << writeLine;
496	mmeineke	377	}
497	tim	926
498	tim	934	for(k = 0; k < outFile.size(); k++)
499			outFile[k]->flush();
500
501	gezelter	907	sprintf( checkPointMsg,
502			"Sucessfully took a dump.\n");
503	chuckv	949
504	gezelter	907	MPIcheckPoint();
505	chuckv	949
506	tim	919	delete[] potatoes;
507	chuckv	949
508	gezelter	415	} else {
509	tim	837
510	gezelter	907	// worldRank != 0, so I'm a remote node.
511	gezelter	916
512			// Set my magic potato to 0:
513
514			myPotato = 0;
515	tim	929	currentIndex = 0;
516	gezelter	907
517			for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
518
519			// Am I the node which has this atom?
520
521			if (AtomToProcMap[i] == worldRank) {
522	tim	837
523	gezelter	916	if (myPotato + 3 >= MAXTAG) {
524	chuckv	949
525	gezelter	916	// The potato was going to exceed the maximum value,
526			// so wrap this processor potato back to 0 (and block until
527			// node 0 says we can go:
528	chuckv	949
529	gezelter	916	MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
530
531			}
532	tim	919	which_atom = i;
533	chuckv	949
534			//local_index = -1;
535
536			//for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
537			// if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
538			//}
539
540			//if (local_index != -1) {
541
542			local_index = indexArray[currentIndex].first;
543	tim	920
544	chuckv	949	if (which_atom == indexArray[currentIndex].second) {
545	gezelter	907
546			atomTypeString = atoms[local_index]->getType();
547	chuckv	949
548			atoms[local_index]->getPos(pos);
549			atoms[local_index]->getVel(vel);
550
551	gezelter	916	atomData6[0] = pos[0];
552			atomData6[1] = pos[1];
553			atomData6[2] = pos[2];
554	tim	837
555	gezelter	916	atomData6[3] = vel[0];
556			atomData6[4] = vel[1];
557			atomData6[5] = vel[2];
558	gezelter	907
559			isDirectional = 0;
560	tim	837
561	gezelter	907	if( atoms[local_index]->isDirectional() ){
562	mmeineke	377
563	gezelter	907	isDirectional = 1;
564
565			dAtom = (DirectionalAtom *)atoms[local_index];
566			dAtom->getQ( q );
567
568	gezelter	916	for (int j = 0; j < 6 ; j++)
569			atomData13[j] = atomData6[j];
570
571			atomData13[6] = q[0];
572			atomData13[7] = q[1];
573			atomData13[8] = q[2];
574			atomData13[9] = q[3];
575	tim	934
576	gezelter	916	atomData13[10] = dAtom->getJx();
577			atomData13[11] = dAtom->getJy();
578			atomData13[12] = dAtom->getJz();
579	gezelter	907	}
580	tim	837
581	gezelter	907	} else {
582	chuckv	949	sprintf(painCave.errMsg,
583			"Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
584			which_atom, worldRank, currentIndex, local_index );
585			haveError= 1;
586			simError();
587			}
588
589	gezelter	916	strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
590	tim	837
591	gezelter	916	// null terminate the string before sending (just in case):
592			MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
593	tim	837
594	gezelter	910	MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
595	tim	934	myPotato, MPI_COMM_WORLD);
596	gezelter	907
597	gezelter	916	myPotato++;
598
599	gezelter	907	MPI_Send(&isDirectional, 1, MPI_INT, 0,
600	tim	934	myPotato, MPI_COMM_WORLD);
601	gezelter	907
602	gezelter	916	myPotato++;
603
604	gezelter	907	if (isDirectional) {
605
606	gezelter	916	MPI_Send(atomData13, 13, MPI_DOUBLE, 0,
607			myPotato, MPI_COMM_WORLD);
608	gezelter	907
609	gezelter	916	} else {
610
611			MPI_Send(atomData6, 6, MPI_DOUBLE, 0,
612			myPotato, MPI_COMM_WORLD);
613	gezelter	907	}
614	gezelter	916
615	tim	929	myPotato++;
616			currentIndex++;
617	mmeineke	377	}
618	gezelter	907	}
619	mmeineke	440
620	gezelter	907	sprintf( checkPointMsg,
621	gezelter	916	"Sucessfully took a dump.\n");
622	gezelter	907	MPIcheckPoint();
623
624	gezelter	916	}
625	gezelter	907
626	mmeineke	377	#endif // is_mpi
627			}
628	mmeineke	440
629			#ifdef IS_MPI
630
631			// a couple of functions to let us escape the write loop
632
633	gezelter	907	void dWrite::DieDieDie( void ){
634	tim	837
635	mmeineke	440	MPI_Finalize();
636			exit (0);
637			}
638
639			#endif //is_mpi