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