root/trunk/tsuite/programs/mpi/mpi.cpp

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*main program to call cloudy to generate a large loc data cube */
#include "cddefines.h"
#include "cddrive.h"
#include "version.h"

/* This flag indicates where we are multi-processing.\
 * set with -DMPI on compiler command line
 * if this is true then all output happens at the end,
 * if false then each set of results is printed and flushed when it happens.
 * must be set true on parallel machines, false on pc */
#ifdef MPI 
#       include <mpi.h>
#endif
/*lint -e506  constant value Boolean */


#ifdef _MSC_VER
        /* disable conditional expression is constant */
#       pragma warning( disable : 4127 )
#endif

/* ========================================================================= */
/*
Main output, file ending in .lin 
format, tab separated fields
Cloudy exit condition for that model
number of warnings,
elapsed time, 
hden
flux
par1
par2
*/

#define TITLE "std con Fe2 sol, col 24 " 

/* this is name of resulting output files, FILENAME.lin will be main results */
#define FILENAME "STD_con_sol" 

/* path to directory where we will write, must end in / or \ for current system
 * final location will be PATHNAME + FILENAME */
/*#define PATHNAME "c:\\projects\\cloudy\\tests\\"*/
#define PATHNAME ""

/* this is an additive offset for hden and flux, usually 0 */
#define OFFSET_HDEN (0.)
#define OFFSET_FLUX (0.)

/* these following pairs make 8x8 square grid */
/* lower and upper bounds on hydrogen density - log cm^-3
 * for BLR work usually 7 and 14 */
/*#define HDENINIT (12.5+OFFSET_HDEN)*/
/*#define HDENLIMIT (12.5+OFFSET_HDEN)*/
#define HDENINIT (7.+OFFSET_HDEN)
#define HDENLIMIT (14.+OFFSET_HDEN)

/* low and up bounds on flux of hydrogen ionizing photons - log cm^-2 s-1 
 * for BLR work usually 17 and 24 */
#define FLUXINIT (17.+OFFSET_FLUX) 
#define FLUXLIMIT (24.+OFFSET_FLUX)
/*#define FLUXINIT (20.5+OFFSET_FLUX)
#define FLUXLIMIT (20.5+OFFSET_FLUX)*/

/* the increment on both axis */
/*#define INCREMENT 0.5*/
#define INCREMENT 1./**/
/*#define INCREMENT 2.*/
/*#define INCREMENT 0.05*/

/* this is the number of iterations to do - <0 means to convergence, 1 will do 2 iterations,
 * 0 will do 1  */
#define ITERATIONS  2  
/*#define ITERATIONS  -1  */

/* the column density, set to false if not used */
#define COLDEN 24

/* this says whether or not to include the FeII atom */
#define FEII true

/* quick test - 3 zones, constant temperature, 2 iterations */
#define QUICK_MODEL false

/* if true the only call cdNoExec */
#define VERY_QUICK_MODEL false

/* says whether to make the large cloudy output, if true then output created,
 * if false then not output, only the emission line intensities. */
#define DO_PRINT false

/* print last iteration */
#define PRINT_LAST false

/* set no buffering on output - if true then print as it happens*/
#define NO_BUFFERING false

/* ============================================================================ */

/* this structure will hold results for final printout */
#define NLINTOT 1500 /* the maximum number of lines */

typedef struct
{
        /* save these parameters for each grid point,
         * density, flux, and full spectrum */
        double hden , flux , pred[NLINTOT] ;
        /* number of cautions and warnings for this calc, all should be zero */
        int nWarnings , nCautions , nTFail;
        /* flag returned by Cloudy, 0 == OK, 1 for failure */
        int lgBadEnd ;
        /* the execution time for this model */
        double etime;
        /* two extra parameters */
        double par1 , par2;
} GRIDTAG;

