grid_xspec.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*gridXspec handles all grid calculations, called by griddo */
/*gridFunc */
/*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */
#include "cddefines.h"
#include "cddrive.h"
#include "continuum.h"
#include "rfield.h"
#include "grid.h"
#include "ipoint.h"
#include "called.h"
#include "prt.h"

/*gridXspec handles all grid calculations, called by grid_do */
void gridXspec(realnum xc[], long int nInterpVars)
{
        long int i;

        DEBUG_ENTRY( "gridXspec()" );

        if( nInterpVars > LIMPAR )
        {
                fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" );
                cdEXIT(EXIT_FAILURE);
        }

        optimize.nOptimiz = 0;
        grid.nintparm = nInterpVars;

        /* if this is changed there must be some change made to actually
         * stuff the additional parameter information. */
        grid.naddparm = 0;

        ASSERT( grid.nintparm + grid.naddparm <= LIMPAR );

        grid.totNumModels = 1;
        /* >>chng 06 aug 21, allow the number of parameter values to be different for different parameters. */
        for( i=0; i<nInterpVars; i++ )
        {
                /* >>chng 06 sep 4, use grid variable instead of passing to routine. */
                grid.totNumModels *= grid.numParamValues[i];
        
                if( grid.paramValuesFromList[i].size() != 0 && grid.lgSaveXspec )
                {
                        fprintf( ioQQQ, " PROBLEM The 'save XSPEC' and 'grid list' options do not work together.\n" );
                        fprintf( ioQQQ, " If any 'save XSPEC' command is given, all grid commands must follow this syntax:\n" );
                        fprintf( ioQQQ, "       grid <p1> <p2> <p3> \n" );
                        fprintf( ioQQQ, " where p1 and p2 are the limits and p3 is a regular increment.\n Sorry.\n" );
                        cdEXIT( EXIT_FAILURE );
                }
        }

        // option to cycle through grid multiple times, default is 1
        grid.totNumModels *= grid.nCycle;
        ASSERT( grid.totNumModels > 1 );

        grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.paramRange = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.paramData = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.interpParameters = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) );

        for( i=0; i<nInterpVars+grid.naddparm; i++ )
        {
                grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) );
                grid.paramRange[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(6) );
                grid.paramData[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numParamValues[i]) );

                sprintf( grid.paramNames[i],    "%s%ld", "PARAM", i+1 );
                /* Method is 0 for linear, 1 for logarithmic */
                grid.paramMethods[i] = 0;
                /* Initial */
                grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f;
                /* Delta */
                grid.paramRange[i][1] = grid.paramIncrements[i]/10.f;
                /* Minimum */
                grid.paramRange[i][2] = xc[i];
                /* Bottom */
                grid.paramRange[i][3] = xc[i]+grid.paramIncrements[i]/10.f;
                /* Top */
                grid.paramRange[i][4] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)-grid.paramIncrements[i]/10.f;
                /* Maximum */
                grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f);

                for( long j=0; j<grid.numParamValues[i]; j++ )
                {
                        grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j;
                }
        }

        for( i=0; i<grid.totNumModels; i++ )
        {
                grid.interpParameters[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nInterpVars) );
        }

        for( i=0; i< grid.totNumModels; i++ )
        {
                long j;
                realnum variableVector[LIMPAR];

                for( j=0; j<nInterpVars; j++ )
                {
                        int index;
                        long volumeOtherDimensions = 1;

                        /* by "volume", we simply mean the product of the parameter dimensions 
                         * AFTER the present index, i.e.:
                         * first "volume" is product of grid.numParamValues[1]*grid.numParamValues[2]*....grid.numParamValues[n]
                         * second "volume" is product of grid.numParamValues[2]*grid.numParamValues[3]*....grid.numParamValues[n]
                         * last "volume" is unity.  */
                        for( long k=j+1; k<nInterpVars; k++ )
                        {
                                volumeOtherDimensions *= grid.numParamValues[k];
                        }

                        /* For each successive parameter, the "volume" is less than the previous one.
                         * So the left-hand side of this modulus operation increases for each parameter,
                         * which causes the index of each parameter to be incremented more often than the
                         * index of the previous parameter.  Thus, the last dimension is incremented
                         * every time, then the second to last dimension is incremented with each repeat
                         * of the last dimension.  This repeats until, finally, the first dimension is
                         * incremented.  For example, the indices of the parameter vectors for a 2x2x3
                         * box would be ordered as such:
                         * [0][0][0]
                         * [0][0][1]
                         * [0][0][2]
                         * [0][1][0]
                         * [0][1][1]
                         * [0][1][2]
                         * [1][0][0]
                         * [1][0][1]
                         * [1][0][2]
                         * [1][1][0]
                         * [1][1][1]
                         * [1][1][2]
                         */
                        index = (int)( (i/volumeOtherDimensions)%(grid.numParamValues[j]) );

                        /* this prevents parameter incrementation for debugging purposes.  */
                        if( grid.lgStrictRepeat )
                                variableVector[j] = xc[j];
                        else if( grid.paramValuesFromList[j].size() != 0U )
                                variableVector[j] = grid.paramValuesFromList[j][index];
                        else
                                variableVector[j] = xc[j] + grid.paramIncrements[j]*index;

                        grid.interpParameters[i][j] = variableVector[j];

                        if( grid.lgLinearIncrements[j] && !optimize.lgOptimizeAsLinear[j] )
                                variableVector[j] = log10(variableVector[j]);
                }

                for( j=nInterpVars; j<LIMPAR; j++ )
                {
                        variableVector[j] = xc[j];
                }

                if( i == grid.totNumModels - 1 )
                {
                        fixit("is this needed ??");
                        called.lgTalk = cpu.i().lgMPI_talk();
                        called.lgTalkIsOK = cpu.i().lgMPI_talk();
                        prt.lgFaintOn = true;
                        grid.lgGridDone = true;
                }

                (void)optimize_func(variableVector,i);
        }
        return;
}

