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
#	Content
1	#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