/* this is where we will save results */
#ifdef MPI
void Build_derived_type(GRIDTAG gridtag, MPI_Datatype* message_type_ptr)
{
  int block_lengths[10];
  MPI_Aint displacements[10];
  MPI_Aint addresses[11];
  MPI_Datatype typelist[10];

  /*First specify the types */
  typelist[0] = MPI_DOUBLE;    /*hden*/
  typelist[1] = MPI_DOUBLE;    /*flux*/
  typelist[2] = MPI_DOUBLE;    /*pred*/
  typelist[3] = MPI_INT;       /*nWarnings*/
  typelist[4] = MPI_INT;       /*nCautions*/
  typelist[5] = MPI_INT;       /*nTFail*/
  typelist[6] = MPI_INT;       /*lgBadEnd*/
  typelist[7] = MPI_DOUBLE;    /*etime*/
  typelist[8] = MPI_DOUBLE;    /*par1*/
  typelist[9] = MPI_DOUBLE;    /*par2*/

  /* Specify the number of elements of each type */
  block_lengths[0] = block_lengths[1] = block_lengths[3] = 1;
  block_lengths[4] = block_lengths[5] = block_lengths[6] = 1;
  block_lengths[7] = 1;
  /* this one is special since it is NLINTOT long */
  block_lengths[2] = NLINTOT;
  /* two new parameters are each one element long */
  block_lengths[8] = block_lengths[9] = 1;

  /* Calculate the displacement of the members relative to indata */
  MPI_Address( &gridtag,            &addresses[0]);
  MPI_Address(&(gridtag.hden),      &addresses[1]);
  MPI_Address(&(gridtag.flux),      &addresses[2]);
  MPI_Address(&(gridtag.pred),      &addresses[3]);
  MPI_Address(&(gridtag.nWarnings), &addresses[4]);
  MPI_Address(&(gridtag.nCautions), &addresses[5]);
  MPI_Address(&(gridtag.nTFail),    &addresses[6]);
  MPI_Address(&(gridtag.lgBadEnd),  &addresses[7]);
  MPI_Address(&(gridtag.etime),     &addresses[8]);
  MPI_Address(&(gridtag.par1),      &addresses[9]);
  MPI_Address(&(gridtag.par2),      &addresses[10]);
 
  /* now far each is from the beginning of the structure */
  displacements[0] = addresses[1] - addresses[0];
  displacements[1] = addresses[2] - addresses[0];
  displacements[2] = addresses[3] - addresses[0];
  displacements[3] = addresses[4] - addresses[0];
  displacements[4] = addresses[5] - addresses[0];
  displacements[5] = addresses[6] - addresses[0];
  displacements[6] = addresses[7] - addresses[0];
  displacements[7] = addresses[8] - addresses[0];
  displacements[8] = addresses[9] - addresses[0];
  displacements[9] = addresses[10]- addresses[0];

  /*Create the derived type */
  MPI_Type_struct(10, block_lengths, displacements, typelist,message_type_ptr);
  
  /*Commit it so that it can be used */
  MPI_Type_commit(message_type_ptr);

} /*Build_derived_type */

# endif

