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