cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
dynamics.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /* DynaIterEnd called at end of iteration when advection is turned on */
4 /* DynaStartZone called at start of zone calculation when advection is turned on */
5 /* DynaEndZone called at end of zone calculation when advection is turned on */
6 /* DynaIonize, called from ionize to evaluate advective terms for current conditions */
7 /* DynaPresChngFactor, called from PressureChange to evaluate new density needed for
8  * current conditions and wind solution, returns ratio of new to old density */
9 /* timestep_next - find next time step in time dependent case */
10 /* DynaZero zero some dynamics variables, called from zero.c */
11 /* DynaCreateArrays allocate some space needed to save the dynamics structure variables,
12  * called from DynaCreateArrays */
13 /* DynaPrtZone - called to print zone results */
14 /* DynaSave save info related to advection */
15 /* DynaSave, save output for dynamics solutions */
16 /* ParseDynaTime parse the time command, called from ParseCommands */
17 /* ParseDynaWind parse the wind command, called from ParseCommands */
18 #include "cddefines.h"
19 #include "cddrive.h"
20 #include "struc.h"
21 #include "colden.h"
22 #include "radius.h"
23 #include "stopcalc.h"
24 #include "hextra.h"
25 #include "iterations.h"
26 #include "conv.h"
27 #include "timesc.h"
28 #include "dense.h"
29 #include "mole.h"
30 #include "thermal.h"
31 #include "pressure.h"
32 #include "phycon.h"
33 #include "wind.h"
34 #include "iso.h"
35 #include "dynamics.h"
36 #include "cosmology.h"
37 #include "parser.h"
38 #include "rfield.h"
39 
40 static const bool lgPrintDynamics = false;
41 
43 static int ipUpstream=-1,iphUpstream=-1,ipyUpstream=-1;
44 
45 /*
46  * >>chng 01 mar 16, incorporate advection within dynamical solutions
47  * this file contains the routines that incorporate effects of dynamics and advection upon the
48  * thermal and ionization solutions.
49  *
50  * This code was originally developed in March 2001 during
51  * a visit to UNAM Morelia, in collaboration with Robin Williams, Jane Arthur, and Will Henney.
52  * Development was continued by email, and in a meeting in July/Aug 2001 at U Cardiff
53  * Further development June 2002, UNAM Morelia (Cloudy and the Duendes Verdes)
54  */
55 
56 /* the interpolated upstream densities of all ionization stages,
57  * [element][ion] */
58 static double **UpstreamIon;
59 static double ***UpstreamStatesElem;
60 /* total abundance of each element per hydrogen */
61 static vector<double> UpstreamElem;
62 
63 /* hydrogen molecules */
64 static vector<double> Upstream_molecules;
65 
66 /* space to save continuum when time command is used
67 static realnum *dyna_flux_save;*/
68 
69 /* array of times and continuum values we will interpolate upon */
70 static vector<double> time_elapsed_time ,
72  time_dt,
74 static bool lgtime_dt_specified;
75 static vector<int> lgtime_Recom;
76 #define NTIME 200
77 
78 /* number of time steps actually read in */
79 static long int nTime_flux=0;
80 
81 /* routine called at end of iteration to determine new step sizes */
82 STATIC void DynaNewStep(void);
83 
84 /* routine called at end of iteration to save values in previous iteration */
85 STATIC void DynaSaveLast(void);
86 
87 /* routine called to determine mass flux at given distance */
88 /* static realnum DynaFlux(double depth); */
89 
90 /* lookahead distance, as derived in DynaIterEnd */
91 static double Dyn_dr;
92 
93 /* advected work */
94 static double AdvecSpecificEnthalpy;
95 
96 /* depth of position in previous iteration */
97 static vector<double> Old_depth/*[NZLIM]*/;
98 
99 /* HI ionization structure from previous iteration */
100 static vector<realnum> Old_histr/*[NZLIM]*/ ,
101  /* Lyman continuum optical depth from previous iteration */
102  Old_xLyman_depth/*[NZLIM]*/,
103  /* old n_p density from previous iteration */
104  Old_hiistr/*[NZLIM]*/,
105  /* old pressure from previous iteration */
106  Old_pressure/*[NZLIM]*/,
107  /* H density - particles per unit vol */
108  Old_density/*[NZLIM]*/ ,
109  /* density - total grams per unit vol */
110  Old_DenMass/*[NZLIM]*/ ,
111  /* sum of enthalpy and kinetic energy per gram */
112  EnthalpyDensity/*[NZLIM]*/,
113  /* old electron density from previous iteration */
114  Old_ednstr/*[NZLIM]*/,
115  /* sum of enthalpy and kinetic energy per gram */
116  Old_EnthalpyDensity/*[NZLIM]*/;
117 
119 
120 /* the ionization fractions from the previous iteration */
122 
123 /* the gas phase abundances from the previous iteration */
125 
126 /* the iso levels from the previous iteration */
128 
129 /* the number of zones that were saved in the previous iteration */
130 static long int nOld_zone;
131 
132 STATIC bool lgNeedTimestep ();
133 STATIC void InitDynaTimestep();
134 
135 
136 /*timestep_next - find next time step in time dependent case */
137 STATIC double timestep_next( void )
138 {
139  DEBUG_ENTRY( "timestep_next()" );
140 
141  double timestep_return = dynamics.timestep;
142 
143  if( dynamics.lgRecom )
144  {
145  timestep_return = 0.04 * timesc.time_therm_short;
146  }
147  else
148  {
149  timestep_return = dynamics.timestep_init;
150  }
151  if( lgPrintDynamics )
152  fprintf( ioQQQ, "DEBUG timestep_next returns %.3e\n", timestep_return );
153 
154  return( timestep_return );
155 }
156 
157 /* ============================================================================== */
158 /* DynaIonize, called from ConvBase to evaluate advective terms for current conditions,
159  * calculates terms to add to ionization balance equation */
160 void DynaIonize(void)
161 {
162  DEBUG_ENTRY( "DynaIonize()" );
163 
164  /* the time (s) needed for gas to move dynamics.Dyn_dr */
165  /* >>chng 02 dec 11 rjrw increase dynamics.dynamics.timestep when beyond end of previous zone, to allow -> eqm */
166  if( !dynamics.lgTimeDependentStatic )
167  {
168  /* in time dependent model dynamics.timestep only changes once at end of iteration
169  * and is constant across a model */
170  /* usual case - full dynamics */
171  dynamics.timestep = -Dyn_dr/wind.windv;
172  }
173  /* printf("%d %g %g %g %g\n",nzone,radius.depth,Dyn_dr,radius.depth-Old_depth[nOld_zone-1],dynamics.timestep); */
174 
176  if( nzone > 0 )
178 
179  /* do nothing on first iteration or when looking beyond previous iteration */
180  /* >>chng 05 jan 27, from hardwired "iteration < 2" to more general case,
181  * this is set with SET DYNAMICS RELAX command with the default of 2 */
182  //For advection cases, switch to local equilibrium when depth is outside
183  //the region of previous iteration.
184  //Possibly should limit range of dynamical sources further by adding Dyn_dr???
185  double depth = radius.depth;
186  if( iteration < dynamics.n_initial_relax+1 ||
187  ( ! dynamics.lgTimeDependentStatic &&
188  ( depth < 0 || depth > dynamics.oldFullDepth ) ) )
189  {
190  /* first iteration, return zero */
191  dynamics.Cool_r = 0.;
192  dynamics.Heat_v = 0.;
193  dynamics.dHeatdT = 0.;
194 
195  dynamics.Rate = 0.;
196 
197  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
198  {
199  for( long ion=0; ion<nelem+2; ++ion )
200  {
201  dynamics.Source[nelem][ion] = 0.;
202  }
203  }
204  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
205  {
206  for( long nelem=ipISO; nelem<LIMELM; ++nelem)
207  {
208  if( dense.lgElmtOn[nelem] )
209  {
210  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; ++level )
211  {
212  dynamics.StatesElem[nelem][nelem-ipISO][level] = 0.;
213  }
214  }
215  }
216  }
217  for( long mol=0;mol<mole_global.num_calc;mol++)
218  {
219  dynamics.molecules[mol] = 0.;
220  }
221  return;
222  }
223 
224  if( dynamics.lgTracePrint )
225  {
226  fprintf(ioQQQ,"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
227  nzone,
230  5./2.*pressure.PresGasCurr
231  );
232  }
233 
234  /* net cooling due to advection */
235  /* >>chng 01 aug 01, removed hden from dilution, put back into here */
236  /* >>chng 01 sep 15, do not update this variable the very first time this routine
237  * is called at the new zone. */
238  dynamics.Cool_r = 1./dynamics.timestep*dynamics.lgCoolHeat;
239  dynamics.Heat_v = AdvecSpecificEnthalpy/dynamics.timestep*dynamics.lgCoolHeat;
240  dynamics.dHeatdT = 0.*dynamics.lgCoolHeat;
241 
242 # if 0
243  if( dynamics.lgTracePrint || nzone>17 && iteration == 10)
244  {
245  fprintf(ioQQQ,
246  "dynamics cool-heat\t%li\t%.3e\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
247  nzone,
248  phycon.te,
249  dynamics.lgCoolHeat,
250  thermal.htot,
255  scalingDensity(),
256  dynamics.timestep);
257  }
258 # endif
259 
260  /* second or greater iteration, have advective terms */
261  /* this will evaluate advective terms for current physical conditions */
262 
263  /* the rate (s^-1) atoms drift in from upstream, a source term for the ground */
264 
265  /* dynamics.Hatom/dynamics.timestep is the source (cm^-3 s^-1) of neutrals,
266  normalized to (s^-1) by the next higher ionization state as for all
267  other recombination terms */
268 
269  /* dynamics.xIonDense[ipHYDROGEN][1]/dynamics.timestep is the sink (cm^-3 s^-1) of
270  ions, normalized to (s^-1) by the ionization state as for all other
271  ionization terms */
272 
273  dynamics.Rate = 1./dynamics.timestep;
274 
275  for( long mol=0;mol<mole_global.num_calc;mol++)
276  {
277  dynamics.molecules[mol] = Upstream_molecules[mol]*scalingDensity() / dynamics.timestep;
278  }
279 
280  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
281  {
282  if( dense.lgElmtOn[nelem] )
283  {
284  /* check that the total number of each element is conserved along the flow */
285  if(fabs(UpstreamElem[nelem]*scalingDensity() -dense.gas_phase[nelem])/dense.gas_phase[nelem]>=1e-3)
286  {
287  fprintf(ioQQQ,
288  "PROBLEM conservation error: zn %li elem %li upstream %.8e abund %.8e (up-ab)/up %.2e\n",
289  nzone ,
290  nelem,
291  UpstreamElem[nelem]*scalingDensity(),
292  dense.gas_phase[nelem] ,
293  (UpstreamElem[nelem]*scalingDensity()-dense.gas_phase[nelem]) /
294  (UpstreamElem[nelem]*scalingDensity()) );
295  }
296  /* ASSERT( fabs(UpstreamElem[nelem]*scalingDensity() -dense.gas_phase[nelem])/dense.gas_phase[nelem]<1e-3); */
297  for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
298  {
299  dynamics.Source[nelem][ion] = 0.;
300  }
301  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
302  {
303  /* Normalize to next higher state in current zone, except at first iteration
304  where upstream version may be a better estimate (required for
305  convergence in the small dynamics.timestep limit) */
306 
307  dynamics.Source[nelem][ion] =
308  /* UpstreamIon is ion number per unit hydrogen because dilution is 1/hden
309  * this forms the ratio of upstream atom over current ion, per dynamics.timestep,
310  * so Source has units cm-3 s-1 */
311  UpstreamIon[nelem][ion] * scalingDensity() / dynamics.timestep;
312 
313  }
314  for( long ion=dense.IonHigh[nelem]+1;ion<nelem+2; ++ion )
315  {
316  dynamics.Source[nelem][ion] = 0.;
317  dynamics.Source[nelem][dense.IonHigh[nelem]] +=
318  UpstreamIon[nelem][ion] * scalingDensity() / dynamics.timestep;
319  }
320  }
321  }
322 
323  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
324  {
325  for( long nelem=ipISO; nelem<LIMELM; ++nelem)
326  {
327  if( dense.lgElmtOn[nelem] )
328  {
329  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; ++level )
330  {
331  dynamics.StatesElem[nelem][nelem-ipISO][level] =
332  UpstreamStatesElem[nelem][nelem-ipISO][level] * scalingDensity() / dynamics.timestep;
333  }
334  }
335  }
336  }
337 
338 # if 0
339  fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
340  nzone,
341  dynamics.Rate,
342  dynamics.Source[ipHYDROGEN][0],
343  dynamics.Rate,
344  dynamics.Source[ipCARBON][3]);
345 # endif
346 #if 0
347  long nelem = ipCARBON;
348  long ion = 3;
349  /*if( nzone > 160 && iteration > 1 )*/
350  fprintf(ioQQQ,"dynaionizeeee\t%li\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
351  nzone,
352  ipUpstream,
353  radius.depth ,
355  dense.xIonDense[nelem][ion],
356  UpstreamIon[nelem][ion]* scalingDensity(),
357  Old_xIonDense[ipUpstream][nelem][ion] ,
358  dense.xIonDense[nelem][ion+1],
359  UpstreamIon[nelem][ion+1]* scalingDensity(),
360  Old_xIonDense[ipUpstream][nelem][ion+1] ,
361  dynamics.timestep,
362  dynamics.Source[nelem][ion]
363  );
364 #endif
365  if( dynamics.lgTracePrint )
366  {
367  fprintf(ioQQQ," DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
368  nzone,dynamics.Rate , dynamics.Source[0][0] );
369  }
370  return;
371 }
372 
373 /* ============================================================================== */
374 /* DynaStartZone called at start of zone calculation when advection is turned on */
375 void DynaStartZone(void)
376 {
377  /* this routine is called at the start of a zone calculation, by ZoneStart:
378  *
379  * it sets deduced variables to zero if first iteration,
380  *
381  * if upstream depth is is outside the computed structure on previous iteration,
382  * return value at shielded face
383  *
384  * Also calculates discretization_error, an estimate of the accuracy of the source terms.
385  *
386  * */
387 
388  /* this is index of upstream cell in structure stored from previous iteration */
389  double upstream, dilution, dilutionleft, dilutionright, frac_next;
390 
391  /* Properties for cell half as far upstream, used to converge dynamics.timestep */
392  double hupstream, hnextfrac=-BIGFLOAT, hion, hmol, hiso;
393 
394  /* Properties for cell at present depth, used to converge dynamics.timestep */
395  double ynextfrac=-BIGFLOAT, yion, ymol, yiso;
396 
397  long int nelem , ion, mol;
398 
399  DEBUG_ENTRY( "DynaStartZone()" );
400 
401  /* do nothing on first iteration */
402  if( iteration < 2 )
403  {
404  dynamics.Upstream_density = 0.;
406  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
407  {
408  for( ion=0; ion<nelem+2; ++ion )
409  {
410  UpstreamIon[nelem][ion] = 0.;
411  }
412  }
413  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
414  {
415  for( nelem=ipISO; nelem<LIMELM; ++nelem)
416  {
417  if( dense.lgElmtOn[nelem] )
418  {
419  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
420  {
421  UpstreamStatesElem[nelem][nelem-ipISO][level] = 0.;
422  }
423  }
424  }
425  }
426  /* hydrogen molecules */
427  for(mol=0; mol<mole_global.num_calc;mol++)
428  {
429  Upstream_molecules[mol] = 0;
430  }
431 
432  ipUpstream = -1;
433  iphUpstream = -1;
434  ipyUpstream = -1;
435  return;
436  }
437 
438  /* radius.depth is distance from illuminated face of cloud to outer edge of
439  * current zone, which has thickness radius.drad */
440 
441  /* find where the down stream point is, in previous iteration: we
442  * are looking for point where gas in current cell was in previous
443  * iteration */
444 
445  /* true, both advection and wind solution */
446  /* don't interpolate to the illuminated side of the first cell */
447  upstream = MAX2(Old_depth[0] , radius.depth + Dyn_dr);
448  hupstream = 0.5*(radius.depth + upstream);
449 
450  /* now locate upstream point in previous stored structure,
451  * will be at the same point or higher than we found previously */
452  while( (Old_depth[ipUpstream+1] < upstream ) &&
453  ( ipUpstream < nOld_zone-1 ) )
454  {
455  ++ipUpstream;
456  }
457  ASSERT( ipUpstream <= nOld_zone-1 );
458 
459  /* above loop will give ipUpstream == nOld_zone-1 if computed structure has been overrun */
460 
462  {
463  /* we have not overrun radius scale of previous iteration */
464  ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] );
465 
466  frac_next = ( upstream - Old_depth[ipUpstream])/
467  (Old_depth[ipUpstream+1] - Old_depth[ipUpstream]);
468 
469  if (! (frac_next >= 0.0 && frac_next <= 1.0))
470  {
471  fprintf(ioQQQ,"PROBLEM frac_next out of range %.16e %.16e %.16e %.16e\n",
472  Old_depth[ipUpstream],upstream,Old_depth[ipUpstream+1],frac_next);
473  }
474 
475  ASSERT (frac_next >= 0.0 && frac_next <= 1.0);
476  dynamics.Upstream_density = (realnum)(Old_density[ipUpstream] +
477  (Old_density[ipUpstream+1] - Old_density[ipUpstream])*
478  frac_next);
479  dilutionleft = 1./Old_density[ipUpstream];
480  dilutionright = 1./Old_density[ipUpstream+1];
481 
482  /* fractional changes in density from passive advection */
483  /* >>chng 01 May 02 rjrw: use hden for dilution */
484  /* >>chng 01 aug 01, remove hden here, put back into vars when used in DynaIonize */
485  dilution = 1./dynamics.Upstream_density;
486 
487  // we have a mix of realnum and double here, so make sure double is used for the test...
488  double lo = double(Old_EnthalpyDensity[ipUpstream])*dilutionleft,
489  hi = double(Old_EnthalpyDensity[ipUpstream+1])*dilutionright;
490  /* the advected enthalphy */
491  AdvecSpecificEnthalpy = lo + (hi-lo)*frac_next;
492  double x = AdvecSpecificEnthalpy;
493 
494  if( ! fp_bound(lo,x,hi) )
495  {
496  fprintf(ioQQQ,"PROBLEM interpolated enthalpy density is not within range %.16e\t%.16e\t%.16e\t%e\t%e\n",
497  lo, x, hi, (hi-x)/(hi-lo), (x-lo)/(hi-lo));
498  }
499 
500  ASSERT( fp_bound(lo,x,hi) );
501 
502  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
503  {
504  UpstreamElem[nelem] = 0.;
505  for( ion=0; ion<nelem+2; ++ion )
506  {
507  /* Easier to bring out the division by the old hydrogen
508  * density rather than putting in dilution and then
509  * looking for how dilution is defined. The code is
510  * essentially equivalent, but slower */
511  /* UpstreamIon is ion number per unit material (measured
512  * as defined in scalingdensity()), at the upstream
513  * position */
514  UpstreamIon[nelem][ion] =
516  (Old_xIonDense[ipUpstream+1][nelem][ion]/Old_density[ipUpstream+1] -
518  frac_next;
519 
520  UpstreamElem[nelem] += UpstreamIon[nelem][ion];
521  }
522  }
523 
524  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
525  {
526  for( nelem=ipISO; nelem<LIMELM; ++nelem)
527  {
528  if( dense.lgElmtOn[nelem] )
529  {
530  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
531  {
532  UpstreamStatesElem[nelem][nelem-ipISO][level] =
533  Old_StatesElem[ipUpstream][nelem][nelem-ipISO][level]/Old_density[ipUpstream] +
534  (Old_StatesElem[ipUpstream+1][nelem][nelem-ipISO][level]/Old_density[ipUpstream+1] -
535  Old_StatesElem[ipUpstream][nelem][nelem-ipISO][level]/Old_density[ipUpstream])*
536  frac_next;
537  /* check for NaN */
538  ASSERT( !isnan( UpstreamStatesElem[nelem][nelem-ipISO][level] ) );
539  }
540  }
541  }
542  }
543 
544  for(mol=0;mol<mole_global.num_calc;mol++)
545  {
547  (Old_molecules[ipUpstream+1][mol]/Old_density[ipUpstream+1] -
549  frac_next;
550  /* Externally represented species are already counted in UpstreamElem */
551  if(mole.species[mol].location == NULL && mole_global.list[mol]->isIsotopicTotalSpecies())
552  {
553  for(molecule::nNucsMap::iterator atom= mole_global.list[mol]->nNuclide.begin();
554  atom != mole_global.list[mol]->nNuclide.end(); ++atom)
555  {
556  UpstreamElem[atom->first->el->Z-1] +=
557  Upstream_molecules[mol]*atom->second;
558  }
559  }
560  }
561  }
562  else
563  {
564  /* SPECIAL CASE - we have overrun the previous iteration's radius */
565  long ipBound = ipUpstream;
566  if (ipBound == -1)
567  ipBound = 0;
568  dynamics.Upstream_density = Old_density[ipBound];
569  /* fractional changes in density from passive advection */
570  dilution = 1./dynamics.Upstream_density;
571  /* AdvecSpecificEnthalpy enters as heat term */
573  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
574  {
575  UpstreamElem[nelem] = 0.;
576  for( ion=0; ion<nelem+2; ++ion )
577  {
578  /* UpstreamIon is ion number per unit hydrogen */
579  UpstreamIon[nelem][ion] =
580  Old_xIonDense[ipBound][nelem][ion]/Old_density[ipBound];
581  UpstreamElem[nelem] += UpstreamIon[nelem][ion];
582  }
583  }
584 
585  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
586  {
587  for( nelem=ipISO; nelem<LIMELM; ++nelem)
588  {
589  if( dense.lgElmtOn[nelem] )
590  {
591  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
592  {
593  UpstreamStatesElem[nelem][nelem-ipISO][level] =
594  Old_StatesElem[ipBound][nelem][nelem-ipISO][level]/Old_density[ipBound];
595  /* check for NaN */
596  ASSERT( !isnan( UpstreamStatesElem[nelem][nelem-ipISO][level] ) );
597  }
598  }
599  }
600  }
601 
602  for(mol=0;mol<mole_global.num_calc;mol++)
603  {
604  Upstream_molecules[mol] = Old_molecules[ipBound][mol]/Old_density[ipBound];
605  if(mole.species[mol].location == NULL && mole_global.list[mol]->isIsotopicTotalSpecies())
606  {
607  for(molecule::nNucsMap::iterator atom=mole_global.list[mol]->nNuclide.begin();
608  atom != mole_global.list[mol]->nNuclide.end(); ++atom)
609  {
610  UpstreamElem[atom->first->el->Z-1] +=
611  Upstream_molecules[mol]*atom->second;
612  }
613  }
614  }
615  }
616 
617  /* Repeat enough of the above for half-step and no-step to judge convergence:
618  * the result of this code is the increment of discretization_error */
619 
620  while( (Old_depth[iphUpstream+1] < hupstream ) &&
621  ( iphUpstream < nOld_zone-1 ) )
622  {
623  ++iphUpstream;
624  }
625  ASSERT( iphUpstream <= nOld_zone-1 );
626 
627  while( (Old_depth[ipyUpstream+1] < radius.depth ) &&
628  ( ipyUpstream < nOld_zone-1 ) )
629  {
630  ++ipyUpstream;
631  }
632  ASSERT( ipyUpstream <= nOld_zone-1 );
633 
634  dynamics.dRad = BIGFLOAT;
635 
637  hnextfrac = ( hupstream - Old_depth[iphUpstream])/
639  else
640  hnextfrac = 0.;
641 
643  ynextfrac = ( radius.depth - Old_depth[ipyUpstream])/
644  (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]);
645  else
646  ynextfrac = 0.;
647 
648  // Shouldn't be jumping over large numbers of upstream cells
649  if(ipUpstream != -1 && ipUpstream < nOld_zone-1)
650  dynamics.dRad = MIN2(dynamics.dRad,
651  5*(Old_depth[ipUpstream+1] - Old_depth[ipUpstream]));
652 
653 
654  // Value for scaling zonal changes to set zone width
655  const double STEP_FACTOR=0.05;
656 
657  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
658  {
659  for( ion=0; ion<nelem+2; ++ion )
660  {
661  double f1;
662  if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
663  {
664  yion =
666  (Old_xIonDense[ipyUpstream+1][nelem][ion]/Old_density[ipyUpstream+1] -
668  ynextfrac;
669  }
670  else
671  {
672  long ipBound = ipyUpstream;
673  if (-1 == ipBound)
674  ipBound = 0;
675  yion = Old_xIonDense[ipBound][nelem][ion]/Old_density[ipBound];
676  }
677  if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
678  {
679  hion =
681  (Old_xIonDense[iphUpstream+1][nelem][ion]/Old_density[iphUpstream+1] -
683  hnextfrac;
684  }
685  else
686  {
687  long ipBound = iphUpstream;
688  if (-1 == ipBound)
689  ipBound = 0;
690  hion = Old_xIonDense[ipBound][nelem][ion]/Old_density[ipBound];
691  }
692 
693  /* the proposed thickness of the next zone */
694  f1 = fabs(yion - UpstreamIon[nelem][ion] );
695  if( f1 > SMALLFLOAT )
696  {
697  dynamics.dRad = MIN2(dynamics.dRad,STEP_FACTOR*fabs(Dyn_dr) *
698  /* don't pay attention to species with abundance relative to H less than 1e-6 */
699  MAX2(yion + UpstreamIon[nelem][ion],1e-6 ) / f1);
700  }
701 
702  /* Must be consistent with convergence_error below */
703  /* this error is error due to the advection length not being zero - a finite
704  * advection length. no need to bring convergence error to below
705  * discretization error. when convergece error is lower than a fraction of this
706  * errror we reduce the advection length. */
707  dynamics.discretization_error += POW2(yion-2.*hion+UpstreamIon[nelem][ion]);
708  dynamics.error_scale2 += POW2(UpstreamIon[nelem][ion]-yion);
709  }
710  }
711 
712  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
713  {
714  for( nelem=ipISO; nelem<LIMELM; ++nelem)
715  {
716  if( dense.lgElmtOn[nelem] )
717  {
718  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
719  {
720  double f1;
721  if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 &&
722  (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
723  {
724  yiso =
725  Old_StatesElem[ipyUpstream][nelem][nelem-ipISO][level]/Old_density[ipyUpstream] +
726  (Old_StatesElem[ipyUpstream+1][nelem][nelem-ipISO][level]/Old_density[ipyUpstream+1] -
727  Old_StatesElem[ipyUpstream][nelem][nelem-ipISO][level]/Old_density[ipyUpstream])*
728  ynextfrac;
729  }
730  else
731  {
732  long ipBound = ipyUpstream;
733  if (-1 == ipBound)
734  ipBound = 0;
735  yiso = Old_StatesElem[ipBound][nelem][nelem-ipISO][level]/Old_density[ipBound];
736  }
737  if(iphUpstream != -1 && iphUpstream != nOld_zone-1 &&
738  (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
739  {
740  hiso =
741  Old_StatesElem[iphUpstream][nelem][nelem-ipISO][level]/Old_density[iphUpstream] +
742  (Old_StatesElem[iphUpstream+1][nelem][nelem-ipISO][level]/Old_density[iphUpstream+1] -
743  Old_StatesElem[iphUpstream][nelem][nelem-ipISO][level]/Old_density[iphUpstream])*
744  hnextfrac;
745  }
746  else
747  {
748  long ipBound = iphUpstream;
749  if (-1 == ipBound)
750  ipBound = 0;
751  hiso = Old_StatesElem[ipBound][nelem][nelem-ipISO][level]/Old_density[ipBound];
752  }
753 
754  /* the proposed thickness of the next zone */
755  f1 = fabs(yiso - UpstreamStatesElem[nelem][nelem-ipISO][level] );
756  if( f1 > SMALLFLOAT )
757  {
758  dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
759  /* don't pay attention to species with abundance relative to H less tghan 1e-6 */
760  MAX2(yiso + UpstreamStatesElem[nelem][nelem-ipISO][level],1e-6 ) / f1);
761  }
762  /* Must be consistent with convergence_error below */
763  /* this error is error due to the advection length not being zero - a finite
764  * advection length. no need to bring convergence error to below
765  * discretization error. when convergece error is lower than a fraction of this
766  * error we reduce the advection length. */
767  dynamics.discretization_error += POW2(yiso-2.*hiso+UpstreamStatesElem[nelem][nelem-ipISO][level]);
768  dynamics.error_scale2 += POW2(UpstreamStatesElem[nelem][nelem-ipISO][level]);
769  }
770  }
771  }
772  }
773 
774  for(mol=0; mol < mole_global.num_calc; mol++)
775  {
776  double f1;
777  if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
778  {
779  ymol =
781  (Old_molecules[ipyUpstream+1][mol]/Old_density[ipyUpstream+1] -
783  ynextfrac;
784  }
785  else
786  {
787  long ipBound = ipyUpstream;
788  if (-1 == ipBound)
789  ipBound = 0;
790  ymol = Old_molecules[ipBound][mol]/Old_density[ipBound];
791  }
792  if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
793  {
794  hmol =
796  (Old_molecules[iphUpstream+1][mol]/Old_density[iphUpstream+1] -
798  hnextfrac;
799  }
800  else
801  {
802  long ipBound = iphUpstream;
803  if (-1 == ipBound)
804  ipBound = 0;
805  hmol = Old_molecules[ipBound][mol]/Old_density[ipBound];
806  }
807 
808  /* the proposed thickness of the next zone */
809  f1 = fabs(ymol - Upstream_molecules[mol] );
810  if( f1 > SMALLFLOAT )
811  {
812  dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
813  /* don't pay attention to species with abundance relative to H less than 1e-6 */
814  MAX2(ymol + Upstream_molecules[mol],1e-6 ) / f1 );
815  }
816 
817  /* Must be consistent with convergence_error below */
818  /* >>chngf 01 aug 01, remove scalingDensity() from HAtom */
819  /* >>chngf 02 aug 01, multiply by cell width */
820  dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_molecules[mol]);
821  dynamics.error_scale2 += POW2(Upstream_molecules[mol]-ymol);
822  }
823 
824  if( dynamics.lgTracePrint )
825  {
826  fprintf(ioQQQ," DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
827  nzone,dynamics.Rate , dynamics.Source[0][0] , dilution*scalingDensity() );
828  }
829  return;
830 }
831 
832 /* DynaEndZone called at end of zone calculation when advection is turned on */
833 void DynaEndZone(void)
834 {
835  DEBUG_ENTRY( "DynaEndZone()" );
836 
837  /* this routine is called at the end of a zone calculation, by ZoneEnd */
838 
840 
841  if(dynamics.lgTracePrint)
842  fprintf(ioQQQ,"Check dp: %g %g mom %g %g mas %g\n",
847  DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth));
848  return;
849 }
850 
851 
852 /* ============================================================================== */
853 /* DynaIterEnd called at end of iteration when advection is turned on */
854 void DynaIterEnd(void)
855 {
856  /* this routine is called by IterRestart at the end of an iteration
857  * when advection is turned on. currently it only derives a
858  * dynamics.timestep by looking at the spatial derivative of
859  * some stored quantities */
860  long int i;
861  static long int nTime_dt_array_element = 0;
862 
863  DEBUG_ENTRY( "DynaIterEnd()" );
864 
865  /* set stopping outer radius to current outer radius once we have let
866  * solution relax dynamics.n_initial_relax times
867  * note the off by one confusion
868  * want to do two static iterations then start dynamics
869  * iteration was incremented before this call so iteration == 2 at
870  * end of first iteration */
871  if( iteration == dynamics.n_initial_relax+1)
872  {
873  fprintf(ioQQQ,"DYNAMICS DynaIterEnd sets stop radius to %.2e after "
874  "dynamics.n_initial_relax=%li iterations.\n",
875  radius.depth,
876  dynamics.n_initial_relax);
877  for( i=iteration-1; i<iterations.iter_malloc; ++i )
878  /* set stopping radius to current radius, this stops
879  * dynamical solutions from overrunning the upstream scale
880  * and extending to very large radius due to unphysical heat
881  * appearing to advect into region */
883 
884  /* in the event time dependent cooling is required,
885  * but no timestep has been given,
886  * make a guess of the initial timestep based
887  * on the shortest cooling time in the model */
888  if( lgNeedTimestep() )
889  {
891  }
892  }
893 
894  dynamics.DivergePresInteg = 0.;
895 
896  /* This routine is only called if advection is turned on at the end of
897  * each iteration. The test
898  * is to check whether wind velocity is also set by dynamics code */
899 
900  /* !dynamics.lgStatic true - a true dynamical model */
901  if( !dynamics.lgTimeDependentStatic )
902  {
903  if(iteration == dynamics.n_initial_relax+1 )
904  {
905  /* let model settle down for n_initial_relax iterations before we begin
906  * dynamics.n_initial_relax set with "set dynamics relax"
907  * command - it gives the first iteration where we do dynamics
908  * note the off by one confusion - relax is 2 by default,
909  * want to do two static iterations then start dynamics
910  * iteration was incremented before this call so iteration == 2
911  * at end of first iteration */
912  if( dynamics.AdvecLengthInit> 0. )
913  {
914  Dyn_dr = dynamics.AdvecLengthInit;
915  }
916  else
917  {
918  /* -ve dynamics.adveclengthlimit sets length as fraction of first iter total depth */
919  Dyn_dr = -dynamics.AdvecLengthInit*radius.depth;
920  }
921 
922  if (wind.windv0 > 0)
923  Dyn_dr = -Dyn_dr;
924 
925  if( dynamics.lgTracePrint )
926  {
927  fprintf(ioQQQ," DynaIterEnd, dr=%.2e \n",
928  Dyn_dr );
929  }
930  }
931  else if(iteration > dynamics.n_initial_relax+1 )
932  {
933  /* evaluate errors and update Dyn_dr */
934  DynaNewStep();
935  }
936  }
937  else
938  {
939  /* this is time-dependent static model */
940  static double HeatInitial=-1. , HeatRadiated=-1. ,
941  DensityInitial=-1;
942  /* n_initial_relax is number of time-steady models before we start
943  * to evolve, set with "set dynamics relax" command */
944  Dyn_dr = 0.;
945  if( lgPrintDynamics )
946  fprintf(ioQQQ,
947  "DEBUG times enter dynamics.timestep %.2e elapsed_time %.2e iteration %li relax %li \n",
948  dynamics.timestep ,
949  dynamics.time_elapsed,
950  iteration , dynamics.n_initial_relax);
951  if( iteration > dynamics.n_initial_relax )
952  {
953  /* evaluate errors */
954  DynaNewStep();
955 
956  /* this is set true on CORONAL XXX TIME INIT command, says to use constant
957  * temperature for first n_initial_relax iterations, then let run free */
958  if( dynamics.lg_coronal_time_init )
959  {
961  thermal.ConstTemp = 0.;
962  }
963 
964  /* time variable branch, now without dynamics */
965  /* elapsed time - don't increase dynamics.time_elapsed during first two
966  * two iterations since this sets static model */
967  dynamics.time_elapsed += dynamics.timestep;
968  /* dynamics.timestep_stop is -1 if not set with explicit stop time */
969  if( dynamics.timestep_stop > 0 && dynamics.time_elapsed > dynamics.timestep_stop )
970  {
971  dynamics.lgStatic_completed = true;
972  }
973 
974  /* stop lowest temperature time command */
977  dynamics.lgStatic_completed = true;
978 
979  /* this is heat radiated, after correction for change of H density in constant
980  * pressure cloud */
981  HeatRadiated += (thermal.ctot-dynamics.Cool()) * dynamics.timestep *
982  (DensityInitial / scalingDensity());
983  }
984  else
985  {
986  /* this branch, during initial relaxation of solution */
987  HeatInitial = 1.5 * pressure.PresGasCurr;
988  HeatRadiated = 0.;
989  DensityInitial = scalingDensity();
990  if( lgPrintDynamics )
991  fprintf(ioQQQ,"DEBUG relaxing times requested %li this is step %li\n",
992  dynamics.n_initial_relax , iteration);
993  }
994  if( lgPrintDynamics )
995  fprintf(ioQQQ,"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n",
996  HeatInitial , HeatRadiated );
997 
998  if( dynamics.lgStatic_completed )
999  {
1000  // Do nothing but talk.
1001  if( lgPrintDynamics )
1002  fprintf(ioQQQ,"DEBUG lgtimes -- stop time reached.\n" );
1003  }
1004  /* at this point dynamics.time_elapsed is the time at the end of the
1005  * previous iteration. We need dt for the next iteration */
1006  else if( time_elapsed_time.size() > 0 && dynamics.time_elapsed > time_elapsed_time[nTime_dt_array_element] )
1007  {
1008  /* time has advanced to next table point,
1009  * set dynamics.timestep to specified value */
1010  ++nTime_dt_array_element;
1011  /* this is an assert since we already qualified the array
1012  * element above */
1013  ASSERT( nTime_dt_array_element < nTime_flux );
1014 
1015  /* option to set flag for recombination logic */
1016  if( lgtime_Recom[nTime_dt_array_element] )
1017  {
1018  if( lgPrintDynamics )
1019  fprintf(ioQQQ,"DEBUG dynamics turn on recombination logic\n");
1020  dynamics.lgRecom = true;
1021  /* set largest possible zone thickness to value on previous
1022  * iteration when light was on - during recombination conditions
1023  * become much more homogeneous and dr can get too large,
1024  * crashing into H i-front */
1026  radius.lgSdrmaxRel = false;
1027  }
1028 
1029  if( lgtime_dt_specified )
1030  {
1031  if( lgPrintDynamics )
1032  fprintf(ioQQQ,"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element,
1033  dynamics.timestep);
1034  /* this is the new dynamics.timestep */
1035  dynamics.timestep = time_dt[nTime_dt_array_element];
1036  /* option to change time step factor - default is 1.2 set in DynaZero */
1037  if( time_dt_scale_factor[nTime_dt_array_element] > 0. )
1038  dynamics.timestep_factor = time_dt_scale_factor[nTime_dt_array_element];
1039  }
1040  }
1041  else if( lgtime_dt_specified )
1042  {
1043  /* we are between two points in the table, increase dynamics.timestep */
1044  dynamics.timestep *= dynamics.timestep_factor;
1045  if( lgPrintDynamics )
1046  fprintf(ioQQQ,"DEBUG lgtimes increment Timeby dynamics.timestep_factor to %li %.2e\n" ,
1047  nTime_dt_array_element,
1048  dynamics.timestep );
1049  if( time_elapsed_time.size() > 0 && dynamics.time_elapsed+dynamics.timestep > time_elapsed_time[nTime_dt_array_element] )
1050  {
1051  if( lgPrintDynamics )
1052  fprintf(ioQQQ,"DEBUG lgtimes but reset to %.2e\n" ,dynamics.timestep );
1053  dynamics.timestep = 1.0001*(time_elapsed_time[nTime_dt_array_element]-dynamics.time_elapsed);
1054  }
1055  }
1056  else
1057  {
1058  /* time not specified in third column, so use initial */
1059  dynamics.timestep = timestep_next();
1060  }
1061 
1062  if( cosmology.lgDo )
1063  {
1067  }
1068 
1069  if( lgPrintDynamics )
1070  fprintf(ioQQQ,"DEBUG times exit dynamics.timestep %.2e elapsed_time %.2e scale %.2e ",
1071  dynamics.timestep ,
1072  dynamics.time_elapsed,
1074 
1075  if( cosmology.lgDo )
1076  {
1077  fprintf(ioQQQ,"redshift %.3e ", cosmology.redshift_current );
1078  }
1079 
1080  fprintf(ioQQQ,"\n" );
1081  }
1082 
1083  /* Dyn_dr == 0 is for static time dependent cloud */
1084  ASSERT( (iteration < dynamics.n_initial_relax+1) ||
1085  Dyn_dr != 0. || (Dyn_dr == 0. && wind.lgStatic()) );
1086 
1087  /* reset the upstream counters */
1089  dynamics.discretization_error = 0.;
1090  dynamics.error_scale2 = 0.;
1091 
1092  /* save results from previous iteration */
1093  DynaSaveLast();
1094  return;
1095 }
1096 
1097 /*DynaNewStep work out convergence errors */
1099 {
1100  long int ilast = 0,
1101  i,
1102  nelem,
1103  ion,
1104  mol;
1105 
1106  double frac_next=-BIGFLOAT,
1107  Oldi_density,
1108  Oldi_ion,
1109  Oldi_iso,
1110  Oldi_mol;
1111 
1112  DEBUG_ENTRY( "DynaNewStep()" );
1113 
1114  /*n = MIN2(nzone, NZLIM-1);*/
1115  dynamics.convergence_error = 0;
1116  dynamics.error_scale1 = 0.;
1117 
1118  ASSERT( nzone < struc.nzlim);
1119  for(i=0;i<nzone;++i)
1120  {
1121  /* Interpolate for present position in previous solution */
1122  while( (Old_depth[ilast] < struc.depth[i] ) &&
1123  ( ilast < nOld_zone-1 ) )
1124  {
1125  ++ilast;
1126  }
1127  ASSERT( ilast <= nOld_zone-1 );
1128 
1129  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1130  {
1131  frac_next = ( struc.depth[i] - Old_depth[ilast])/
1132  (Old_depth[ilast+1] - Old_depth[ilast]);
1133  Oldi_density = Old_density[ilast] +
1134  (Old_density[ilast+1] - Old_density[ilast])*
1135  frac_next;
1136  }
1137  else
1138  {
1139  Oldi_density = Old_density[ilast];
1140  }
1141  /* Must be consistent with discretization_error above */
1142  /* >>chngf 02 aug 01, multiply by cell width */
1143  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1144  {
1145  for( ion=0; ion<nelem+2; ++ion )
1146  {
1147  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1148  {
1149  Oldi_ion = (Old_xIonDense[ilast][nelem][ion] +
1150  (Old_xIonDense[ilast+1][nelem][ion]-Old_xIonDense[ilast][nelem][ion])*
1151  frac_next);
1152  }
1153  else
1154  {
1155  Oldi_ion = Old_xIonDense[ilast][nelem][ion];
1156  }
1157  dynamics.convergence_error += POW2((double)Oldi_ion/Oldi_density-struc.xIonDense[i][nelem][ion]/scalingZoneDensity(i)) /* *struc.dr[i] */;
1158 
1159  /* >>chng 02 nov 11, add first error scale estimate from Robin */
1160  //fprintf(ioQQQ,"%g %g %g\n",dynamics.error_scale1,
1161  // struc.xIonDense[nelem][ion][i],scalingZoneDensity(i));
1162  dynamics.error_scale1 += POW2((double)struc.xIonDense[i][nelem][ion]/scalingZoneDensity(i));
1163  }
1164  }
1165  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1166  {
1167  for( nelem=ipISO; nelem<LIMELM; ++nelem)
1168  {
1169  if( dense.lgElmtOn[nelem] )
1170  {
1171  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; ++level )
1172  {
1173  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1174  {
1175  Oldi_iso = (Old_StatesElem[ilast][nelem][nelem-ipISO][level] +
1176  (Old_StatesElem[ilast+1][nelem][nelem-ipISO][level]-Old_StatesElem[ilast][nelem][nelem-ipISO][level])*
1177  frac_next);
1178  }
1179  else
1180  {
1181  Oldi_iso = Old_StatesElem[ilast][nelem][nelem-ipISO][level];
1182  }
1183  dynamics.convergence_error += POW2(Oldi_iso/Oldi_density-struc.StatesElem[i][nelem][nelem-ipISO][level]/struc.hden[i]) /* *struc.dr[i] */;
1184 
1185  /* >>chng 02 nov 11, add first error scale estimate from Robin */
1186  dynamics.error_scale1 += POW2(struc.StatesElem[i][nelem][nelem-ipISO][level]/struc.hden[i]);
1187  }
1188  }
1189  }
1190  }
1191 
1192  for( mol=0; mol < mole_global.num_calc; mol++)
1193  {
1194  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1195  {
1196  Oldi_mol = (Old_molecules[ilast][mol] +
1197  (Old_molecules[ilast+1][mol]-Old_molecules[ilast][mol])*
1198  frac_next);
1199  }
1200  else
1201  {
1202  Oldi_mol = Old_molecules[ilast][mol];
1203  }
1204  dynamics.convergence_error += POW2((double)Oldi_mol/Oldi_density-struc.molecules[i][mol]/scalingZoneDensity(i)) /* *struc.dr[i] */;
1205 
1206  /* >>chng 02 nov 11, add first error scale estimate from Robin
1207  * used to normalize the above convergence_error */
1208  dynamics.error_scale1 += POW2((double)struc.molecules[i][mol]/scalingZoneDensity(i));
1209  }
1210  }
1211 
1212  /* convergence_error is an estimate of the convergence of the solution from its change during the last iteration,
1213  discretization_error is an estimate of the accuracy of the advective terms, calculated in DynaStartZone above:
1214  if dominant error is from the advective terms, need to make them more accurate.
1215  */
1216 
1217  /* report properties of previous iteration */
1218  fprintf(ioQQQ,"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
1219  Dyn_dr, dynamics.convergence_error , dynamics.discretization_error ,
1220  dynamics.error_scale1 , dynamics.error_scale2
1221  );
1222 
1223  /* >>chng 02 nov 29, dynamics.convergence_tolerance is now set to 0.1 in init routine */
1224  if( dynamics.convergence_error < dynamics.convergence_tolerance*dynamics.discretization_error )
1225  Dyn_dr /= 1.5;
1226  return;
1227 }
1228 
1229 /*DynaSaveLast save results from previous iteration */
1231 {
1232  long int i,
1233  ion,
1234  nelem,
1235  mol;
1236 
1237  DEBUG_ENTRY( "DynaSaveLast()" );
1238 
1239  /* Save results from previous iteration */
1240  nOld_zone = nzone;
1241  dynamics.oldFullDepth = struc.depth[nzone-1];
1242  ASSERT( nzone < struc.nzlim );
1243  for( i=0; i<nzone; ++i )
1244  {
1245  Old_histr[i] = struc.histr[i];
1246  Old_depth[i] = struc.depth[i];
1248  /* old n_p density from previous iteration */
1249  Old_hiistr[i] = struc.hiistr[i];
1250  /* old pressure from previous iteration */
1251  Old_pressure[i] = struc.pressure[i];
1252  /* old electron density from previous iteration */
1253  Old_ednstr[i] = struc.ednstr[i];
1254  /* energy term */
1257  Old_DenMass[i] = struc.DenMass[i];
1258 
1259  for(mol=0;mol<mole_global.num_calc;mol++)
1260  {
1261  Old_molecules[i][mol] = struc.molecules[i][mol];
1262  }
1263 
1264  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1265  {
1266  Old_gas_phase[i][nelem] = dense.gas_phase[nelem];
1267  for( ion=0; ion<nelem+2; ++ion )
1268  {
1269  Old_xIonDense[i][nelem][ion] = struc.xIonDense[i][nelem][ion];
1270  }
1271  }
1272  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1273  {
1274  for( nelem=ipISO; nelem<LIMELM; ++nelem)
1275  {
1276  if( dense.lgElmtOn[nelem] )
1277  {
1278  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
1279  {
1280  Old_StatesElem[i][nelem][nelem-ipISO][level] = struc.StatesElem[i][nelem][nelem-ipISO][level];
1281  ASSERT( !isnan( Old_StatesElem[i][nelem][nelem-ipISO][level] ) );
1282  }
1283  }
1284  }
1285  }
1286  }
1287  return;
1288 }
1289 
1290 realnum DynaFlux(double depth)
1291 {
1292  realnum flux;
1293 
1294  DEBUG_ENTRY( "DynaFlux()" );
1295 
1296  if(dynamics.FluxIndex == 0)
1297  {
1298  flux = (realnum)dynamics.FluxScale;
1299  }
1300  else
1301  {
1302  flux = (realnum)(dynamics.FluxScale*pow(fabs(depth-dynamics.FluxCenter),dynamics.FluxIndex));
1303  if(depth < dynamics.FluxCenter)
1304  flux = -flux;
1305  }
1306  if(dynamics.lgFluxDScale)
1307  {
1308  /*flux *= struc.DenMass[0]; */
1309  /* WJH 21 may 04, changed to use dense.xMassDensity0, which should be strictly constant */
1310  flux *= dense.xMassDensity0;
1311  }
1312  return flux;
1313 }
1314 
1315 /* ============================================================================== */
1316 /*DynaZero zero some dynamics variables, called from zero.c,
1317  * before parsing commands */
1318 void t_dynamics::zero( void )
1319 {
1320  int ipISO;
1321 
1322  DEBUG_ENTRY( "t_dynamics::zero()" );
1323 
1324  /* the number of zones in the previous iteration */
1325  nOld_zone = 0;
1326 
1327  /* by default advection is turned off */
1328  lgAdvection = false;
1329  /*Velocity = 0.;*/
1330  AdvecSpecificEnthalpy = 0.;
1331  Cool_r = 0.;
1332  Heat_v = 0.;
1333  dHeatdT = 0.;
1334  HeatMax = 0.;
1335  CoolMax = 0.;
1336  Rate = 0.;
1337 
1338  /* sets recombination logic, keyword RECOMBINATION on a time step line */
1339  lgRecom = false;
1340 
1341  /* don't force populations to equilibrium levels */
1342  lgEquilibrium = false;
1343 
1344  /* set true if time dependent calculation is finished */
1345  lgStatic_completed = false;
1346 
1347  /* vars that determine whether time dependent soln only - set with time command */
1348  lgTimeDependentStatic = false;
1349  timestep_init = -1.;
1350  /* this factor multiplies the time step */
1351  timestep_factor = 1.2;
1352  time_elapsed = 0.;
1353 
1354  /* set the first iteration to include dynamics rather than constant pressure */
1355  /* iteration number, initial iteration is 1, default is 3 - changed with SET DYNAMICS FIRST command */
1356  n_initial_relax = 3;
1357 
1358  /* set initial value of the advection length,
1359  * neg => fraction of depth of init model, + length cm */
1360  AdvecLengthInit = -0.1;
1361 
1362  /* this is a tolerance for determining whether dynamics has converged */
1363  convergence_tolerance = 0.1;
1364 
1365  /* this says that set dynamics pressure mode was set */
1366  lgSetPresMode = false;
1367 
1368  /* set default values for uniform mass flux */
1369  FluxScale = 0.;
1370  lgFluxDScale = false;
1371  FluxCenter = 0.;
1372  FluxIndex = 0.;
1373  dRad = BIGFLOAT;
1374 
1375  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1376  {
1377  /* factor to allow turning off advection for one of the iso seq,
1378  * this is done with command "no advection h-like" or "he-like"
1379  * only for testing */
1380  lgISO[ipISO] = true;
1381  }
1382  /* turn off advection for rest of ions, command "no advection metals" */
1383  lgMETALS = true;
1384  /* turn off thermal effects of advection, command "no advection cooling" */
1385  lgCoolHeat = true;
1386  DivergePresInteg = 0.;
1387 
1388  discretization_error = 0.;
1389  error_scale2 = 0.;
1390  Upstream_density = 0.;
1391  return;
1392 }
1393 
1394 
1395 /* ============================================================================== */
1396 /* DynaCreateArrays allocate some space needed to save the dynamics structure variables,
1397  * called from DynaCreateArrays */
1398 void DynaCreateArrays( void )
1399 {
1400  long int nelem,
1401  ns,
1402  i,
1403  ion,
1404  mol;
1405 
1406  DEBUG_ENTRY( "DynaCreateArrays()" );
1407 
1409  dynamics.molecules.resize(mole_global.num_calc);
1410 
1411  UpstreamElem.resize(LIMELM);
1412 
1413  dynamics.Source = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
1414  UpstreamIon = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
1415  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1416  {
1417  dynamics.Source[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
1418  UpstreamIon[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
1419  for( ion=0; ion<nelem+2; ++ion )
1420  {
1421  dynamics.Source[nelem][ion] = 0.;
1422  }
1423  }
1424 
1425  UpstreamStatesElem = ((double***)MALLOC( (size_t)LIMELM*sizeof(double **) ));
1426  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1427  {
1428  if( dense.lgElmtOn[nelem] )
1429  {
1430  UpstreamStatesElem[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
1431  for( long ion=0; ion<nelem+1; ion++ )
1432  {
1433  long ipISO = nelem-ion;
1434  if( ipISO < NISO )
1435  {
1436  UpstreamStatesElem[nelem][nelem-ipISO] = (double*)MALLOC(sizeof(double)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
1437  }
1438  else
1439  {
1440  fixit("for now, point non-iso ions to NULL");
1441  UpstreamStatesElem[nelem][nelem-ipISO] = NULL;
1442  }
1443  }
1444  }
1445  }
1446 
1447 
1448  dynamics.StatesElem = ((double***)MALLOC( (size_t)LIMELM*sizeof(double **) ));
1449  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1450  {
1451  if( dense.lgElmtOn[nelem] )
1452  {
1453  dynamics.StatesElem[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
1454  for( long ion=0; ion<nelem+1; ion++ )
1455  {
1456  long ipISO = nelem-ion;
1457  if( ipISO < NISO )
1458  {
1459  dynamics.StatesElem[nelem][nelem-ipISO] = (double*)MALLOC(sizeof(double)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
1460  }
1461  else
1462  {
1463  fixit("for now, point non-iso ions to NULL");
1464  dynamics.StatesElem[nelem][nelem-ipISO] = NULL;
1465  }
1466  }
1467  }
1468  }
1469 
1470  dynamics.Rate = 0.;
1471 
1473 
1474  Old_ednstr.resize(struc.nzlim);
1475 
1476  EnthalpyDensity.resize(struc.nzlim);
1477 
1478  Old_DenMass.resize(struc.nzlim);
1479 
1480  Old_density.resize(struc.nzlim);
1481 
1482  Old_pressure.resize(struc.nzlim);
1483 
1484  Old_histr.resize(struc.nzlim);
1485 
1486  Old_hiistr.resize(struc.nzlim);
1487 
1488  Old_depth.resize(struc.nzlim);
1489 
1490  Old_xLyman_depth.resize(struc.nzlim);
1491 
1492  Old_xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) );
1493 
1494  Old_StatesElem = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(struc.nzlim) );
1495 
1496  Old_gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
1497 
1498  Old_molecules = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
1499 
1500  /* now create diagonal of space for ionization arrays */
1501  for( ns=0; ns < struc.nzlim; ++ns )
1502  {
1503  Old_xIonDense[ns] =
1504  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM) );
1505 
1506  Old_StatesElem[ns] =
1507  (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(LIMELM) );
1508 
1509  Old_gas_phase[ns] =
1510  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM) );
1511 
1512  if (mole_global.num_calc != 0)
1513  {
1514  Old_molecules[ns] =
1515  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(mole_global.num_calc) );
1516  }
1517  else
1518  {
1519  Old_molecules[ns] = NULL;
1520  }
1521 
1522  for( nelem=0; nelem<LIMELM; ++nelem )
1523  {
1524  Old_xIonDense[ns][nelem] =
1525  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+1) );
1526  }
1527 
1528  for( nelem=0; nelem< LIMELM; ++nelem )
1529  {
1530  if( dense.lgElmtOn[nelem] )
1531  {
1532  Old_StatesElem[ns][nelem] =
1533  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+1) );
1534  for( ion=0; ion<nelem+1; ion++ )
1535  {
1536  long ipISO = nelem-ion;
1537  if( ipISO < NISO )
1538  {
1539  Old_StatesElem[ns][nelem][ion] =
1540  (realnum*)MALLOC(sizeof(realnum)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
1541  }
1542  else
1543  {
1544  fixit("for now, point non-iso ions to NULL");
1545  Old_StatesElem[ns][nelem][ion] = NULL;
1546  }
1547  }
1548  }
1549  }
1550  }
1551 
1552  for( i=0; i < struc.nzlim; i++ )
1553  {
1554  /* these are values if H0 and tau_912 from previous iteration */
1555  Old_histr[i] = 0.;
1556  Old_xLyman_depth[i] = 0.;
1557  Old_depth[i] = 0.;
1558  dynamics.oldFullDepth = 0.;
1559  /* old n_p density from previous iteration */
1560  Old_hiistr[i] = 0.;
1561  /* old pressure from previous iteration */
1562  Old_pressure[i] = 0.;
1563  /* old electron density from previous iteration */
1564  Old_ednstr[i] = 0.;
1565  Old_density[i] = 0.;
1566  Old_DenMass[i] = 0.;
1567  Old_EnthalpyDensity[i] = 0.;
1568  for( nelem=0; nelem<LIMELM; ++nelem )
1569  {
1570  for( ion=0; ion<LIMELM+1; ++ion )
1571  {
1572  Old_xIonDense[i][nelem][ion] = 0.;
1573  }
1574  }
1575  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1576  {
1577  for( nelem=ipISO; nelem<LIMELM; ++nelem)
1578  {
1579  if( dense.lgElmtOn[nelem] )
1580  {
1581  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
1582  {
1583  Old_StatesElem[i][nelem][nelem-ipISO][level] = 0.;
1584  }
1585  }
1586  }
1587  }
1588 
1589  for(mol=0;mol<mole_global.num_calc;mol++)
1590  {
1591  Old_molecules[i][mol] = 0.;
1592  }
1593  }
1594  return;
1595 }
1596 
1597 /*advection_set_default - called to set default conditions
1598  * when time and wind commands are parsed,
1599  * lgWind is true if dynamics, false if time dependent */
1600 STATIC void advection_set_default( bool lgWind )
1601 {
1602 
1603  DEBUG_ENTRY( "advection_set_default()" );
1604 
1605  /* turn off prediction of next zone's temperature, as guessed in ZoneStart,
1606  * also set with no tepredictor */
1607  thermal.lgPredNextTe = false;
1608 
1609  /* use the new temperature solver
1610  strcpy( conv.chSolverEden , "new" ); */
1611 
1612  /* constant total pressure, gas+rad+incident continuum
1613  * turn on radiation pressure */
1616  pressure.lgPres_ram_ON = true;
1617 
1618  /* we need to make the solvers much more exact when advection is in place */
1619  if( lgWind )
1620  {
1621  /* turn on advection */
1622  dynamics.lgAdvection = true;
1623 
1624  /* increase precision of solution */
1625  conv.EdenErrorAllowed = 1e-3;
1626  /* the actual relative error is relative to the total heating and cooling,
1627  * which include the dynamics.heat() and .cool(), which are the advected heating/cooling.
1628  * the two terms can be large and nearly cancel, what is written to the .heat() and cool files
1629  * by save files has had the smaller of the two subtracted, leaving only the net advected
1630  * heating and cooling */
1631  conv.HeatCoolRelErrorAllowed = 3e-4f;
1632  conv.PressureErrorAllowed = 1e-3f;
1633  }
1634 
1635  return;
1636 }
1637 
1638 /* ============================================================================== */
1639 /* lgNeedTimestep check if an initial timestep is required */
1641 {
1642  DEBUG_ENTRY( "lgNeedTimestep()" );
1643 
1644  return (dynamics.lgTimeDependentStatic && dynamics.lgRecom && dynamics.timestep_init <= 0.);
1645 }
1646 
1647 /* ============================================================================== */
1648 /* InitDynaTimestep - initialize timestep based on shortest zone cooling time; used
1649  * with 'coronal init time' without 'time first timestep' */
1651 {
1652  DEBUG_ENTRY( "InitDynaTimestep()" );
1653 
1654  /* set default flags - false says that time dependent, not dynamical solution */
1655  advection_set_default(false);
1656 
1657  wind.windv0 = 0.;
1658  wind.setStatic();
1659  wind.windv = wind.windv0;
1660 
1662  dynamics.timestep_init = dynamics.timestep = 1e-4 * timesc.time_therm_short;
1663 
1664  return;
1665 }
1666 
1667 /* ============================================================================== */
1668 /* ParseDynaTime parse the time command, called from ParseCommands */
1670 {
1671  DEBUG_ENTRY( "ParseDynaTime()" );
1672 
1673  /*flag set true when time dependent only */
1674  dynamics.lgTimeDependentStatic = true;
1675 
1676  dynamics.timestep_init = p.getNumberCheckAlwaysLogLim("dynamics.timestep",30.);
1677 
1678  dynamics.timestep = dynamics.timestep_init;
1679  if( p.nMatch( "TRAC" ) )
1680  dynamics.lgTracePrint = true;
1681 
1682  /* this is the stop time and is optional */
1683  dynamics.timestep_stop = p.getNumberDefaultAlwaysLog("stop time", -1.);
1684 
1685  /* set default flags - false says that time dependent, not dynamical solution */
1686  advection_set_default(false);
1687 
1688  wind.windv0 = 0.;
1689  wind.setStatic();
1690  wind.windv = wind.windv0;
1691 
1692  /* create time step and flux arrays */
1693  time_elapsed_time.resize(NTIME);
1694  time_flux_ratio.resize(NTIME);
1695  time_dt.resize(NTIME);
1696  time_dt_scale_factor.resize(NTIME);
1697  lgtime_Recom.resize(NTIME);
1698 
1699  /* number of lines we will save */
1700  nTime_flux = 0;
1701 
1702  /* get the next line, and check for eof */
1703  p.getline();
1704  if( p.m_lgEOF )
1705  {
1706  fprintf( ioQQQ,
1707  " Hit EOF while reading time-continuum list; use END to end list.\n" );
1709  }
1710 
1711  /* third column might set dt - if any third column is missing then
1712  * this is set false and only time on command line is used */
1713  lgtime_dt_specified = true;
1714 
1715  while( ! p.hasCommand("END") )
1716  {
1717  if( nTime_flux >= NTIME )
1718  {
1719  fprintf( ioQQQ,
1720  " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n",
1721  NTIME );
1723  }
1724 
1725  if( p.nMatch("CYCLE") )
1726  {
1727  double period = p.getNumberCheckAlwaysLog("log time");
1728  ASSERT( period > time_elapsed_time[nTime_flux-1] );
1729  long pointsPerPeriod = nTime_flux;
1730  while( nTime_flux < NTIME - 1 )
1731  {
1732  time_elapsed_time[nTime_flux] = period + time_elapsed_time[nTime_flux-pointsPerPeriod];
1734  time_dt[nTime_flux] = time_dt[nTime_flux-pointsPerPeriod];
1736  nTime_flux++;
1737  }
1738  //Tell the code to continue cyclically by equating two named time points
1739  fprintf( ioQQQ, " Adding cycles with period = %e s.\n", period );
1740 
1741  /* get next line and check for eof */
1742  p.getline();
1743  if( p.m_lgEOF )
1744  {
1745  fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
1747  }
1748  continue;
1749  }
1750 
1752  if( nTime_flux >= 1 )
1754  time_flux_ratio[nTime_flux] = p.getNumberCheckAlwaysLog("log flux ratio");
1755 
1756  /* this is optional dt to set time step - if not given then initial
1757  * time step is always used */
1758  time_dt[nTime_flux] = p.getNumberDefaultAlwaysLog("log time step",-1.);
1759 
1760  /* if any of these are not specified then do not use times array */
1761  if( time_dt[nTime_flux] < 0.0 )
1762  lgtime_dt_specified = false;
1763 
1764  /* this is optional scale factor to increase time */
1766  "scale factor to increase time",-1.);
1767 
1768  /* turn on recombination front logic */
1769  if( p.nMatch("RECOMBIN") )
1770  {
1771  /* this sets flag dynamics.lgRecom true so that all of code knows recombination
1772  * is taking place */
1773  lgtime_Recom[nTime_flux] = true;
1774  }
1775  else
1776  {
1777  lgtime_Recom[nTime_flux] = false;
1778  }
1779 
1780  /* this is total number stored so far */
1781  ++nTime_flux;
1782 
1783  /* get next line and check for eof */
1784  p.getline();
1785  if( p.m_lgEOF )
1786  {
1787  fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
1789  }
1790 
1791  }
1792 
1793  if( nTime_flux < 2 )
1794  {
1795  fprintf( ioQQQ, " At least two instances of time must be specified. There is an implicit instance at t=0.\n"
1796  " The user must specify at least one additional time. Sorry.\n" );
1798  }
1799 
1800  if( lgPrintDynamics )
1801  {
1802  for( long i=0; i < nTime_flux; i++ )
1803  {
1804  fprintf( ioQQQ, "DEBUG time dep %.2e %.2e %.2e %.2e\n",
1805  time_elapsed_time[i],
1806  time_flux_ratio[i] ,
1807  time_dt[i],
1809  }
1810  fprintf( ioQQQ, "\n" );
1811  }
1812  return;
1813 }
1814 /* ============================================================================== */
1815 /* ParseDynaWind parse the wind command, called from ParseCommands */
1817 {
1818  int iVelocity_Type;
1819  bool lgModeSet=false;
1820  /* compiler flagged possible paths where dfdr used but not set -
1821  * this is for safety/keep it happy */
1822  double dfdr=-BIGDOUBLE;
1823 
1824  DEBUG_ENTRY( "ParseDynaWind()" );
1825 
1826  if( p.nMatch( "TRAC" ) )
1827  dynamics.lgTracePrint = true;
1828 
1829  /* Flag for type of velocity law:
1830  * 1 is original, give initial velocity at illuminated face
1831  * 2 is face flux gradient (useful if face velocity is zero),
1832  * set to zero, but will be reset if velocity specified */
1833  iVelocity_Type = 0;
1834  /* wind structure, parameters are initial velocity and optional mass
1835  * v read in in km s-1 and convert to cm s-1, mass in solar masses */
1836  if( p.nMatch( "VELO" ) )
1837  {
1838  wind.windv0 = (realnum)(p.getNumberPlain("velocity")*1e5);
1839  wind.windv = wind.windv0;
1840  wind.setDefault();
1841  iVelocity_Type = 1;
1842  }
1843 
1844  if( p.nMatch( "BALL" ) )
1845  {
1846  wind.setBallistic();
1847  lgModeSet = true;
1848  }
1849 
1850  if( p.nMatch( "STAT" ) )
1851  {
1852  wind.windv0 = 0.;
1853  wind.setStatic();
1854  lgModeSet = true;
1855  iVelocity_Type = 1;
1856  }
1857 
1858  if ( 1 == iVelocity_Type && !lgModeSet)
1859  {
1860  if (wind.windv0 > 0.)
1861  {
1862  fprintf(ioQQQ,"Warning, BALListic option needed to switch off pressure gradient terms\n");
1863  }
1864  else if (wind.windv0 == 0.)
1865  {
1866  fprintf(ioQQQ,"Warning, STATic option needed for zero speed solutions\n");
1867  }
1868  }
1869 
1870  if( p.nMatch("DFDR") )
1871  {
1872  /* velocity not specified, rather mass flux gradient */
1873  dfdr = p.getNumberPlain("flux gradient");
1874  iVelocity_Type = 2;
1875  }
1876 
1877  /* center option, gives xxx */
1878  if( p.nMatch("CENT") )
1879  {
1880  /* physical length in cm, can be either sign */
1881  dynamics.FluxCenter = p.getNumberPlain(
1882  "centre of mass flux distribution");
1883  }
1884 
1885  /* flux index */
1886  if( p.nMatch("INDE") )
1887  {
1888  /* power law index */
1889  dynamics.FluxIndex = p.getNumberPlain(
1890  "power law index of mass flux distribution");
1891  }
1892 
1893  /* the case where velocity was set */
1894  if(iVelocity_Type == 1)
1895  {
1896  /* was flux index also set? */
1897  if(dynamics.FluxIndex == 0)
1898  {
1899  dynamics.FluxScale = wind.windv0;
1900  dynamics.lgFluxDScale = true;
1901  /* Center doesn't mean much in this case -- make sure it's
1902  * in front of grid so DynaFlux doesn't swap signs where
1903  * it shouldn't */
1904  dynamics.FluxCenter = -1.;
1905  }
1906  else
1907  {
1910  /* velocity was set but flux index was not set - estimate it */
1911  dynamics.FluxScale = wind.windv0*
1912  pow(fabs(dynamics.FluxCenter),-dynamics.FluxIndex);
1913 
1914  dynamics.lgFluxDScale = true;
1915  if(dynamics.FluxCenter > 0)
1916  {
1917  dynamics.FluxScale = -dynamics.FluxScale;
1918  }
1919  }
1920  }
1921  /* the case where flux gradient is set */
1922  else if(iVelocity_Type == 2)
1923  {
1924  if(dynamics.FluxIndex == 0)
1925  {
1926  fprintf(ioQQQ,"Can't specify gradient when flux is constant!\n");
1927  /* use this exit handler, which closes out MPI when multiprocessing */
1929  }
1932  /* Can't specify FluxScale from dvdr rather than dfdr, as
1933  * d(rho)/dr != 0 */
1934  dynamics.FluxScale = dfdr/dynamics.FluxIndex*
1935  pow(fabs(dynamics.FluxCenter),1.-dynamics.FluxIndex);
1936  if(dynamics.FluxCenter > 0)
1937  {
1938  dynamics.FluxScale = -dynamics.FluxScale;
1939  }
1940  dynamics.lgFluxDScale = false;
1941 
1942  /* put in bogus value simply as flag -- assume that surface velocity
1943  * is small or we wouldn't be using this to specify. */
1944  wind.windv0 = -0.01f;
1945  wind.setDefault();
1946  }
1947  else
1948  {
1949  /* assume old-style velocity-only specification */
1950  /* wind structure, parameters are initial velocity and optional mass
1951  * v in km/sec, mass in solar masses */
1952  wind.windv0 = (realnum)(p.getNumberCheck("wind velocity")*1e5);
1953  if (wind.windv0 < 0.)
1954  {
1955  wind.setDefault();
1956  }
1957  else if (wind.windv0 > 0.)
1958  {
1959  wind.setBallistic();
1960  }
1961  else
1962  {
1963  wind.setStatic();
1964  }
1965 
1966  dynamics.FluxScale = wind.windv0;
1967  dynamics.FluxIndex = 0.;
1968  dynamics.lgFluxDScale = true;
1969  /* Center doesn't mean much in this case -- make sure it's
1970  * in front of grid so DynaFlux doesn't swap signs where
1971  * it shouldn't */
1972  dynamics.FluxCenter = -1.;
1973  }
1974 
1975  wind.windv = wind.windv0;
1976 
1977 # ifdef FOO
1978  fprintf(ioQQQ,"Scale %g (*%c) Index %g Center %g\n",
1979  dynamics.FluxScale,(dynamics.lgFluxDScale)?'D':'1',
1980  dynamics.FluxIndex,dynamics.FluxCenter);
1981 # endif
1982 
1983  /* option to include advection */
1984  if( p.nMatch( "ADVE" ) )
1985  {
1986  /* set default flags - true says dynamical solution */
1987  advection_set_default(true);
1988  strcpy( dense.chDenseLaw, "DYNA" );
1989  }
1990 
1991  else
1992  {
1993  /* this is usual hypersonic outflow */
1994  if( wind.windv0 <= 1.e6 )
1995  {
1996  /* speed of sound roughly 10 km/s */
1997  fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
1998  wind.lgWindOK = false;
1999  }
2000 
2001  /* set the central object mass, in solar masses */
2002  wind.comass = (realnum)p.getNumberDefault("central object mass",1.);
2003  /* default is one solar mass */
2004 
2005  /* option for rotating disk, keyword is disk */
2006  wind.lgDisk = false;
2007  if( p.nMatch( "DISK") )
2008  wind.lgDisk = true;
2009 
2010  strcpy( dense.chDenseLaw, "WIND" );
2011  }
2012 
2013 
2014  /* option to turn off continuum radiative acceleration */
2015  if( p.nMatch("NO CO") )
2016  {
2017  pressure.lgContRadPresOn = false;
2018  }
2019  else
2020  {
2021  pressure.lgContRadPresOn = true;
2022  }
2023  return;
2024 }
2025 
2026 /*DynaPrtZone called to print zone results */
2027 void DynaPrtZone( void )
2028 {
2029 
2030  DEBUG_ENTRY( "DynaPrtZone()" );
2031 
2032  ASSERT( nzone>0 && nzone<struc.nzlim );
2033 
2034  if( nzone > 0 )
2035  {
2036  fprintf(ioQQQ," DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
2038  wind.windv/1e5 ,
2039  dynamics.Cool()/thermal.ctot,
2040  dynamics.Heat()/thermal.ctot);
2041  }
2042 
2043  ASSERT( EnthalpyDensity[nzone-1] > 0. );
2044 
2045  fprintf(ioQQQ," DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
2050  5./2.*pressure.PresGasCurr ,
2052  );
2053  return;
2054 }
2055 
2056 /*DynaPunchTimeDep - save info about time dependent solution */
2057 void DynaPunchTimeDep( FILE* ipPnunit , const char *chJob )
2058 {
2059 
2060  DEBUG_ENTRY( "DynaPunchTimeDep()" );
2061 
2062  if( strcmp( chJob , "END" ) == 0 )
2063  {
2064  double te_mean,
2065  H2_mean,
2066  H0_mean,
2067  Hp_mean,
2068  Hep_mean;
2069  /* save info at end */
2070  if( cdTemp(
2071  /* four char string, null terminated, giving the element name */
2072  "HYDR",
2073  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2074  * 0 means that chLabel is a special case */
2075  2,
2076  /* will be temperature */
2077  &te_mean,
2078  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2079  "RADIUS" ) )
2080  {
2081  TotalInsanity();
2082  }
2083  if( cdIonFrac(
2084  /* four char string, null terminated, giving the element name */
2085  "HYDR",
2086  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2087  * 0 says special case */
2088  2,
2089  /* will be fractional ionization */
2090  &Hp_mean,
2091  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2092  "RADIUS" ,
2093  /* if true then weighting also has electron density, if false then only volume or radius */
2094  false ) )
2095  {
2096  TotalInsanity();
2097  }
2098  if( cdIonFrac(
2099  /* four char string, null terminated, giving the element name */
2100  "HYDR",
2101  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2102  * 0 says special case */
2103  1,
2104  /* will be fractional ionization */
2105  &H0_mean,
2106  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2107  "RADIUS" ,
2108  /* if true then weighting also has electron density, if false then only volume or radius */
2109  false ) )
2110  {
2111  TotalInsanity();
2112  }
2113  if( cdIonFrac(
2114  /* four char string, null terminated, giving the element name */
2115  "H2 ",
2116  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2117  * 0 says special case */
2118  0,
2119  /* will be fractional ionization */
2120  &H2_mean,
2121  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2122  "RADIUS" ,
2123  /* if true then weighting also has electron density, if false then only volume or radius */
2124  false ) )
2125  {
2126  TotalInsanity();
2127  }
2128  if( cdIonFrac(
2129  /* four char string, null terminated, giving the element name */
2130  "HELI",
2131  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2132  * 0 says special case */
2133  2,
2134  /* will be fractional ionization */
2135  &Hep_mean,
2136  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2137  "RADIUS" ,
2138  /* if true then weighting also has electron density, if false then only volume or radius */
2139  false ) )
2140  {
2141  TotalInsanity();
2142  }
2143  fprintf( ipPnunit ,
2144  "%.12e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" ,
2145  dynamics.time_elapsed ,
2146  dynamics.timestep ,
2148  scalingDensity(),
2149  te_mean ,
2150  Hp_mean ,
2151  H0_mean ,
2152  H2_mean ,
2153  Hep_mean ,
2154  /* ratio of CO to total H column densities */
2158  );
2159  }
2160  else
2161  TotalInsanity();
2162  return;
2163 }
2164 
2165 /*DynaSave save dynamics - info related to advection */
2166 void DynaSave(FILE* ipPnunit , char chJob )
2167 {
2168  DEBUG_ENTRY( "DynaSave()" );
2169 
2170  if( chJob=='a' )
2171  {
2172  /* this is save dynamics advection, the only save dynamics */
2173  fprintf( ipPnunit , "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2175  thermal.htot ,
2176  dynamics.Cool() ,
2177  dynamics.Heat() ,
2178  dynamics.dCooldT() ,
2179  dynamics.Source[ipHYDROGEN][ipHYDROGEN],
2180  dynamics.Rate,
2183  );
2184  }
2185  else
2186  TotalInsanity();
2187  return;
2188 }
2189 
2190 #define MERGE 0
2192 {
2193  double heat = Heat_v*scalingDensity();
2194  if (MERGE)
2195  {
2196  double cool = Cool_r*phycon.EnthalpyDensity;
2197  if (heat > cool)
2198  return heat-cool;
2199  else
2200  return 0.;
2201  }
2202  return heat;
2203 }
2204 
2206 {
2207  double cool = Cool_r*phycon.EnthalpyDensity;
2208  if (MERGE)
2209  {
2210  double heat = Heat_v*scalingDensity();
2211  if (heat < cool)
2212  return cool-heat;
2213  else
2214  return 0.;
2215  }
2216  return cool;
2217 }
2218 #undef MERGE
2219 
2221 {
2222  return Cool_r*5./2.*pressure.PresGasCurr/phycon.te;
2223 }
2224 
2225 void DynaIterStart(void)
2226 {
2227  DEBUG_ENTRY( "DynaIterStart()" );
2228 
2229  if( 0 == nTime_flux || iteration <= dynamics.n_initial_relax )
2230  {
2232  return;
2233  }
2234  else if( dynamics.time_elapsed <= time_elapsed_time[0] )
2235  {
2236  /* if very early times not specified assume no flux variation yet */
2238  }
2239  else if( dynamics.time_elapsed > time_elapsed_time[nTime_flux-1] )
2240  {
2241  if( dynamics.lgStatic_completed )
2243  else
2244  {
2245  fprintf( ioQQQ, " PROBLEM - DynaIterStart - I need the continuum at time %.2e but the table ends at %.2e.\n" ,
2248  }
2249  }
2250  else
2251  {
2253  /* the times in seconds */
2254  &time_elapsed_time[0],
2255  /* the rfield.time_continuum_scale factors */
2256  &time_flux_ratio[0],
2257  /* the number of rfield.time_continuum_scale factors */
2258  nTime_flux,
2259  /* the desired time */
2260  dynamics.time_elapsed);
2261  }
2262 
2263  if( lgPrintDynamics )
2264  fprintf(ioQQQ,"DEBUG time dep reset continuum iter %ld dynamics.timestep %.2e elapsed time %.2e scale %.2e",
2265  iteration,
2266  dynamics.timestep ,
2267  dynamics.time_elapsed,
2269  if( dynamics.lgRecom )
2270  {
2271  fprintf(ioQQQ," recom");
2272  }
2273  fprintf(ioQQQ,"\n");
2274 
2275  /* make sure that at least one continuum source is variable */
2276  long int nTimeVary = 0;
2277  for( long int i=0; i < rfield.nShape; i++ )
2278  {
2279  /* this is set true if particular continuum source can vary with time
2280  * set true if TIME appears on intensity / luminosity command line */
2281  if( rfield.lgTimeVary[i] )
2282  ++nTimeVary;
2283  }
2284 
2286  {
2287  /* vary extra heating */
2289  if( lgPrintDynamics )
2290  fprintf(ioQQQ,"DEBUG TurbHeat vary new heat %.2e\n",
2291  hextra.TurbHeat);
2292  }
2293 }
realnum * hden
Definition: struc.h:25
bool nMatch(const char *chKey) const
Definition: parser.h:140
bool lgContRadPresOn
Definition: pressure.h:65
void DynaPrtZone(void)
Definition: dynamics.cpp:2027
t_mole_global mole_global
Definition: mole.cpp:7
long int iter_malloc
Definition: iterations.h:40
void calc_therm_timesc(long izone)
Definition: timesc.cpp:26
bool hasCommand(const char *s2)
Definition: parser.cpp:705
double depth
Definition: radius.h:31
void setDefault(void)
Definition: wind.h:82
double htot
Definition: thermal.h:169
static realnum *** Old_xIonDense
Definition: dynamics.cpp:121
t_thermal thermal
Definition: thermal.cpp:6
static vector< realnum > Old_DenMass
Definition: dynamics.cpp:100
static vector< double > time_elapsed_time
Definition: dynamics.cpp:70
STATIC void advection_set_default(bool lgWind)
Definition: dynamics.cpp:1600
t_colden colden
Definition: colden.cpp:5
static vector< realnum > Old_density
Definition: dynamics.cpp:100
double FluxCenter
Definition: dynamics.h:117
double Cool()
Definition: dynamics.cpp:2205
double EdenErrorAllowed
Definition: conv.h:265
void ParseDynaWind(Parser &p)
Definition: dynamics.cpp:1816
static vector< realnum > Old_xLyman_depth
Definition: dynamics.cpp:100
double convergence_tolerance
Definition: dynamics.h:162
static vector< realnum > Old_pressure
Definition: dynamics.cpp:100
void DynaIterEnd(void)
Definition: dynamics.cpp:854
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
const double BIGDOUBLE
Definition: cpu.h:249
bool lgWindOK
Definition: wind.h:42
static double *** UpstreamStatesElem
Definition: dynamics.cpp:59
int num_calc
Definition: mole.h:362
STATIC void InitDynaTimestep()
Definition: dynamics.cpp:1650
t_struc struc
Definition: struc.cpp:6
static vector< realnum > Old_histr
Definition: dynamics.cpp:100
vector< double > molecules
Definition: dynamics.h:86
void DynaCreateArrays(void)
Definition: dynamics.cpp:1398
static realnum ** Old_molecules
Definition: dynamics.cpp:118
bool lgDisk
Definition: wind.h:74
const realnum SMALLFLOAT
Definition: cpu.h:246
static vector< double > time_dt_scale_factor
Definition: dynamics.cpp:70
const int NISO
Definition: cddefines.h:310
double EnergyIonization
Definition: phycon.h:41
long int IonHigh[LIMELM+1]
Definition: dense.h:130
bool lgFluxDScale
Definition: dynamics.h:138
realnum windv0
Definition: wind.h:11
static vector< double > Upstream_molecules
Definition: dynamics.cpp:64
void DynaIterStart(void)
Definition: dynamics.cpp:2225
static vector< double > UpstreamElem
Definition: dynamics.cpp:61
realnum redshift_step
Definition: cosmology.h:26
void DynaSave(FILE *ipPnunit, char chJob)
Definition: dynamics.cpp:2166
bool lgMETALS
Definition: dynamics.h:92
bool lgAdvection
Definition: dynamics.h:66
bool lgTimeDependentStatic
Definition: dynamics.h:102
bool lgTracePrint
Definition: dynamics.h:183
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
Definition: cddrive.cpp:912
t_StopCalc StopCalc
Definition: stopcalc.cpp:7
void DynaStartZone(void)
Definition: dynamics.cpp:375
t_conv conv
Definition: conv.cpp:5
realnum * ednstr
Definition: struc.h:25
t_hextra hextra
Definition: hextra.cpp:5
static double ** UpstreamIon
Definition: dynamics.cpp:58
realnum ** molecules
Definition: struc.h:71
t_phycon phycon
Definition: phycon.cpp:6
double EnthalpyDensity
Definition: phycon.h:50
double time_therm_short
Definition: timesc.h:29
realnum TurbHeatSave
Definition: hextra.h:42
double getNumberDefaultAlwaysLog(const char *chDesc, double fdef)
Definition: parser.cpp:412
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition: cddrive.cpp:1296
double sdrmax
Definition: radius.h:158
static const bool lgPrintDynamics
Definition: dynamics.cpp:40
bool lgTemperatureConstant
Definition: thermal.h:44
vector< double > StopThickness
Definition: iterations.h:77
FILE * ioQQQ
Definition: cddefines.cpp:7
molezone * findspecieslocal(const char buf[])
long int nzone
Definition: cddefines.cpp:14
static int ipUpstream
Definition: dynamics.cpp:43
#define NTIME
Definition: dynamics.cpp:76
t_dynamics dynamics
Definition: dynamics.cpp:42
realnum time_continuum_scale
Definition: rfield.h:205
#define MIN2(a, b)
Definition: cddefines.h:807
bool lgDo
Definition: cosmology.h:44
Definition: parser.h:42
bool lgISO[NISO]
Definition: dynamics.h:89
realnum PressureErrorAllowed
Definition: conv.h:270
double FluxScale
Definition: dynamics.h:135
bool lgStatic_completed
Definition: dynamics.h:111
double getNumberPlain(const char *chDesc)
Definition: parser.cpp:354
t_dense dense
Definition: global.cpp:15
double error_scale2
Definition: dynamics.h:168
double Heat_v
Definition: dynamics.h:69
static vector< int > lgtime_Recom
Definition: dynamics.cpp:75
bool lgTimeVary[LIMSPC]
Definition: rfield.h:292
realnum xMassDensity0
Definition: dense.h:105
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
bool lgPredNextTe
Definition: thermal.h:40
realnum * DenMass
Definition: struc.h:25
void DynaEndZone(void)
Definition: dynamics.cpp:833
static vector< realnum > Old_hiistr
Definition: dynamics.cpp:100
double Heat()
Definition: dynamics.cpp:2191
static long int nOld_zone
Definition: dynamics.cpp:130
STATIC void DynaNewStep(void)
Definition: dynamics.cpp:1098
Wind wind
Definition: wind.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
long int iteration
Definition: cddefines.cpp:16
STATIC double timestep_next(void)
Definition: dynamics.cpp:137
#define MALLOC(exp)
Definition: cddefines.h:556
double drad
Definition: radius.h:31
realnum * depth
Definition: struc.h:25
long int IonLow[LIMELM+1]
Definition: dense.h:129
realnum redshift_current
Definition: cosmology.h:26
double timestep_stop
Definition: dynamics.h:187
void zero()
Definition: dynamics.cpp:1318
#define POW2
Definition: cddefines.h:983
#define STATIC
Definition: cddefines.h:118
static double Dyn_dr
Definition: dynamics.cpp:91
void setStatic(void)
Definition: wind.h:88
realnum DynaFlux(double depth)
Definition: dynamics.cpp:1290
realnum TurbHeat
Definition: hextra.h:42
t_mole_local mole
Definition: mole.cpp:8
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
void DynaIonize(void)
Definition: dynamics.cpp:160
double timestep
Definition: dynamics.h:187
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
double dHeatdT
Definition: dynamics.h:69
#define EXIT_FAILURE
Definition: cddefines.h:168
realnum **** StatesElem
Definition: struc.h:67
const realnum BIGFLOAT
Definition: cpu.h:244
bool lgPres_magnetic_ON
Definition: pressure.h:91
bool lgSdrmaxRel
Definition: radius.h:166
bool lgEquilibrium
Definition: dynamics.h:180
double dr_max_last_iter
Definition: radius.h:182
double discretization_error
Definition: dynamics.h:165
bool lgTurbHeatVaryTime
Definition: hextra.h:76
#define cdEXIT(FAIL)
Definition: cddefines.h:484
STATIC bool lgNeedTimestep()
Definition: dynamics.cpp:1640
STATIC void DynaSaveLast(void)
Definition: dynamics.cpp:1230
bool fp_bound(sys_float lo, sys_float x, sys_float hi, int n=3)
Definition: cddefines.h:931
double depth_mid_zone
Definition: radius.h:31
t_iterations iterations
Definition: iterations.cpp:6
double column(const genericState &gs)
void ParseDynaTime(Parser &p)
Definition: dynamics.cpp:1669
double PresGasCurr
Definition: pressure.h:46
realnum HeatCoolRelErrorAllowed
Definition: conv.h:276
realnum scalingZoneDensity(long i)
Definition: dense.cpp:416
t_radius radius
Definition: radius.cpp:5
t_timesc timesc
Definition: timesc.cpp:7
double dRad
Definition: dynamics.h:144
void setBallistic(void)
Definition: wind.h:94
bool lgElmtOn[LIMELM]
Definition: dense.h:160
static long int nTime_flux
Definition: dynamics.cpp:79
double CoolMax
Definition: dynamics.h:74
double Cool_r
Definition: dynamics.h:69
long int nzlim
Definition: struc.h:19
realnum gas_phase[LIMELM]
Definition: dense.h:76
realnum TempLoStopIteration
Definition: stopcalc.h:45
double timestep_init
Definition: dynamics.h:187
double getNumberCheckAlwaysLogLim(const char *chDesc, double flim)
Definition: parser.cpp:399
#define ASSERT(exp)
Definition: cddefines.h:617
double sound_speed_adiabatic
Definition: timesc.h:58
static realnum ** Old_gas_phase
Definition: dynamics.cpp:124
char chDenseLaw[5]
Definition: dense.h:176
bool lgSetPresMode
Definition: dynamics.h:172
double error_scale1
Definition: dynamics.h:168
double EnergyBinding
Definition: phycon.h:54
double Rate
Definition: dynamics.h:77
realnum ConstTemp
Definition: thermal.h:56
bool getline()
Definition: parser.cpp:249
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
#define isnan
Definition: cddefines.h:663
t_cosmology cosmology
Definition: cosmology.cpp:8
double AdvecLengthInit
Definition: dynamics.h:114
realnum scalingDensity(void)
Definition: dense.cpp:409
realnum xMassDensity
Definition: dense.h:101
double linint(const double x[], const double y[], long n, double xval)
double getNumberCheckAlwaysLog(const char *chDesc)
Definition: parser.cpp:393
double dCooldT()
Definition: dynamics.cpp:2220
static vector< realnum > EnthalpyDensity
Definition: dynamics.cpp:100
double getNumberDefault(const char *chDesc, double fdef)
Definition: parser.cpp:367
double *** StatesElem
Definition: dynamics.h:83
realnum * histr
Definition: struc.h:25
double eden
Definition: dense.h:201
bool lgCoolHeat
Definition: dynamics.h:95
double EnergyExcitation
Definition: phycon.h:47
realnum * xLyman_depth
Definition: struc.h:25
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
static int iphUpstream
Definition: dynamics.cpp:43
MoleculeList list
Definition: mole.h:365
bool lgRecom
Definition: dynamics.h:108
double ** Source
Definition: dynamics.h:80
realnum GetHubbleFactor(realnum z)
Definition: cosmology.cpp:10
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
double pow(double x, int i)
Definition: cddefines.h:782
double time_elapsed
Definition: dynamics.h:105
long int n_initial_relax
Definition: dynamics.h:132
const int ipCARBON
Definition: cddefines.h:353
bool lgPres_ram_ON
Definition: pressure.h:92
static vector< realnum > Old_EnthalpyDensity
Definition: dynamics.cpp:100
long int nShape
Definition: rfield.h:308
void DynaPunchTimeDep(FILE *ipPnunit, const char *chJob)
Definition: dynamics.cpp:2057
long int numLevels_max
Definition: iso.h:524
bool lgStatic(void) const
Definition: wind.h:24
realnum Upstream_density
Definition: dynamics.h:175
static realnum **** Old_StatesElem
Definition: dynamics.cpp:127
#define MERGE
Definition: dynamics.cpp:2190
double timestep_factor
Definition: dynamics.h:187
#define fixit(a)
Definition: cddefines.h:416
bool m_lgEOF
Definition: parser.h:53
double convergence_error
Definition: dynamics.h:159
double te
Definition: phycon.h:21
static vector< double > time_dt
Definition: dynamics.cpp:70
static vector< double > Old_depth
Definition: dynamics.cpp:97
const int ipHYDROGEN
Definition: cddefines.h:348
realnum * hiistr
Definition: struc.h:25
static double AdvecSpecificEnthalpy
Definition: dynamics.cpp:94
realnum colden[NCOLD]
Definition: colden.h:32
double HeatMax
Definition: dynamics.h:74
double oldFullDepth
Definition: dynamics.h:147
realnum *** xIonDense
Definition: struc.h:64
bool lg_coronal_time_init
Definition: dynamics.h:99
bool lgPres_radiation_ON
Definition: pressure.h:90
long int numLevels_local
Definition: iso.h:529
static vector< realnum > Old_ednstr
Definition: dynamics.cpp:100
realnum windv
Definition: wind.h:18
realnum comass
Definition: wind.h:14
double FluxIndex
Definition: dynamics.h:141
static int ipyUpstream
Definition: dynamics.cpp:43
static bool lgtime_dt_specified
Definition: dynamics.cpp:74
realnum DivergePresInteg
Definition: dynamics.h:177
realnum TempHiStopIteration
Definition: stopcalc.h:38
realnum * pressure
Definition: struc.h:25
static vector< double > time_flux_ratio
Definition: dynamics.cpp:70
double ctot
Definition: thermal.h:130
double getNumberCheck(const char *chDesc)
Definition: parser.cpp:358