cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cloudy.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*cloudy the main routine, this IS Cloudy, return 0 normal exit, 1 error exit,
4  * called by maincl when used as standalone program */
5 /*BadStart announce that things are so bad the calculation cannot even start */
6 #include "cddefines.h"
7 #include "punch.h"
8 #include "noexec.h"
9 #include "lines.h"
10 #include "abund.h"
11 #include "continuum.h"
12 #include "warnings.h"
13 #include "atmdat.h"
14 #include "prt.h"
15 #include "conv.h"
16 #include "parse.h"
17 #include "dynamics.h"
18 #include "init.h"
19 #include "opacity.h"
20 #include "state.h"
21 #include "rt.h"
22 #include "assertresults.h"
23 #include "zones.h"
24 #include "iterations.h"
25 #include "plot.h"
26 #include "radius.h"
27 #include "grid.h"
28 #include "cloudy.h"
29 #include "ionbal.h"
30 
31 /*BadStart announce that things are so bad the calculation cannot even start */
32 STATIC void BadStart(void);
33 
34 /* returns 1 if disaster strikes, 0 if everything appears ok */
35 bool cloudy(void)
36 {
37  bool lgOK,
38  /* will be used to remember why we stopped,
39  * return error exit in some cases */
40  lgBadEnd;
41 
42  /*
43  * this is the main routine to drive the modules that compute a model
44  * when cloudy is used as a stand-alone program
45  * the main program (maincl) calls cdInit then cdDrive
46  * this sub is called by cdDrive which returns upon exiting
47  *
48  * this routine uses the following variables:
49  *
50  * nzone
51  * this is the zone number, and is incremented here
52  * is zero during search phase, 1 for first zone at illuminated face
53  * logical function iter_end_check returns a true condition if NZONE reaches
54  * NEND(ITER), the limit to the number of zones in this iteration
55  * nzone is totally controlled in this subroutine
56  *
57  * iteration
58  * this is the iteration number, it is 1 for the first iteration
59  * and is incremented in this subroutine at the end of each iteration
60  *
61  * iterations.itermx
62  * this is the limit to the number of iterations, and is entered
63  * by the user
64  * This routine returns when iteration > iterations.itermx
65  */
66 
67  /* nzone is zero while initial search for conditions takes place
68  * and is incremented to 1 for start of first zone */
69  nzone = 0;
70  fnzone = 0.;
71 
72  /* iteration is iteration number, calculation is complete when iteration > iterations.itermx */
73  iteration = 1;
74 
75  /* flag for error exit */
76  lgBadEnd = false;
77 
78  /* this initializes variables at the start of each simulation
79  * in a grid, before the parser is called - this must set any values
80  * that may be changed by the command parser */
82 
83  /* scan in and parse input data */
84  ParseCommands();
85 
86  /* fix abundances of elements, in abundances.cpp */
87  AbundancesSet();
88 
89  /* one time creation of some variables */
91 
92  /* initialize vars that can only be done after parsing the commands
93  * sets up CO network since this depends on which elements are on
94  * and whether chemistry is enabled */
96 
97  /* ContCreateMesh calls fill to set up continuum energy mesh if first call,
98  * otherwise reset to original mesh.
99  * This is AFTER ParseCommands so that
100  * path and mesh size can be set with commands before mesh is set */
101  ContCreateMesh();
102 
103  /* create several data sets by allocating memory and reading in data files,
104  * but only if this is first call */
105  atmdat_readin();
106 
107  /* fix pointers to ionization edges and frequency array
108  * calls iso_create
109  * this routine only returns if this is a later call of code */
111 
112  /* Badnell_rec_init This code is written by Terry Yun, 2005 *
113  * It reads dielectronic recombination rate coefficient fits into 3D arrays */
115 
116  /* set continuum normalization, continuum itself, and ionization stage limits */
117  if( ContSetIntensity() )
118  {
119  /* this happens when disaster strikes in the initial setup of the continuum intensity array,
120  * BadStart is in this file, below */
121  BadStart();
122  return(true);
123  }
124 
125  /* print header */
126  PrtHeader();
127 
128  /* this is an option to stop after initial setup */
129  if( noexec.lgNoExec )
130  return(false);
131 
132  /* guess some outward optical depths, set inward optical depths,
133  * also calls RT_line_all so escape probabilities ok before printout of trace */
134  RT_tau_init();
135 
136  /* generate initial set of opacities, but only if this is the first call
137  * in this core load
138  * grains only exist AFTER this routine exits */
140 
141  /* this checks that various parts of the code still work properly */
142  SanityCheck("begin");
143 
144  /* option to recover state from previous calculation */
145  if( state.lgGet_state )
146  state_get_put( "get" );
147 
148  /* find the initial temperature, punt if initial conditions outside range of code,
149  * abort condition set by flag lgAbort */
150  if( ConvInitSolution() )
151  {
152  BadStart();
153  return(true);
154  }
155 
156  /* set thickness of first zone */
157  radius_first();
158 
159  /* find thickness of next zone */
160  if( radius_next() )
161  {
162  BadStart();
163  return(true);
164  }
165 
166  /* set up some zone variables, correct continuum for sphericity,
167  * after return, radius is equal to the inner radius, not outer radius
168  * of this zone */
169  ZoneStart("init");
170 
171  /* print all abundances, gas phase and grains, in abundances.c */
172  /* >>chng 06 mar 05, move AbundancesPrt() after ZoneStart("init") so that
173  * GrnVryDpth() gives correct values for the inner face of the cloud, PvH */
174  AbundancesPrt();
175 
176  /* this is an option to stop after printing header only */
177  if( prt.lgOnlyHead )
178  return(false);
179 
180  plot("FIRST");
181 
182  /* outer loop is over number of iterations
183  * >>chng 05 mar 12, chng from test on true to not aborting */
184  while( !lgAbort )
185  {
186  IterStart();
187  nzone = 0;
188  fnzone = 0.;
189 
190  /* loop over zones across cloud for this iteration,
191  * iter_end_check checks whether model is complete and this iteration is done
192  * returns true is current iteration is complete */
193  while( !iter_end_check() )
194  {
195  /* the zone number, 0 during search phase, first zone is 1 */
196  ++nzone;
197  /* this is the zone number plus the number of calls to bottom solvers
198  * from top pressure solver, divided by 100 */
199  fnzone = (double)nzone;
200 
201  /* use changes in opacity, temperature, ionization, to fix new dr for next zone */
202  /* >>chng 03 oct 29, move radius_next up to here from below, so that
203  * precise correct zone thickness will be used in current zone, so that
204  * column density and thickness will be exact
205  * abort condition is possible */
206  if( radius_next() )
207  break;
208 
209  /* following sets zone thickness, drad, to drnext */
210  /* set variable dealing with new radius, in zones.c */
211  ZoneStart("incr");
212 
213  /* converge the pressure-temperature-ionization solution for this zone
214  * NB ignoring return value - should be ok (return 1 for disaster) */
216 
217  /* generate diffuse emission from this zone, add to outward & reflected beams */
218  RT_diffuse();
219 
220  /* do work associated with incrementing this radius,
221  * total attenuation and dilution of radiation fields takes place here,
222  * reflected continuum incremented here
223  * various mean and extremes of solution are also remembered here */
225 
226  /* increment optical depths */
227  RT_tau_inc();
228 
229  /* fill in emission line array, adds outward lines */
230  /* >>chng 99 dec 29, moved to here from below RT_tau_inc,
231  * lines adds lines to outward beam,
232  * and these are attenuated in radius_increment */
233  lines();
234 
235  /* possibly punch out some results from this zone */
236  PunchDo("MIDL");
237 
238  /* do some things to finish out this zone */
239  ZoneEnd();
240  }
241  /* end loop over zones */
242 
243  /* close out this iteration, in startenditer.c */
244  IterEnd();
245 
246  /* print out some comments, generate warning and cautions*/
247  PrtComment();
248 
249  /* punch stuff only needed on completion of this iteration */
250  PunchDo("LAST" );
251 
252  /* second call to plot routine, to complete plots for this iteration */
253  plot("SECND");
254 
255  /* print out the results */
256  PrtFinal();
257 
258  /*ConvIterCheck check whether model has converged or whether more iterations
259  * are needed - implements the iter to convergence command */
260  ConvIterCheck();
261 
262  /* reset limiting and initial optical depth variables */
263  RT_tau_reset();
264 
265  /* option to save state */
266  if( state.lgPut_state )
267  state_get_put( "put" );
268 
269  /* >>chng 06 mar 03 move block to here, had been after PrtFinal,
270  * so that save state will save reset optical depths */
271  /* this is the normal exit, occurs if we reached limit to number of iterations,
272  * or if code has set busted */
273  /* >>chng 06 mar 18, add flag for time dependent simulation completed */
275  break;
276 
277  /* increment iteration counter */
278  ++iteration;
279 
280  /* reinitialize some variables to initial conditions at previous first zone
281  * routine in startenditer.c */
282  IterRestart();
283 
284  /* reset zone number to 0 - done here since this is only routine that sets nzone */
285  nzone = 0;
286  fnzone = 0.;
287 
288  /* now redo escape probabilities */
289  RT_line_all(true , false);
290 
291  ZoneStart("init");
292 
293  /* find new initial temperature, punt if initial conditions outside range of code,
294  * abort condition set by flag lgAbort */
295  if( ConvInitSolution() )
296  {
297  lgBadEnd = true;
298  break;
299  }
300  }
301 
302  ClosePunchFiles( false );
303 
304  /* this checks that various parts of the code worked properly */
305  SanityCheck("final");
306 
307  /* check whether any asserts were present and failed.
308  * return is true if ok, false if not. routine also checks
309  * number of warnings and returns false if not ok */
310  lgOK = lgCheckAsserts(ioQQQ);
311 
312  if( lgOK && !warnings.lgWarngs && !lgBadEnd)
313  {
314  /* no failed asserts or warnings */
315  return(false);
316  }
317  else
318  {
319  /* there were failed asserts or warnings */
320  return(true);
321  }
322 }
323 
324 /*BadStart announce that things are so bad the calculation cannot even start */
325 STATIC void BadStart(void)
326 {
327  char chLine[INPUT_LINE_LENGTH];
328 
329  DEBUG_ENTRY( "BadStart()" );
330 
331  /* initialize the line saver */
332  wcnint();
333  sprintf( warnings.chRgcln[0], " Calculation stopped because initial conditions out of bounds." );
334  sprintf( chLine, " W-Calculation could not begin." );
335  warnin(chLine);
336 
337  if( grid.lgGrid )
338  {
339  /* possibly punch out some results from this zone when doing grid
340  * since grid punch files expect something at this grid point */
341  PunchDo("MIDL");
342  PunchDo("LAST" );
343  }
344  return;
345 }

Generated for cloudy by doxygen 1.8.3.1