source: proiecte/NBody/Tree codes/treeio.c @ 170

/****************************************************************************/
/* TREEIO.C: I/O routines for hierarchical N-body code. Public routines:    */
/* inputdata(), startoutput(), output(), savestate(), restorestate().       */
/* Copyright (c) 2001 by Joshua E. Barnes, Honolulu, Hawai`i.               */
/****************************************************************************/

#include "stdinc.h"
#include "mathfns.h"
#include "vectmath.h"
#include "getparam.h"
#include "treecode.h"

#include <sys/types.h>
#include <sys/stat.h>
#include <strings.h>

/*
 * Prototypes for local routines.
 */

local void outputdata(void);                    /* write N-body data        */
local void diagnostics(void);                   /* eval N-body diagnostics  */
local void in_int(stream, int *);               /* input integer value      */
local void in_real(stream, real *);             /* input real value         */
local void in_vector(stream, vector);           /* input vector of reals    */
local void out_int(stream, int);                /* output integer value     */
local void out_real(stream, real);              /* output real value        */
local void out_vector(stream, vector);          /* output vector of reals   */

/*
 * Diagnositc output variables.
 */

local real mtot;                                /* total mass of system     */
local real etot[3];                             /* Etot, KE, PE of system   */
local matrix keten;                             /* kinetic energy tensor    */
local matrix peten;                             /* potential energy tensor  */
local vector cmpos;                             /* center of mass position  */
local vector cmvel;                             /* center of mass velocity  */
local vector amvec;                             /* angular momentum vector  */

#define IFMT  "%d "                             /* output format for ints   */
#define RFMT  "%f "                         /* output format for reals  */


/*
 * INPUTDATA: read initial conditions from input file.
 */

void inputdata(void)
{
    stream instr;
    int ndim;
    bodyptr p;

    instr = stropen(infile, "r");               /* open input stream        */
    in_int(instr, &nbody);                      /* read number of bodies    */
    if (nbody < 1)
        error("inputdata: nbody = %d is absurd\n", nbody);
    in_int(instr, &ndim);                       /* read number of dims      */
    if (ndim != NDIM)
        error("inputdata: ndim = %d; expected %d\n", ndim, NDIM);
    in_real(instr, &tnow);                      /* read starting time       */
    bodytab = (bodyptr) allocate(nbody * sizeof(body));
                                                /* allocate body array      */
    for (p = bodytab; p < bodytab+nbody; p++)   /* loop over all bodies     */
        in_real(instr, &Mass(p));               /* read mass of each        */
    for (p = bodytab; p < bodytab+nbody; p++)
        in_vector(instr, Pos(p));               /* read position of each    */
    for (p = bodytab; p < bodytab+nbody; p++)
        in_vector(instr, Vel(p));               /* read velocity of each    */
    fclose(instr);                              /* close input stream       */
    if (scanopt(options, "reset-time"))         /* reset starting time?     */
        tnow = 0.0;                             /* then set it to zero      */
    for (p = bodytab; p < bodytab+nbody; p++)   /* loop over new bodies     */
        Type(p) = BODY;                         /* initialize type field    */
}
void inputdata_2(void)
{
   stream instr;
    int ndim;
    bodyptr p;
   
    
    instr = stropen(infile, "r");               /* open input stream        */
    in_int(instr, &nbody);                      /* read number of bodies    */
    if (nbody < 1)
        error("inputdata: nbody = %d is absurd\n", nbody);
    bodytab = (bodyptr) allocate(nbody * sizeof(body));
                                                /* allocate body array      */
    in_real(instr, &tnow);
    for (p = bodytab; p < bodytab+nbody; p++)   /* loop over all bodies     */
    {
        in_real(instr, &Mass(p));               /* read mass of each        */  
        in_vector(instr, Pos(p));               /* read position of each    */   
        in_vector(instr, Vel(p));               /* read velocity of each    */  
    }
    fclose(instr);                              /* close input stream       */
    if (scanopt(options, "reset-time"))         /* reset starting time?     */
        tnow = 0.0;                             /* then set it to zero      */
    for (p = bodytab; p < bodytab+nbody; p++)   /* loop over new bodies     */
        Type(p) = BODY;                         /* initialize type field    */
}
/*
 * STARTOUTPUT: begin output to log file.
 */

