1 /*************************************************************************/
3 /* Copyright (c) 1994 Stanford University */
5 /* All rights reserved. */
7 /* Permission is given to use, copy, and modify this software for any */
8 /* non-commercial purpose as long as this copyright notice is not */
9 /* removed. All other uses, including redistribution in whole or in */
10 /* part, are forbidden without prior written permission. */
12 /* This software is provided with absolutely no warranty and no */
15 /*************************************************************************/
25 void in_int (), in_real (), in_vector ();
26 void out_int (), out_real (), out_vector ();
27 void diagnostics (unsigned int ProcessId);
30 * INPUTDATA: read initial conditions from input file.
36 permanent char headbuf[128];
42 fprintf(stderr,"reading input file : %s\n",infile);
44 instr = fopen(infile, "r");
46 error("inputdata: cannot find file %s\n", infile);
47 sprintf(headbuf, "Hack code: input file %s\n", infile);
49 in_int(instr, &nbody);
51 error("inputdata: nbody = %d is absurd\n", nbody);
54 error("inputdata: NDIM = %d ndim = %d is absurd\n", NDIM,ndim);
55 in_real(instr, &tnow);
56 for (i = 0; i < MAX_PROC; i++) {
59 bodytab = (bodyptr) G_MALLOC(nbody * sizeof(body));
61 error("inputdata: not enuf memory\n");
62 for (p = bodytab; p < bodytab+nbody; p++) {
68 for (p = bodytab; p < bodytab+nbody; p++)
69 in_real(instr, &Mass(p));
70 for (p = bodytab; p < bodytab+nbody; p++)
71 in_vector(instr, Pos(p));
72 for (p = bodytab; p < bodytab+nbody; p++)
73 in_vector(instr, Vel(p));
78 * INITOUTPUT: initialize output routines.
84 printf("\n\t\t%s\n\n", headline);
85 printf("%10s%10s%10s%10s%10s%10s%10s%10s\n",
86 "nbody", "dtime", "eps", "tol", "dtout", "tstop","fcells","NPROC");
87 printf("%10d%10.5f%10.4f%10.2f%10.3f%10.3f%10.2f%10d\n\n",
88 nbody, dtime, eps, tol, dtout, tstop, fcells, NPROC);
92 * STOPOUTPUT: finish up after a run.
97 * OUTPUT: compute diagnostics and output data.
102 unsigned int ProcessId;
104 int nttot, nbavg, ncavg,k;
107 vector tempv1,tempv2;
109 if ((Local[ProcessId].tout - 0.01 * dtime) <= Local[ProcessId].tnow) {
110 Local[ProcessId].tout += dtout;
113 diagnostics(ProcessId);
115 if (Local[ProcessId].mymtot!=0) {
116 LOCK(Global->CountLock);
117 Global->n2bcalc += Local[ProcessId].myn2bcalc;
118 Global->nbccalc += Local[ProcessId].mynbccalc;
119 Global->selfint += Local[ProcessId].myselfint;
120 ADDM(Global->keten, Global-> keten, Local[ProcessId].myketen);
121 ADDM(Global->peten, Global-> peten, Local[ProcessId].mypeten);
122 for (k=0;k<3;k++) Global->etot[k] += Local[ProcessId].myetot[k];
123 ADDV(Global->amvec, Global-> amvec, Local[ProcessId].myamvec);
125 MULVS(tempv1, Global->cmphase[0],Global->mtot);
126 MULVS(tempv2, Local[ProcessId].mycmphase[0], Local[ProcessId].mymtot);
127 ADDV(tempv1, tempv1, tempv2);
128 DIVVS(Global->cmphase[0], tempv1, Global->mtot+Local[ProcessId].mymtot);
130 MULVS(tempv1, Global->cmphase[1],Global->mtot);
131 MULVS(tempv2, Local[ProcessId].mycmphase[1], Local[ProcessId].mymtot);
132 ADDV(tempv1, tempv1, tempv2);
133 DIVVS(Global->cmphase[1], tempv1, Global->mtot+Local[ProcessId].mymtot);
134 Global->mtot +=Local[ProcessId].mymtot;
135 UNLOCK(Global->CountLock);
138 BARRIER(Global->Baraccel,NPROC);
141 nttot = Global->n2bcalc + Global->nbccalc;
142 nbavg = (int) ((real) Global->n2bcalc / (real) nbody);
143 ncavg = (int) ((real) Global->nbccalc / (real) nbody);
150 * DIAGNOSTICS: compute set of dynamical diagnostics.
154 diagnostics (ProcessId)
155 unsigned int ProcessId;
157 register bodyptr p,*pp;
162 Local[ProcessId].mymtot = 0.0;
163 Local[ProcessId].myetot[1] = Local[ProcessId].myetot[2] = 0.0;
164 CLRM(Local[ProcessId].myketen);
165 CLRM(Local[ProcessId].mypeten);
166 CLRV(Local[ProcessId].mycmphase[0]);
167 CLRV(Local[ProcessId].mycmphase[1]);
168 CLRV(Local[ProcessId].myamvec);
169 for (pp = Local[ProcessId].mybodytab+Local[ProcessId].mynbody -1;
170 pp >= Local[ProcessId].mybodytab; pp--) {
172 Local[ProcessId].mymtot += Mass(p);
173 DOTVP(velsq, Vel(p), Vel(p));
174 Local[ProcessId].myetot[1] += 0.5 * Mass(p) * velsq;
175 Local[ProcessId].myetot[2] += 0.5 * Mass(p) * Phi(p);
176 MULVS(tmpv, Vel(p), 0.5 * Mass(p));
177 OUTVP(tmpt, tmpv, Vel(p));
178 ADDM(Local[ProcessId].myketen, Local[ProcessId].myketen, tmpt);
179 MULVS(tmpv, Pos(p), Mass(p));
180 OUTVP(tmpt, tmpv, Acc(p));
181 ADDM(Local[ProcessId].mypeten, Local[ProcessId].mypeten, tmpt);
182 MULVS(tmpv, Pos(p), Mass(p));
183 ADDV(Local[ProcessId].mycmphase[0], Local[ProcessId].mycmphase[0], tmpv);
184 MULVS(tmpv, Vel(p), Mass(p));
185 ADDV(Local[ProcessId].mycmphase[1], Local[ProcessId].mycmphase[1], tmpv);
186 CROSSVP(tmpv, Pos(p), Vel(p));
187 MULVS(tmpv, tmpv, Mass(p));
188 ADDV(Local[ProcessId].myamvec, Local[ProcessId].myamvec, tmpv);
190 Local[ProcessId].myetot[0] = Local[ProcessId].myetot[1]
191 + Local[ProcessId].myetot[2];
192 if (Local[ProcessId].mymtot!=0){
193 DIVVS(Local[ProcessId].mycmphase[0], Local[ProcessId].mycmphase[0],
194 Local[ProcessId].mymtot);
195 DIVVS(Local[ProcessId].mycmphase[1], Local[ProcessId].mycmphase[1],
196 Local[ProcessId].mymtot);
203 * Low-level input and output operations.
206 void in_int(str, iptr)
210 if (fscanf(str, "%d", iptr) != 1)
211 error("in_int: input conversion error\n");
214 void in_real(str, rptr)
220 if (fscanf(str, "%lf", &tmp) != 1)
221 error("in_real: input conversion error\n");
225 void in_vector(str, vec)
229 double tmpx, tmpy, tmpz;
231 if (fscanf(str, "%lf%lf%lf", &tmpx, &tmpy, &tmpz) != 3)
232 error("in_vector: input conversion error\n");
233 vec[0] = tmpx; vec[1] = tmpy; vec[2] = tmpz;
236 void out_int(str, ival)
240 fprintf(str, " %d\n", ival);
243 void out_real(str, rval)
247 fprintf(str, " %21.14E\n", rval);
250 void out_vector(str, vec)
254 fprintf(str, " %21.14E %21.14E", vec[0], vec[1]);
255 fprintf(str, " %21.14E\n",vec[2]);