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