void startoutput(void)
{
    printf("\n%s\n", headline);                 /* print headline, params   */
#if defined(USEFREQ)
    printf("\n%8s%10s%10s", "nbody", "freq", "eps");
#else
    printf("\n%8s%10s%10s", "nbody", "dtime", "eps");
#endif
#if !defined(QUICKSCAN)
    printf("%10s", "theta");
#endif
#if defined(USEFREQ)
    printf("%10s%10s%10s\n", "usequad", "freqout", "tstop");
    printf("%8d%10.2f%10.4f", nbody, freq, eps);
#else
    printf("%10s%10s%10s\n", "usequad", "dtout", "tstop");
    printf("%8d%10.5f%10.4f", nbody, dtime, eps);
#endif
#if !defined(QUICKSCAN)
    printf("%10.2f", theta);
#endif
#if defined(USEFREQ)
    printf("%10s%10.2f%10.4f\n", usequad ? "true" : "false", freqout, tstop);
#else
    printf("%10s%10.5f%10.4f\n", usequad ? "true" : "false", dtout, tstop);
#endif
    if (! strnull(options))                     /* print options, if any    */
        printf("\n\toptions: %s\n", options);
    if (! strnull(savefile))                    /* was state file given?    */
        savestate(savefile);                    /* save initial data        */
}

/*
 * FORCEREPORT: print staristics on tree construction and force calculation.
 */

void forcereport(void)
{
    printf("\n\t%8s%8s%8s%8s%10s%10s%8s\n",
           "rsize", "tdepth", "ftree",
           "actmax", "nbbtot", "nbctot", "CPUfc");
    printf("\t%8.1f%8d%8.3f%8d%10d%10d%8.3f\n",
           rsize, tdepth, (nbody + ncell - 1) / ((real) ncell),
           actmax, nbbcalc, nbccalc, cpuforce);
}

/*
 * OUTPUT: compute diagnostics and output body data.
 */

void output(void)
{
    real cmabs, amabs, teff;

    diagnostics();                              /* compute std diagnostics  */
    ABSV(cmabs, cmvel);                         /* find magnitude of cm vel */
    ABSV(amabs, amvec);                         /* find magnitude of J vect */
    printf("\n    %8s%8s%8s%8s%8s%8s%8s%8s\n",
           "time", "|T+U|", "T", "-U", "-T/U", "|Vcom|", "|Jtot|", "CPUtot");
    printf("    %8.3f%8.5f%8.5f%8.5f%8.5f%8.5f%8.5f%8.3f\n",
           tnow, ABS(etot[0]), etot[1], -etot[2], -etot[1]/etot[2],
           cmabs, amabs, cputime());
#if defined(USEFREQ)
    teff = tnow + (freq > 0 ? 0.125/freq : 0);  /* anticipate slightly...   */
#else
    teff = tnow + dtime/8;                      /* anticipate slightly...   */
#endif
    if (! strnull(outfile) && teff >= tout)     /* time for data output?    */
        //outputdata();
        {}
    if (! strnull(savefile))                    /* was state file given?    */
        savestate(savefile);                    /* save data for restart    */
}

/*
 * OUTPUTDATA: output body data.
 */

/*void outputdata(void)
{
    char namebuf[256];
    struct stat buf;
    stream outstr;
    bodyptr p;

    sprintf(namebuf, outfile, nstep);           // construct output name    
    if (stat(namebuf, &buf) != 0)               // no output file exists?   
        outstr = stropen(namebuf, "w");         // create & open for output 
    else                                        // else file already exists 
        outstr = stropen(namebuf, "a");         // reopen and append output 
    out_int(outstr, nbody);                     // write number of bodies   
    out_int(outstr, NDIM);                      // number of dimensions     
    out_real(outstr, tnow);                     // and current time value   
    for (p = bodytab; p < bodytab+nbody; p++)   // loop over all bodies     
        out_real(outstr, Mass(p));              // output mass of each      
    for (p = bodytab; p < bodytab+nbody; p++)
        out_vector(outstr, Pos(p));             // output positions         
    for (p = bodytab; p < bodytab+nbody; p++)
        out_vector(outstr, Vel(p));             // output velocities        
    if (scanopt(options, "out-phi"))            // potentials requested?    
        for (p = bodytab; p < bodytab+nbody; p++)
            out_real(outstr, Phi(p));           // output potentials        
    if (scanopt(options, "out-acc"))            // accelerations requested? 
        for (p = bodytab; p < bodytab+nbody; p++)
            out_vector(outstr, Acc(p));         // output accelerations     
    fclose(outstr);                             // close up output file     
    printf("\n\tdata output to file %s at time %f\n", namebuf, tnow);
#if defined(USEFREQ)
    tout += 1.0 / freqout;                      // schedule next output     
#else
    tout += dtout;                              // schedule next output     
#endif
}
*/

