src/io/DumpWriter.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"

namespace dWrite {
    void DieDieDie(void);

}

using namespace dWrite;

#endif //is_mpi

#include "io/ReadWrite.hpp"
#include "utils/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.c_str(), ios::out | ios::trunc);

        if (!dumpFile) {
            sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
                    entry_plug->sampleName.c_str());
            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() {
    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(), indexSortingCriterion);
}

#endif

void DumpWriter::writeDump(double currentTime) {
    ofstream finalOut;
    vector<ofstream *>fileStreams;

#ifdef IS_MPI

    if (worldRank == 0) {
#endif

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

        if (!finalOut) {
            sprintf(painCave.errMsg,
                    "Could not open \"%s\" for final dump output.\n",
                    entry_plug->finalName.c_str());
            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.c_str(), ios::out | ios::trunc);

        if (!finalOut) {
            sprintf(painCave.errMsg,
                    "Could not open \"%s\" for final dump output.\n",
                    entry_plug->finalName.c_str());
            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;
    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 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;
    int which_node;
    int done;
    int which_atom;
    int local_index;
    int currentIndex;
    double atomData[13];
    int isDirectional;
    char * atomTypeString;
    char MPIatomTypeString[MINIBUFFERSIZE];
    int nObjects;
    int msgLen; // the length of message actually recieved at master nodes

#endif          //is_mpi

    Quat4d q;
    Vector3d ji;
    DirectionalAtom * dAtom;
    Vector3d pos;
    Vector3d vel;

    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] << nTotObjects << "\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
        //another circular reference nightmare
        *outFile[k] << entry_plug->the_integrator->getAdditionalParameters()
             << endl;
    }

    for(i = 0; i < entry_plug->n_mol; i++) {
        integrableObjects = entry_plug->molecules[i].getIntegrableObjects();
        molName
            = (entry_plug->compStamps[entry_plug->molecules[i].getStampID()])->getID();

        for(iter = integrableObjects.begin();
            iter != integrableObjects.end(); ++iter) {
            sd = *iter;
            pos = sd->getPos();
            vel = sd->getVel();

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

            strcpy(writeLine, tempBuffer);

            if (sd->isDirectional()) {
                q = sd->getQ();
                ji = sd->getJ();

                sprintf(tempBuffer, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", q[0],
                        q[1], q[2], q[3],
                        ji[0], ji[1], ji[2]);
                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 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] << nTotObjects << "\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->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] + 2 >= 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(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR,
                             which_node, myPotato, MPI_COMM_WORLD,
                             &istatus);

                    atomTypeString = MPIatomTypeString;

                    myPotato++;

                    MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato,
                             MPI_COMM_WORLD, &istatus);
                    myPotato++;

                    MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);

                    if (msgLen == 13)
                        isDirectional = 1;
                    else
                        isDirectional = 0;

                    // 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, atomData[0],
                                atomData[1], atomData[2],
                                atomData[3], atomData[4],
                                atomData[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,
                                atomData[0],
                                atomData[1],
                                atomData[2],
                                atomData[3],
                                atomData[4],
                                atomData[5],
                                atomData[6],
                                atomData[7],
                                atomData[8],
                                atomData[9],
                                atomData[10],
                                atomData[11],
                                atomData[12]);
                    }

                    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;
                    atomTypeString = sd->getType();

                    pos = sd->getPos();
                    vel = sd->getVel();

                    atomData[0] = pos[0];
                    atomData[1] = pos[1];
                    atomData[2] = pos[2];

                    atomData[3] = vel[0];
                    atomData[4] = vel[1];
                    atomData[5] = vel[2];

                    isDirectional = 0;

                    if (sd->isDirectional()) {
                        isDirectional = 1;

                        q = sd->getQ();
                        ji = sd->getJ();

                        for(int j = 0; j < 6; j++) {
                            atomData[j] = atomData[j];
                        }

                        atomData[6] = q[0];
                        atomData[7] = q[1];
                        atomData[8] = q[2];
                        atomData[9] = q[3];

                        atomData[10] = ji[0];
                        atomData[11] = ji[1];
                        atomData[12] = ji[2];
                    }

                    // 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, atomData[0],
                                atomData[1], atomData[2],
                                atomData[3], atomData[4],
                                atomData[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,
                                atomData[0],
                                atomData[1],
                                atomData[2],
                                atomData[3],
                                atomData[4],
                                atomData[5],
                                atomData[6],
                                atomData[7],
                                atomData[8],
                                atomData[9],
                                atomData[10],
                                atomData[11],
                                atomData[12]);
                    }

                    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, "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->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 + 2 >= 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;

                    atomTypeString = sd->getType();

                    pos = sd->getPos();
                    vel = sd->getVel();

                    atomData[0] = pos[0];
                    atomData[1] = pos[1];
                    atomData[2] = pos[2];

                    atomData[3] = vel[0];
                    atomData[4] = vel[1];
                    atomData[5] = vel[2];

                    isDirectional = 0;

                    if (sd->isDirectional()) {
                        isDirectional = 1;

                        q = sd->getQ();
                        ji = sd->getJ();

                        atomData[6] = q[0];
                        atomData[7] = q[1];
                        atomData[8] = q[2];
                        atomData[9] = q[3];

                        atomData[10] = ji[0];
                        atomData[11] = ji[1];
                        atomData[12] = ji[2];
                    }

                    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++;

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