00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "save.h"
00008 #include "warnings.h"
00009 #include "optimize.h"
00010 #include "cddrive.h"
00011 #include "continuum.h"
00012 #include "rfield.h"
00013 #include "grid.h"
00014 #include "ipoint.h"
00015 #include "called.h"
00016 #include "physconst.h"
00017 #include "prt.h"
00018 #include "mpi_utilities.h"
00019
00020
00021 void gridXspec(realnum xc[], long int nInterpVars)
00022 {
00023 long int i;
00024
00025 DEBUG_ENTRY( "gridXspec()" );
00026
00027 if( nInterpVars > LIMPAR )
00028 {
00029 fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" );
00030 cdEXIT(EXIT_FAILURE);
00031 }
00032
00033 optimize.nOptimiz = 0;
00034 grid.nintparm = nInterpVars;
00035
00036
00037
00038 grid.naddparm = 0;
00039
00040 ASSERT( grid.nintparm + grid.naddparm <= LIMPAR );
00041
00042 grid.totNumModels = 1;
00043
00044 for( i=0; i<nInterpVars; i++ )
00045 {
00046
00047 grid.totNumModels *= grid.numParamValues[i];
00048 }
00049
00050 ASSERT( grid.totNumModels > 1 );
00051
00052 grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) );
00053 grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) );
00054 grid.paramRange = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
00055 grid.paramData = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
00056 grid.interpParameters = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) );
00057
00058 for( i=0; i<nInterpVars+grid.naddparm; i++ )
00059 {
00060 grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) );
00061 grid.paramRange[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(6) );
00062 grid.paramData[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numParamValues[i]) );
00063
00064 sprintf( grid.paramNames[i], "%s%ld", "PARAM", i+1 );
00065
00066 grid.paramMethods[i] = grid.lgLinearIncrements[i] ? 0 : 1;
00067
00068 grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f;
00069
00070 grid.paramRange[i][1] = grid.paramIncrements[i]/10.f;
00071
00072 grid.paramRange[i][2] = xc[i];
00073
00074 grid.paramRange[i][3] = xc[i]+grid.paramIncrements[i]/10.f;
00075
00076 grid.paramRange[i][4] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)-grid.paramIncrements[i]/10.f;
00077
00078 grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f);
00079
00080 for( long j=0; j<grid.numParamValues[i]; j++ )
00081 {
00082 grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j;
00083 }
00084 }
00085
00086 for( i=0; i<grid.totNumModels; i++ )
00087 {
00088 grid.interpParameters[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nInterpVars) );
00089 }
00090
00091 for( i=0; i< grid.totNumModels; i++ )
00092 {
00093 long j;
00094 realnum variableVector[LIMPAR];
00095
00096 for( j=0; j<nInterpVars; j++ )
00097 {
00098 int index;
00099 long volumeOtherDimensions = 1;
00100
00101
00102
00103
00104
00105
00106 for( long k=j+1; k<nInterpVars; k++ )
00107 {
00108 volumeOtherDimensions *= grid.numParamValues[k];
00109 }
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 index = (int)( (i/volumeOtherDimensions)%(grid.numParamValues[j]) );
00133
00134
00135 if( grid.lgStrictRepeat )
00136 variableVector[j] = xc[j];
00137 else
00138 variableVector[j] = xc[j] + grid.paramIncrements[j]*index;
00139
00140 grid.interpParameters[i][j] = variableVector[j];
00141
00142 if( grid.lgLinearIncrements[j] && !optimize.lgOptimizeAsLinear[j] )
00143 variableVector[j] = log10(variableVector[j]);
00144 }
00145
00146 for( j=nInterpVars; j<LIMPAR; j++ )
00147 {
00148 variableVector[j] = xc[j];
00149 }
00150
00151 if( i == grid.totNumModels - 1 )
00152 {
00153 fixit();
00154 called.lgTalk = cpu.lgMPI_talk();
00155 called.lgTalkIsOK = cpu.lgMPI_talk();
00156 prt.lgFaintOn = true;
00157 grid.lgGridDone = true;
00158 }
00159
00160 (void)optimize_func(variableVector,i);
00161 }
00162 return;
00163 }
00164
00165
00166 void GridGatherInCloudy( void )
00167 {
00168 long i;
00169
00170 DEBUG_ENTRY( "GridGatherInCloudy()" );
00171
00172 ASSERT( grid.lgGrid );
00173
00174
00175 if( grid.Energies.empty() )
00176 {
00177 long i1, i2;
00178
00179
00180
00181 if( rfield.nupper <= 0 )
00182 ContCreateMesh();
00183
00184 if( grid.LoEnergy_keV == 0. )
00185 grid.ipLoEnergy = 0;
00186 else
00187 grid.ipLoEnergy = ipoint( grid.LoEnergy_keV * 1000. / EVRYD );
00188
00189 if( grid.HiEnergy_keV == 0. || grid.HiEnergy_keV >= continuum.filbnd[continuum.nrange] )
00190 grid.ipHiEnergy = rfield.nupper - 1;
00191 else
00192 grid.ipHiEnergy = ipoint( grid.HiEnergy_keV * 1000. / EVRYD );
00193
00194 grid.numEnergies = grid.ipHiEnergy - grid.ipLoEnergy + 1;
00195 ASSERT( grid.numEnergies > 0 );
00196 grid.Energies.resize( grid.numEnergies );
00197 grid.Spectra.reserve(NUM_OUTPUT_TYPES);
00198 for( i1=0; i1 < NUM_OUTPUT_TYPES; i1++ )
00199 {
00200 if( grid.lgOutputTypeOn[i1] )
00201 {
00202 grid.Spectra.reserve(i1,grid.totNumModels);
00203 for( i2=0; i2 < grid.totNumModels; i2++ )
00204 {
00205 grid.Spectra.reserve(i1,i2,grid.numEnergies);
00206 }
00207 }
00208 }
00209 grid.Spectra.alloc();
00210
00211 grid.Spectra.zero();
00212
00213 for( i1=0; i1<grid.numEnergies; i1++ )
00214 {
00215 long j = grid.ipLoEnergy + i1;
00216 grid.Energies[i1] = rfield.AnuOrg[j];
00217 }
00218 }
00219
00220 if( optimize.nOptimiz < grid.totNumModels )
00221 {
00222 ASSERT( optimize.nOptimiz >= 0 );
00223
00224 for( i=0; i < NUM_OUTPUT_TYPES; i++ )
00225 {
00226
00227
00228 if( grid.lgOutputTypeOn[i] )
00229 cdSPEC2( i, grid.numEnergies, grid.ipLoEnergy, grid.ipHiEnergy,
00230 &grid.Spectra[i][optimize.nOptimiz][0]);
00231 }
00232 }
00233 else if( optimize.nOptimiz == grid.totNumModels )
00234 {
00235 ASSERT( cpu.lgMPI() );
00236
00237 multi_arr<realnum,3> Spectra_Copy = grid.Spectra;
00238
00239
00240
00241
00242 for( int i=0; i < NUM_OUTPUT_TYPES; ++i )
00243 {
00244 if( grid.lgOutputTypeOn[i] )
00245 {
00246 for( int j=0; j < grid.totNumModels; ++j )
00247 {
00248 MPI::COMM_WORLD.Reduce( &Spectra_Copy[i][j][0], &grid.Spectra[i][j][0],
00249 grid.numEnergies, MPI::type(grid.Spectra[i][j][0]),
00250 MPI::SUM, 0 );
00251 }
00252 }
00253 }
00254 MPI::COMM_WORLD.Barrier();
00255 }
00256 else
00257 {
00258 TotalInsanity();
00259 }
00260 return;
00261 }