void outputdata(void)
{
    char namebuf[256];
    struct stat buf;
    stream outstr;
    bodyptr p;
    int tnow_int;

    sprintf(namebuf, outfile, nstep);           // construct output name    
    if (stat(namebuf, &buf) != 0)               // no output file exists?   
        outstr = stropen(namebuf, "w");         // create & open for output 
    else                                        // else file already exists 
        outstr = stropen(namebuf, "a");         // reopen and append output 
    tnow_int = tnow;
    out_int(outstr, nbody);                     // write number of bodies         
    out_real(outstr, tnow);                        
    //out_int(outstr, tnow_int);                        // and current time value
    for (p = bodytab; p < bodytab+nbody; p++)   // loop over all bodies     
    {
      //out_real(outstr, Mass(p));              // output mass of each   
      fprintf(outstr, RFMT, Mass(p));
      out_vector(outstr, Pos(p));             // output positions 
      out_vector(outstr, Vel(p));             // output velocities
      fprintf(outstr,"\n");
    }  
    //fprintf(outstr,"\n");
    fclose(outstr);                             // close up output file     
    printf("\n\tdata output to file %s at time %f\n", namebuf, tnow);
#if defined(USEFREQ)
    tout += 1.0 / freqout;                      // schedule next output     
#else
    tout += dtout;                              // schedule next output     
#endif
}

/*
 * DIAGNOSTICS: compute set of dynamical diagnostics.
 */

local void diagnostics(void)
{
    register bodyptr p;
    real velsq;
    vector tmpv;
    matrix tmpt;

    mtot = 0.0;                                 /* zero total mass          */
    etot[1] = etot[2] = 0.0;                    /* zero total KE and PE     */
    CLRM(keten);                                /* zero ke tensor           */
    CLRM(peten);                                /* zero pe tensor           */
    CLRV(amvec);                                /* zero am vector           */
    CLRV(cmpos);                                /* zero c. of m. position   */
    CLRV(cmvel);                                /* zero c. of m. velocity   */
    for (p = bodytab; p < bodytab+nbody; p++) { /* loop over all particles  */
        mtot += Mass(p);                        /* sum particle masses      */
        DOTVP(velsq, Vel(p), Vel(p));           /* square vel vector        */
        etot[1] += 0.5 * Mass(p) * velsq;       /* sum current KE           */
        etot[2] += 0.5 * Mass(p) * Phi(p);      /* and current PE           */
        MULVS(tmpv, Vel(p), 0.5 * Mass(p));     /* sum 0.5 m v_i v_j        */
        OUTVP(tmpt, tmpv, Vel(p));
        ADDM(keten, keten, tmpt);
        MULVS(tmpv, Pos(p), Mass(p));           /* sum m r_i a_j            */
        OUTVP(tmpt, tmpv, Acc(p));
        ADDM(peten, peten, tmpt);
        CROSSVP(tmpv, Vel(p), Pos(p));          /* sum angular momentum     */
        MULVS(tmpv, tmpv, Mass(p));
        ADDV(amvec, amvec, tmpv);
        MULVS(tmpv, Pos(p), Mass(p));           /* sum cm position          */
        ADDV(cmpos, cmpos, tmpv);
        MULVS(tmpv, Vel(p), Mass(p));           /* sum cm momentum          */
        ADDV(cmvel, cmvel, tmpv);
    }
    etot[0] = etot[1] + etot[2];                /* sum KE and PE            */
    DIVVS(cmpos, cmpos, mtot);                  /* normalize cm coords      */
    DIVVS(cmvel, cmvel, mtot);
}

/*
 * IN_INT, IN_REAL, IN_VECTOR: low level input routines.
 */

local void in_int(stream str, int *iptr)
{
#if !defined(BINARYIO)
    if (fscanf(str, "%d", iptr) != 1)
        error("in_int: input conversion error\n");
#else
    if (fread((void *) iptr, sizeof(int), 1, str) != 1)
        error("in_int: fread failed\n");
#endif
}

local void in_real(stream str, real *rptr)
{
    double tmp;

#if !defined(BINARYIO)
    if (fscanf(str, "%lf", &tmp) != 1)
        error("in_real: input conversion error\n");
    *rptr = tmp;
#else
    if (fread((void *) rptr, sizeof(real), 1, str) != 1)
        error("in_real: fread failed\n");
#endif
}

local void in_vector(stream str, vector vec)
{
    double tmpx, tmpy, tmpz;

#if !defined(BINARYIO)
    if (fscanf(str, "%lf%lf%lf", &tmpx, &tmpy, &tmpz) != 3)
        error("in_vector: input conversion error\n");
    vec[0] = tmpx;
    vec[1] = tmpy;
    vec[2] = tmpz;
#else
    if (fread((void *) vec, sizeof(real), NDIM, str) != NDIM)
        error("in_vector: fread failed\n");
#endif
}

/*
 * OUT_INT, OUT_REAL, OUT_VECTOR: low level output routines. */


