cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
init_defaults_preparse.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 /*InitDefaultsPreparse initialization at start of simulation, called from cloudy
4 * before parser, sets initial values of quantities changed by parser */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "phycon.h"
8 #include "trace.h"
9 #include "noexec.h"
10 #include "stopcalc.h"
11 #include "plot.h"
12 #include "rt.h"
13 #include "fudgec.h"
14 #include "abund.h"
15 #include "extinc.h"
16 #include "h2.h"
17 #include "hextra.h"
18 #include "mole.h"
19 #include "wind.h"
20 #include "pressure.h"
21 #include "thermal.h"
22 #include "optimize.h"
23 #include "punch.h"
24 #include "dense.h"
25 #include "hcmap.h"
26 #include "prt.h"
27 #include "init.h"
28 #include "taulines.h"
29 #include "lines_service.h"
30 
31 /*InitDefaultsPreparse initialization at start of simulation, called from cloudy
32  * before parser, will now be called one time per sim in grid but long term
33  * goal is to call parser only one time in a grid. This would be called
34  * first to set defaults that may be changed by the parser. This routine
35  * should only set defaults for variables that change in the parser.
36  * It does not initialize variables for the start of a calculation. That
37  * is done by InitSimPostparse */
39 {
40  long int i,
41  ipISO,
42  nelem;
43 
44  DEBUG_ENTRY( "InitDefaultsPreparse()" );
45 
46  /* init vars before parsing commands for each sim */
47 
48  /* option to read in all input quantiles and NOT execute the actual model
49  * only check on input parameters - set by calling cdNoExec */
50  noexec.lgNoExec = false;
51 
52  /* constant for the extinguish command */
53  extinc.excolm = 0.;
54  extinc.exleak = 0.;
55  extinc.exlow = 1.;
56 
57  /* parameters having to do with thermal map */
58  hcmap.RangeMap[0] = 10.f;
59  hcmap.RangeMap[1] = .99e10f;
60  /* zone where map is to be done */
61  hcmap.MapZone = -1;
62  hcmap.nMapStep = 20;
63 
67  thermal.ConstTemp = 0.;
68  thermal.lgTeHigh = false;
69  thermal.lgTeBD96 = false;
70  thermal.lgTLaw = false;
71  thermal.lgTeSN99 = false;
72  /* try toe predict next zone's temperature in constant density models,
73  * as done in ZoneStart. Turned off with no tepred command */
74  thermal.lgPredNextTe = true;
75 
76  /* turbulent heating - set with hextra command */
77  hextra.TurbHeat = 0.;
78  /* set true with TIME option on hextra command for time dependent sims */
79  hextra.lgTurbHeatVaryTime = false;
80  /* options for for extra heating to depend on scale radius */
81  hextra.lgHextraDepth = false;
82  /* options for for extra heating to depend on density */
83  hextra.lgHextraDensity = false;
84 
85  /* options set by cosmic ray command */
86  hextra.cryden = 0.;
89  hextra.crtemp = 0.;
90  hextra.crpowr = 0.;
92 
93  /* this is number of electronic levels to include in the print and punch output
94  * changed with the PRINT LINE H2 ELECTRONIC and PUNCH H2 commands
95  * by default only include X */
96  h2.nElecLevelOutput = 1;
97 
98  /* holds options for trace set with atom h2 command */
99  mole.nH2_TRACE = false;
100 
101  /* option to scramble collision data */
102  mole.lgH2_NOISE = false;
103  mole.lgH2_NOISECOSMIC = false;
104 
105  /* option to turn off or on gbar collisions of the collision rate,
106  * default is to have it on */
107  /* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
108  mole.lgColl_gbar = true;
109 
110  /* include collision rates that come from real calculations,
111  * off with atom h2 collisions off command */
112  mole.lgColl_deexec_Calc = true;
113 
114  mole.lgColl_dissoc_coll = true;
115 
116  /* option to turn off H2 - grain collision & deexcitation,
117  * atom h2 grain collision on/off */
119 
120  /* option to turn off ortho-para collisions, command ATOM H2 COLLISIONS ORTHO PARA OFF */
122 
123  /* which set of H2 - He collision data to use? default is Le Bourlot Meudon set,
124  * set to ORNL with command
125  * atom H2 He collisions ORNL */
126  mole.lgH2_He_ORNL = true;
127 
128  /* use 1999 or 2007 H2 - H collision set atomic data mole H2 h */
129  h2.lgH2_H_coll_07 = true;
130 
131  /* parameters to do with wind */
132  wind.lgWindOK = true;
133  wind.DiskRadius = 0;
134  wind.lgDisk = false;
135  wind.windv0 = 0.;
136  wind.comass = 0.;
137  wind.windv = 0.;
138  wind.emdot = 0.;
139  wind.AccelAver = 0.;
140  wind.acldr = 0.;
141  wind.AccelGravity = 0.;
142  wind.AccelTot = 0.;
143  wind.AccelCont = 0.;
144  wind.AccelLine = 0.;
145  wind.AccelPres = 0.;
146  wind.AccelMax = 0.;
147  wind.fmul = 0.;
148  wind.lgVelPos = true;
149 
150  /* argument on ELEMENT LIMIT OFF XXX command, lowest abundance
151  * element to include in the calculation */
152  dense.AbundanceLimit = 0.;
153 
154  /* controls density fluctuations, when true variations are due to density changing,
155  * when false are due to abundances changing if changes are on */
156  dense.lgDenFlucOn = true;
157  dense.lgDenFlucRadius = true;
158  dense.flong = 0.;
159  dense.flcPhase = 0.;
160 
161  /* this says keep initial density constant,
162  * so pressure from iter to iter not really const */
163  dense.lgDenseInitConstant = true;
164 
165  /* individual terms for the pressure equation of state */
166  /* >>chng 05 jan 01, all three are set true at start *
167  * code default is constant density, so all this is ignored
168  * is turned to constant pressure then these flags must be adjusted so
169  * that we get the pressure we expect. with all true, all three
170  * contributors to pressure will be counted - with constant gas pressure
171  * these are turned off */
174  pressure.lgPres_ram_ON = true;
175  /* flag for constant pressure, include continuum too */
176  pressure.lgContRadPresOn = true;
177  /* constant density is the default */
178  strcpy( dense.chDenseLaw, "CDEN" );
179  /* number on line is log of nT - option to specify initial pressure */
181  /* this is log of nT product - if not present then set zero */
183 
184  /* select certain atomic data, changed with set atomic data command
185  * this says to use Zeippen 1982 [OII] transition probabilities */
186  dense.lgAsChoose[ipOXYGEN][1] = false;
187 
188  abund.lgAbnSolar = false;
189 
190  /* the ipSolar array is a set of pointers used for reading
191  * in abundances with the "abundances" command
192  * the hydrogen abundance of unity is not read in */
193  abund.npSolar = LIMELM - 1;
194  for( i=0; i < abund.npSolar; i++ )
195  {
196  abund.ipSolar[i] = i + 2;
197  }
198 
199  /* option to turn off an element */
200  for( nelem=0; nelem < LIMELM; nelem++ )
201  {
202  /* set of scale factors for changing abundances with elements command */
203  abund.ScaleElement[nelem] = 1.;
204  abund.solar[nelem] = abund.SolarSave[nelem];
205  }
206 
207  abund.lgAbTaON = false;
208 
209  /* option to turn off an element */
210  for( nelem=0; nelem < LIMELM; nelem++ )
211  {
212  /* option to have abundances from table */
213  abund.lgAbunTabl[nelem] = false;
214  }
215 
216  /* threshold for faintest heating cooling to punch with punch heating or
217  * punch cooling commands, reset with PUNCH WEAKHEATCOOL command */
218  punch.WeakHeatCool = 0.05f;
219 
220  /* set of variables used to control punch command */
221  punch.npunch = 0;
222  punch.cp_npun = 0;
223  punch.lgPunContinuum = false;
224  punch.lgDRPLst = false;
225  punch.lgDRHash = true;
227  /* this is the string that will appear after each model in the punch output,
228  * reset with the "set punch hash" command */
229  strcpy( punch.chHashString , " ###########################" );
230  /* this is the string that will be added as a prefix to filenames for punch output,
231  * reset with the "set punch prefix" command */
232  strcpy( punch.chFilenamePrefix , "" );
233  /* punch every one continuum point, set skip punch command */
234  punch.ncPunchSkip = 1;
235  /* flush file after every iteration - set punch flush */
236  punch.lgFLUSH = false;
237  for( i=0; i<LIMPUN; ++i )
238  {
239  punch.lgHashEndIter[i] = true;
240  /* set false for time dependent calculations*/
241  punch.lg_separate_iterations[i] = true;
242  punch.lgPunchEveryZone[i] = false;
243  punch.nPunchEveryZone[i] = -1;
244  }
245 
246  /* default line width for plots is 1e3 km/s, reset with
247  * set punch line width command */
248  punch.PunchLWidth = (realnum)(SPEEDLIGHT/1000.e5);
249 
250  /* default no printing of optical depths, TooFaint is .1 */
251  prt.lgPrtTau = false;
252  prt.PrtTauFnt = 0.1f;
253  prt.lgPrtShort = false;
254  prt.TooFaint = 1e-3f;
255  prt.lgFaintOn = true;
256  prt.lgFntSet = false;
257  prt.lgPrnLineCell = false;
258  prt.nPrnLineCell = -1000;
259 
260  /* if true then print main block of lines as array,
261  * set false with print lines column, will then
262  * do a single column of lines */
263  prt.lgPrtLineArray = true;
264 
265  /* when printing a column this is option to print linear rather than log */
266  prt.lgPrtLineLog = true;
267 
268  /* print emergent intensities even when sphere, set with print line emergenet */
269  prt.lgPrtLineEmergent = false;
270 
271  /* print ages */
272  prt.lgPrnAges = false;
273 
274  /* print column densities */
275  prt.lgPrintColumns = true;
276 
277  /* option to sort lines by wavelength, print sort command */
278  prt.lgSortLines = false;
279 
280  prt.lgPrtMaser = false;
281  prt.lgPrintTime = true;
282 
283  prt.lgPrtContIndices = false;
284  prt.lgPrtCont = false;
285  prt.lgPrnDiff = false;
286  prt.lgPrnPump = false;
287  prt.lgPrnInwd = false;
288  prt.lgPrnColl = false;
289  prt.lgPrnHeat = false;
290  /* >>chng 00 dec 08, these determine the standard items included in "nFnu", PvH */
291  prt.lgSourceReflected = true;
292  prt.lgSourceTransmitted = false;
293  prt.lgDiffuseInward = true;
294  prt.lgDiffuseOutward = true;
295  prt.lgPrtLastIt = false;
296  prt.lgOnlyZone = false;
297  prt.lgOnlyHead = false;
298  prt.lgPrtStart = false;
299  prt.nstart = 0;
300 
301  /* turn off printing of heating agents */
302  prt.lgPrintHeating = false;
303 
304  /* flag saying to print all matrix elements for ionization balance
305  * set with PRINT ARRAYS command */
306  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
307  {
308  prt.lgPrtArry[nelem] = false;
309  }
310 
311  /* print line flux at earth */
312  prt.lgPrintFluxEarth = false;
313 
314  /* print line surface brightness, def sr, option arcsec */
315  prt.lgSurfaceBrightness = false;
317 
318  prt.nzdump = -100;
319 
320  trace.lgSecIon = false;
321  trace.lgTrOvrd = true;
322  trace.lgOpacBug = false;
323  trace.nTrConvg = 0;
324  trace.lgTr8446 = false;
325  trace.lgTrLevN = false;
326  trace.lgTrGant = false;
327  trace.lgOptcBug = false;
328  trace.lgTrace3Bod = false;
329  trace.lgOTSBug = false;
330  trace.lgESOURCE = false;
331  trace.lgTr_CO_Mole = false;
332  trace.lgTr_H2_Mole = false;
333  trace.lgHeatBug = false;
334  trace.lgHeavyBug = false;
335  trace.lgBug2nu = false;
336  trace.lgDrBug = false;
337  trace.lgWind = false;
338  trace.lgPtrace = false;
339  trace.lgDrv_cdLine = false;
340  trace.lgDustBug = false;
341  trace.lgComBug = false;
342  trace.lgHeBug = false;
343  trace.lgCarBug = false;
344  trace.lgCalBug = false;
345  trace.lgConBug = false;
346  trace.lgNeBug = false;
347  trace.lgFeBug = false;
348  trace.lgHBug = false;
349  trace.lgTrLine = false;
350  trace.nznbug = 10000;
351  trace.npsbug = 10000;
352  trace.lgTrace = false;
353  trace.lgPointBug = false;
354  trace.lgNeonBug = false;
355  trace.lgCoolTr = false;
356  trace.lgTrDiff = false;
357  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
358  trace.lgIsoTraceFull[ipISO] = false;
359 
360  /* variables used in stop ... command */
361 
362  /* various criteria for stopping model */
363  /* >>chng 04 dec 21, remove from here and init to 1e30 in zero */
364  /*StopCalc.tauend = 0.;*/
365  StopCalc.tauend = 1e30f;
366 
367  /* >>chng 05 nov 22 - NPA. Stop calculation when fraction of oxygen frozen
368  * out on grains ices gets too high - formation of ices */
369  /*StopCalc.StopDepleteFrac = 0.99f;*/
370  /* >>chng 05 dec 16, with revised ion solver logic, code should be able to
371  * converge away from situation where ices have disturbed the chemistry and
372  * net negative atomic abundances result. now we say solution not converged and
373  * soldier on
374  * this test should not be necessary */
375  StopCalc.StopDepleteFrac = 1.02f;
376 
377  StopCalc.xMass = 0.;
378  StopCalc.taunu = 0.;
379  StopCalc.iptnu = -1;
380  /* stopping AV */
381  StopCalc.AV_extended = 1e30f;
382  StopCalc.AV_point = 1e30f;
383  /* highest allowed temperature */
385 
386  /* the floor sets a limit to the temperature in the calculation -
387  * if te falls below this, we do a constant temperature cloud at
388  * this temperature */
389  StopCalc.TeFloor = 0.;
390 
391  /* this is in cddefines.h and is 4000 */
393  /* ending column densities */
401 
402  /* stopping electron density */
404 
405  /* stopping electron and molecular fractions */
406  StopCalc.StopElecFrac = -FLT_MAX;
407  StopCalc.StopHPlusFrac = -FLT_MAX;
408  /* stopping molecular fraction has opposite sign - want to stop when 2H_2/NH gt this */
409  StopCalc.StopH2MoleFrac = FLT_MAX;
410  /* this flag says that 21cm line optical depth is the stop quantity */
411  StopCalc.lgStop21cm = false;
412  /* debugging aid - stop when conv.nTotalIoniz reaches this value */
414  /* stop when absolute value of velocity falls below this */
415  StopCalc.StopVelocity = 0.;
416  /* number of stop line commands entered */
417  StopCalc.nstpl = 0;
418 
419  /* initialize some variables for the optimizer */
420  optimize.nIterOptim = 20;
421  optimize.ioOptim = NULL;
422  optimize.OptGlobalErr = 0.10f;
423  optimize.nlobs = 0;
424  optimize.nTempObs = 0;
425  optimize.ncobs = 0;
426  optimize.nRangeSet = 0;
427  strcpy( optimize.chOptRtn, "SUBP" );
428 
429  /* flags says what is to be matched */
430  optimize.lgOptLin = false;
431  optimize.lgOptLum = false;
432  optimize.lgOptCol = false;
433  optimize.lgOptTemp = false;
434  optimize.lgOptimize = false;
435 
436  /* trace flag for optimization process */
437  optimize.lgTrOpt = false;
438 
439  optimize.lgOptimFlow = false;
440  optimize.optint = 0.;
441  optimize.optier = 0.;
442 # ifdef __unix
443  optimize.lgParallel = true;
444 # else
445  optimize.lgParallel = false;
446 # endif
447  optimize.lgOptCont = false;
448 
449  /* the fudge factors command */
450  fudgec.nfudge = 0;
451  fudgec.lgFudgeUsed = false;
452  for( i=0; i < NFUDGC; i++ )
453  fudgec.fudgea[i] = 0.;
454 
455  EmLineZero( &DummyEmis );
456  DummyEmis.iRedisFun = 0;
457  DummyEmis.ipFine = -1;
458  DummyEmis.gf = 0.;
459  DummyEmis.damp = 0.;
460  DummyEmis.dampXvel = 0.;
461  DummyEmis.opacity = 0.;
462  DummyEmis.Aul = 1e-30f;
463  DummyEmis.tran = NULL;
464  DummyEmis.next = NULL;
465 
466  /* parameters dealing with printer plots */
467  for( i=0; i < NDPLOT; i++ )
468  {
469  plotCom.lgPltTrace[i] = false;
470  }
471 
472  /* this says what types of printer plots we will make */
473  for( i=0; i < NDPLOT; i++ )
474  {
475  strcpy( plotCom.chPType[i], "NONE" );
476  }
477  plotCom.lgPlotON = false;
478 
479  /* following were block data logic */
480  rt.lgStarkON = true;
481 
482  /* by default use Federman form of shielding function */
484 
485  return;
486 }

Generated for cloudy by doxygen 1.8.1.1