/*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */
void GridGatherInCloudy()
{
        long i;

        DEBUG_ENTRY( "GridGatherInCloudy()" );

        ASSERT( grid.lgGrid );

        /* malloc some arrays if first call and save continuum energies. */
        if( grid.Energies.empty() )
        {
                long i1, i2;

                // this will happen if we have more MPI ranks than grid points
                // the highest ranks will not have executed any model
                if( rfield.nflux_with_check <= 0 )
                        ContCreateMesh();

                if( grid.LoEnergy_keV == 0. )
                        grid.ipLoEnergy = 0;
                else
                        grid.ipLoEnergy = ipoint( grid.LoEnergy_keV * 1000. / EVRYD );

                if( grid.HiEnergy_keV == 0. || grid.HiEnergy_keV * 1000. / EVRYD > rfield.egamry() )
                        grid.ipHiEnergy = rfield.nflux_with_check - 1;
                else
                        grid.ipHiEnergy = ipoint( grid.HiEnergy_keV * 1000. / EVRYD );

                grid.numEnergies = grid.ipHiEnergy - grid.ipLoEnergy + 1;
                ASSERT( grid.numEnergies > 0 );
                grid.Energies.resize( grid.numEnergies );
                grid.Spectra.reserve(NUM_OUTPUT_TYPES);
                for( i1=0; i1 < NUM_OUTPUT_TYPES; i1++ )
                {
                        if( grid.lgOutputTypeOn[i1] )
                        {
                                grid.Spectra.reserve(i1,grid.totNumModels);
                                for( i2=0; i2 < grid.totNumModels; i2++ )
                                {
                                        grid.Spectra.reserve(i1,i2,grid.numEnergies);
                                }
                        }
                }
                grid.Spectra.alloc();
                // this needs to be zeroed out for MPI runs
                grid.Spectra.zero();

                for( i1=0; i1<grid.numEnergies; i1++ )
                {
                        long j = grid.ipLoEnergy + i1;
                        grid.Energies[i1] = rfield.anu(j);
                }
        }

        if( optimize.nOptimiz < grid.totNumModels )
        {
                ASSERT( optimize.nOptimiz >= 0 );

                for( i=0; i < NUM_OUTPUT_TYPES; i++ )
                {
                        /* Grab spectrum for xspec printout 
                         * at this point nOptimiz has already been incremented for first model */
                        if( grid.lgOutputTypeOn[i] )
                                cdSPEC2( i, grid.numEnergies, grid.ipLoEnergy, grid.ipHiEnergy,
                                         &grid.Spectra[i][optimize.nOptimiz][0]);
                }
        }
        else if( optimize.nOptimiz == grid.totNumModels )
        {
                if( cpu.i().lgMPI() )
                {
                        multi_arr<realnum,3> Spectra_Copy = grid.Spectra;

                        // combine the grid.Spectra data from all ranks. This is done by adding up
                        // the results from all ranks. All but one should be zero. This is needed
                        // because we do not know which rank calculated which grid point...
                        for( int i=0; i < NUM_OUTPUT_TYPES; ++i )
                        {
                                if( grid.lgOutputTypeOn[i] )
                                {
                                        for( int j=0; j < grid.totNumModels; ++j )
                                        {
                                                MPI_Reduce( &Spectra_Copy[i][j][0], &grid.Spectra[i][j][0], 
                                                            grid.numEnergies, MPI_type(grid.Spectra[i][j][0]),
                                                            MPI_SUM, 0, MPI_COMM_WORLD );
                                        }
                                }
                        }
                        MPI_Barrier( MPI_COMM_WORLD );
                }
        }
        else
        {
                TotalInsanity();
        }
        return;
}