cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*grid_do called by cdDrive, this returns 0 if things went ok, 1 for disaster */ 00004 #include "cddefines.h" 00005 #include "conv.h" 00006 #include "input.h" 00007 #include "called.h" 00008 #include "version.h" 00009 #include "init.h" 00010 #include "prt.h" 00011 #include "trace.h" 00012 #include "grains.h" 00013 #include "parse.h" 00014 #include "punch.h" 00015 #include "optimize.h" 00016 #include "grid.h" 00017 00018 /*grid_do called by cdDrive, calls gridXspec, returns false if ok, true for disaster */ 00019 bool grid_do(void) 00020 { 00021 char chLine[INPUT_LINE_LENGTH], 00022 chNote[8]; 00023 long int i, 00024 ii, 00025 j; 00026 realnum ptem[LIMPAR]; 00027 00028 DEBUG_ENTRY( "grid_do()" ); 00029 00030 /* main driver for optimization runs 00031 * Drives cloudy to grid variables;*/ 00032 00033 /* code originally written by R.F. Carswell, IOA Cambridge */ 00034 00035 /* variables with optimizer */ 00036 for( i=0; i < LIMPAR; i++ ) 00037 { 00038 optimize.OptIncrm[i] = 0.; 00039 optimize.varang[i][0] = -FLT_MAX; 00040 optimize.varang[i][1] = FLT_MAX; 00041 /* this should be overwritten by format of vary line */ 00042 strcpy( optimize.chVarFmt[i], "error - no optimizer line image was set" ); 00043 } 00044 00045 /* necessary to do this to keep all lines in */ 00046 prt.lgFaintOn = false; 00047 conv.LimFail = 1000; 00048 00049 /* this initializes variables at the start of each simulation 00050 * in a grid, before the parser is called - this must set any values 00051 * that may be changed by the command parser */ 00052 InitDefaultsPreparse(); 00053 00054 /* call READR the first time to scan off all variable options */ 00055 /* this is just an initial parsing to get the number of iterations and 00056 * the number of varied parameters. The other Init* routines are not 00057 * called after this because this is all done again later for each grid point */ 00058 ParseCommands(); 00059 00060 /* done parsing, now turn off punch headers. */ 00061 punch.lgPunHeader = false; 00062 00063 /* >>chng 06 may 20, do not turn off printout if grid calculation is in place 00064 * rather than optimization */ 00065 if( !grid.lgGrid ) 00066 { 00067 called.lgTalk = false; 00068 /* this flag is needed to turn print on to have effect */ 00069 called.lgTalkIsOK = false; 00070 } 00071 00072 /* >>chng 00 aug 09, return memory allocated for grains, they are not used, PvH */ 00073 ReturnGrainBins(); 00074 00075 optimize.nvary = optimize.nparm; 00076 00077 /* check that more than 1 observed intensities or column densities were entered */ 00078 if( ((optimize.nlobs + (long)optimize.lgOptLum + optimize.nTempObs + optimize.ncobs) < 1 ) && 00079 !grid.lgGrid ) 00080 { 00081 fprintf( ioQQQ, " The input stream has vary commands, but\n" ); 00082 fprintf( ioQQQ, " no observed quantities were entered. Whats up?\n" ); 00083 fprintf( ioQQQ, " Use the NO VARY command to input vary options but not try to perform this.\n" ); 00084 cdEXIT(EXIT_FAILURE); 00085 } 00086 00087 /* check that the total number of parameters to vary is greater than 1 */ 00088 if( optimize.nvary < 1 ) 00089 { 00090 fprintf( ioQQQ, " No parameters to vary were entered. Whats up?\n" ); 00091 cdEXIT(EXIT_FAILURE); 00092 } 00093 00094 if( ( strcmp(optimize.chOptRtn,"XSPE") == 0 ) && ( optimize.nRangeSet != optimize.nvary ) ) 00095 { 00096 fprintf( ioQQQ, " Every parameter with a VARY option must have a GRID specified,\n" ); 00097 fprintf( ioQQQ, " and the GRID must be specified after the VARY option.\n" ); 00098 fprintf( ioQQQ, " These requirements were not satisfied for %ld parameter(s).\n", abs(optimize.nvary - optimize.nRangeSet) ); 00099 cdEXIT(EXIT_FAILURE); 00100 } 00101 00102 /* lgTrOptm set with trace grid command */ 00103 if( trace.lgTrOptm ) 00104 { 00105 for( i=0; i < optimize.nvary; i++ ) 00106 { 00107 /*print the command format as debugging aid */ 00108 fprintf( ioQQQ, "%s\n", optimize.chVarFmt[i]); 00109 00110 /* now generate the actual command with parameter, 00111 * there will be from 1 to 3 numbers on the line */ 00112 if( optimize.nvarxt[i] == 1 ) 00113 { 00114 /* case with 1 parameter */ 00115 sprintf( chLine , optimize.chVarFmt[i], optimize.vparm[0][i] ); 00116 } 00117 00118 else if( optimize.nvarxt[i] == 2 ) 00119 { 00120 /* case with 2 parameter */ 00121 sprintf( chLine , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]); 00122 } 00123 00124 else if( optimize.nvarxt[i] == 3 ) 00125 { 00126 /* case with 3 parameter */ 00127 sprintf( chLine , optimize.chVarFmt[i], 00128 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] ); 00129 } 00130 00131 else if( optimize.nvarxt[i] == 4 ) 00132 { 00133 /* case with 4 parameter */ 00134 sprintf( chLine , optimize.chVarFmt[i], 00135 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] ); 00136 } 00137 00138 else if( optimize.nvarxt[i] == 5 ) 00139 { 00140 /* case with 5 parameter */ 00141 sprintf( chLine , optimize.chVarFmt[i], 00142 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], 00143 optimize.vparm[3][i], optimize.vparm[4][i]); 00144 } 00145 00146 else 00147 { 00148 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me1\n"); 00149 cdEXIT(EXIT_FAILURE); 00150 } 00151 00152 /* print the resulting command line*/ 00153 fprintf( ioQQQ, "%s\n", chLine ); 00154 } 00155 } 00156 00157 /* option to change default increments; if zero then leave as is */ 00158 for( i=0; i < LIMPAR; i++ ) 00159 { 00160 if( optimize.OptIncrm[i] != 0. ) 00161 { 00162 optimize.vincr[i] = optimize.OptIncrm[i]; 00163 } 00164 } 00165 00166 /* say who we are */ 00167 if( strcmp(optimize.chOptRtn,"XSPE") == 0 ) 00168 { 00169 fprintf( ioQQQ, " Grid Driver\n" ); 00170 } 00171 else 00172 { 00173 fprintf( ioQQQ, " Optimization Driver\n" ); 00174 } 00175 int indent = (int)((122 - strlen(t_version::Inst().chVersion))/2); 00176 fprintf( ioQQQ, "%*cCloudy %s\n\n",indent,' ',t_version::Inst().chVersion); 00177 fprintf( ioQQQ, "%23c**************************************%7.7s**************************************\n", ' ', t_version::Inst().chDate ); 00178 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' ); 00179 00180 /* now echo initial input quantities with flag for vary */ 00181 /* first loop steps over all command lines entered */ 00182 for( i=0; i <= input.nSave; i++ ) 00183 { 00184 /* put space to start line, overwrite if vary found */ 00185 strcpy( chNote, " " ); 00186 /* loop over all vary commands, see if this is one */ 00187 for( j=0; j < optimize.nvary; j++ ) 00188 { 00189 /* input.nSave is on C array counting, rest are on fortran */ 00190 if( i == optimize.nvfpnt[j] ) 00191 { 00192 /* this is a vary command, put keyword at start */ 00193 strcpy( chNote, "VARY>>>" ); 00194 } 00195 } 00196 00197 fprintf( ioQQQ, "%22.7s * %-80s*\n", chNote, input.chCardSav[i] ); 00198 } 00199 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' ); 00200 fprintf( ioQQQ, "%23c***********************************************************************************\n\n\n", ' ' ); 00201 00202 /* option to trace logical flow within this sub */ 00203 if( optimize.lgOptimFlow ) 00204 { 00205 for( j=0; j < optimize.nvary; j++ ) 00206 { 00207 i = optimize.nvfpnt[j]; 00208 fprintf( ioQQQ, " trace:%80.80s\n", input.chCardSav[i]); 00209 fprintf( ioQQQ, "%80.80s\n", optimize.chVarFmt[j]); 00210 fprintf( ioQQQ, " number of variables on line:%4ld\n", 00211 optimize.nvarxt[j] ); 00212 fprintf( ioQQQ, " Values:" ); 00213 for( ii=1; ii <= optimize.nvarxt[j]; ii++ ) 00214 { 00215 fprintf( ioQQQ, "%10.2e", optimize.vparm[ii-1][j] ); 00216 } 00217 fprintf( ioQQQ, "\n" ); 00218 } 00219 } 00220 00221 00222 if( strcmp(optimize.chOptRtn,"PHYM") == 0 ) 00223 { 00224 fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 00225 optimize.nIterOptim ); 00226 fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 00227 chOptimFileName ); 00228 00229 fprintf( ioQQQ, " The optimize_phymir method will be used" ); 00230 if( optimize.lgParallel ) { 00231 fprintf( ioQQQ, " in parallel mode.\n The maximum no. of CPU's to be used is %1ld.\n",optimize.useCPU ); 00232 } 00233 else { 00234 fprintf( ioQQQ, " in sequential mode.\n" ); 00235 } 00236 } 00237 00238 else if( strcmp(optimize.chOptRtn,"SUBP") == 0 ) 00239 { 00240 fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 00241 optimize.nIterOptim ); 00242 fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 00243 chOptimFileName ); 00244 00245 fprintf( ioQQQ, " Subplex method will be used.\n" ); 00246 } 00247 00248 else if( strcmp(optimize.chOptRtn,"XSPE") == 0 ) 00249 { 00250 fprintf( ioQQQ, " Producing grid output.\n" ); 00251 } 00252 00253 else 00254 { 00255 fprintf( ioQQQ, " I do not understand what method to use.\n" ); 00256 fprintf( ioQQQ, " Sorry.\n" ); 00257 cdEXIT(EXIT_FAILURE); 00258 } 00259 00260 fprintf( ioQQQ, "\n%4ld parameter(s) will be varied. The first lines, and the increments are:\n", 00261 optimize.nvary ); 00262 00263 for( i=0; i < optimize.nvary; i++ ) 00264 { 00265 optimize.varmax[i] = -FLT_MAX; 00266 optimize.varmin[i] = FLT_MAX; 00267 /* write formatted to output using the format held in chVarFmt(np) */ 00268 00269 /* now generate the actual command with parameter, 00270 * there will be from 1 to 3 numbers on the line */ 00271 if( optimize.nvarxt[i] == 1 ) 00272 { 00273 /* case with 1 parameter */ 00274 sprintf( chLine , optimize.chVarFmt[i], optimize.vparm[0][i] ); 00275 } 00276 00277 else if( optimize.nvarxt[i] == 2 ) 00278 { 00279 /* case with 2 parameter */ 00280 sprintf( chLine , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]); 00281 } 00282 00283 else if( optimize.nvarxt[i] == 3 ) 00284 { 00285 /* case with 3 parameter */ 00286 sprintf( chLine , optimize.chVarFmt[i], 00287 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] ); 00288 } 00289 00290 else if( optimize.nvarxt[i] == 4 ) 00291 { 00292 /* case with 4 parameter */ 00293 sprintf( chLine , optimize.chVarFmt[i], 00294 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] , optimize.vparm[3][i] ); 00295 } 00296 00297 else if( optimize.nvarxt[i] == 5 ) 00298 { 00299 /* case with 5 parameter */ 00300 sprintf( chLine , optimize.chVarFmt[i], 00301 optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], 00302 optimize.vparm[3][i] , optimize.vparm[4][i]); 00303 } 00304 00305 else 00306 { 00307 fprintf(ioQQQ,"The number of variable options on this line makes no sense to me2\n"); 00308 cdEXIT(EXIT_FAILURE); 00309 } 00310 00311 fprintf( ioQQQ, "\n %s\n", chLine ); 00312 fprintf( ioQQQ, " Initial increment is%6.3f, the limits are%10.2e to %10.2e\n", 00313 optimize.vincr[i], optimize.varang[i][0], optimize.varang[i][1] ); 00314 } 00315 00316 /* this will be number of times grid calls cloudy */ 00317 optimize.nOptimiz = 0; 00318 /* >>chng 06 jan 26, moved this to an option on the command line */ 00319 /* grid.numParamValues = 10; */ 00320 00321 if( strcmp(optimize.chOptRtn,"XSPE") == 0 ) 00322 { 00323 for( j=0; j < optimize.nvary; j++ ) 00324 { 00325 /* ptem[j] = optimize.vparm[0][j]; */ 00326 ptem[j] = optimize.varang[j][0]; 00327 /* >> chng 06 sep 4, delta is now read at command line. */ 00328 /* delta[j] = optimize.vincr[j]; */ 00329 /* delta[j] = ( optimize.varang[j][1] - optimize.varang[j][0] )/ 00330 ((realnum)grid.numParamValues[j] - 1.f); */ 00331 } 00332 /* >>chng 06 aug 23, fill rest of array with zeros. */ 00333 for( j=optimize.nvary; j < LIMPAR; j++ ) 00334 { 00335 ptem[j] = 0.f; 00336 grid.paramIncrements[j] = 0.f; 00337 } 00338 00339 gridXspec(ptem,optimize.nvary); 00340 for( j=0; j < optimize.nvary; j++ ) 00341 { 00342 optimize.vparm[0][j] = ptem[j]; 00343 } 00344 00345 fprintf( ioQQQ, " **************************************************\n" ); 00346 fprintf( ioQQQ, " **************************************************\n" ); 00347 fprintf( ioQQQ, " **************************************************\n" ); 00348 fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", optimize.nOptimiz ); 00349 00350 /* though a page eject */ 00351 fprintf( ioQQQ, "\f" ); 00352 } 00353 else 00354 { 00355 lgAbort = lgOptimize_do(); 00356 } 00357 00358 if( lgAbort ) 00359 { 00360 /* busted set means there were serious problems somewhere */ 00361 return 1; 00362 } 00363 else 00364 { 00365 /* return 0 if everything is ok */ 00366 return 0; 00367 } 00368 }