int main( int argc, char *argv[] )
{
        DEBUG_ENTRY( "main()" );

        try {
                bool lgAbort;
                int IntegerMod;
                int  myrank ;
                double hden , absolute , relative , flux;
                char chVer[10] , chPath[1000] , chFilename[1000] , chFile[2000] ;

                double *xpar , *ypar, *par1, *par2;

                /* following will become array of line wavelengths and the number of lines 
                 * in this array */
                long int mod, nLines, LimModels, nModels;

                /* these will be passed to cdGetLineList and will become arrays of
                 * labels and wavelengths */
                char **chLabel ;
                realnum *wl;

                long int 
                        NumberWarnings1, NumberCautions1, NumberNotes1, NumberSurprises1, NumberTempFailures1,
                        NumberPresFailures1, NumberIonFailures1, NumberNeFailures1 ;
                /* number of time Cloudy returned with error condition */
                int nErrorExits;

                FILE *ioOUT , *ioDATA ;
                char chLine[100];

                long int n;
                /* number of processors, =1 for single non-MPI */
                int numprocs = 1;

                /* this will hold a single grid point */
                GRIDTAG grid;

                /* start MPI if -DMPI on command line */
#               ifdef MPI
                GRIDTAG gridtag;
                GRIDTAG *grids;
                
                MPI_Datatype message_type; /* Arguments to */

                MPI_Init(&argc, &argv);
                MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
                MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

                Build_derived_type(gridtag, &message_type);
#               else
                myrank = 0;
#               endif

                /* the following is the path for the directory where the output should appear,
                 * it needs to end with the OS standard directory delimiter, a "/" on unix */
                strcpy( chPath , PATHNAME );

                /* this is the first part of the name of the resulting data file.  The final
                 * file will have this name with a ".lin" after it. */
                strcpy( chFilename , FILENAME );

                /* ====================================================================== */
                /*
                 * make arrays of parameters for this grid */

                /* total number of models */
                LimModels = (long)((FLUXLIMIT - FLUXINIT)/INCREMENT + 1.1);
                LimModels *= (long)((HDENLIMIT - HDENINIT)/INCREMENT+1.1);
                /* add number of processors since this is most we would possibly have to 
                 * add on to make integer number of mpi calls */
                LimModels += numprocs;

                xpar = (double*)malloc(LimModels*sizeof(double) );
                ypar = (double*)malloc(LimModels*sizeof(double) );
                par1 = (double*)malloc(LimModels*sizeof(double) );
                par2 = (double*)malloc(LimModels*sizeof(double) );

                if( xpar==NULL || ypar==NULL || par1==NULL || par2==NULL )
                {
                        fprintf(stderr,"malloc failed\n");
#                       ifdef MPI 
                        MPI_Finalize(); 
#                       endif
                        cdEXIT(EXIT_FAILURE);
                }
                /* now define set of parameters */
                hden = HDENINIT;
                flux = FLUXINIT;
                nModels = 0;
                nErrorExits = 0;
#               ifdef MPI 
                /* loop along constant U diagonals, starting from
                 * bottom right of loc plane, increasing density, 
                 * and ending on the major diagonal.
                 * this is for load leveling on parallel machines */
                for( hden = HDENLIMIT; hden>= HDENINIT; hden-=INCREMENT )
                {
                        for( flux=FLUXINIT; flux<=FLUXLIMIT;  flux+= INCREMENT )
                        {
                                if( hden + (flux-FLUXINIT) > HDENLIMIT ) break;
                                assert( nModels < LimModels );
                                xpar[nModels] = hden + (flux-FLUXINIT);
                                ypar[nModels] = flux;
                                par1[nModels] = COLDEN;
                                par2[nModels] = 0.;
                                ++nModels;
                        }
                }
                /* loop increasing flux, starting at lower left plus increment */
                for( flux=FLUXINIT+INCREMENT; flux<=FLUXLIMIT;  flux+= INCREMENT )
                {
                        for( hden = HDENINIT; hden<= HDENLIMIT;  hden+=INCREMENT )
                        {
                                if( flux + (hden-HDENINIT) > FLUXLIMIT ) break;
                                assert( nModels < LimModels );
                                xpar[nModels] = hden ;
                                ypar[nModels] = flux + (hden-HDENINIT);
                                par1[nModels] = COLDEN;
                                par2[nModels] = 0.;
                                ++nModels;
                        }
                }
#               else
                /* on scalar machines do models in xy order */
                while( hden < 1.00001 * HDENLIMIT )
                {
                        while( flux < 1.00001*FLUXLIMIT )
                        {
                                xpar[nModels] = hden;
                                ypar[nModels] = flux;
                                par1[nModels] = COLDEN;
                                par2[nModels] = 0.;
                                ++nModels;
                                flux += INCREMENT;
                        }
                        flux = FLUXINIT;
                        hden += INCREMENT;
                }
#               endif
                /* increase total number to an integer multiple of number of processors */
                if( nModels%numprocs)
                {
                        /* add extra bit */
                        IntegerMod = nModels + (numprocs - nModels%numprocs);
                }
                else
                {
                        IntegerMod = nModels;
                }

                /* make up data for the "extra" models */
                for( mod=nModels; mod<IntegerMod; ++mod)
                {
                        xpar[mod] = xpar[0];
                        ypar[mod] = ypar[0];
                        par1[mod] = par1[0];
                        par2[mod] = par2[0];
                }

                /* ====================================================================== */

#               ifdef MPI
                /* allocate memory for grids */ 
                if( (grids = ((GRIDTAG *)malloc( (size_t)(numprocs)*sizeof(GRIDTAG )))) == NULL )
                {
                        fprintf(stderr," main could not malloc grid\n");
#                       ifdef MPI 
                        MPI_Finalize(); 
#                       endif
                        cdEXIT(EXIT_FAILURE);
                }
#               endif

                /* initialize the code so that cdGetLineList can be called */
                cdInit();
                /* get list of lines from standard data file.  With no arguments this reads
                 * the file BLRLineList, part of the standard data distribution.  There could
                 * have been another file name within the quotes - that file would have been
                 * used instead */
                nLines = cdGetLineList("" , &chLabel , &wl);
                /*nLines = cdGetLineList("NLRLineList.dat" , &chLabel , &wl);*/

                /* now copy version into string to print */
                cdVersion(chVer);

                /* open file for large output if this is desired */
                if( DO_PRINT )
                {
                        /* now create final filename */
                        strcpy( chFile , chPath );
                        strcat( chFile , chFilename );
                        strcat( chFile , ".out" );
                        ioOUT = fopen(chFile,"w");
                        if( ioOUT == NULL )
                        {
                                printf(" main could not open %s for writing.\n", chFile);
#                               ifdef MPI 
                                MPI_Finalize(); 
#                               endif
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                else
                {
                        /* there should be nothing coming here, but something might */
                        ioOUT = stderr;
                }

                /* calculation's results, all the lines in file ending .lin */
                strcpy( chFile , chPath );
                strcat( chFile , chFilename );
                strcat( chFile , ".lin" );
                ioDATA = fopen(chFile,"w");
                if( ioDATA == NULL )
                {
                        printf(" could not open %s for writing.\n" , chFile );
#                       ifdef MPI 
                        MPI_Finalize(); 
#                       endif
                        cdEXIT(EXIT_FAILURE);
                }

                /* print the header information before we begin the calculation */
                /* print version number on both headers */
                fprintf(ioDATA,"Cloudy version %s, %s\n",chVer,TITLE);

                for( mod=myrank; mod<IntegerMod; mod+=numprocs )
                {
                        hden = xpar[mod];
                        flux = ypar[mod];
                        /* zero out grid, which will hold all info for a single grid point */
                        memset(&grid, 0, sizeof(GRIDTAG));

                        /* initialize the code for this run */
                        cdInit();

                        /* direct the standard output to this file */
                        cdOutp(ioOUT);

                        /* send flag saying whether to create large output */
                        cdTalk( DO_PRINT );

                        /* title for grid */
                        sprintf(chLine,"title %s", TITLE);
                        n = cdRead( chLine  );

                        /* print only last iteration? */
#                       if PRINT_LAST
                        n = cdRead("print last iteration ");
#                       endif

                        /* option for fast model, or only header */
#                       if QUICK_MODEL
                        n = cdRead( "stop zone 3 "  );/**/
                        n = cdRead( "constant temper 4 "  );/**/
#                       endif
#                       if VERY_QUICK_MODEL
                        cdNoExec();
#                       endif

                        /* iterate if not fast model */
#                       if ( (!QUICK_MODEL) && (!VERY_QUICK_MODEL) &&(ITERATIONS!=0) )
#                       if ITERATIONS < 0
                        n = cdRead( "iterate convergence "  );/**/
#                       else 
                        sprintf(chLine,"iterate %i", ITERATIONS);
                        n = cdRead( chLine  );
#                       endif
#                       endif

#                       ifdef MPI 
                        /* set flag so that exit handler will call MPI_Finalize */
                        cdMPI();
#                       endif

                        /* option to turn off file buffering if code might crash */
#                       ifdef NO_BUFFERING
                        n = cdRead("no buffering ");
#                       endif

                        /* if we have run off end of array of models say not to execute,
                         * still must do something on this processor so that the gather can occur
                         * at the bottom of the loop */
                        if( mod >= nModels )
                                cdNoExec();

                        /* flux and density */
                        sprintf(chLine,"phi(h) %f", flux);
                        n = cdRead( chLine  );

                        sprintf(chLine,"hden %f", hden);
                        n = cdRead( chLine  );

                        /* turn on large FeII atom, VERY slow */
#                       if FEII
                        n = cdRead( "atom feii  "  );
#                       endif

                        /* this is the continuum we used in the large apjs paper */
                        n = cdRead( "agn 6.00 -1.40 -0.50 -1.0  "  );/**/

                        /* a bare powerlaw with a break at 1 microns */
                        /*n = cdRead( "table powerlaw -1.0 1 micron break "  );*/

                        /* the much maligned Mathews&Ferland continuum */
                        /*n = cdRead( "table agn  "  );*/

                        /*n = cdRead( "blackbody 142,000  "  );*/

                        /* broken power law used in Fred's paper */
                        /*n = cdRead( "interpolate (0.00000001, -14.899), (0.0091126732, 0.000)  "  );
                          n = cdRead( "continue (1.00, -1.2242), (73.54, -4.2106)  "  );
                          n = cdRead( "continue (3675, -5.7395), (7354000, -11.2526)  "  );*/

                        /* broken power law used in Mark's paper */
                        /*n = cdRead( "interpolate (10.517, -30.850) (13.383, -23.684)   "  );
                          n = cdRead( "continue (16.684, -26.986) (17.383, -27.9996)    "  );
                          n = cdRead( "continue (18.383, -28.8899) (19.124, -29.268)    "  );
                          n = cdRead( "continue (22.383, -35.788)  "  );*/

                        /* add cosmic IR background as well */
                        n = cdRead( "background z=0  "  );

                        sprintf(chLine,"stop column density %f ", par1[mod]);
                        n = cdRead( chLine  );/**/

                        n = cdRead( "failures 3 "  );

                        /* options to change abundances */
#                       if 0
                        n = cdRead( "metals 10 "  );/**/
                        n = cdRead( "element nitrogen scale 10 "  );/**/
                        n = cdRead( "element helium scale 1.66 "  );/**/
                        n = cdRead( "metals 20 "  );/**/
                        n = cdRead(