cloudy  trunk
 All Data Structures 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-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /* DynaEndIter 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 /* DynaPunch punch info related to advection */
15 /* DynaPunch, punch 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 "input.h"
22 #include "colden.h"
23 #include "radius.h"
24 #include "thirdparty.h"
25 #include "hextra.h"
26 #include "rfield.h"
27 #include "iterations.h"
28 #include "trace.h"
29 #include "conv.h"
30 #include "timesc.h"
31 #include "dense.h"
32 #include "mole.h"
33 #include "thermal.h"
34 #include "pressure.h"
35 #include "phycon.h"
36 #include "wind.h"
37 #include "hmi.h"
38 #include "dynamics.h"
40 
41 /*
42  * >>chng 01 mar 16, incorporate advection within dynamical solutions
43  * this file contains the routines that incorporate effects of dynamics and advection upon the
44  * thermal and ionization solutions.
45  *
46  * This code was originally developed in March 2001 during
47  * a visit to UNAM Morelia, in collaboration with Robin Williams, Jane Arthur, and Will Henney.
48  * Development was continued by email, and in a meeting in July/Aug 2001 at U Cardiff
49  * Further development June 2002, UNAM Morelia (Cloudy and the Duendes Verdes)
50  */
51 
52 /* turn this on for dynamics printout */
53 #define DIAG_PRINT false
54 #define MAINPRINT false
55 
56 /* the interpolated upstream densities of all ionization stages,
57  * [element][ion] */
58 static double **UpstreamIon;
59 /* total abundance of each element per hydrogen */
60 static double *UpstreamElem;
61 
62 /* hydrogen molecules */
63 static double *Upstream_H2_molec;
64 
65 /* CO molecules */
66 static double *Upstream_CO_molec;
67 
68 /* space to save continuum when time command is used
69 static realnum *dyna_flux_save;*/
70 
71 /* initial timestep (seconds) read off time command,
72  * each iteration accounts for this much time */
73 static double timestep_init,
74  timestep,
77 
78 /* array of times and continuum values we will interpolate upon */
79 static double *time_elapsed_time ,
81  *time_dt,
85 #define NTIME 200
86 
87 /* number of time steps actually read in */
88 static long int nTime_flux;
89 
90 /* routine called at end of iteration to determine new step sizes */
91 STATIC void DynaNewStep(void);
92 
93 /* routine called at end of iteration to save values in previous iteration */
94 STATIC void DynaSaveLast(void);
95 
96 /* routine called to determine mass flux at given distance */
97 /* static realnum DynaFlux(double depth); */
98 
99 /* lookahead distance, as derived in DynaEndIter */
100 static double Dyn_dr;
101 
102 /* advected work */
103 static double AdvecSpecificEnthalpy;
104 
105 /* the upstream hydrogen atoms per unit vol*/
106 static double Upstream_hden;
107 
108 /* HI ionization structure from previous iteration */
109 static realnum *Old_histr/*[NZLIM]*/ ,
110  /* Lyman continuum optical depth from previous iteration */
111  *Old_xLyman_depth/*[NZLIM]*/,
112  /* depth of position in previous iteration */
113  *Old_depth/*[NZLIM]*/,
114  /* old n_p density from previous iteration */
115  *Old_hiistr/*[NZLIM]*/,
116  /* old pressure from previous iteration */
117  *Old_pressure/*[NZLIM]*/,
118  /* H density - particles per unit vol */
119  *Old_hden/*[NZLIM]*/ ,
120  /* density - total grams per unit vol */
121  *Old_DenMass/*[NZLIM]*/ ,
122  /* sum of enthalpy and kinetic energy per gram */
123  *EnthalpyDensity/*[NZLIM]*/,
124  /* old electron density from previous iteration */
125  *Old_ednstr/*[NZLIM]*/,
126  /* sum of enthalpy and kinetic energy per gram */
127  *Old_EnthalpyDensity/*[NZLIM]*/;
128 
131 
132 /* the ionization fractions from the previous iteration */
134 
135 /* the gas phase abundances from the previous iteration */
137 
138 /* the number of zones that were saved in the previous iteration */
139 static long int nOld_zone;
140 
141 /* the number of zones that were saved in the previous iteration */
143 
144 /*timestep_next - find next time step in time dependent case */
145 STATIC double timestep_next( void )
146 {
147  static double te_old=-1;
148  double timestep_Hp_temp , timestep_return;
149 
150  DEBUG_ENTRY( "timestep_next()" );
151 
152  timestep_return = timestep;
153 
154  if( dynamics.lgRecom )
155  {
156  double te_new;
157  if( cdTemp(
158  /* four char string, null terminzed, giving the element name */
159  "hydr",
160  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
161  2 ,
162  /* will be temperature */
163  &te_new,
164  /* how to weight the average, must be "VOLUME" or "RADIUS" */
165  "VOLUME" ) )
166  TotalInsanity();
167 
168  if( te_old>0 )
169  {
170  double dTdStep = fabs(te_new-te_old)/te_new;
171  /* this is the change factor to get 0.1 frac change in mean temp */
172  double dT_factor = 0.04/SDIV(dTdStep);
173  dT_factor = MIN2( 2.0 , dT_factor );
174  dT_factor = MAX2( 0.01 , dT_factor );
175  timestep_Hp_temp = timestep * dT_factor;
176  }
177  else
178  {
179  timestep_Hp_temp = -1.;
180  }
181  te_old = te_new;
182  if( timestep_Hp_temp > 0. )
183  timestep_return = timestep_Hp_temp;
184  }
185  else
186  {
187  timestep_return = timestep_init;
188  }
189  fprintf(ioQQQ,"DEBUG timestep_next returns %.3e, old temp %.2e\n" , timestep_return, te_old );
190 
191  return( timestep_return );
192 }
193 
194 /* ============================================================================== */
195 /* DynaPresChngFactor, called from PressureChange to evaluate factor needed
196  * to find new density needed for
197  * current conditions and wind solution, returns ratio of new to old density,
198  * called when wind velocity is negative */
199 
200 /* object is to get the local ram pressure
201  * RamNeeded = pressure.PresTotlInit + pressure.PresGasCurr + pressure.PresInteg;
202  * to equal the sum of the inital pressur at the illuminated face, the local gas pressure,
203  * and the integrate radiative acceleration from the incident continuum
204  *
205  * the local gas pressure is linear in density if the temperature is constant,
206  *
207  * the local ram pressure is inversely linear in density because of the relationship
208  * between velocity and density introduced by the mass flux conservation
209  */
210 #define SUBSONIC 1
211 #define SUPERSONIC 2
212 /*#define FREE 3*/
213 #define STRONGD 4
214 #define ORIGINAL 5
215 #define SHOCK 6
216 #define ANTISHOCK 7
217 #define ANTISHOCK2 8
218 
219 double DynaPresChngFactor(void)
220 {
221  double factor,
222  er,
223  width;
224  static double dp = -1.,
225  dpp = -1.,
226  erp = -1.,
227  erpp = -1.;
230  static int lastzone = -1,
231  zonepresmode,
232  globalpresmode;
233  int ipPRE;
234 
235  DEBUG_ENTRY( "DynaPresChngFactor()" );
236 
237  /* update the current pressure and radiative acceleration */
238  PresTotCurrent();
239 
240  /* update the current desired pressure */
243  if( trace.lgTrace )
244  {
245  fprintf( ioQQQ,
246  " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n",
249  }
250 
251  if( DIAG_PRINT)
252  fprintf(stdout,"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
255  /*fprintf(ioQQQ,"DEBUG\t%.2f\thden\t%.4e\tPtot\t%.4e\tPgas\t%.4e\n",
256  fnzone, dense.gas_phase[ipHYDROGEN],pressure.PresTotlCurr,pressure.PresGasCurr );*/
257 
258  /* dynamics.lgSetPresMode is flag to indicate sane value of dynamics.chPresMode.
259  * If set true with SET DYNAMICS PRESSURE MODE
260  * then do not override that choice */
261  if( !dynamics.lgSetPresMode )
262  {
263  /* above set true if pressure mode was set with
264  * SET DYNAMICS PRESSURE MODE - if we got here
265  * it was not set, and must make a guess */
267  strcpy( dynamics.chPresMode , "supersonic" );
268  else
269  strcpy( dynamics.chPresMode , "subsonic" );
270  /* clear the flag - pressure mode has been set */
271  dynamics.lgSetPresMode = true;
272  }
273 
274  /* if globally looking for transonic solution, then locally sometimes
275  * supersonic, sometimes subsonic - this branch sets global flag,
276  * which can also be set with SET DYNAMICS PRESSURE MODE.
277  * Under default conditions, ChPresMode was set in previous branch
278  * to sub or super sonic depending on current velocity on first time*/
279  if( strcmp( dynamics.chPresMode , "original" ) == 0 )
280  {
281  globalpresmode = ORIGINAL;
283  }
284  else if( strcmp( dynamics.chPresMode , "subsonic" ) == 0 )
285  {
286  globalpresmode = SUBSONIC;
288  }
289  /* supersonic */
290  else if( strcmp( dynamics.chPresMode , "supersonic" ) == 0 )
291  {
292  globalpresmode = SUPERSONIC;
294  }
295  /* strong d */
296  else if( strcmp( dynamics.chPresMode , "strongd" ) == 0 )
297  {
298  globalpresmode = STRONGD;
300  }
301  else if( strcmp( dynamics.chPresMode , "shock" ) == 0 )
302  {
303  globalpresmode = SHOCK;
305  }
306  else if( strcmp( dynamics.chPresMode , "antishock-by-mach" ) == 0 )
307  {
308  globalpresmode = ANTISHOCK2;
310  }
311  else if( strcmp( dynamics.chPresMode , "antishock" ) == 0 )
312  {
313  globalpresmode = ANTISHOCK;
315  }
316 
317  /* this sets pressure mode for the current zone based on global mode
318  * and local conditions */
319  if(globalpresmode == ORIGINAL)
320  {
321  /* in this mode use comparison between ram and gas pressure to determine
322  * whether sub or super sonic */
324  {
325  zonepresmode = SUBSONIC;
326  }
327  else
328  {
329  zonepresmode = SUPERSONIC;
330  }
331  }
332  else if(globalpresmode == STRONGD)
333  {
334  if(nzone <= 1)
335  zonepresmode = SUPERSONIC;
336  }
337  else if(globalpresmode == SUBSONIC)
338  {
339  zonepresmode = SUBSONIC;
340  }
341  else if(globalpresmode == SUPERSONIC)
342  {
343  zonepresmode = SUPERSONIC;
344  }
345  else if(globalpresmode == SHOCK)
346  {
348  {
349  zonepresmode = SUBSONIC;
350  }
351  else
352  {
353  zonepresmode = SUPERSONIC;
354  }
355  }
356  else if(globalpresmode == ANTISHOCK)
357  {
359  {
360  zonepresmode = SUPERSONIC;
361  }
362  else
363  {
364  zonepresmode = SUBSONIC;
365  }
366  }
367  else if(globalpresmode == ANTISHOCK2)
368  {
369  /* WJH 19 May 2004: This version of the antishock mode will
370  insert the antishock at the point where the isothermal Mach
371  number falls below a certain value, dynamics.ShockMach */
373  {
374  zonepresmode = SUPERSONIC;
375  }
376  else
377  {
378  zonepresmode = SUBSONIC;
379  }
380  }
381  else
382  {
383  printf("Need to set global pressure mode\n");
384  cdEXIT( EXIT_FAILURE );
385  }
386 
388  /* fprintf(ioQQQ,"Ds %ld: %ld; %g %g %g %g %g %g %g %g\n",iteration,nzone,dense.gas_phase[ipHYDROGEN],er,pressure.PresTotlCorrect,pressure.PresTotlCurr,phycon.te,thermal.ctot,thermal.htot,phycon.EnthalpyDensity); */
389  /* fprintf(ioQQQ,"Ds %ld: %ld; %g %g %g %g\n",iteration,nzone,dense.gas_phase[ipHYDROGEN],er,pressure.PresTotlCurr,phycon.te); */
390  if(globalpresmode == ORIGINAL || lastzone != nzone || fabs(er-erp) < SMALLFLOAT)
391  {
392  /* First time in this zone, or when last step made no change,
393  * take step hopefully in the right direction...
394  * ...or at least somewhere */
395  if( zonepresmode == SUBSONIC )
396  {
397  /* when gas pressure dominated need to increase density to increase pressure */
399  ipPRE = 0;
400  }
401  else
402  {
403  /* when ram pressure dominated need to decrease density to increase pressure */
405  ipPRE = 1;
406  }
407  if(fabs(factor-1.) > 0.01)
408  {
409  factor = 1.+sign(0.01,factor-1.);
410  }
411  erp = er;
412  dp = dense.gas_phase[ipHYDROGEN];
413  erpp = -1.;
414  dpp = -1.;
415  }
416  else
417  {
418 #if 0
419  printf("Ds: %d; %g %g %g; %g %g %g tot %g\n",nzone,dense.gas_phase[ipHYDROGEN],dp,dpp,er,erp,erpp,
421 #endif
422  if(1 || dpp < 0. || fabs((dense.gas_phase[ipHYDROGEN]-dp)*(dp-dpp)*(dpp-dense.gas_phase[ipHYDROGEN])) < SMALLFLOAT)
423  {
424  /* second iteration on this zone, do linerar fit to previous Pres - rho curve
425  * Linear approximation to guess root with two independent samples */
426  factor = (dense.gas_phase[ipHYDROGEN]*erp-dp*er)/(erp-er)/dense.gas_phase[ipHYDROGEN];
427  /* Limit step length to `reasonable' extrapolation */
428  width = fabs(1.-dp/dense.gas_phase[ipHYDROGEN]);
429  if(width > 1e-2)
430  width = 1e-2;
431 
432  /* Subsonic case: pressure ^ with density ^ => increase density further */
433  /* Super " case: pressure ^ with density v => decrease density further */
434 
435  /* printf("Presmode %d flag %g factor %g\n",zonepresmode,(er-erp)*(dense.gas_phase[ipHYDROGEN]-dp),factor); */
436  /* WJH 21 May 2004: I think that this part is to force the solution
437  * towards the required branch when it appears to have the "wrong" value
438  * of dP/drho */
439  if(zonepresmode == SUBSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) < 0)
440  {
441  factor = 1+3*width;
442  }
443  else if(zonepresmode == SUPERSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) > 0)
444  {
445  factor = 1-3*width;
446  }
447 
448  if(fabs(factor-1.) > 3*width)
449  {
450  factor = 1.+sign(3*width,factor-1.);
451  }
452  ipPRE = 2;
453  if(fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT)
454  {
455  dpp = dp;
456  erpp = erp;
457  dp = dense.gas_phase[ipHYDROGEN];
458  erp = er;
459  }
460  }
461  else
462  {
463  /* 3rd or more solutions in this zone so have extensive data
464  * on pressure - density relation
465  * Use quadratic fit to last three errors to estimate optimum */
466  double a, b, c, q, dmin, emin, dsol, f1 , f2;
467  int iroot;
468  a = er/(dense.gas_phase[ipHYDROGEN]-dp)/(dense.gas_phase[ipHYDROGEN]-dpp) +
469  erp/(dp-dense.gas_phase[ipHYDROGEN])/(dp-dpp)+
470  erpp/(dpp-dense.gas_phase[ipHYDROGEN])/(dpp-dp);
471  b = (erp-erpp)/(dp-dpp) - a * (dp+dpp);
472  c = erp - dp*(a*dp+b);
473  dmin = (-0.5*b/a);
474  emin = (a*dmin+b)*dmin+c;
475  if(a < 0)
476  {
477  a = -a;
478  b = -b;
479  c = -c;
480  }
481 #if 0
482  fprintf(ioQQQ,"Check 1: %g %g\n",er,(a*dense.gas_phase[ipHYDROGEN]+b)*dense.gas_phase[ipHYDROGEN]+c);
483  fprintf(ioQQQ,"Check 2: %g %g\n",erp,(a*dp+b)*dp+c);
484  fprintf(ioQQQ,"Check 3: %g %g\n",erpp,(a*dpp+b)*dpp+c);
485 #endif
486  q = b*b-4*a*c;
487  if(q < 0)
488  {
489  /* Imaginary root, search for local minimum */
490  /* printf("No root at %d (%g cf %g) => %g\n",nzone,q,b*b,dmin); */
491  factor = dmin/dense.gas_phase[ipHYDROGEN];
492 
496  if(globalpresmode == STRONGD && -q > 1e-3*b*b)
497  {
498  zonepresmode = SUBSONIC;
499  }
500  }
501  else
502  {
503  /* Look for nearest root */
504  /* if(zonepresmode == SUPERSONIC || (zonepresmode != SUBSONIC && (dense.gas_phase[ipHYDROGEN]-dmin) < 0)) */
505 
506  if(zonepresmode == SUPERSONIC)
507  iroot = 1;
508  else
509  iroot = 0;
510 
511  if(emin > 0)
512  iroot = ! iroot;
513 
514  if(iroot)
515  {
516  /* Look for smaller root */
517  if(b > 0)
518  {
519  dsol = -(b+sqrt(q))/(2*a);
520  }
521  else
522  {
523  dsol = 2*c/(-b+sqrt(q));
524  }
525  }
526  else
527  {
528  if(b < 0)
529  {
530  dsol = (-b+sqrt(q))/(2*a);
531  }
532  else
533  {
534  dsol = -2*c/(b+sqrt(q));
535  }
536  }
537  factor = dsol/dense.gas_phase[ipHYDROGEN];
538 #if 0
539  fprintf(ioQQQ,"Check 4: %g %g %d %g %g\n",dsol,(a*dsol+b)*dsol+c,iroot,dmin,emin);
540 #endif
541  }
542  /* Limit step length */
543  f1 = fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]);
544  f2 = fabs(1.- dp/dense.gas_phase[ipHYDROGEN]);
545  /*width = MAX2(fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]),fabs(1.-dp/dense.gas_phase[ipHYDROGEN]));*/
546  width = MAX2(f1,f2);
547  /* width = MAX2(width,1e-2); */
548  if(fabs(factor-1.) > 3*width)
549  {
550  factor = 1.+sign(3*width,factor-1.);
551  }
552  ipPRE = 3;
553  if(fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT)
554  {
555  dpp = dp;
556  erpp = erp;
557  dp = dense.gas_phase[ipHYDROGEN];
558  erp = er;
559  }
560  }
561  }
562 
563 #if 0
564  printf("Out: %d; %g; %g %g; %g %g\n",nzone,factor*dense.gas_phase[ipHYDROGEN],dp,dpp,erp,erpp);
565 #endif
566  lastzone = nzone;
567 
568  if( DIAG_PRINT )
569  fprintf(ioQQQ,"windv %li r %g v %g f %g\n",
571 
572  /* convergence trace at this level */
573  if( trace.nTrConvg>=1 )
574  {
575  fprintf( ioQQQ,
576  " DynaPresChngFactor mode %s scale factor %.3f vel %.3e\n",
577  dynamics.chPresMode , factor , wind.windv );
578  }
579 
580  {
581  /*@-redef@*/
582  enum{DEBUG_LOC=false};
583  /*@+redef@*/
584  if( DEBUG_LOC )
585  {
586  char chPRE[][4] = {"gas" , "ram", "sec", "par" };
587  fprintf(ioQQQ,
588  "pre %s\tfac\t%.5f\n",
589  chPRE[ipPRE],
590  factor
591  );
592  }
593  }
594 
595  /* identically zero velocities cannot occur */
596  ASSERT( wind.windv != 0. || (wind.windv == 0. && dynamics.lgStatic) );
597 
598  return factor;
599 }
600 
601 /* ============================================================================== */
602 /* DynaIonize, called from ConvBase to evaluate advective terms for current conditions,
603  * calculates terms to add to ionization balance equation */
604 void DynaIonize(void)
605 {
606  long int nelem,
607  ion,
608  mol ,
609  i;
610 
611  DEBUG_ENTRY( "DynaIonize()" );
612 
613  /* the time (s) needed for gas to move Dyn_dr */
614  /* >>chng 02 dec 11 rjrw increase timestep when beyond end of previous zone, to allow -> eqm */
615  if( !dynamics.lgStatic )
616  {
617  /* in time dependent model timestep only changes once at end of iteration
618  * and is constant across a model */
619  /* usual case - full dynamics */
621  }
622  /* printf("%d %g %g %g %g\n",nzone,radius.depth,Dyn_dr,radius.depth-Old_depth[nOld_zone-1],timestep); */
623 
625  if( nzone > 0 )
627 
628  /* do nothing on first iteration or when looking beyond previous iteration */
629  /* >>chng 05 jan 27, from hardwired "iteration < 2" to more general case,
630  * this is set with SET DYNAMICS RELAX command with the default of 2 */
631  /*if( iteration < 2 || radius.depth > dynamics.oldFullDepth)*/
633  {
634  /* first iteration, return zero */
635  dynamics.Cool = 0.;
636  dynamics.Heat = 0.;
637  dynamics.dCooldT = 0.;
638  dynamics.dHeatdT = 0.;
639 
640  dynamics.Rate = 0.;
641 
642  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
643  {
644  for( ion=0; ion<nelem+2; ++ion )
645  {
646  dynamics.Source[nelem][ion] = 0.;
647  }
648  }
649 
650  for(mol=0;mol<N_H_MOLEC;mol++)
651  {
652  dynamics.H2_molec[mol] = 0.;
653  }
654  for(mol=0;mol<mole.num_comole_calc;mol++)
655  {
656  dynamics.CO_molec[mol] = 0.;
657  }
658  return;
659  }
660 
661  if( MAINPRINT )
662  {
663  fprintf(ioQQQ,"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
664  nzone,
667  5./2.*pressure.PresGasCurr
668  );
669  }
670 
671  /* net cooling due to advection */
672  /* >>chng 01 aug 01, removed hden from dilution, put back into here */
673  /* >>chng 01 sep 15, do not update this variable the very first time this routine
674  * is called at the new zone. */
677  /* temp deriv of cooling minus heating */
680 
681 #if 0
682  /* printf("DynaCool %g %g %g\n", dynamics.Cool, phycon.EnthalpyDensity/timestep,AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep); */
683  if(dynamics.Cool > 0) {
684  dynamics.Heat = 0.;
685  /* temp deriv of cooling minus heating */
688  } else {
690  dynamics.Cool = 0.;
691  /* temp deriv of cooling minus heating */
694  }
695 #endif
696 
697 # if 0
698  if( MAINPRINT || nzone>17 && iteration == 10)
699  {
700  fprintf(ioQQQ,
701  "dynamics cool-heat\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t %.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
702  nzone,
703  phycon.te,
704  dynamics.CoolHeat,
705  thermal.htot,
710  dense.gas_phase[ipHYDROGEN],
711  timestep);
712  }
713 # endif
714 
715  /* second or greater iteration, have advective terms */
716  /* this will evaluate advective terms for current physical conditions */
717 
718  /* the rate (s^-1) atoms drift in from upstream, a source term for the ground */
719 
720  /* dynamics.Hatom/timestep is the source (cm^-3 s^-1) of neutrals,
721  normalized to (s^-1) by the next higher ionization state as for all
722  other recombination terms */
723 
724  /* dynamics.xIonDense[ipHYDROGEN][1]/timestep is the sink (cm^-3 s^-1) of
725  ions, normalized to (s^-1) by the ionization state as for all other
726  ionization terms */
727 
728  dynamics.Rate = 1./timestep;
729 
730  for(mol=0;mol<N_H_MOLEC;mol++)
731  {
733  }
734 
735  for(mol=0;mol<mole.num_comole_calc;mol++)
736  {
738  }
739 
740  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
741  {
742  if( dense.lgElmtOn[nelem] )
743  {
744  /* check that the total number of each element is conserved along the flow */
745  if(fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]>=1e-3)
746  {
747  /* sum of all H in standard H chem */
748  realnum sumh = 0.;
749  for(i=0;i<N_H_MOLEC;i++)
750  {
751  sumh += hmi.Hmolec[i]*hmi.nProton[i];
752  }
753  fprintf(ioQQQ,
754  "PROBLEM conservation error: zn %li elem %li upstream %.8e abund %.8e (up-ab)/up %.2e f1(H n CO) %.2e f2(H n CO) %.2e\n",
755  nzone ,
756  nelem,
758  dense.gas_phase[nelem] ,
759  (UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN]-dense.gas_phase[nelem]) /
760  (UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN]) ,
761  (dense.gas_phase[ipHYDROGEN]-sumh) / dense.gas_phase[ipHYDROGEN] ,
762  dense.H_sum_in_CO / dense.gas_phase[ipHYDROGEN] );
763  }
764  /* ASSERT( fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]<1e-3); */
765  for( ion=0; ion<dense.IonLow[nelem]; ++ion )
766  {
767  dynamics.Source[nelem][ion] = 0.;
768  }
769  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
770  {
771  /* Normalize to next higher state in current zone, except at first iteration
772  where upstream version may be a better estimate (required for
773  convergence in the small timestep limit) */
774 
775  dynamics.Source[nelem][ion] =
776  /* UpstreamIon is ion number per unit hydrogen because dilution is 1/hden
777  * this forms the ratio of upstream atom over current ion, per timestep,
778  * so Recomb has units s-1 */
779  UpstreamIon[nelem][ion] * dense.gas_phase[ipHYDROGEN] / timestep;
780 
781  }
782  for( ion=dense.IonHigh[nelem]+1;ion<nelem+2; ++ion )
783  {
784  dynamics.Source[nelem][ion] = 0.;
785  }
786  }
787  }
788 # if 0
789  fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
790  nzone,
791  dynamics.Rate,
793  dynamics.Rate,
794  dynamics.Source[ipCARBON][3]);
795 # endif
796 #if 0
797  nelem = ipCARBON;
798  ion = 3;
799  /*if( nzone > 160 && iteration > 1 )*/
800  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",
801  nzone,
802  ipUpstream,
803  radius.depth ,
805  dense.xIonDense[nelem][ion],
806  UpstreamIon[nelem][ion]* dense.gas_phase[ipHYDROGEN],
807  Old_xIonDense[ipUpstream][nelem][ion] ,
808  dense.xIonDense[nelem][ion+1],
809  UpstreamIon[nelem][ion+1]* dense.gas_phase[ipHYDROGEN],
810  Old_xIonDense[ipUpstream][nelem][ion+1] ,
811  timestep,
812  dynamics.Source[nelem][ion]
813  );
814 #endif
815  if( MAINPRINT )
816  {
817  fprintf(ioQQQ," DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
818  nzone,dynamics.Rate , dynamics.Source[0][0] );
819  }
820  return;
821 }
822 
823 /* ============================================================================== */
824 /* DynaStartZone called at start of zone calculation when advection is turned on */
825 void DynaStartZone(void)
826 {
827  /* this routine is called at the start of a zone calculation, by ZoneStart:
828  *
829  * it sets deduced variables to zero if first iteration,
830  *
831  * if upstream depth is is outside the computed structure on previous iteration,
832  * return value at shielded face
833  *
834  * Also calculates discretization_error, an estimate of the accuracy of the source terms.
835  *
836  * */
837 
838  /* this is index of upstream cell in structure stored from previous iteration */
839  double upstream, dilution, dilutionleft, dilutionright, frac_next;
840 
841  /* Properties for cell half as far upstream, used to converge timestep */
842  double hupstream, hnextfrac=-BIGFLOAT, hion, hmol;
843 
844  /* Properties for cell at present depth, used to converge timestep */
845  double ynextfrac=-BIGFLOAT, yion, ymol;
846 
847  long int nelem , ion, mol;
848 
849  DEBUG_ENTRY( "DynaStartZone()" );
850 
851  /* do nothing on first iteration */
852  if( iteration < 2 )
853  {
854  Upstream_hden = 0.;
856  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
857  {
858  for( ion=0; ion<nelem+2; ++ion )
859  {
860  UpstreamIon[nelem][ion] = 0.;
861  }
862  }
863  /* hydrogen molecules */
864  for(mol=0; mol<N_H_MOLEC;mol++)
865  {
866  Upstream_H2_molec[mol] = 0;
867  }
868  for(mol=0; mol<mole.num_comole_calc;mol++)
869  {
870  Upstream_CO_molec[mol] = 0;
871  }
872 
873  ipUpstream = 0;
874  iphUpstream = 0;
875  ipyUpstream = 0;
876  return;
877  }
878 
879  /* radius.depth is distance from illuminated face of cloud to outer edge of
880  * current zone, which has thickness radius.drad */
881 
882  /* find where the down stream point is, in previous iteration,
883  * NB neg sign since velocity is negative, we are looking for point
884  * where current gas cell was, in previous iteration */
885 
886  /* true, both advection and wind solution */
887  /* don't interpolate to the illuminated side of the first cell */
888  upstream = MAX2(Old_depth[0] , radius.depth + Dyn_dr);
889  hupstream = 0.5*(radius.depth + upstream);
890 
891  /* now locate upstream point in previous stored structure,
892  * will be at the same point or higher than we found previously */
893  while( (Old_depth[ipUpstream+1] < upstream ) &&
894  ( ipUpstream < nOld_zone-1 ) )
895  {
896  ++ipUpstream;
897  }
898  ASSERT( ipUpstream <= nOld_zone-1 );
899 
900  /* above loop will give ipUpstream == nOld_zone-1 if computed structure has been overrun */
901 
903  {
904  /* we have not overrun radius scale of previous iteration */
905  frac_next = ( upstream - Old_depth[ipUpstream])/
906  (Old_depth[ipUpstream+1] - Old_depth[ipUpstream]);
908  (Old_hden[ipUpstream+1] - Old_hden[ipUpstream])*
909  frac_next;
910  dilutionleft = 1./Old_hden[ipUpstream];
911  dilutionright = 1./Old_hden[ipUpstream+1];
912 
913  /* fractional changes in density from passive advection */
914  /* >>chng 01 May 02 rjrw: use hden for dilution */
915  /* >>chng 01 aug 01, remove hden here, put back into vars when used in DynaIonize */
916  dilution = 1./Upstream_hden;
917 
918  /* the advected enthalphy */
920  (Old_EnthalpyDensity[ipUpstream+1]*dilutionright - Old_EnthalpyDensity[ipUpstream]*dilutionleft)*
921  frac_next);
922 
923  ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] );
924  if( (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) *
925  (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) < -SMALLFLOAT )
926  {
927  fprintf(ioQQQ,"PROBLEM advected enthalpy density is not conserved %e\t%e\t%e\t%e\n",
928  (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) *
929  (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) ,
930  Old_EnthalpyDensity[ipUpstream]*dilutionleft,
932  Old_EnthalpyDensity[ipUpstream+1]*dilutionright);
933  }
934 
935  ASSERT( (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) *
936  (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) > -SMALLFLOAT );
937  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
938  {
939  UpstreamElem[nelem] = 0.;
940  for( ion=0; ion<nelem+2; ++ion )
941  {
942  /* Robin - I made several changes like the following - it seems easier to
943  * bring out the division by the old hydrogen density rather than putting in
944  * dilution and then looking for how diluation is defined. I think the code
945  * is equivalent */
946  /* UpstreamIon is ion number per unit hydrogen, both at the upstream position */
947  UpstreamIon[nelem][ion] =
949  (Old_xIonDense[ipUpstream+1][nelem][ion]/Old_hden[ipUpstream+1] -
951  frac_next;
952 
953  UpstreamElem[nelem] += UpstreamIon[nelem][ion];
954  }
955  }
956 
957  for(mol=0;mol<N_H_MOLEC;mol++)
958  {
960  (Old_H2_molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] -
962  frac_next;
963  /* test on mol > 1, first two elements are H0 and H+, which are alread
964  * counted in upstreamElem */
965  if(mol > 1)
967  }
968  /* loop is only up to mole.num_heavy_molec, not mole.num_comole_calc, since do not want
969  * to consider atoms and ions that have also been done */
970  for(mol=0;mol<mole.num_comole_calc;mol++)
971  {
973  (Old_CO_molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] -
975  frac_next;
976 
977  if(COmole[mol]->n_nuclei > 1)
978  {
979  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
980  {
981  if( mole.lgElem_in_chemistry[nelem] )
982  {
983  UpstreamElem[nelem] +=
984  COmole[mol]->nElem[nelem]*Upstream_CO_molec[mol];
985  }
986  }
987  }
988  }
989  }
990  else
991  {
992  /* SPECIAL CASE - we have overrun the previous iteration's radius */
994  /* fractional changes in density from passive advection */
995  dilution = 1./Upstream_hden;
996  /* AdvecSpecificEnthalpy enters as heat term */
998  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
999  {
1000  UpstreamElem[nelem] = 0.;
1001  for( ion=0; ion<nelem+2; ++ion )
1002  {
1003  /* UpstreamIon is ion number per unit hydrogen */
1004  UpstreamIon[nelem][ion] =
1006  UpstreamElem[nelem] += UpstreamIon[nelem][ion];
1007  }
1008  }
1009 
1010  for(mol=0;mol<N_H_MOLEC;mol++)
1011  {
1013  if(mol > 1)
1015  }
1016  for(mol=0;mol<mole.num_comole_calc;mol++)
1017  {
1019  if(COmole[mol]->n_nuclei > 1)
1020  {
1021  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1022  {
1023  if( mole.lgElem_in_chemistry[nelem] )
1024  {
1025  UpstreamElem[nelem] +=
1026  Upstream_CO_molec[mol]*COmole[mol]->nElem[nelem];
1027  }
1028  }
1029  }
1030  }
1031  }
1032 
1033  /* Repeat enough of the above for half-step and no-step to judge convergence:
1034  the result of this code is the increment of discretization_error */
1035 
1036  while( (Old_depth[iphUpstream+1] < hupstream ) &&
1037  ( iphUpstream < nOld_zone-1 ) )
1038  {
1039  ++iphUpstream;
1040  }
1041  ASSERT( iphUpstream <= nOld_zone-1 );
1042 
1043  while( (Old_depth[ipyUpstream+1] < radius.depth ) &&
1044  ( ipyUpstream < nOld_zone-1 ) )
1045  {
1046  ++ipyUpstream;
1047  }
1048  ASSERT( ipyUpstream <= nOld_zone-1 );
1049 
1051 
1053  {
1054  hnextfrac = ( hupstream - Old_depth[iphUpstream])/
1056  /* >>chng 02 nov 11, hhden never appears on rhs of = */
1057  /*hhden = Old_hden[iphUpstream] +
1058  (Old_hden[iphUpstream+1] - Old_hden[iphUpstream])*
1059  hnextfrac;*/
1060  }
1061  else
1062  {
1063  /* >>chng 06 mar 17, set this to zero */
1064  hnextfrac = 0.;
1065  /* >>chng 02 nov 11, hhden never appears on rhs of = */
1066  /*hhden = Old_hden[iphUpstream];*/
1067  }
1069  {
1070  ynextfrac = ( radius.depth - Old_depth[ipyUpstream])/
1072  /* >>chng 02 nov 11, yhden never appears on rhsof = */
1073  /*yhden = Old_hden[ipyUpstream] +
1074  (Old_hden[ipyUpstream+1] - Old_hden[ipyUpstream])*
1075  ynextfrac;*/
1076  }
1077  else
1078  {
1079  /* >>chng 06 mar 17, set this to zero */
1080  ynextfrac = 0.;
1081  /* >>chng 02 nov 11, yhden never appears on rhsof = */
1082  /*yhden = Old_hden[ipyUpstream];*/
1083  }
1084 
1085  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1086  {
1087  for( ion=0; ion<nelem+2; ++ion )
1088  {
1089  double f1;
1091  {
1092  yion =
1094  (Old_xIonDense[ipyUpstream+1][nelem][ion]/Old_hden[ipyUpstream+1] -
1096  ynextfrac;
1097  }
1098  else
1099  {
1100  yion = Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream];
1101  }
1103  {
1104  hion =
1106  (Old_xIonDense[iphUpstream+1][nelem][ion]/Old_hden[iphUpstream+1] -
1108  hnextfrac;
1109  }
1110  else
1111  {
1112  hion = Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream];
1113  }
1114 
1115  /* the proposed thickness of the next zone, there will be a scale factor, something
1116  * like 1/500, that will be applied when this is used in nextdr */
1117  f1 = fabs(yion - UpstreamIon[nelem][ion] );
1118  f1 = SDIV( f1 );
1120  /* don't pay attention to species with abundance relative to H less tghan 1e-6 */
1121  MAX2(yion + UpstreamIon[nelem][ion],1e-6 ) / f1);
1122  /*MAX2(fabs(yion - UpstreamIon[nelem][ion] ),SMALLFLOAT));*/
1123  /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */
1124 
1125  /* Must be consistent with convergence_error below */
1126  /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */
1127  /* >>chngf 02 aug 01, multiply by cell width */
1128  /* this error is error due to the advection length not being zero - a finite
1129  * advection length. no need to bring convergence error to below
1130  * discretization error. when convergece error is lower than a fraction of this
1131  * errror we reduce the advection length. */
1132  dynamics.discretization_error += POW2(yion-2.*hion+UpstreamIon[nelem][ion]);
1133  dynamics.error_scale2 += POW2(UpstreamIon[nelem][ion]-yion);
1134  }
1135  }
1136 
1137  for(mol=0; mol < N_H_MOLEC; mol++)
1138  {
1139  double f1;
1141  {
1142  ymol =
1146  ynextfrac;
1147  }
1148  else
1149  {
1151  }
1153  {
1154  hmol =
1158  hnextfrac;
1159  }
1160  else
1161  {
1163  }
1164 
1165  /* the proposed thickness of the next zone, there will be a scale factor, something
1166  * like 1/500, that will be applied when this is used in nextdr */
1167  f1 = fabs(ymol - Upstream_H2_molec[mol] );
1168  f1 = SDIV( f1 );
1170  /* don't pay attention to species with abundance relative to H less tghan 1e-6 */
1171  MAX2(ymol + Upstream_H2_molec[mol],1e-6 ) / f1 );
1172  /*MAX2(fabs(ymol - Upstream_H2_molec[mol] ),SMALLFLOAT));*/
1173  /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */
1174 
1175  /* Must be consistent with convergence_error below */
1176  /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */
1177  /* >>chngf 02 aug 01, multiply by cell width */
1178  dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_H2_molec[mol]);
1180  }
1181 
1182  for(mol=0; mol < mole.num_comole_calc; mol++)
1183  {
1184  double f1;
1186  {
1187  ymol =
1191  ynextfrac;
1192  }
1193  else
1194  {
1196  }
1198  {
1199  hmol =
1203  hnextfrac;
1204  }
1205  else
1206  {
1208  }
1209 
1210  /* the proposed thickness of the next zone, there will be a scale factor, something
1211  * like 1/500, that will be applied when this is used in nextdr */
1212  f1 = fabs(ymol - Upstream_CO_molec[mol] );
1213  f1 = SDIV( f1 );
1215  /* don't pay attention to species with abundance relative to H less tghan 1e-6 */
1216  MAX2(ymol + Upstream_CO_molec[mol],1e-6 ) / f1 );
1217  /*MAX2(fabs(ymol - Upstream_H2_molec[mol] ),SMALLFLOAT));*/
1218  /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */
1219 
1220  /* Must be consistent with convergence_error below */
1221  /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */
1222  /* >>chngf 02 aug 01, multiply by cell width */
1223  dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_CO_molec[mol]);
1225  }
1226 
1227  if( MAINPRINT )
1228  {
1229  fprintf(ioQQQ," DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
1230  nzone,dynamics.Rate , dynamics.Source[0][0] , dilution*dense.gas_phase[ipHYDROGEN] );
1231  }
1232  return;
1233 }
1234 
1235 /* DynaEndZone called at end of zone calculation when advection is turned on */
1236 void DynaEndZone(void)
1237 {
1238  DEBUG_ENTRY( "DynaEndZone()" );
1239 
1240  /* this routine is called at the end of a zone calculation, by ZoneEnd */
1241 
1243  if(DIAG_PRINT)
1244  fprintf(ioQQQ,"Check dp: %g %g mom %g %g mas %g\n",
1249  DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth));
1250  return;
1251 }
1252 
1253 
1254 /* ============================================================================== */
1255 /* DynaEndIter called at end of iteration when advection is turned on */
1256 void DynaEndIter(void)
1257 {
1258  /* this routine is called by IterRestart at the end of an iteration
1259  * when advection is turned on. currently it only derives a
1260  * timestep by looking at the spatial derivative of
1261  * some stored quantities */
1262  long int i;
1263  static long int nTime_dt_array_element = 0;
1264 
1265  DEBUG_ENTRY( "DynaEndIter()" );
1266 
1267  /* set stopping outer radius to current outer radius once we have let
1268  * solution relax dynamics.n_initial_relax times
1269  * note the off by one confusion - relax is 2 by default,
1270  * want to do two static iterations then start dynamics
1271  * iteration was incremented before this call so iteration == 2 at
1272  * end of first iteration */
1274  {
1275  fprintf(ioQQQ,"DYNAMICS DynaEndIter sets stop radius to %.2e after"
1276  "dynamics.n_initial_relax=%li iterations.\n",
1277  radius.depth,
1279  for( i=iteration-1; i<iterations.iter_malloc; ++i )
1280  /* set stopping radius to current radius, this stops
1281  * dynamical solutions from overrunning the upstream scale
1282  * and extending to very large radius due to unphysical heat
1283  * appearing to advect into region */
1284  radius.router[i] = radius.depth;
1285  }
1286 
1287  DivergePresInteg = 0.;
1288 
1289  /* This routine is only called if advection is turned on at the end of
1290  * each iteration. The test
1291  * is to check whether wind velocity is also set by dynamics code */
1292 
1293  /* !dynamics.lgStatic true - a true dynamical model */
1294  if( !dynamics.lgStatic )
1295  {
1297  {
1298  /* let model settle down for n_initial_relax iterations before we begin
1299  * dynamics.n_initial_relax set with "set dynamics relax"
1300  * command - it gives the first iteration where we do dynamics
1301  * note the off by one confusion - relax is 2 by default,
1302  * want to do two static iterations then start dynamics
1303  * iteration was incremented before this call so iteration == 2
1304  * at end of first iteration */
1305  if( dynamics.AdvecLengthInit> 0. )
1306  {
1308  }
1309  else
1310  {
1311  /* -ve dynamics.adveclengthlimit sets length as fraction of first iter total depth */
1313  }
1314 
1315  if( MAINPRINT )
1316  {
1317  fprintf(ioQQQ," DynaEndIter, dr=%.2e \n",
1318  Dyn_dr );
1319  }
1320  }
1321  else if(iteration > dynamics.n_initial_relax+1 )
1322  {
1323  /* evaluate errors and update Dyn_dr */
1324  DynaNewStep();
1325  }
1326  }
1327  else
1328  {
1329  /* this is time-dependent static model */
1330  static double HeatInitial=-1. , HeatRadiated=-1. ,
1331  DensityInitial=-1;
1332  /* n_initial_relax is number of time-steady models before we start
1333  * to evolve, set with "set dynamics relax" command */
1334  Dyn_dr = 0.;
1335  fprintf(ioQQQ,
1336  "DEBUG times enter timestep %.2e elapsed_time %.2e iteration %li relax %li \n",
1337  timestep ,
1341  {
1342  /* this is into the changing part of the simulation */
1343  long int nTimeVary;
1344 
1345  /* evaluate errors */
1346  DynaNewStep();
1347 
1348  /* this is set true on CORONAL XXX TIME INIT command, says to use constant
1349  * temperature for first n_initial_relax iterations, then let run free */
1351  {
1353  thermal.ConstTemp = 0.;
1354  }
1355 
1356  /* time variable branch, now without dynamics */
1357  /* elapsed time - don't increase dynamics.time_elapsed during first two
1358  * two iterations since this sets static model */
1360  /* timestep_stop is -1 if not set with explicit stop time */
1362  {
1364  }
1365 
1367  {
1368  /* if very early times not specified assume no flux variation yet */
1370  }
1372  {
1373  fprintf( ioQQQ,
1374  " PROBLEM - DynaEnditer - I need the continuum at time %.2e but the table ends at %.2e.\n" ,
1377  cdEXIT(EXIT_FAILURE);
1378  }
1379  else
1380  {
1382  /* the times in seconds */
1384  /* the rfield.time_continuum_scale factors */
1385  time_flux_ratio,
1386  /* the number of rfield.time_continuum_scale factors */
1387  nTime_flux,
1388  /* the desired time */
1390  }
1391 
1392  fprintf(ioQQQ,"DEBUG time dep reset continuum zero timestep %.2e elapsed time %.2e scale %.2e",
1393  timestep ,
1396  if( dynamics.lgRecom )
1397  {
1398  fprintf(ioQQQ," recom");
1399  }
1400  fprintf(ioQQQ,"\n");
1401 
1402  /* make sure that at least one continuum source is variable */
1403  nTimeVary = 0;
1404  for( i=0; i < rfield.nspec; i++ )
1405  {
1406  /* this is set true if particular continuum source can vary with time
1407  * set true if TIME appears on intensity / luminosity command line */
1408  if( rfield.lgTimeVary[i] )
1409  ++nTimeVary;
1410  }
1411 
1413  {
1414  /* vary extra heating */
1416  fprintf(ioQQQ,"DEBUG TurbHeat vary new heat %.2e\n",
1417  hextra.TurbHeat);
1418  }
1419 # if 0
1420  /* >>chng 07 may 23, rm this time dependent logic from here
1421  * time dependent continua implemented when continuum is reset
1422  * in inter_startend */
1423  else if( nTimeVary )
1424  {
1425  for( i=0; i<rfield.nupper; ++i )
1426  {
1427  /* >>chng 06 mar 22, break into constant and time dependent parts */
1428  rfield. flux[i] = rfield.flux_beam_const[i] +
1432  }
1433  }
1434 # endif
1435  else if( !nTimeVary )
1436  {
1437  fprintf(ioQQQ," DISASTER - there were no variable continua "
1438  "or heat sources - put TIME option on at least one "
1439  "luminosity or hextra command.\n");
1440  cdEXIT(EXIT_FAILURE);
1441  }
1442  /* this is heat radiated, after correction for change of H density in constant
1443  * pressure cloud */
1444  HeatRadiated += (thermal.ctot-dynamics.Cool) * timestep *
1445  (DensityInitial / dense.gas_phase[ipHYDROGEN]);
1446  }
1447  else
1448  {
1449  /* this branch, during initial relaxation of solution */
1450  HeatInitial = 1.5 * pressure.PresGasCurr;
1451  HeatRadiated = 0.;
1452  DensityInitial = dense.gas_phase[ipHYDROGEN];
1453  fprintf(ioQQQ,"DEBUG relaxing times requested %li this is step %li\n",
1455  }
1456  fprintf(ioQQQ,"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n",
1457  HeatInitial , HeatRadiated );
1458 
1459  /* at this point dynamics.time_elapsed is the time at the end of the
1460  * previous iteration. We need dt for the next iteration */
1461  if( dynamics.time_elapsed > time_elapsed_time[nTime_dt_array_element] )
1462  {
1463  /* time has advanced to next table point,
1464  * set timestep to specified value */
1465  ++nTime_dt_array_element;
1466  /* this is an assert since we already qualified the array
1467  * element above */
1468  ASSERT( nTime_dt_array_element < nTime_flux );
1469 
1470  /* option to set flag for recombination logic */
1471  if( lgtime_Recom[nTime_dt_array_element] )
1472  {
1473  fprintf(ioQQQ,"DEBUG dynamics turn on recombination logic\n");
1474  dynamics.lgRecom = true;
1475  /* set largest possible zone thickness to value on previous
1476  * iteration when light was on - during recombination conditions
1477  * become much more homogeneous and dr can get too large,
1478  * crashing into H i-front */
1480  }
1481 
1482  if( lgtime_dt_specified )
1483  {
1484  /* this is the new timestep */
1485  fprintf(ioQQQ,"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element,
1486  timestep);
1487  timestep = time_dt[nTime_dt_array_element];
1488  /* option to change time step factor - default is 1.2 set in DynaZero */
1489  if( time_dt_scale_factor[nTime_dt_array_element] > 0. )
1490  timestep_factor = time_dt_scale_factor[nTime_dt_array_element];
1491  }
1492  }
1493  else if( lgtime_dt_specified )
1494  {
1495  /* we are between two points in the table, increase timestep */
1497  fprintf(ioQQQ,"DEBUG lgtimes increment Timeby timestep_factor to %li %.2e\n" ,
1498  nTime_dt_array_element,
1499  timestep );
1500  if(dynamics.time_elapsed+timestep > time_elapsed_time[nTime_dt_array_element] )
1501  {
1502  fprintf(ioQQQ,"DEBUG lgtimes but reset to %.2e\n" ,timestep );
1503  timestep = 1.0001*(time_elapsed_time[nTime_dt_array_element]-dynamics.time_elapsed);
1504  }
1505  }
1506  else
1507  {
1508  /* time not specified in third column, so use initial */
1509  timestep = timestep_next();
1510  }
1511  fprintf(ioQQQ,"DEBUG times exit timestep %.2e elapsed_time %.2e scale %.2e \n",
1512  timestep ,
1515  }
1516 
1517  /* Dyn_dr == 0 is for static time dependent cloud */
1519  Dyn_dr > 0. ||
1520  (Dyn_dr == 0. && wind.windv==0.) );
1521 
1522  /* reset the upstream counters */
1525  dynamics.error_scale2 = 0.;
1526 
1527  /* save results from previous iteration */
1528  DynaSaveLast();
1529  return;
1530 }
1531 
1532 /*DynaNewStep work out convergence errors */
1534 {
1535  long int ilast = 0,
1536  i,
1537  nelem,
1538  ion,
1539  mol;
1540 
1541  double frac_next=-BIGFLOAT,
1542  Oldi_hden,
1543  Oldi_ion,
1544  Oldi_mol;
1545 
1546  DEBUG_ENTRY( "DynaNewStep()" );
1547 
1548  /*n = MIN2(nzone, NZLIM-1);*/
1550  dynamics.error_scale1 = 0.;
1551 
1552  ASSERT( nzone < struc.nzlim);
1553  for(i=0;i<nzone;++i)
1554  {
1555  /* Interpolate for present position in previous solution */
1556  while( (Old_depth[ilast] < struc.depth[i] ) &&
1557  ( ilast < nOld_zone-1 ) )
1558  {
1559  ++ilast;
1560  }
1561  ASSERT( ilast <= nOld_zone-1 );
1562 
1563  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1564  {
1565  frac_next = ( struc.depth[i] - Old_depth[ilast])/
1566  (Old_depth[ilast+1] - Old_depth[ilast]);
1567  Oldi_hden = Old_hden[ilast] +
1568  (Old_hden[ilast+1] - Old_hden[ilast])*
1569  frac_next;
1570  }
1571  else
1572  {
1573  Oldi_hden = Old_hden[ilast];
1574  }
1575  /* Must be consistent with discretization_error above */
1576  /* >>chngf 02 aug 01, multiply by cell width */
1577  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1578  {
1579  for( ion=0; ion<nelem+2; ++ion )
1580  {
1581  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1582  {
1583  Oldi_ion = (Old_xIonDense[ilast][nelem][ion] +
1584  (Old_xIonDense[ilast+1][nelem][ion]-Old_xIonDense[ilast][nelem][ion])*
1585  frac_next);
1586  }
1587  else
1588  {
1589  Oldi_ion = Old_xIonDense[ilast][nelem][ion];
1590  }
1591  dynamics.convergence_error += POW2(Oldi_ion/Oldi_hden-struc.xIonDense[nelem][ion][i]/struc.hden[i]) /* *struc.dr[i] */;
1592 
1593  /* >>chng 02 nov 11, add first error scale estimate from Robin */
1594  dynamics.error_scale1 += POW2(struc.xIonDense[nelem][ion][i]/struc.hden[i]);
1595  }
1596  }
1597  for( mol=0; mol < N_H_MOLEC; mol++)
1598  {
1599  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1600  {
1601  Oldi_mol = (Old_H2_molec[ilast][mol] +
1602  (Old_H2_molec[ilast+1][mol]-Old_H2_molec[ilast][mol])*
1603  frac_next);
1604  }
1605  else
1606  {
1607  Oldi_mol = Old_H2_molec[ilast][mol];
1608  }
1609  dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.H2_molec[mol][i]/struc.hden[i]) /* *struc.dr[i] */;
1610 
1611  /* >>chng 02 nov 11, add first error scale estimate from Robin
1612  * used to normalize the above convergence_error */
1613  dynamics.error_scale1 += POW2(struc.H2_molec[mol][i]/struc.hden[i]);
1614  }
1615  for( mol=0; mol < mole.num_comole_calc; mol++)
1616  {
1617  if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
1618  {
1619  Oldi_mol = (Old_CO_molec[ilast][mol] +
1620  (Old_CO_molec[ilast+1][mol]-Old_CO_molec[ilast][mol])*
1621  frac_next);
1622  }
1623  else
1624  {
1625  Oldi_mol = Old_CO_molec[ilast][mol];
1626  }
1627  dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.CO_molec[mol][i]/struc.hden[i]);
1628 
1629  /* >>chng 02 nov 11, add first error scale estimate from Robin
1630  * used to normalize the above convergence_error */
1632  }
1633  }
1634 
1635  /* convergence_error is an estimate of the convergence of the solution from its change during the last iteration,
1636  discretization_error is an estimate of the accuracy of the advective terms, calculated in DynaStartZone above:
1637  if dominant error is from the advective terms, need to make them more accurate.
1638  */
1639 
1640  /* report properties of previous iteration */
1641  fprintf(ioQQQ,"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
1644  );
1645 
1646  /* >>chng 02 nov 29, dynamics.convergence_tolerance is now set to 0.1 in init routine */
1648  Dyn_dr /= 1.5;
1649  return;
1650 }
1651 
1652 /*DynaSaveLast save results from previous iteration */
1654 {
1655  long int i,
1656  ion,
1657  nelem,
1658  mol;
1659 
1660  DEBUG_ENTRY( "DynaSaveLast()" );
1661 
1662  /* Save results from previous iteration */
1663  nOld_zone = nzone;
1665  ASSERT( nzone < struc.nzlim );
1666  for( i=0; i<nzone; ++i )
1667  {
1668  Old_histr[i] = struc.histr[i];
1669  Old_depth[i] = struc.depth[i];
1671  /* old n_p density from previous iteration */
1672  Old_hiistr[i] = struc.hiistr[i];
1673  /* old pressure from previous iteration */
1674  Old_pressure[i] = struc.pressure[i];
1675  /* old electron density from previous iteration */
1676  Old_ednstr[i] = struc.ednstr[i];
1677  /* energy term */
1679  /* >>chng 02 May 2001 rjrw: add hden for dilution */
1680  Old_hden[i] = struc.hden[i];
1681  Old_DenMass[i] = struc.DenMass[i];
1682 
1683  for(mol=0;mol<N_H_MOLEC;mol++)
1684  {
1685  Old_H2_molec[i][mol] = struc.H2_molec[mol][i];
1686  }
1687  for(mol=0;mol<mole.num_comole_calc;mol++)
1688  {
1689  Old_CO_molec[i][mol] = struc.CO_molec[mol][i];
1690  }
1691 
1692  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
1693  {
1694  Old_gas_phase[i][nelem] = dense.gas_phase[nelem];
1695  for( ion=0; ion<nelem+2; ++ion )
1696  {
1697  Old_xIonDense[i][nelem][ion] = struc.xIonDense[nelem][ion][i];
1698  }
1699  }
1700  }
1701  return;
1702 }
1703 
1704 realnum DynaFlux(double depth)
1705 
1706 {
1707  realnum flux;
1708 
1709  DEBUG_ENTRY( "DynaFlux()" );
1710 
1711  if(dynamics.FluxIndex == 0)
1712  {
1713  flux = (realnum)dynamics.FluxScale;
1714  }
1715  else
1716  {
1717  flux = (realnum)(dynamics.FluxScale*pow(fabs(depth-dynamics.FluxCenter),dynamics.FluxIndex));
1718  if(depth < dynamics.FluxCenter)
1719  flux = -flux;
1720  }
1721  if(dynamics.lgFluxDScale)
1722  {
1723  /*flux *= struc.DenMass[0]; */
1724  /* WJH 21 may 04, changed to use dense.xMassDensity0, which should be strictly constant */
1725  flux *= dense.xMassDensity0;
1726  }
1727  return flux;
1728 }
1729 
1730 /* ============================================================================== */
1731 /*DynaZero zero some dynamics variables, called from zero.c,
1732  * before parsing commands */
1733 void DynaZero( void )
1734 {
1735  int ipISO;
1736 
1737  DEBUG_ENTRY( "DynaZero()" );
1738 
1739  /* the number of zones in the previous iteration */
1740  nOld_zone = 0;
1741 
1742  /* by default advection is turned off */
1743  dynamics.lgAdvection = false;
1744  /*dynamics.Velocity = 0.;*/
1745  AdvecSpecificEnthalpy = 0.;
1746  dynamics.Cool = 0.;
1747  dynamics.Heat = 0.;
1748  dynamics.dCooldT = 0.;
1749  dynamics.dHeatdT = 0.;
1750  dynamics.HeatMax = 0.;
1751  dynamics.CoolMax = 0.;
1752  dynamics.Rate = 0.;
1753 
1754  /* sets recombination logic, keyword RECOMBINATION on a time step line */
1755  dynamics.lgRecom = false;
1756 
1757  /* set true if time dependent calculation is finished */
1758  dynamics.lgStatic_completed = false;
1759 
1760  /* vars that determine whether time dependent soln only - set with time command */
1761  dynamics.lgStatic = false;
1762  timestep_init = -1.;
1763  /* this factor multiplies the time step */
1764  timestep_factor = 1.2;
1765  dynamics.time_elapsed = 0.;
1766 
1767  /* set the first iteration to include dynamics rather than constant pressure */
1768  /* iteration number, initial iteration is 1, default is 2 - changed with SET DYNAMICS FIRST command */
1770 
1771  /* set initial value of the advection length,
1772  * neg => fraction of depth of init model, + length cm */
1773  dynamics.AdvecLengthInit = -0.1;
1774 
1775  /* this is a tolerance for determining whether dynamics has converged */
1777 
1778  /* this says that set dynamics pressure mode was set */
1779  dynamics.lgSetPresMode = false;
1780 
1781  /* set default values for uniform mass flux */
1782  dynamics.FluxScale = 0.;
1783  dynamics.lgFluxDScale = false;
1784  dynamics.FluxCenter = 0.;
1785  dynamics.FluxIndex = 0.;
1787 
1788  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1789  {
1790  /* factor to allow turning off advection for one of the iso seq,
1791  * this is done with command "no advection h-like" or "he-like"
1792  * only for testing */
1793  dynamics.lgISO[ipISO] = true;
1794  }
1795  /* turn off advection for rest of ions, command "no advection metals" */
1796  dynamics.lgMETALS = true;
1797  /* turn off thermal effects of advection, command "no advection cooling" */
1798  dynamics.lgCoolHeat = true;
1799  DivergePresInteg = 0.;
1800 
1802  dynamics.error_scale2 = 0.;
1803  return;
1804 }
1805 
1806 
1807 /* ============================================================================== */
1808 /* DynaCreateArrays allocate some space needed to save the dynamics structure variables,
1809  * called from DynaCreateArrays */
1810 void DynaCreateArrays( void )
1811 {
1812  long int nelem,
1813  ns,
1814  i,
1815  ion,
1816  mol;
1817 
1818  DEBUG_ENTRY( "DynaCreateArrays()" );
1819 
1820  Upstream_H2_molec = (double*)MALLOC((size_t)N_H_MOLEC*sizeof(double) );
1821  Upstream_CO_molec = (double*)MALLOC((size_t)mole.num_comole_calc*sizeof(double) );
1822 
1823  dynamics.H2_molec = (double*)MALLOC((size_t)N_H_MOLEC*sizeof(double) );
1824  dynamics.CO_molec = (double*)MALLOC((size_t)mole.num_comole_calc*sizeof(double) );
1825 
1826  UpstreamElem = (double*)MALLOC((size_t)LIMELM*sizeof(double) );
1827 
1828  dynamics.Source = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
1829  UpstreamIon = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
1830  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
1831  {
1832  dynamics.Source[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
1833  UpstreamIon[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
1834  for( ion=0; ion<nelem+2; ++ion )
1835  {
1836  dynamics.Source[nelem][ion] = 0.;
1837  }
1838  }
1839  dynamics.Rate = 0.;
1840 
1841  Old_EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1842 
1843  Old_ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1844 
1845  EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1846 
1847  Old_DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1848 
1849  Old_hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1850 
1851  Old_pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1852 
1853  Old_histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1854 
1855  Old_hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1856 
1857  Old_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1858 
1859  Old_xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
1860 
1861  Old_xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) );
1862 
1863  Old_gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
1864 
1865  Old_H2_molec = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
1866 
1867  Old_CO_molec = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
1868 
1869  /* now create diagonal of space for ionization arrays */
1870  for( ns=0; ns < struc.nzlim; ++ns )
1871  {
1872  Old_xIonDense[ns] =
1873  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+3) );
1874 
1875  Old_gas_phase[ns] =
1876  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+3) );
1877 
1878  Old_H2_molec[ns] =
1879  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(N_H_MOLEC) );
1880 
1881  Old_CO_molec[ns] =
1882  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(mole.num_comole_calc) );
1883 
1884  for( nelem=0; nelem< (LIMELM+3);++nelem )
1885  {
1886  Old_xIonDense[ns][nelem] =
1887  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+1) );
1888  }
1889  }
1890 
1891  for( i=0; i < struc.nzlim; i++ )
1892  {
1893  /* these are values if H0 and tau_912 from previous iteration */
1894  Old_histr[i] = 0.;
1895  Old_xLyman_depth[i] = 0.;
1896  Old_depth[i] = 0.;
1897  dynamics.oldFullDepth = 0.;
1898  /* old n_p density from previous iteration */
1899  Old_hiistr[i] = 0.;
1900  /* old pressure from previous iteration */
1901  Old_pressure[i] = 0.;
1902  /* old electron density from previous iteration */
1903  Old_ednstr[i] = 0.;
1904  Old_hden[i] = 0.;
1905  Old_DenMass[i] = 0.;
1906  Old_EnthalpyDensity[i] = 0.;
1907  for( nelem=0; nelem< (LIMELM+3);++nelem )
1908  {
1909  for( ion=0; ion<LIMELM+1; ++ion )
1910  {
1911  Old_xIonDense[i][nelem][ion] = 0.;
1912  }
1913  }
1914  for(mol=0;mol<N_H_MOLEC;mol++)
1915  {
1916  Old_H2_molec[i][mol] = 0.;
1917  }
1918  for(mol=0;mol<mole.num_comole_calc;mol++)
1919  {
1920  Old_CO_molec[i][mol] = 0.;
1921  }
1922  }
1923  return;
1924 }
1925 
1926 /*advection_set_detault - called to set default conditions
1927  * when time and wind commands are parsed,
1928  * lgWind is true if dynamics, false if time dependent */
1929 STATIC void advection_set_detault( bool lgWind )
1930 {
1931 
1932  DEBUG_ENTRY( "advection_set_detault()" );
1933 
1934  /* turn on advection */
1935  dynamics.lgAdvection = true;
1936 
1937  /* turn off prediction of next zone's temperature, as guessed in ZoneStart,
1938  * also set with no tepredictor */
1939  thermal.lgPredNextTe = false;
1940 
1941  /* turn off both CO and H2 networks since advection not included, also say physical
1942  * conditions are not ok*/
1943  /* the heavy element molecules are turned off for the moment,
1944  * until I am confident that the H-only models work robustly */
1945  /* co.lgNoH2Mole = true; */
1947  /* >>chng 06 jun 29, add advective terms to CO solver */
1948  co.lgNoCOMole = true;
1949 # if 0
1950  co.lgNoCOMole = true; /* >>chng 04 apr 23, tried to rm this line - problems */
1951  phycon.lgPhysOK = false;/* >>chng 04 apr 23, rm this line */
1952 # endif
1953  /* >>chng 06 nov 29, there is a conservation problem in the ionization -
1954  * molecular solvers that is demonstrated by the */
1957  /* use the new temperature solver
1958  strcpy( conv.chSolverEden , "new" ); */
1959 
1960  /* constant total pressure, gas+rad+incident continuum
1961  * turn on radiation pressure */
1964  pressure.lgPres_ram_ON = true;
1965 
1966  /* we need to make the solvers much more exact when advection is in place */
1967  if( lgWind )
1968  {
1969  /* increse precision of solution */
1970  conv.EdenErrorAllowed = 1e-3;
1971  /* the actual relative error is relative to the total heating and cooling,
1972  * which include the dynamics.heat and .cool, which are the advected heating/cooling.
1973  * the two terms can be large and nearly cancel, what is written to the .heat and cool files
1974  * by punch files has had the smaller of the two subtracted, leaving only the net advected
1975  * heating and cooling */
1976  conv.HeatCoolRelErrorAllowed = 3e-4f;
1977  conv.PressureErrorAllowed = 1e-3f;
1978  }
1979  return;
1980 }
1981 
1982 /* ============================================================================== */
1983 /* ParseDynaTime parse the time command, called from ParseCommands */
1984 void ParseDynaTime( char *chCard )
1985 {
1986  long int i;
1987  bool lgEOL;
1988  char chCap[INPUT_LINE_LENGTH];
1989 
1990  DEBUG_ENTRY( "ParseDynaTime()" );
1991 
1992  /*flag set true when time dependent only */
1993  dynamics.lgStatic = true;
1994 
1995  i = 5;
1996  timestep_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1997  /* this is log of timestep */
1998  if( timestep_init > 30. )
1999  {
2000  /* this number is large and will almost certainly crash */
2001  fprintf(ioQQQ,"WARNING - the log of the initial time step is too "
2002  "large, I shall probably crash. The value was log t = %.2e\n",
2003  timestep_init );
2004  fflush(ioQQQ);
2005  }
2006  timestep_init = pow( 10. , timestep_init );
2008  if( lgEOL )
2009  {
2010  NoNumb(chCard);
2011  }
2012 
2013  /* this is the stop time and is optional */
2014  timestep_stop = pow( 10. , FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) );
2015  if( lgEOL )
2016  {
2017  timestep_stop = -1.;
2018  }
2019 
2020  /* set default flags - false says that time dependent, not dynamical solution */
2021  advection_set_detault(false);
2022 
2023  wind.windv0 = 0.;
2024  wind.windv = wind.windv0;
2025 
2026  /* this is used in convpres to say wind solution - both cases use this*/
2027  strcpy( dense.chDenseLaw, "WIND" );
2028 
2029  /* create time step and flux arrays */
2030  time_elapsed_time = (double*)MALLOC((size_t)NTIME*sizeof(double));
2031  time_flux_ratio = (double*)MALLOC((size_t)NTIME*sizeof(double));
2032  time_dt = (double*)MALLOC((size_t)NTIME*sizeof(double));
2033  time_dt_scale_factor = (double*)MALLOC((size_t)NTIME*sizeof(double));
2034  lgtime_Recom = (int*)MALLOC((size_t)NTIME*sizeof(int));
2035 
2036  /* number of lines we will save */
2037  nTime_flux = 0;
2038 
2039  /* get the next line, and check for eof */
2040  input_readarray(chCard,&lgEOL);
2041  if( lgEOL )
2042  {
2043  fprintf( ioQQQ,
2044  " Hit EOF while reading time-continuum list; use END to end list.\n" );
2045  cdEXIT(EXIT_FAILURE);
2046  }
2047 
2048  /* convert line to caps */
2049  strcpy( chCap, chCard );
2050  caps(chCap);
2051 
2052  /* third column might set dt - if any third column is missing then
2053  * this is set false and only time on command line is used */
2054  lgtime_dt_specified = true;
2055 
2056  while( strncmp(chCap, "END" ,3 ) != 0 )
2057  {
2058  if( nTime_flux >= NTIME )
2059  {
2060  fprintf( ioQQQ,
2061  " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n",
2062  NTIME );
2063  cdEXIT(EXIT_FAILURE);
2064  }
2065 
2066  i = 1;
2067  time_elapsed_time[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
2068  time_flux_ratio[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
2069  if( lgEOL )
2070  {
2071  fprintf( ioQQQ,
2072  " each line must have two numbers, log of time (s0 and flux ratio.\n");
2073  cdEXIT(EXIT_FAILURE);
2074  }
2075  /* this is optional dt to set time step - if not given then initial
2076  * time step is always used */
2077  time_dt[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
2078 
2079  /* if any of these are not specified then do not use times array */
2080  if( lgEOL )
2081  lgtime_dt_specified = false;
2082 
2083  /* this is optional scale factor to increase time */
2084  time_dt_scale_factor[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL));
2085  if( lgEOL )
2087 
2088  /* turn on recombination front logic */
2089  if( nMatch("RECO" , chCap ) && nMatch("MBIN" , chCap ) )
2090  {
2091  /* this sets flag dynamics.lgRecom true so that all of code knows recombination
2092  * is taking place */
2093  lgtime_Recom[nTime_flux] = true;
2094  }
2095  else
2096  {
2097  lgtime_Recom[nTime_flux] = false;
2098  }
2099 
2100  /* this is total number stored so far */
2101  ++nTime_flux;
2102 
2103  /* get next line and check for eof */
2104  input_readarray(chCard,&lgEOL);
2105  if( lgEOL )
2106  {
2107  fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
2108  cdEXIT(EXIT_FAILURE);
2109  }
2110 
2111  /* convert it to caps */
2112  strcpy( chCap, chCard );
2113  caps(chCap);
2114  }
2115 
2116  for( i=0; i < nTime_flux; i++ )
2117  {
2118  fprintf( ioQQQ, "DEBUG time dep %.2e %.2e %.2e %.2e\n",
2119  time_elapsed_time[i],
2120  time_flux_ratio[i] ,
2121  time_dt[i],
2123  }
2124  fprintf( ioQQQ, "\n" );
2125  return;
2126 }
2127 /* ============================================================================== */
2128 /* ParseDynaWind parse the wind command, called from ParseCommands */
2129 void ParseDynaWind( char *chCard )
2130 {
2131  long int i;
2132  int iVelocity_Type;
2133  bool lgEOL;
2134  /* compiler flagged possible paths where dfdr used but not set -
2135  * this is for safety/keep it happy */
2136  double dfdr=-BIGDOUBLE;
2137 
2138  DEBUG_ENTRY( "ParseDynaWind()" );
2139 
2140  /* Flag for type of velocity law:
2141  * 1 is original, give initial velocity at illuminated face
2142  * 2 is face flux gradient (useful if face velocity is zero),
2143  * set to zero, but will be reset if velocity specified */
2144  iVelocity_Type = 0;
2145  /* wind structure, parameters are initial velocity and optional mass
2146  * v read in in km s-1 and convert to cm s-1, mass in solar masses */
2147  if( (i = nMatch("VELO",chCard))>0 )
2148  {
2149  i += 5;
2150  wind.windv0 = (realnum)(FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)*1e5);
2151  wind.windv = wind.windv0;
2152  iVelocity_Type = 1;
2153  }
2154 
2155  if( (i = nMatch("DFDR",chCard))>0 )
2156  {
2157  /* velocity not specified, rather mass flux gradient */
2158  i += 5;
2159  dfdr = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
2160  iVelocity_Type = 2;
2161  }
2162 
2163  /* center option, gives xxx */
2164  if( (i = nMatch("CENT",chCard))>0 )
2165  {
2166  /* physical length in cm, can be either sign */
2167  i += 5;
2168  dynamics.FluxCenter = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
2169  }
2170 
2171  /* flux index */
2172  if( (i = nMatch("INDE",chCard))>0 )
2173  {
2174  /* power law index */
2175  i += 5;
2176  dynamics.FluxIndex = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
2177  }
2178 
2179  /* the case where velocity was set */
2180  if(iVelocity_Type == 1)
2181  {
2182  /* was flux index also set? */
2183  if(dynamics.FluxIndex == 0)
2184  {
2186  dynamics.lgFluxDScale = true;
2187  /* Center doesn't mean much in this case -- make sure it's
2188  * in front of grid so DynaFlux doesn't swap signs where
2189  * it shouldn't */
2190  dynamics.FluxCenter = -1.;
2191  }
2192  else
2193  {
2196  /* velocity was set but flux index was not set - estimate it */
2198  pow(fabs(dynamics.FluxCenter),-dynamics.FluxIndex);
2199 
2200  dynamics.lgFluxDScale = true;
2201  if(dynamics.FluxCenter > 0)
2202  {
2204  }
2205  }
2206  }
2207  /* the case where flux gradient is set */
2208  else if(iVelocity_Type == 2)
2209  {
2210  if(dynamics.FluxIndex == 0)
2211  {
2212  fprintf(ioQQQ,"Can't specify gradient when flux is constant!\n");
2213  /* use this exit handler, which closes out MPI when multiprocessing */
2214  cdEXIT(EXIT_FAILURE);
2215  }
2218  /* Can't specify FluxScale from dvdr rather than dfdr, as
2219  * d(rho)/dr != 0 */
2221  pow(fabs(dynamics.FluxCenter),1.-dynamics.FluxIndex);
2222  if(dynamics.FluxCenter > 0)
2223  {
2225  }
2226  dynamics.lgFluxDScale = false;
2227 
2228  /* put in bogus value simply as flag -- assume that surface velocity
2229  * is small or we wouldn't be using this to specify. */
2230  wind.windv0 = -0.01f;
2231  }
2232  else
2233  {
2234  /* assume old-style velocity-only specification */
2235  /* wind structure, parameters are initial velocity and optional mass
2236  * v in km/sec, mass in solar masses */
2237  i = 5;
2238  wind.windv0 = (realnum)(FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)*1e5);
2240  dynamics.FluxIndex = 0.;
2241  dynamics.lgFluxDScale = true;
2242  /* Center doesn't mean much in this case -- make sure it's
2243  * in front of grid so DynaFlux doesn't swap signs where
2244  * it shouldn't */
2245  dynamics.FluxCenter = -1.;
2246  if( lgEOL )
2247  {
2248  NoNumb(chCard);
2249  }
2250  }
2251 
2252  wind.windv = wind.windv0;
2253 
2254 # ifdef FOO
2255  fprintf(ioQQQ,"Scale %g (*%c) Index %g Center %g\n",
2258 # endif
2259 
2260  /* option to include advection */
2261  if( nMatch( "ADVE" , chCard ) )
2262  {
2263  /* set default flags - true says dynamical solution */
2264  advection_set_detault(true);
2265  }
2266 
2267  else
2268  {
2269  /* this is usual hypersonic outflow */
2270  if( wind.windv0 <= 1.e6 )
2271  {
2272  /* speed of sound roughly 10 km/s */
2273  fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
2274  wind.lgWindOK = false;
2275  }
2276 
2277  /* set the central object mass, in solar masses */
2278  wind.comass = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
2279  /* default is one solar mass */
2280  if( lgEOL )
2281  wind.comass = 1.;
2282 
2283  /* option for rotating disk, keyword is disk */
2284  wind.lgDisk = false;
2285  if( nMatch( "DISK" , chCard ) )
2286  wind.lgDisk = true;
2287 
2288  }
2289 
2290  /* this is used in convpres to say wind solution - both cases use this*/
2291  strcpy( dense.chDenseLaw, "WIND" );
2292 
2293  /* option to turn off continuum radiative acceleration */
2294  if( nMatch("O CO",chCard) )
2295  {
2296  pressure.lgContRadPresOn = false;
2297  }
2298  else
2299  {
2300  pressure.lgContRadPresOn = true;
2301  }
2302  return;
2303 }
2304 
2305 /*DynaPrtZone called to print zone results */
2306 void DynaPrtZone( void )
2307 {
2308 
2309  DEBUG_ENTRY( "DynaPrtZone()" );
2310 
2311  ASSERT( nzone>0 && nzone<struc.nzlim );
2312 
2313  if( nzone > 0 )
2314  {
2315  fprintf(ioQQQ," DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
2317  wind.windv/1e5 ,
2320  }
2321 
2322  ASSERT( EnthalpyDensity[nzone-1] > 0. );
2323 
2324  fprintf(ioQQQ," DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
2329  5./2.*pressure.PresGasCurr ,
2331  );
2332  return;
2333 }
2334 
2335 /*DynaPunchTimeDep - punch info about time dependent solution */
2336 void DynaPunchTimeDep( FILE* ipPnunit , const char *chJob )
2337 {
2338 
2339  DEBUG_ENTRY( "DynaPunchTimeDep()" );
2340 
2341  if( strcmp( chJob , "END" ) == 0 )
2342  {
2343  double te_mean,
2344  H2_mean,
2345  H0_mean,
2346  Hp_mean,
2347  Hep_mean;
2348  /* punch info at end */
2349  if( cdTemp(
2350  /* four char string, null terminated, giving the element name */
2351  "HYDR",
2352  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2353  * 0 means that chLabel is a special case */
2354  2,
2355  /* will be temperature */
2356  &te_mean,
2357  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2358  "RADIUS" ) )
2359  {
2360  TotalInsanity();
2361  }
2362  if( cdIonFrac(
2363  /* four char string, null terminated, giving the element name */
2364  "HYDR",
2365  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2366  * 0 says special case */
2367  2,
2368  /* will be fractional ionization */
2369  &Hp_mean,
2370  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2371  "RADIUS" ,
2372  /* if true then weighting also has electron density, if false then only volume or radius */
2373  false ) )
2374  {
2375  TotalInsanity();
2376  }
2377  if( cdIonFrac(
2378  /* four char string, null terminated, giving the element name */
2379  "HYDR",
2380  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2381  * 0 says special case */
2382  1,
2383  /* will be fractional ionization */
2384  &H0_mean,
2385  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2386  "RADIUS" ,
2387  /* if true then weighting also has electron density, if false then only volume or radius */
2388  false ) )
2389  {
2390  TotalInsanity();
2391  }
2392  if( cdIonFrac(
2393  /* four char string, null terminated, giving the element name */
2394  "H2 ",
2395  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2396  * 0 says special case */
2397  0,
2398  /* will be fractional ionization */
2399  &H2_mean,
2400  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2401  "RADIUS" ,
2402  /* if true then weighting also has electron density, if false then only volume or radius */
2403  false ) )
2404  {
2405  TotalInsanity();
2406  }
2407  if( cdIonFrac(
2408  /* four char string, null terminated, giving the element name */
2409  "HELI",
2410  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2411  * 0 says special case */
2412  2,
2413  /* will be fractional ionization */
2414  &Hep_mean,
2415  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2416  "RADIUS" ,
2417  /* if true then weighting also has electron density, if false then only volume or radius */
2418  false ) )
2419  {
2420  TotalInsanity();
2421  }
2422  fprintf( ipPnunit ,
2423  "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" ,
2425  timestep ,
2428  te_mean ,
2429  Hp_mean ,
2430  H0_mean ,
2431  H2_mean ,
2432  Hep_mean ,
2433  /* ratio of CO to total H column densities */
2434  findspecies("CO")->hevcol / SDIV( colden.colden[ipCOL_HTOT] ));
2435  }
2436  else
2437  TotalInsanity();
2438  return;
2439 }
2440 
2441 /*DynaPunch punch dynamics - info related to advection */
2442 void DynaPunch(FILE* ipPnunit , char chJob )
2443 {
2444  DEBUG_ENTRY( "DynaPunch()" );
2445 
2446  if( chJob=='a' )
2447  {
2448  /* this is punch dynamics advection, the only punch dynamics */
2449  fprintf( ipPnunit , "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2451  thermal.htot ,
2452  dynamics.Cool ,
2453  dynamics.Heat ,
2454  dynamics.dCooldT ,
2456  dynamics.Rate,
2457  phycon.EnthalpyDensity/dense.gas_phase[ipHYDROGEN] ,
2459  );
2460  }
2461  else
2462  TotalInsanity();
2463  return;
2464 }

Generated for cloudy by doxygen 1.8.3.1