00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*cloudy the main routine, this IS Cloudy, return 0 normal exit, 1 error exit, 00004 * called by maincl when used as standalone program */ 00005 /*BadStart announce that things are so bad the calculation cannot even start */ 00006 #include "cddefines.h" 00007 #include "save.h" 00008 #include "noexec.h" 00009 #include "lines.h" 00010 #include "abund.h" 00011 #include "continuum.h" 00012 #include "warnings.h" 00013 #include "atmdat.h" 00014 #include "prt.h" 00015 #include "conv.h" 00016 #include "parse.h" 00017 #include "dynamics.h" 00018 #include "init.h" 00019 #include "opacity.h" 00020 #include "state.h" 00021 #include "rt.h" 00022 #include "monitor_results.h" 00023 #include "zones.h" 00024 #include "iterations.h" 00025 #include "plot.h" 00026 #include "radius.h" 00027 #include "grid.h" 00028 #include "cloudy.h" 00029 #include "ionbal.h" 00030 #include "dense.h" 00031 #include "energy.h" 00032 #include "called.h" 00033 00034 /*BadStart announce that things are so bad the calculation cannot even start */ 00035 STATIC void BadStart(); 00036 00037 /* returns 1 if disaster strikes, 0 if everything appears ok */ 00038 bool cloudy() 00039 { 00040 DEBUG_ENTRY( "cloudy()" ); 00041 00042 bool lgOK; 00043 00044 /* 00045 * this is the main routine to drive the modules that compute a model 00046 * when cloudy is used as a stand-alone program 00047 * the main program (maincl) calls cdInit then cdDrive 00048 * this sub is called by cdDrive which returns upon exiting 00049 * 00050 * this routine uses the following variables: 00051 * 00052 * nzone 00053 * this is the zone number, and is incremented here 00054 * is zero during search phase, 1 for first zone at illuminated face 00055 * logical function iter_end_check returns a true condition if NZONE reaches 00056 * NEND(ITER), the limit to the number of zones in this iteration 00057 * nzone is totally controlled in this subroutine 00058 * 00059 * iteration 00060 * this is the iteration number, it is 1 for the first iteration 00061 * and is incremented in this subroutine at the end of each iteration 00062 * 00063 * iterations.itermx 00064 * this is the limit to the number of iterations, and is entered 00065 * by the user 00066 * This routine returns when iteration > iterations.itermx 00067 */ 00068 00069 /* nzone is zero while initial search for conditions takes place 00070 * and is incremented to 1 for start of first zone */ 00071 nzone = 0; 00072 fnzone = 0.; 00073 00074 /* iteration is iteration number, calculation is complete when iteration > iterations.itermx */ 00075 iteration = 1; 00076 00077 /* this initializes variables at the start of each simulation 00078 * in a grid, before the parser is called - this must set any values 00079 * that may be changed by the command parser */ 00080 InitDefaultsPreparse(); 00081 00082 /* scan in and parse input data */ 00083 ParseCommands(); 00084 00085 /* fix abundances of elements, in abundances.cpp */ 00086 AbundancesSet(); 00087 00088 ASSERT(lgElemsConserved()); 00089 00090 /* one time creation of some variables */ 00091 InitCoreloadPostparse(); 00092 00093 /* initialize vars that can only be done after parsing the commands 00094 * sets up CO network since this depends on which elements are on 00095 * and whether chemistry is enabled */ 00096 InitSimPostparse(); 00097 00098 /* ContCreateMesh calls fill to set up continuum energy mesh if first call, 00099 * otherwise reset to original mesh. 00100 * This is AFTER ParseCommands so that 00101 * path and mesh size can be set with commands before mesh is set */ 00102 ContCreateMesh(); 00103 00104 /* create several data sets by allocating memory and reading in data files, 00105 * but only if this is first call */ 00106 atmdat_readin(); 00107 00108 /* fix pointers to ionization edges and frequency array 00109 * calls iso_create 00110 * this routine only returns if this is a later call of code */ 00111 ContCreatePointers(); 00112 00113 /* Badnell_rec_init This code is written by Terry Yun, 2005 * 00114 * It reads dielectronic recombination rate coefficient fits into 3D arrays */ 00115 Badnell_rec_init(); 00116 00117 ASSERT(lgElemsConserved()); 00118 00119 /* set continuum normalization, continuum itself, and ionization stage limits */ 00120 ContSetIntensity(); 00121 00122 ASSERT(lgElemsConserved()); 00123 00124 /* print header */ 00125 PrtHeader(); 00126 00127 /* this is an option to stop after initial setup */ 00128 if( noexec.lgNoExec ) 00129 return false; 00130 00131 /* guess some outward optical depths, set inward optical depths, 00132 * also calls RT_line_all so escape probabilities ok before printout of trace */ 00133 RT_tau_init(); 00134 00135 /* generate initial set of opacities, but only if this is the first call 00136 * in this core load 00137 * grains only exist AFTER this routine exits */ 00138 OpacityCreateAll(); 00139 00140 ASSERT(lgElemsConserved()); 00141 00142 /* this checks that various parts of the code still work properly */ 00143 SanityCheck("begin"); 00144 00145 /* option to recover state from previous calculation */ 00146 if( state.lgGet_state ) 00147 state_get_put( "get" ); 00148 00149 ASSERT(lgElemsConserved()); 00150 00151 /* find the initial temperature, punt if initial conditions outside range of code, 00152 * abort condition set by flag lgAbort */ 00153 if( ConvInitSolution() ) 00154 { 00155 // create line stacks for possible printout before bailing out 00156 LineStackCreate(); 00157 BadStart(); 00158 return true; 00159 } 00160 // create line stacks ... 00161 LineStackCreate(); 00162 00163 /* set thickness of first zone */ 00164 radius_first(); 00165 00166 /* find thickness of next zone */ 00167 if( radius_next() ) 00168 { 00169 BadStart(); 00170 return true; 00171 } 00172 00173 /* set up some zone variables, correct continuum for sphericity, 00174 * after return, radius is equal to the inner radius, not outer radius 00175 * of this zone */ 00176 ZoneStart("init"); 00177 00178 /* print all abundances, gas phase and grains, in abundances.c */ 00179 /* >>chng 06 mar 05, move AbundancesPrt() after ZoneStart("init") so that 00180 * GrnVryDpth() gives correct values for the inner face of the cloud, PvH */ 00181 AbundancesPrt(); 00182 00183 /* this is an option to stop after printing header only */ 00184 if( prt.lgOnlyHead ) 00185 return false; 00186 00187 plot("FIRST"); 00188 00189 /* outer loop is over number of iterations 00190 * >>chng 05 mar 12, chng from test on true to not aborting */ 00191 while( !lgAbort ) 00192 { 00193 IterStart(); 00194 nzone = 0; 00195 fnzone = 0.; 00196 00197 /* loop over zones across cloud for this iteration, 00198 * iter_end_check checks whether model is complete and this iteration is done 00199 * returns true is current iteration is complete 00200 * calls PrtZone to print zone results */ 00201 while( !iter_end_check() ) 00202 { 00203 /* the zone number, 0 during search phase, first zone is 1 */ 00204 ++nzone; 00205 /* this is the zone number plus the number of calls to bottom solvers 00206 * from top pressure solver, divided by 100 */ 00207 fnzone = (double)nzone; 00208 00209 /* use changes in opacity, temperature, ionization, to fix new dr for next zone */ 00210 /* >>chng 03 oct 29, move radius_next up to here from below, so that 00211 * precise correct zone thickness will be used in current zone, so that 00212 * column density and thickness will be exact 00213 * abort condition is possible */ 00214 if( radius_next() ) 00215 break; 00216 00217 /* following sets zone thickness, drad, to drnext */ 00218 /* set variable dealing with new radius, in zones.c */ 00219 ZoneStart("incr"); 00220 00221 /* converge the pressure-temperature-ionization solution for this zone 00222 * NB ignoring return value - should be ok (return 1 for abort) 00223 * we can't break here in case of abort since there is still housekeeping 00224 * that must be done in following routines */ 00225 ConvPresTempEdenIoniz(); 00226 00227 /* generate diffuse emission from this zone, add to outward & reflected beams */ 00228 RT_diffuse(); 00229 00230 /* do work associated with incrementing this radius, 00231 * total attenuation and dilution of radiation fields takes place here, 00232 * reflected continuum incremented here 00233 * various mean and extremes of solution are also remembered here */ 00234 radius_increment(); 00235 00236 /* attenuation of diffuse and beamed continua */ 00237 RT_continuum(); 00238 00239 /* increment optical depths 00240 * final evaluation of line rates and cooling */ 00241 RT_tau_inc(); 00242 00243 /* fill in emission line array, adds outward lines */ 00244 /* >>chng 99 dec 29, moved to here from below RT_tau_inc, 00245 * lines adds lines to outward beam, 00246 * and these are attenuated in radius_increment */ 00247 lines(); 00248 00249 /* possibly save out some results from this zone */ 00250 SaveDo("MIDL"); 00251 00252 /* do some things to finish out this zone */ 00253 ZoneEnd(); 00254 00255 // default false due to slow time - set true with set check energy every zone 00256 // to allow per zone check of energy conservation, relatively slow 00257 // when chianti is fully enabled 00258 if( continuum.lgCheckEnergyEveryZone ) 00259 { 00260 /* Ensure that energy has been conserved. Complain and punt if not */ 00261 if( !lgConserveEnergy() ) 00262 { 00263 fprintf( ioQQQ, " PROBLEM DISASTER Energy was not conserved at zone %li\n", nzone ); 00264 ShowMe(); 00265 lgAbort = true; 00266 } 00267 } 00268 } 00269 /* end loop over zones */ 00270 00271 /* close out this iteration, in startenditer.c */ 00272 IterEnd(); 00273 00274 /* print out some comments, generate warning and cautions*/ 00275 PrtComment(); 00276 00277 /* save stuff only needed on completion of this iteration */ 00278 SaveDo("LAST" ); 00279 00280 /* second call to plot routine, to complete plots for this iteration */ 00281 plot("SECND"); 00282 00283 /* print out the results */ 00284 PrtFinal(); 00285 00286 /*ConvIterCheck check whether model has converged or more iterations 00287 * are needed - implements the iter to convergence command 00288 * acts by setting iterations.itermx to iteration if converged*/ 00289 ConvIterCheck(); 00290 00291 /* option to save state */ 00292 if( state.lgPut_state ) 00293 state_get_put( "put" ); 00294 00295 /* >>chng 06 mar 03 move block to here, had been after PrtFinal, 00296 * so that save state will save reset optical depths */ 00297 /* this is the normal exit, occurs if we reached limit to number of iterations, 00298 * or if code has set busted */ 00299 /* >>chng 06 mar 18, add flag for time dependent simulation completed */ 00300 if( iteration > iterations.itermx || lgAbort || dynamics.lgStatic_completed ) 00301 break; 00302 00303 /* reset limiting and initial optical depth variables */ 00304 RT_tau_reset(); 00305 00306 /* increment iteration counter */ 00307 ++iteration; 00308 00309 /* reinitialize some variables to initial conditions at previous first zone 00310 * routine in startenditer.c */ 00311 IterRestart(); 00312 00313 /* reset zone number to 0 - done here since this is only routine that sets nzone */ 00314 nzone = 0; 00315 fnzone = 0.; 00316 00317 ZoneStart("init"); 00318 00319 /* find new initial temperature, punt if initial conditions outside range of code, 00320 * if ConvInitSolution() fails, lgAbort will always be true, so we check that */ 00321 if( ConvInitSolution() ) 00322 break; 00323 } 00324 00325 CloseSaveFiles( false ); 00326 00327 /* this checks that various parts of the code worked properly */ 00328 SanityCheck("final"); 00329 00330 if( called.lgTalk ) 00331 { 00332 fprintf(ioQQQ,"---------------Convergence statistics---------------\n"); 00333 fprintf(ioQQQ,"%10.3g mean iterations/state convergence\n",((double)conv.getCounter(CONV_BASE_LOOPS))/ 00334 (MAX2((double)conv.getCounter(CONV_BASE_CALLS),1.0))); 00335 fprintf(ioQQQ,"%10.3g mean cx acceleration loops/iteration\n",((double)conv.getCounter(CONV_BASE_ACCELS))/ 00336 (MAX2((double)conv.getCounter(CONV_BASE_LOOPS),1.0))); 00337 fprintf(ioQQQ,"%10.3g mean iso convergence loops/ion solve\n",((double)conv.getCounter(ISO_LOOPS))/ 00338 (MAX2((double)conv.getCounter(ION_SOLVES),1.0))); 00339 fprintf(ioQQQ,"%10.3g mean steps/chemistry solve\n",((double)conv.getCounter(MOLE_SOLVE_STEPS))/ 00340 (MAX2((double)conv.getCounter(MOLE_SOLVE),1.0))); 00341 fprintf(ioQQQ,"%10.3g mean step length searches/chemistry step\n",((double)conv.getCounter(NEWTON_LOOP))/ 00342 (MAX2((double)conv.getCounter(NEWTON),1.0))); 00343 fprintf(ioQQQ,"----------------------------------------------------\n\n"); 00344 } 00345 00346 /* check whether any asserts were present and failed. 00347 * return is true if ok, false if not. routine also checks 00348 * number of warnings and returns false if not ok */ 00349 lgOK = lgCheckMonitors(ioQQQ); 00350 00351 if( lgOK && !warnings.lgWarngs && !lgAbort ) 00352 { 00353 /* no failed asserts or warnings */ 00354 return false; 00355 } 00356 else 00357 { 00358 /* there were failed asserts or warnings */ 00359 return true; 00360 } 00361 } 00362 00363 /*BadStart announce that things are so bad the calculation cannot even start */ 00364 STATIC void BadStart() 00365 { 00366 char chLine[INPUT_LINE_LENGTH]; 00367 00368 DEBUG_ENTRY( "BadStart()" ); 00369 00370 /* initialize the line saver */ 00371 wcnint(); 00372 sprintf( warnings.chRgcln[0], " Calculation stopped because initial conditions out of bounds." ); 00373 sprintf( chLine, " W-Calculation could not begin." ); 00374 warnin(chLine); 00375 00376 if( grid.lgGrid ) 00377 { 00378 /* possibly save out some results from this zone when doing grid 00379 * since grid save files expect something at this grid point */ 00380 SaveDo("MIDL"); 00381 SaveDo("LAST" ); 00382 } 00383 return; 00384 }