cloudy trunk
Loading...
Searching...
No Matches
cloudy.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 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 "save.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 "monitor_results.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#include "dense.h"
31#include "energy.h"
32#include "called.h"
33
34/*BadStart announce that things are so bad the calculation cannot even start */
35STATIC void BadStart();
36
37/* returns 1 if disaster strikes, 0 if everything appears ok */
38bool cloudy()
39{
40 DEBUG_ENTRY( "cloudy()" );
41
42 bool lgOK;
43
44 /*
45 * this is the main routine to drive the modules that compute a model
46 * when cloudy is used as a stand-alone program
47 * the main program (maincl) calls cdInit then cdDrive
48 * this sub is called by cdDrive which returns upon exiting
49 *
50 * this routine uses the following variables:
51 *
52 * nzone
53 * this is the zone number, and is incremented here
54 * is zero during search phase, 1 for first zone at illuminated face
55 * logical function iter_end_check returns a true condition if NZONE reaches
56 * NEND(ITER), the limit to the number of zones in this iteration
57 * nzone is totally controlled in this subroutine
58 *
59 * iteration
60 * this is the iteration number, it is 1 for the first iteration
61 * and is incremented in this subroutine at the end of each iteration
62 *
63 * iterations.itermx
64 * this is the limit to the number of iterations, and is entered
65 * by the user
66 * This routine returns when iteration > iterations.itermx
67 */
68
69 /* nzone is zero while initial search for conditions takes place
70 * and is incremented to 1 for start of first zone */
71 nzone = 0;
72 fnzone = 0.;
73
74 /* iteration is iteration number, calculation is complete when iteration > iterations.itermx */
75 iteration = 1;
76
77 /* this initializes variables at the start of each simulation
78 * in a grid, before the parser is called - this must set any values
79 * that may be changed by the command parser */
81
82 /* scan in and parse input data */
84
85 /* fix abundances of elements, in abundances.cpp */
87
89
90 /* one time creation of some variables */
92
93 /* initialize vars that can only be done after parsing the commands
94 * sets up CO network since this depends on which elements are on
95 * and whether chemistry is enabled */
97
98 /* ContCreateMesh calls fill to set up continuum energy mesh if first call,
99 * otherwise reset to original mesh.
100 * This is AFTER ParseCommands so that
101 * path and mesh size can be set with commands before mesh is set */
103
104 /* create several data sets by allocating memory and reading in data files,
105 * but only if this is first call */
107
108 /* fix pointers to ionization edges and frequency array
109 * calls iso_create
110 * this routine only returns if this is a later call of code */
112
113 /* Badnell_rec_init This code is written by Terry Yun, 2005 *
114 * It reads dielectronic recombination rate coefficient fits into 3D arrays */
116
118
119 /* set continuum normalization, continuum itself, and ionization stage limits */
121
123
124 /* print header */
125 PrtHeader();
126
127 /* this is an option to stop after initial setup */
128 if( noexec.lgNoExec )
129 return false;
130
131 /* guess some outward optical depths, set inward optical depths,
132 * also calls RT_line_all so escape probabilities ok before printout of trace */
133 RT_tau_init();
134
135 /* generate initial set of opacities, but only if this is the first call
136 * in this core load
137 * grains only exist AFTER this routine exits */
139
141
142 /* this checks that various parts of the code still work properly */
143 SanityCheck("begin");
144
145 /* option to recover state from previous calculation */
146 if( state.lgGet_state )
147 state_get_put( "get" );
148
150
151 /* find the initial temperature, punt if initial conditions outside range of code,
152 * abort condition set by flag lgAbort */
153 if( ConvInitSolution() )
154 {
155 // create line stacks for possible printout before bailing out
157 BadStart();
158 return true;
159 }
160 // create line stacks ...
162
163 /* set thickness of first zone */
164 radius_first();
165
166 /* find thickness of next zone */
167 if( radius_next() )
168 {
169 BadStart();
170 return true;
171 }
172
173 /* set up some zone variables, correct continuum for sphericity,
174 * after return, radius is equal to the inner radius, not outer radius
175 * of this zone */
176 ZoneStart("init");
177
178 /* print all abundances, gas phase and grains, in abundances.c */
179 /* >>chng 06 mar 05, move AbundancesPrt() after ZoneStart("init") so that
180 * GrnVryDpth() gives correct values for the inner face of the cloud, PvH */
182
183 /* this is an option to stop after printing header only */
184 if( prt.lgOnlyHead )
185 return false;
186
187 plot("FIRST");
188
189 /* outer loop is over number of iterations
190 * >>chng 05 mar 12, chng from test on true to not aborting */
191 while( !lgAbort )
192 {
193 IterStart();
194 nzone = 0;
195 fnzone = 0.;
196
197 /* loop over zones across cloud for this iteration,
198 * iter_end_check checks whether model is complete and this iteration is done
199 * returns true is current iteration is complete
200 * calls PrtZone to print zone results */
201 while( !iter_end_check() )
202 {
203 /* the zone number, 0 during search phase, first zone is 1 */
204 ++nzone;
205 /* this is the zone number plus the number of calls to bottom solvers
206 * from top pressure solver, divided by 100 */
207 fnzone = (double)nzone;
208
209 /* use changes in opacity, temperature, ionization, to fix new dr for next zone */
210 /* >>chng 03 oct 29, move radius_next up to here from below, so that
211 * precise correct zone thickness will be used in current zone, so that
212 * column density and thickness will be exact
213 * abort condition is possible */
214 if( radius_next() )
215 break;
216
217 /* following sets zone thickness, drad, to drnext */
218 /* set variable dealing with new radius, in zones.c */
219 ZoneStart("incr");
220
221 /* converge the pressure-temperature-ionization solution for this zone
222 * NB ignoring return value - should be ok (return 1 for abort)
223 * we can't break here in case of abort since there is still housekeeping
224 * that must be done in following routines */
226
227 /* generate diffuse emission from this zone, add to outward & reflected beams */
228 RT_diffuse();
229
230 /* do work associated with incrementing this radius,
231 * total attenuation and dilution of radiation fields takes place here,
232 * reflected continuum incremented here
233 * various mean and extremes of solution are also remembered here */
235
236 /* attenuation of diffuse and beamed continua */
237 RT_continuum();
238
239 /* increment optical depths
240 * final evaluation of line rates and cooling */
241 RT_tau_inc();
242
243 /* fill in emission line array, adds outward lines */
244 /* >>chng 99 dec 29, moved to here from below RT_tau_inc,
245 * lines adds lines to outward beam,
246 * and these are attenuated in radius_increment */
247 lines();
248
249 /* possibly save out some results from this zone */
250 SaveDo("MIDL");
251
252 /* do some things to finish out this zone */
253 ZoneEnd();
254
255 // default false due to slow time - set true with set check energy every zone
256 // to allow per zone check of energy conservation, relatively slow
257 // when chianti is fully enabled
258 if( continuum.lgCheckEnergyEveryZone )
259 {
260 /* Ensure that energy has been conserved. Complain and punt if not */
261 if( !lgConserveEnergy() )
262 {
263 fprintf( ioQQQ, " PROBLEM DISASTER Energy was not conserved at zone %li\n", nzone );
264 ShowMe();
265 lgAbort = true;
266 }
267 }
268 }
269 /* end loop over zones */
270
271 /* close out this iteration, in startenditer.c */
272 IterEnd();
273
274 /* print out some comments, generate warning and cautions*/
275 PrtComment();
276
277 /* save stuff only needed on completion of this iteration */
278 SaveDo("LAST" );
279
280 /* second call to plot routine, to complete plots for this iteration */
281 plot("SECND");
282
283 /* print out the results */
284 PrtFinal();
285
286 /*ConvIterCheck check whether model has converged or more iterations
287 * are needed - implements the iter to convergence command
288 * acts by setting iterations.itermx to iteration if converged*/
290
291 /* option to save state */
292 if( state.lgPut_state )
293 state_get_put( "put" );
294
295 /* >>chng 06 mar 03 move block to here, had been after PrtFinal,
296 * so that save state will save reset optical depths */
297 /* this is the normal exit, occurs if we reached limit to number of iterations,
298 * or if code has set busted */
299 /* >>chng 06 mar 18, add flag for time dependent simulation completed */
300 if( iteration > iterations.itermx || lgAbort || dynamics.lgStatic_completed )
301 break;
302
303 /* reset limiting and initial optical depth variables */
304 RT_tau_reset();
305
306 /* increment iteration counter */
307 ++iteration;
308
309 /* reinitialize some variables to initial conditions at previous first zone
310 * routine in startenditer.c */
311 IterRestart();
312
313 /* reset zone number to 0 - done here since this is only routine that sets nzone */
314 nzone = 0;
315 fnzone = 0.;
316
317 ZoneStart("init");
318
319 /* find new initial temperature, punt if initial conditions outside range of code,
320 * if ConvInitSolution() fails, lgAbort will always be true, so we check that */
321 if( ConvInitSolution() )
322 break;
323 }
324
325 CloseSaveFiles( false );
326
327 /* this checks that various parts of the code worked properly */
328 SanityCheck("final");
329
330 if( called.lgTalk )
331 {
332 fprintf(ioQQQ,"---------------Convergence statistics---------------\n");
333 fprintf(ioQQQ,"%10.3g mean iterations/state convergence\n",((double)conv.getCounter(CONV_BASE_LOOPS))/
334 (MAX2((double)conv.getCounter(CONV_BASE_CALLS),1.0)));
335 fprintf(ioQQQ,"%10.3g mean cx acceleration loops/iteration\n",((double)conv.getCounter(CONV_BASE_ACCELS))/
336 (MAX2((double)conv.getCounter(CONV_BASE_LOOPS),1.0)));
337 fprintf(ioQQQ,"%10.3g mean iso convergence loops/ion solve\n",((double)conv.getCounter(ISO_LOOPS))/
338 (MAX2((double)conv.getCounter(ION_SOLVES),1.0)));
339 fprintf(ioQQQ,"%10.3g mean steps/chemistry solve\n",((double)conv.getCounter(MOLE_SOLVE_STEPS))/
340 (MAX2((double)conv.getCounter(MOLE_SOLVE),1.0)));
341 fprintf(ioQQQ,"%10.3g mean step length searches/chemistry step\n",((double)conv.getCounter(NEWTON_LOOP))/
342 (MAX2((double)conv.getCounter(NEWTON),1.0)));
343 fprintf(ioQQQ,"----------------------------------------------------\n\n");
344 }
345
346 /* check whether any asserts were present and failed.
347 * return is true if ok, false if not. routine also checks
348 * number of warnings and returns false if not ok */
349 lgOK = lgCheckMonitors(ioQQQ);
350
351 if( lgOK && !warnings.lgWarngs && !lgAbort )
352 {
353 /* no failed asserts or warnings */
354 return false;
355 }
356 else
357 {
358 /* there were failed asserts or warnings */
359 return true;
360 }
361}
362
363/*BadStart announce that things are so bad the calculation cannot even start */
365{
366 char chLine[INPUT_LINE_LENGTH];
367
368 DEBUG_ENTRY( "BadStart()" );
369
370 /* initialize the line saver */
371 wcnint();
372 sprintf( warnings.chRgcln[0], " Calculation stopped because initial conditions out of bounds." );
373 sprintf( chLine, " W-Calculation could not begin." );
374 warnin(chLine);
375
376 if( grid.lgGrid )
377 {
378 /* possibly save out some results from this zone when doing grid
379 * since grid save files expect something at this grid point */
380 SaveDo("MIDL");
381 SaveDo("LAST" );
382 }
383 return;
384}
void AbundancesSet(void)
void AbundancesPrt(void)
void atmdat_readin(void)
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
STATIC void BadStart()
Definition cloudy.cpp:364
bool cloudy()
Definition cloudy.cpp:38
void SanityCheck(const char *chJob)
void ContCreateMesh(void)
void ContCreatePointers(void)
void ContSetIntensity()
t_continuum continuum
Definition continuum.cpp:5
t_conv conv
Definition conv.cpp:5
int ConvPresTempEdenIoniz(void)
void ConvIterCheck(void)
@ CONV_BASE_LOOPS
Definition conv.h:76
@ CONV_BASE_ACCELS
Definition conv.h:77
@ ISO_LOOPS
Definition conv.h:79
@ MOLE_SOLVE
Definition conv.h:71
@ NEWTON_LOOP
Definition conv.h:74
@ MOLE_SOLVE_STEPS
Definition conv.h:72
@ CONV_BASE_CALLS
Definition conv.h:75
@ ION_SOLVES
Definition conv.h:78
@ NEWTON
Definition conv.h:73
bool ConvInitSolution()
bool lgElemsConserved(void)
Definition dense.cpp:99
t_dynamics dynamics
Definition dynamics.cpp:44
bool lgConserveEnergy()
Definition energy.cpp:314
t_grid grid
Definition grid.cpp:5
void InitDefaultsPreparse(void)
void InitSimPostparse(void)
void InitCoreloadPostparse(void)
void Badnell_rec_init(void)
int iter_end_check(void)
void IterRestart(void)
void IterEnd(void)
void IterStart(void)
t_iterations iterations
Definition iterations.cpp:5
void lines(void)
Definition prt_lines.cpp:34
void LineStackCreate(void)
bool lgCheckMonitors(FILE *ioMONITOR)
t_noexec noexec
Definition noexec.cpp:5
void OpacityCreateAll(void)
void ParseCommands(void)
void CloseSaveFiles(bool lgFinal)
void plot(const char *chCall)
Definition plot.cpp:38
t_prt prt
Definition prt.cpp:10
void PrtFinal(void)
Definition prt_final.cpp:75
void PrtHeader(void)
void PrtComment(void)
int radius_next(void)
void radius_increment(void)
void radius_first(void)
void RT_continuum(void)
void RT_tau_inc(void)
void RT_diffuse(void)
void RT_tau_init(void)
void RT_tau_reset(void)
void SaveDo(const char *chTime)
Definition save_do.cpp:573
void state_get_put(const char chJob[])
Definition state.cpp:76
t_state state
Definition state.cpp:19
t_warnings warnings
Definition warnings.cpp:11
void wcnint(void)
Definition warnings.cpp:13
void warnin(char *chLine)
Definition warnings.cpp:27
void ZoneStart(const char *chMode)
void ZoneEnd(void)