examples/charm++/barnes-charm/src/code_io.C

   1 /*************************************************************************/
   2 /*                                                                       */
   3 /*  Copyright (c) 1994 Stanford University                               */
   4 /*                                                                       */
   5 /*  All rights reserved.                                                 */
   6 /*                                                                       */
   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.                */
  11 /*                                                                       */
  12 /*  This software is provided with absolutely no warranty and no         */
  13 /*  support.                                                             */
  14 /*                                                                       */
  15 /*************************************************************************/
  16
  17 /*
  18  * CODE_IO.C:
  19  */
  20 EXTERN_ENV
  21 #define global extern
  22
  23 #include "code.h"
  24
  25 void in_int (), in_real (), in_vector ();
  26 void out_int (), out_real (), out_vector ();
  27 void diagnostics (unsigned int ProcessId);
  28
  29 /*
  30  * INPUTDATA: read initial conditions from input file.
  31  */
  32
  33 inputdata ()
  34 {
  35    stream instr;
  36    permanent char headbuf[128];
  37    int ndim,counter=0;
  38    real tnow;
  39    bodyptr p;
  40    int i;
  41
  42    fprintf(stderr,"reading input file : %s\n",infile);
  43    fflush(stderr);
  44    instr = fopen(infile, "r");
  45    if (instr == NULL)
  46       error("inputdata: cannot find file %s\n", infile);
  47    sprintf(headbuf, "Hack code: input file %s\n", infile);
  48    headline = headbuf;
  49    in_int(instr, &nbody);
  50    if (nbody < 1)
  51       error("inputdata: nbody = %d is absurd\n", nbody);
  52    in_int(instr, &ndim);
  53    if (ndim != NDIM)
  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++) {
  57       Local[i].tnow = tnow;
  58    }
  59    bodytab = (bodyptr) G_MALLOC(nbody * sizeof(body));
  60    if (bodytab == NULL)
  61       error("inputdata: not enuf memory\n");
  62    for (p = bodytab; p < bodytab+nbody; p++) {
  63       Type(p) = BODY;
  64       Cost(p) = 1;
  65       Phi(p) = 0.0;
  66       CLRV(Acc(p));
  67    }
  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));
  74    fclose(instr);
  75 }
  76 \f
  77 /*
  78  * INITOUTPUT: initialize output routines.
  79  */
  80
  81
  82 initoutput()
  83 {
  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);
  89 }
  90
  91 /*
  92  * STOPOUTPUT: finish up after a run.
  93  */
  94
  95 \f
  96 /*
  97  * OUTPUT: compute diagnostics and output data.
  98  */
  99
 100 void
 101 output (ProcessId)
 102    unsigned int ProcessId;
 103 {
 104    int nttot, nbavg, ncavg,k;
 105    double cputime();
 106    bodyptr p, *pp;
 107    vector tempv1,tempv2;
 108
 109    if ((Local[ProcessId].tout - 0.01 * dtime) <= Local[ProcessId].tnow) {
 110       Local[ProcessId].tout += dtout;
 111    }
 112
 113    diagnostics(ProcessId);
 114
 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);
 124
 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);
 129
 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);
 136    }
 137
 138    BARRIER(Global->Baraccel,NPROC);
 139
 140    if (ProcessId==0) {
 141       nttot = Global->n2bcalc + Global->nbccalc;
 142       nbavg = (int) ((real) Global->n2bcalc / (real) nbody);
 143       ncavg = (int) ((real) Global->nbccalc / (real) nbody);
 144    }
 145 }
 146
 147
 148 \f
 149 /*
 150  * DIAGNOSTICS: compute set of dynamical diagnostics.
 151  */
 152
 153 void
 154 diagnostics (ProcessId)
 155    unsigned int ProcessId;
 156 {
 157    register bodyptr p,*pp;
 158    real velsq;
 159    vector tmpv;
 160    matrix tmpt;
 161
 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--) {
 171       p= *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);
 189    }
 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);
 197    }
 198 }
 199
 200 \f
 201
 202 /*
 203  * Low-level input and output operations.
 204  */
 205
 206 void in_int(str, iptr)
 207   stream str;
 208   int *iptr;
 209 {
 210    if (fscanf(str, "%d", iptr) != 1)
 211       error("in_int: input conversion error\n");
 212 }
 213
 214 void in_real(str, rptr)
 215   stream str;
 216   real *rptr;
 217 {
 218    double tmp;
 219
 220    if (fscanf(str, "%lf", &tmp) != 1)
 221       error("in_real: input conversion error\n");
 222    *rptr = tmp;
 223 }
 224
 225 void in_vector(str, vec)
 226   stream str;
 227   vector vec;
 228 {
 229    double tmpx, tmpy, tmpz;
 230
 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;
 234 }
 235
 236 void out_int(str, ival)
 237   stream str;
 238   int ival;
 239 {
 240    fprintf(str, "  %d\n", ival);
 241 }
 242
 243 void out_real(str, rval)
 244   stream str;
 245   real rval;
 246 {
 247    fprintf(str, " %21.14E\n", rval);
 248 }
 249
 250 void out_vector(str, vec)
 251   stream str;
 252   vector vec;
 253 {
 254    fprintf(str, " %21.14E %21.14E", vec[0], vec[1]);
 255    fprintf(str, " %21.14E\n",vec[2]);
 256 }