local void out_int(stream str, int ival)
{
#if !defined(BINARYIO)
    if (fprintf(str, IFMT "\n", ival) < 0)
        error("out_int: fprintf failed\n");
#else
    if (fwrite((void *) &ival, sizeof(int), 1, str) != 1)
        error("out_int: fwrite failed\n");
#endif
}

local void out_real(stream str, real rval)
{
#if !defined(BINARYIO)
    if (fprintf(str, RFMT "\n", rval) < 0)
        error("out_real: fprintf failed\n");
#else
    if (fwrite((void *) &rval, sizeof(real), 1, str) != 1)
        error("out_real: fwrite failed\n");
#endif
}

local void out_vector(stream str, vector vec)
{
#if !defined(BINARYIO)
    if (fprintf(str, RFMT RFMT RFMT, vec[0], vec[1], vec[2]) < 0)
        error("out_vector: fprintf failed\n");
#else
    if (fwrite((void *) vec, sizeof(real), NDIM, str) != NDIM)
        error("out_vector: fwrite failed\n");
#endif
}

/*
 * SAVESTATE: write current state to disk file.
 */

#define safewrite(ptr,len,str)                  \
    if (fwrite((void *) ptr, len, 1, str) != 1) \
        error("savestate: fwrite failed\n")

void savestate(string pattern)
{
    char namebuf[256];
    stream str;
    int nchars;

    sprintf(namebuf, pattern, nstep & 1);       /* construct alternate name */
    str = stropen(namebuf, "w!");
    nchars = strlen(getargv0()) + 1;
    safewrite(&nchars, sizeof(int), str);
    safewrite(getargv0(), nchars * sizeof(char), str);
    nchars = strlen(getversion()) + 1;
    safewrite(&nchars, sizeof(int), str);
    safewrite(getversion(), nchars * sizeof(char), str);
#if defined(USEFREQ)
    safewrite(&freq, sizeof(real), str);
#else
    safewrite(&dtime, sizeof(real), str);
#endif
#if !defined(QUICKSCAN)
    safewrite(&theta, sizeof(real), str);
#endif
    safewrite(&usequad, sizeof(bool), str);
    safewrite(&eps, sizeof(real), str);
    nchars = strlen(options) + 1;
    safewrite(&nchars, sizeof(int), str);
    safewrite(options, nchars * sizeof(char), str);
    safewrite(&tstop, sizeof(real), str);
#if defined(USEFREQ)
    safewrite(&freqout, sizeof(real), str);
#else
    safewrite(&dtout, sizeof(real), str);
#endif
    safewrite(&tnow, sizeof(real), str);
    safewrite(&tout, sizeof(real), str);
    safewrite(&nstep, sizeof(int), str);
    safewrite(&rsize, sizeof(real), str);
    safewrite(&nbody, sizeof(int), str);
    safewrite(bodytab, nbody * sizeof(body), str);
    fclose(str);
}

/*
 * RESTORESTATE: restore state from disk file.
 */

#define saferead(ptr,len,str)                  \
    if (fread((void *) ptr, len, 1, str) != 1) \
        error("restorestate: fread failed\n")

void restorestate(string file)
{
    stream str;
    int nchars;
    string program, version;

    str = stropen(file, "r");
    saferead(&nchars, sizeof(int), str);
    program = (string) allocate(nchars * sizeof(char));
    saferead(program, nchars * sizeof(char), str);
    saferead(&nchars, sizeof(int), str);
    version = (string) allocate(nchars * sizeof(char));
    saferead(version, nchars * sizeof(char), str);
    if (! streq(program, getargv0()) ||         /* check program, version   */
          ! streq(version, getversion()))
        printf("warning: state file may be outdated\n\n");
#if defined(USEFREQ)
    saferead(&freq, sizeof(real), str);
#else
    saferead(&dtime, sizeof(real), str);
#endif
#if !defined(QUICKSCAN)
    saferead(&theta, sizeof(real), str);
#endif
    saferead(&usequad, sizeof(bool), str);
    saferead(&eps, sizeof(real), str);
    saferead(&nchars, sizeof(int), str);
    options = (string) allocate(nchars * sizeof(char));
    saferead(options, nchars * sizeof(char), str);
    saferead(&tstop, sizeof(real), str);
#if defined(USEFREQ)
    saferead(&freqout, sizeof(real), str);
#else
    saferead(&dtout, sizeof(real), str);
#endif
    saferead(&tnow, sizeof(real), str);
    saferead(&tout, sizeof(real), str);
    saferead(&nstep, sizeof(int), str);
    saferead(&rsize, sizeof(real), str);
    saferead(&nbody, sizeof(int), str);
    bodytab = (bodyptr) allocate(nbody * sizeof(body));
    saferead(bodytab, nbody * sizeof(body), str);
    fclose(str);
}