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