00001 /* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 00002 * For conditions of distribution and use see copyright notice in license.txt */ 00003 /*cloudy the main routine, this IS Cloudy, ret 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 "opacity.h" 00019 #include "state.h" 00020 #include "rt.h" 00021 #include "assertresults.h" 00022 #include "zones.h" 00023 #include "iterations.h" 00024 #include "plot.h" 00025 #include "radius.h" 00026 #include "optimize.h" 00027 #include "grid.h" 00028 #include "mole.h" 00029 #include "cloudy.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 /* scan in and parse input data */ 00079 ParseCommands(); 00080 00081 CO_Init(); 00082 CO_zero(); 00083 00084 /* >> chng 06 Nov 24 rjrw: Move reaction definitions here so they can be controlled by parsed commands */ 00085 CO_create_react(); 00086 00087 /* ContCreateMesh calls fill to set up continuum energy mesh if first call, 00088 * otherwise reset to original mesh. 00089 * This is AFTER ParseCommands so that 00090 * path and mesh size can be set with commands before mesh is set */ 00091 ContCreateMesh(); 00092 00093 /* create several data sets by allocating memory and reading in data files, 00094 * but only if this is first call */ 00095 atmdat_readin(); 00096 00097 /* fix pointers to ionization edges and frequency array 00098 * calls iso_create 00099 * this routine only returns if this is a later call of code */ 00100 ContCreatePointers(); 00101 00102 /* fix abundances of elements, in abundances.c */ 00103 AbundancesSet(); 00104 00105 /* set continuum normalization, continuum itself, and ionization stage limits */ 00106 if( ContSetIntensity() ) 00107 { 00108 /* this happens when disaster strikes in the initial setup of the continuum intensity array, 00109 * BadStart is in this file, below */ 00110 BadStart(); 00111 return(true); 00112 } 00113 00114 /* print header */ 00115 PrtHeader(); 00116 00117 /* this is an option to stop after initial setup */ 00118 if( noexec.lgNoExec ) 00119 return(false); 00120 00121 /* guess some outward optical depths, set inward optical depths, 00122 * also calls RT_line_all so escape probabilities ok before printout of trace */ 00123 RT_tau_init(); 00124 00125 /* generate inital set of opacities, but only if this is the first call 00126 * in this core load 00127 * grains only exist AFTER this routine exits */ 00128 OpacityCreateAll(); 00129 00130 /* this checks that various parts of the code still work properly */ 00131 SanityCheck("begin"); 00132 00133 /* option to recover state from previous calculation */ 00134 if( state.lgGet_state ) 00135 state_get_put( "get" ); 00136 00137 /* find the initial temperature, punt if initial conditions outside range of code, 00138 * abort condition set by flag lgAbort */ 00139 if( ConvInitSolution() ) 00140 { 00141 BadStart(); 00142 return(true); 00143 } 00144 00145 /* set thickness of first zone */ 00146 radius_first(); 00147 00148 /* find thickness of next zone */ 00149 if( radius_next() ) 00150 { 00151 BadStart(); 00152 return(true); 00153 } 00154 00155 /* set up some zone variables, correct continuum for sphericity, 00156 * after return, radius is equal to the inner radius, not outer radius 00157 * of this zone */ 00158 ZoneStart("init"); 00159 00160 /* print all abundances, gas phase and grains, in abundances.c */ 00161 /* >>chng 06 mar 05, move AbundancesPrt() after ZoneStart("init") so that 00162 * GrnVryDpth() gives correct values for the inner face of the cloud, PvH */ 00163 AbundancesPrt(); 00164 00165 /* this is an option to stop after printing header only */ 00166 if( prt.lgOnlyHead ) 00167 return(false); 00168 00169 plot("FIRST"); 00170 00171 /* outer loop is over number of iterations 00172 * >>chng 05 mar 12, chng from test on true to not aborting */ 00173 while( !lgAbort ) 00174 { 00175 IterStart(); 00176 nzone = 0; 00177 fnzone = 0.; 00178 00179 /* loop over zones across cloud for this iteration, 00180 * iter_end_check checks whether model is complete and this iteration is done 00181 * returns true is current iteration is complete */ 00182 while( !iter_end_check() ) 00183 { 00184 /* the zone number, 0 during search phase, first zone is 1 */ 00185 ++nzone; 00186 /* this is the zone number plus the number of calls to bottom solvers 00187 * from top pressure solver, divided by 100 */ 00188 fnzone = (double)nzone; 00189 00190 /* use changes in opacity, temperature, ionization, to fix new dr for next zone */ 00191 /* >>chng 03 oct 29, move radius_next up to here from below, so that 00192 * precise correct zone thickness will be used in current zone, so that 00193 * column density and thickness will be exact 00194 * abort condition is possible */ 00195 if( radius_next() ) 00196 continue; 00197 00198 /* following sets zone thickness, drad, to drnext */ 00199 /* set variable dealing with new radius, in zones.c */ 00200 ZoneStart("incr"); 00201 00202 /* converge the pressure-temperature-ionization solution for this zone 00203 * NB ignoring return value - should be ok (ret 1 for disaster) */ 00204 ConvPresTempEdenIoniz(); 00205 00206 /* generate diffuse emission from this zone, add to outward & reflected beams */ 00207 RT_diffuse(); 00208 00209 /* do work associated with incrementing this radius, 00210 * total attenuation and dilution of radiation fields takes place here, 00211 * reflected continuum incremented here 00212 * various mean and extremes of solution are also remembered here */ 00213 radius_increment(); 00214 00215 /* increment optical depths */ 00216 RT_tau_inc(); 00217 00218 /* >>chng 99 may 04, move lines here from after RT_diffuse */ 00219 /* fill in emission line array, adds outward lines */ 00220 /* >>chng 99 dec 29, moved to here from below RT_tau_inc, 00221 * lines adds lines to outward beam, 00222 * and these are attenuated in radius_increment */ 00223 lines(); 00224 00225 /* possibly punch out some results from this zone */ 00226 PunchDo("MIDL"); 00227 00228 /* do some things to finish out this zone */ 00229 ZoneEnd(); 00230 } 00231 /* end loop over zones */ 00232 00233 /* close out this iteration, in startenditer.c */ 00234 IterEnd(); 00235 00236 /* print out some comments, generate warning and cautions*/ 00237 PrtComment(); 00238 00239 GridGather(); 00240 00241 /* punch stuff only needed on completion of this iteration */ 00242 PunchDo("LAST" ); 00243 00244 /* second call to plot routine, to complete plots for this iteration */ 00245 plot("SECND"); 00246 00247 /* print out the results */ 00248 PrtFinal(); 00249 00250 /*ConvIterCheck check whether model has converged or whether more iterations 00251 * are needed - implements the iter to convergence comnd */ 00252 ConvIterCheck(); 00253 00254 /* reset limiting and initial optical depth variables */ 00255 RT_tau_reset(); 00256 00257 /* option to save state */ 00258 if( state.lgPut_state ) 00259 state_get_put( "put" ); 00260 00261 /* >>chng 06 mar 03 move block to here, had been after PrtFinal, 00262 * so that save state will save reset optical depths */ 00263 /* this is the normal exit, occurs if we reached limit to number of iterations, 00264 * or if code has set busted */ 00265 /* >>chng 06 mar 18, add flag for time dependent simulation completed */ 00266 if( iteration > iterations.itermx || lgAbort || dynamics.lgStatic_completed ) 00267 { 00268 break; 00269 } 00270 00271 /* increment iteration counter */ 00272 ++iteration; 00273 00274 /* reinitialize some variables to initial conditions at previous first zone 00275 * routine in startenditer.c */ 00276 IterRestart(); 00277 00278 /* reset zone number to 0 - done here since this is only routine that sets nzone */ 00279 nzone = 0; 00280 fnzone = 0.; 00281 00282 /* now redo escape probabilities */ 00283 RT_line_all(true , false); 00284 00285 ZoneStart("init"); 00286 00287 /* find new initial temperature, punt if initial conditions outside range of code, 00288 * abort condition set by flag lgAbort */ 00289 if( ConvInitSolution() ) 00290 { 00291 lgBadEnd = true; 00292 break; 00293 } 00294 } 00295 00296 ClosePunchFiles(); 00297 00298 /* this checks that various parts of the code worked properly */ 00299 SanityCheck("final"); 00300 00301 /* check whether any asserts were present and failed. 00302 * return is true if ok, false if not. routine also checks 00303 * number of warnings and returns false if not ok */ 00304 lgOK = lgCheckAsserts(ioQQQ); 00305 00306 if( lgOK && !warnings.lgWarngs && !lgBadEnd) 00307 { 00308 /* no failed asserts or warnings */ 00309 return(false); 00310 } 00311 else 00312 { 00313 /* there were failed asserts or warnings */ 00314 return(true); 00315 } 00316 00317 } 00318 00319 /*BadStart announce that things are so bad the calculation cannot even start */ 00320 static void BadStart(void) 00321 { 00322 char chLine[INPUT_LINE_LENGTH]; 00323 00324 DEBUG_ENTRY( "BadStart()" ); 00325 00326 /* initialize the line saver */ 00327 wcnint(); 00328 sprintf( warnings.chRgcln[0], " Calculation stopped because initial conditions out of bounds." ); 00329 sprintf( chLine, " W-Calculation could not begin." ); 00330 warnin(chLine); 00331 00332 DEBUG_EXIT( "BadStart()" ); 00333 return; 00334 } 00335