src/io/RestraintWriter.cpp

#define _LARGEFILE_SOURCE64
#define _FILE_OFFSET_BITS 64

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

#ifdef IS_MPI
#include <mpi.h>
#include "brains/mpiSimulation.hpp"
#endif // is_mpi

#include "io/ReadWrite.hpp"
#include "utils/simError.h"

RestraintWriter::RestraintWriter( SimInfo* the_entry_plug ){

  entry_plug = the_entry_plug;

#ifdef IS_MPI
  //sort the local atoms by global index
  sortByGlobalIndex();
#endif // is_mpi

}

RestraintWriter::~RestraintWriter( ){}

#ifdef IS_MPI

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

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

/**
 * Sorting the local index by global index
 */
 
void RestraintWriter::sortByGlobalIndex(){
  Molecule* mols = entry_plug->molecules;  
  indexArray.clear();
  
  for(int i = 0; i < entry_plug->n_mol;i++) 
    indexArray.push_back(make_pair(i, mols[i].getGlobalIndex()));
  
  sort(indexArray.begin(), indexArray.end(), indexSortingCriterion2);   
}

#endif

void RestraintWriter::writeZangle(double currentTime){

  ofstream angleOut;
  vector<ofstream*> fileStreams;

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

  fileStreams.push_back(&angleOut); 
  writeFrame(fileStreams, currentTime);

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

void RestraintWriter::writeZangle(double currentTime, const char *in_name){

  ofstream finalOut;
  vector<ofstream*> fileStreams;

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

  fileStreams.push_back(&finalOut);
  writeFrame(fileStreams, currentTime);

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

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

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

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

  int i;
  unsigned int k;

#ifdef IS_MPI
  
  /*********************************************************************
   * Documentation?  You want DOCUMENTATION?
   * 
   * Why all the potatoes below?  
   *
   * To make a long story short, the original version of RestraintWriter
   * 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 RestraintWriter 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 RestraintWriter 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 atomData;
  int isDirectional;
  char* atomTypeString;
  char MPIatomTypeString[MINIBUFFERSIZE];
  int nObjects;
  int msgLen; // the length of message actually recieved at master nodes
#endif //is_mpi

  double angle;
  DirectionalAtom* dAtom;
  int nTotObjects;
  StuntDouble* sd;
  char* molName;
  vector<StuntDouble*> integrableObjects;
  vector<StuntDouble*>::iterator iter;
  nTotObjects = entry_plug->getTotIntegrableObjects();
#ifndef IS_MPI
  
  for(k = 0; k < outFile.size(); k++)
    *outFile[k] << currentTime << " : omega values at this time\n";

  for( i=0; i<nTotObjects; i++ ){
    
    integrableObjects = entry_plug->molecules[i].getIntegrableObjects();
    
    for( iter = integrableObjects.begin();iter !=  integrableObjects.end(); ++iter){
      sd = *iter;
      
      sprintf( tempBuffer,
               "%14.10lf\n",
               sd->getZangle());
      strcpy( writeLine, tempBuffer );

      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 nCurObj;
  int *MolToProcMap = mpiSim->getMolToProcMap();

  // 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->getNProcessors();
    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] << currentTime << " fs: omega values at this time\n";
    
    currentIndex = 0;
    
    for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
      
      // Get the Node number which has this atom;
      
      which_node = MolToProcMap[i];
      
      if (which_node != 0) {
        
        if (potatoes[which_node] + 1 >= MAXTAG) {
          // The potato was going to exceed the maximum value, 
          // so wrap this processor potato back to 0:         
          
          potatoes[which_node] = 0;          
          MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, 
                   MPI_COMM_WORLD);
          
        }
        
        myPotato = potatoes[which_node];        
        
        //recieve the number of integrableObject in current molecule
        MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
                 myPotato, MPI_COMM_WORLD, &istatus);
        myPotato++;
        
        for(int l = 0; l < nCurObj; l++){
          
          if (potatoes[which_node] + 1 >= MAXTAG) {
            // The potato was going to exceed the maximum value, 
            // so wrap this processor potato back to 0:         
            
            potatoes[which_node] = 0;          
            MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, 
                     MPI_COMM_WORLD);
            
          }
                  
          MPI_Recv(&atomData, 1, MPI_DOUBLE, which_node, myPotato, MPI_COMM_WORLD, &istatus);
          myPotato++;
          
          // If we've survived to here, format the line:
          sprintf( writeLine,
                   "%14.10lf\n",
                   atomData);

          for(k = 0; k < outFile.size(); k++)
            *outFile[k] << writeLine;            
          
        }// end for(int l =0)
        potatoes[which_node] = myPotato;
        
      }
      else {
        
        haveError = 0;
        
        local_index = indexArray[currentIndex].first;        
        
        integrableObjects = (entry_plug->molecules[local_index]).getIntegrableObjects(); 
        
        for(iter= integrableObjects.begin(); iter != integrableObjects.end(); ++iter){    
          sd = *iter;
          atomData = sd->getZangle();
          
          // If we've survived to here, format the line:
          
          sprintf( writeLine,
                   "%14.10lf\n",
                   atomData);
          
          for(k = 0; k < outFile.size(); k++)
            *outFile[k] << writeLine;

        }//end for(iter = integrableObject.begin())
        
        currentIndex++;
      }

    }//end for(i = 0; i < mpiSim->getNmol())

    for(k = 0; k < outFile.size(); k++)
      outFile[k]->flush();
    
    sprintf( checkPointMsg,
             "Successfully printed a zAngle.\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->getNMolGlobal(); i++ ) {
      
      // Am I the node which has this integrableObject?
      
      if (MolToProcMap[i] == worldRank) {


        if (myPotato + 1 >= 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);
        }
        
        local_index = indexArray[currentIndex].first;        
        integrableObjects = entry_plug->molecules[local_index].getIntegrableObjects();
        
        nCurObj = integrableObjects.size();
        
        MPI_Send(&nCurObj, 1, MPI_INT, 0,
                 myPotato, MPI_COMM_WORLD);
        myPotato++;
        
        for( iter = integrableObjects.begin(); 
             iter != integrableObjects.end(); iter++ ){
          
          if (myPotato + 1 >= 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);
          }
          
          sd = *iter;
          
          atomData = sd->getZangle();
          
          MPI_Send(&atomData, 1, MPI_DOUBLE, 0,
                   myPotato, MPI_COMM_WORLD);
          
          myPotato++;  
        }
        
        currentIndex++;    
      }
    }
    
    sprintf( checkPointMsg,
             "Successfully dropped a zAngle.\n");
    MPIcheckPoint();                
 }
#endif // is_mpi
}

Revision:	1772
Committed:	Tue Nov 23 22:48:31 2004 UTC (19 years, 7 months ago) by chrisfen
File size:	10034 byte(s)
Log Message:	Improvements to restraints
#	User	Rev	Content
1	chrisfen	1772	#define _LARGEFILE_SOURCE64
2			#define _FILE_OFFSET_BITS 64
3
4			#include <string.h>
5			#include <iostream>
6			#include <fstream>
7			#include <algorithm>
8			#include <utility>
9
10			#ifdef IS_MPI
11			#include <mpi.h>
12			#include "brains/mpiSimulation.hpp"
13			#endif // is_mpi
14
15			#include "io/ReadWrite.hpp"
16			#include "utils/simError.h"
17
18			RestraintWriter::RestraintWriter( SimInfo* the_entry_plug ){
19
20			entry_plug = the_entry_plug;
21
22			#ifdef IS_MPI
23			//sort the local atoms by global index
24			sortByGlobalIndex();
25			#endif // is_mpi
26
27			}
28
29			RestraintWriter::~RestraintWriter( ){}
30
31			#ifdef IS_MPI
32
33			/**
34			* A hook function to load balancing
35			*/
36
37			void RestraintWriter::update(){
38			sortByGlobalIndex();
39			}
40
41			/**
42			* Auxiliary sorting function
43			*/
44
45			bool indexSortingCriterion2(const pair<int, int>& p1, const pair<int, int>& p2){
46			return p1.second < p2.second;
47			}
48
49			/**
50			* Sorting the local index by global index
51			*/
52
53			void RestraintWriter::sortByGlobalIndex(){
54			Molecule* mols = entry_plug->molecules;
55			indexArray.clear();
56
57			for(int i = 0; i < entry_plug->n_mol;i++)
58			indexArray.push_back(make_pair(i, mols[i].getGlobalIndex()));
59
60			sort(indexArray.begin(), indexArray.end(), indexSortingCriterion2);
61			}
62
63			#endif
64
65			void RestraintWriter::writeZangle(double currentTime){
66
67			ofstream angleOut;
68			vector<ofstream*> fileStreams;
69
70			#ifdef IS_MPI
71			if(worldRank == 0 ){
72			#endif
73			angleOut.open( entry_plug->zAngleName.c_str(), ios::out \| ios::trunc );
74			if( !angleOut ){
75			sprintf( painCave.errMsg,
76			"Could not open \"%s\" for zAngle output.\n",
77			entry_plug->zAngleName.c_str() );
78			painCave.isFatal = 1;
79			simError();
80			}
81			#ifdef IS_MPI
82			}
83			#endif // is_mpi
84
85			fileStreams.push_back(&angleOut);
86			writeFrame(fileStreams, currentTime);
87
88			#ifdef IS_MPI
89			angleOut.close();
90			#endif
91
92			}
93
94			void RestraintWriter::writeZangle(double currentTime, const char *in_name){
95
96			ofstream finalOut;
97			vector<ofstream*> fileStreams;
98
99			#ifdef IS_MPI
100			if(worldRank == 0 ){
101			#endif
102			finalOut.open( in_name, ios::out \| ios::trunc );
103			if( !finalOut ){
104			sprintf( painCave.errMsg,
105			"Could not open \"%s\" for zAngle output.\n",
106			in_name );
107			painCave.isFatal = 1;
108			simError();
109			}
110			#ifdef IS_MPI
111			}
112			#endif // is_mpi
113
114			fileStreams.push_back(&finalOut);
115			writeFrame(fileStreams, currentTime);
116
117			#ifdef IS_MPI
118			finalOut.close();
119			#endif
120
121			}
122
123			void RestraintWriter::writeFrame( vector<ofstream*>& outFile, double currentTime ){
124
125			const int BUFFERSIZE = 2000;
126			const int MINIBUFFERSIZE = 100;
127
128			char tempBuffer[BUFFERSIZE];
129			char writeLine[BUFFERSIZE];
130
131			int i;
132			unsigned int k;
133
134			#ifdef IS_MPI
135
136			/*********************************************************************
137			* Documentation? You want DOCUMENTATION?
138			*
139			* Why all the potatoes below?
140			*
141			* To make a long story short, the original version of RestraintWriter
142			* worked in the most inefficient way possible. Node 0 would
143			* poke each of the node for an individual atom's formatted data
144			* as node 0 worked its way down the global index. This was particularly
145			* inefficient since the method blocked all processors at every atom
146			* (and did it twice!).
147			*
148			* An intermediate version of RestraintWriter could be described from Node
149			* zero's perspective as follows:
150			*
151			* 1) Have 100 of your friends stand in a circle.
152			* 2) When you say go, have all of them start tossing potatoes at
153			* you (one at a time).
154			* 3) Catch the potatoes.
155			*
156			* It was an improvement, but MPI has buffers and caches that could
157			* best be described in this analogy as "potato nets", so there's no
158			* need to block the processors atom-by-atom.
159			*
160			* This new and improved RestraintWriter works in an even more efficient
161			* way:
162			*
163			* 1) Have 100 of your friend stand in a circle.
164			* 2) When you say go, have them start tossing 5-pound bags of
165			* potatoes at you.
166			* 3) Once you've caught a friend's bag of potatoes,
167			* toss them a spud to let them know they can toss another bag.
168			*
169			* How's THAT for documentation?
170			*
171			*********************************************************************/
172
173			int *potatoes;
174			int myPotato;
175
176			int nProc;
177			int j, which_node, done, which_atom, local_index, currentIndex;
178			double atomData;
179			int isDirectional;
180			char* atomTypeString;
181			char MPIatomTypeString[MINIBUFFERSIZE];
182			int nObjects;
183			int msgLen; // the length of message actually recieved at master nodes
184			#endif //is_mpi
185
186			double angle;
187			DirectionalAtom* dAtom;
188			int nTotObjects;
189			StuntDouble* sd;
190			char* molName;
191			vector<StuntDouble*> integrableObjects;
192			vector<StuntDouble*>::iterator iter;
193			nTotObjects = entry_plug->getTotIntegrableObjects();
194			#ifndef IS_MPI
195
196			for(k = 0; k < outFile.size(); k++)
197			*outFile[k] << currentTime << " : omega values at this time\n";
198
199			for( i=0; i<nTotObjects; i++ ){
200
201			integrableObjects = entry_plug->molecules[i].getIntegrableObjects();
202
203			for( iter = integrableObjects.begin();iter != integrableObjects.end(); ++iter){
204			sd = *iter;
205
206			sprintf( tempBuffer,
207			"%14.10lf\n",
208			sd->getZangle());
209			strcpy( writeLine, tempBuffer );
210
211			for(k = 0; k < outFile.size(); k++)
212			*outFile[k] << writeLine;
213			}
214
215			}
216
217			#else // is_mpi
218
219			/* code to find maximum tag value */
220
221			int *tagub, flag, MAXTAG;
222			MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
223			if (flag) {
224			MAXTAG = *tagub;
225			} else {
226			MAXTAG = 32767;
227			}
228
229			int haveError;
230
231			MPI_Status istatus;
232			int nCurObj;
233			int *MolToProcMap = mpiSim->getMolToProcMap();
234
235			// write out header and node 0's coordinates
236
237			if( worldRank == 0 ){
238
239			// Node 0 needs a list of the magic potatoes for each processor;
240
241			nProc = mpiSim->getNProcessors();
242			potatoes = new int[nProc];
243
244			//write out the comment lines
245			for (i = 0; i < nProc; i++)
246			potatoes[i] = 0;
247
248			for(k = 0; k < outFile.size(); k++)
249			*outFile[k] << currentTime << " fs: omega values at this time\n";
250
251			currentIndex = 0;
252
253			for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
254
255			// Get the Node number which has this atom;
256
257			which_node = MolToProcMap[i];
258
259			if (which_node != 0) {
260
261			if (potatoes[which_node] + 1 >= MAXTAG) {
262			// The potato was going to exceed the maximum value,
263			// so wrap this processor potato back to 0:
264
265			potatoes[which_node] = 0;
266			MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0,
267			MPI_COMM_WORLD);
268
269			}
270
271			myPotato = potatoes[which_node];
272
273			//recieve the number of integrableObject in current molecule
274			MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
275			myPotato, MPI_COMM_WORLD, &istatus);
276			myPotato++;
277
278			for(int l = 0; l < nCurObj; l++){
279
280			if (potatoes[which_node] + 1 >= MAXTAG) {
281			// The potato was going to exceed the maximum value,
282			// so wrap this processor potato back to 0:
283
284			potatoes[which_node] = 0;
285			MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0,
286			MPI_COMM_WORLD);
287
288			}
289
290			MPI_Recv(&atomData, 1, MPI_DOUBLE, which_node, myPotato, MPI_COMM_WORLD, &istatus);
291			myPotato++;
292
293			// If we've survived to here, format the line:
294			sprintf( writeLine,
295			"%14.10lf\n",
296			atomData);
297
298			for(k = 0; k < outFile.size(); k++)
299			*outFile[k] << writeLine;
300
301			}// end for(int l =0)
302			potatoes[which_node] = myPotato;
303
304			}
305			else {
306
307			haveError = 0;
308
309			local_index = indexArray[currentIndex].first;
310
311			integrableObjects = (entry_plug->molecules[local_index]).getIntegrableObjects();
312
313			for(iter= integrableObjects.begin(); iter != integrableObjects.end(); ++iter){
314			sd = *iter;
315			atomData = sd->getZangle();
316
317			// If we've survived to here, format the line:
318
319			sprintf( writeLine,
320			"%14.10lf\n",
321			atomData);
322
323			for(k = 0; k < outFile.size(); k++)
324			*outFile[k] << writeLine;
325
326			}//end for(iter = integrableObject.begin())
327
328			currentIndex++;
329			}
330
331			}//end for(i = 0; i < mpiSim->getNmol())
332
333			for(k = 0; k < outFile.size(); k++)
334			outFile[k]->flush();
335
336			sprintf( checkPointMsg,
337			"Successfully printed a zAngle.\n");
338
339			MPIcheckPoint();
340
341			delete[] potatoes;
342
343			} else {
344
345			// worldRank != 0, so I'm a remote node.
346
347			// Set my magic potato to 0:
348
349			myPotato = 0;
350			currentIndex = 0;
351
352			for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
353
354			// Am I the node which has this integrableObject?
355
356			if (MolToProcMap[i] == worldRank) {
357
358
359			if (myPotato + 1 >= MAXTAG) {
360
361			// The potato was going to exceed the maximum value,
362			// so wrap this processor potato back to 0 (and block until
363			// node 0 says we can go:
364
365			MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
366			}
367
368			local_index = indexArray[currentIndex].first;
369			integrableObjects = entry_plug->molecules[local_index].getIntegrableObjects();
370
371			nCurObj = integrableObjects.size();
372
373			MPI_Send(&nCurObj, 1, MPI_INT, 0,
374			myPotato, MPI_COMM_WORLD);
375			myPotato++;
376
377			for( iter = integrableObjects.begin();
378			iter != integrableObjects.end(); iter++ ){
379
380			if (myPotato + 1 >= MAXTAG) {
381
382			// The potato was going to exceed the maximum value,
383			// so wrap this processor potato back to 0 (and block until
384			// node 0 says we can go:
385
386			MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
387			}
388
389			sd = *iter;
390
391			atomData = sd->getZangle();
392
393			MPI_Send(&atomData, 1, MPI_DOUBLE, 0,
394			myPotato, MPI_COMM_WORLD);
395
396			myPotato++;
397			}
398
399			currentIndex++;
400			}
401			}
402
403			sprintf( checkPointMsg,
404			"Successfully dropped a zAngle.\n");
405			MPIcheckPoint();
406			}
407			#endif // is_mpi
408			}
409