00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "cddefines.h"
00019 #include "cddrive.h"
00020 #include "struc.h"
00021 #include "input.h"
00022 #include "colden.h"
00023 #include "radius.h"
00024 #include "thirdparty.h"
00025 #include "stopcalc.h"
00026 #include "hextra.h"
00027 #include "rfield.h"
00028 #include "iterations.h"
00029 #include "trace.h"
00030 #include "conv.h"
00031 #include "timesc.h"
00032 #include "dense.h"
00033 #include "mole.h"
00034 #include "thermal.h"
00035 #include "pressure.h"
00036 #include "phycon.h"
00037 #include "wind.h"
00038 #include "hmi.h"
00039 #include "iso.h"
00040 #include "dynamics.h"
00041 #include "cosmology.h"
00042 #include "taulines.h"
00043 #include "parser.h"
00044 t_dynamics dynamics;
00045 static int ipUpstream=-1,iphUpstream=-1,ipyUpstream=-1;
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 static double **UpstreamIon;
00061 static double ***UpstreamStatesElem;
00062
00063 static double *UpstreamElem;
00064
00065
00066 static double *Upstream_molecules;
00067
00068
00069
00070
00071
00072 static double *time_elapsed_time ,
00073 *time_flux_ratio ,
00074 *time_dt,
00075 *time_dt_scale_factor;
00076 bool lgtime_dt_specified;
00077 int *lgtime_Recom;
00078 #define NTIME 200
00079
00080
00081 static long int nTime_flux=0;
00082
00083
00084 STATIC void DynaNewStep(void);
00085
00086
00087 STATIC void DynaSaveLast(void);
00088
00089
00090
00091
00092
00093 static double Dyn_dr;
00094
00095
00096 static double AdvecSpecificEnthalpy;
00097
00098
00099 static realnum *Old_histr ,
00100
00101 *Old_xLyman_depth,
00102
00103 *Old_depth,
00104
00105 *Old_hiistr,
00106
00107 *Old_pressure,
00108
00109 *Old_density ,
00110
00111 *Old_DenMass ,
00112
00113 *EnthalpyDensity,
00114
00115 *Old_ednstr,
00116
00117 *Old_EnthalpyDensity;
00118
00119 static realnum **Old_molecules;
00120
00121
00122 static realnum ***Old_xIonDense;
00123
00124
00125 static realnum **Old_gas_phase;
00126
00127
00128 static realnum ****Old_StatesElem;
00129
00130
00131 static long int nOld_zone;
00132
00133
00134 STATIC double timestep_next( void )
00135 {
00136 static double te_old=-1;
00137 double timestep_Hp_temp , timestep_return;
00138
00139 DEBUG_ENTRY( "timestep_next()" );
00140
00141 timestep_return = dynamics.timestep;
00142
00143 if( dynamics.lgRecom )
00144 {
00145 double te_new;
00146 if( cdTemp(
00147
00148 "hydr",
00149
00150 2 ,
00151
00152 &te_new,
00153
00154 "VOLUME" ) )
00155 TotalInsanity();
00156
00157 if( te_old>0 )
00158 {
00159 double dTdStep = fabs(te_new-te_old)/te_new;
00160
00161 double dT_factor = 0.04/SDIV(dTdStep);
00162 dT_factor = MIN2( 2.0 , dT_factor );
00163 dT_factor = MAX2( 0.01 , dT_factor );
00164 timestep_Hp_temp = dynamics.timestep * dT_factor;
00165 }
00166 else
00167 {
00168 timestep_Hp_temp = -1.;
00169 }
00170 te_old = te_new;
00171 if( timestep_Hp_temp > 0. )
00172 timestep_return = timestep_Hp_temp;
00173 }
00174 else
00175 {
00176 timestep_return = dynamics.timestep_init;
00177 }
00178 fprintf(ioQQQ,"DEBUG timestep_next returns %.3e, old temp %.2e\n" , timestep_return, te_old );
00179
00180 return( timestep_return );
00181 }
00182
00183
00184
00185
00186 void DynaIonize(void)
00187 {
00188 DEBUG_ENTRY( "DynaIonize()" );
00189
00190
00191
00192 if( !dynamics.lgTimeDependentStatic )
00193 {
00194
00195
00196
00197 dynamics.timestep = -Dyn_dr/wind.windv;
00198 }
00199
00200
00201 ASSERT(nzone<struc.nzlim );
00202 if( nzone > 0 )
00203 EnthalpyDensity[nzone-1] = (realnum)phycon.EnthalpyDensity;
00204
00205
00206
00207
00208
00209
00210
00211 double depth = radius.depth;
00212 if( iteration < dynamics.n_initial_relax+1 ||
00213 ( ! dynamics.lgTimeDependentStatic &&
00214 ( depth < 0 || depth > dynamics.oldFullDepth ) ) )
00215 {
00216
00217 dynamics.Cool_r = 0.;
00218 dynamics.Heat_v = 0.;
00219 dynamics.dHeatdT = 0.;
00220
00221 dynamics.Rate = 0.;
00222
00223 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00224 {
00225 for( long ion=0; ion<nelem+2; ++ion )
00226 {
00227 dynamics.Source[nelem][ion] = 0.;
00228 }
00229 }
00230 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00231 {
00232 for( long nelem=ipISO; nelem<LIMELM; ++nelem)
00233 {
00234 if( dense.lgElmtOn[nelem] )
00235 {
00236 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; ++level )
00237 {
00238 dynamics.StatesElem[nelem][nelem-ipISO][level] = 0.;
00239 }
00240 }
00241 }
00242 }
00243 for( long mol=0;mol<mole_global.num_calc;mol++)
00244 {
00245 dynamics.molecules[mol] = 0.;
00246 }
00247 return;
00248 }
00249
00250 if( dynamics.lgTracePrint )
00251 {
00252 fprintf(ioQQQ,"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
00253 nzone,
00254 phycon.EnthalpyDensity,
00255 0.5*POW2(wind.windv)*dense.xMassDensity ,
00256 5./2.*pressure.PresGasCurr
00257 );
00258 }
00259
00260
00261
00262
00263
00264 dynamics.Cool_r = 1./dynamics.timestep*dynamics.lgCoolHeat;
00265 dynamics.Heat_v = AdvecSpecificEnthalpy/dynamics.timestep*dynamics.lgCoolHeat;
00266 dynamics.dHeatdT = 0.*dynamics.lgCoolHeat;
00267
00268 # if 0
00269 if( dynamics.lgTracePrint || nzone>17 && iteration == 10)
00270 {
00271 fprintf(ioQQQ,
00272 "dynamics cool-heat\t%li\t%.3e\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
00273 nzone,
00274 phycon.te,
00275 dynamics.lgCoolHeat,
00276 thermal.htot,
00277 phycon.EnthalpyDensity/dynamics.timestep,
00278 AdvecSpecificEnthalpy*scalingDensity()/dynamics.timestep,
00279 phycon.EnthalpyDensity,
00280 AdvecSpecificEnthalpy*scalingDensity(),
00281 scalingDensity(),
00282 dynamics.timestep);
00283 }
00284 # endif
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299 dynamics.Rate = 1./dynamics.timestep;
00300
00301 for( long mol=0;mol<mole_global.num_calc;mol++)
00302 {
00303 dynamics.molecules[mol] = Upstream_molecules[mol]*scalingDensity();
00304 }
00305
00306 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00307 {
00308 if( dense.lgElmtOn[nelem] )
00309 {
00310
00311 if(fabs(UpstreamElem[nelem]*scalingDensity() -dense.gas_phase[nelem])/dense.gas_phase[nelem]>=1e-3)
00312 {
00313 fprintf(ioQQQ,
00314 "PROBLEM conservation error: zn %li elem %li upstream %.8e abund %.8e (up-ab)/up %.2e\n",
00315 nzone ,
00316 nelem,
00317 UpstreamElem[nelem]*scalingDensity(),
00318 dense.gas_phase[nelem] ,
00319 (UpstreamElem[nelem]*scalingDensity()-dense.gas_phase[nelem]) /
00320 (UpstreamElem[nelem]*scalingDensity()) );
00321 }
00322
00323 for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
00324 {
00325 dynamics.Source[nelem][ion] = 0.;
00326 }
00327 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00328 {
00329
00330
00331
00332
00333 dynamics.Source[nelem][ion] =
00334
00335
00336
00337 UpstreamIon[nelem][ion] * scalingDensity() / dynamics.timestep;
00338
00339 }
00340 for( long ion=dense.IonHigh[nelem]+1;ion<nelem+2; ++ion )
00341 {
00342 dynamics.Source[nelem][ion] = 0.;
00343 dynamics.Source[nelem][dense.IonHigh[nelem]] +=
00344 UpstreamIon[nelem][ion] * scalingDensity() / dynamics.timestep;
00345 }
00346 }
00347 }
00348
00349 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00350 {
00351 for( long nelem=ipISO; nelem<LIMELM; ++nelem)
00352 {
00353 if( dense.lgElmtOn[nelem] )
00354 {
00355 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; ++level )
00356 {
00357 dynamics.StatesElem[nelem][nelem-ipISO][level] =
00358 UpstreamStatesElem[nelem][nelem-ipISO][level] * scalingDensity() / dynamics.timestep;
00359 }
00360 }
00361 }
00362 }
00363
00364 # if 0
00365 fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
00366 nzone,
00367 dynamics.Rate,
00368 dynamics.Source[ipHYDROGEN][0],
00369 dynamics.Rate,
00370 dynamics.Source[ipCARBON][3]);
00371 # endif
00372 #if 0
00373 long nelem = ipCARBON;
00374 long ion = 3;
00375
00376 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",
00377 nzone,
00378 ipUpstream,
00379 radius.depth ,
00380 Old_depth[ipUpstream],
00381 dense.xIonDense[nelem][ion],
00382 UpstreamIon[nelem][ion]* scalingDensity(),
00383 Old_xIonDense[ipUpstream][nelem][ion] ,
00384 dense.xIonDense[nelem][ion+1],
00385 UpstreamIon[nelem][ion+1]* scalingDensity(),
00386 Old_xIonDense[ipUpstream][nelem][ion+1] ,
00387 dynamics.timestep,
00388 dynamics.Source[nelem][ion]
00389 );
00390 #endif
00391 if( dynamics.lgTracePrint )
00392 {
00393 fprintf(ioQQQ," DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
00394 nzone,dynamics.Rate , dynamics.Source[0][0] );
00395 }
00396 return;
00397 }
00398
00399
00400
00401 void DynaStartZone(void)
00402 {
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415 double upstream, dilution, dilutionleft, dilutionright, frac_next;
00416
00417
00418 double hupstream, hnextfrac=-BIGFLOAT, hion, hmol, hiso;
00419
00420
00421 double ynextfrac=-BIGFLOAT, yion, ymol, yiso;
00422
00423 long int nelem , ion, mol;
00424
00425 DEBUG_ENTRY( "DynaStartZone()" );
00426
00427
00428 if( iteration < 2 )
00429 {
00430 dynamics.Upstream_density = 0.;
00431 AdvecSpecificEnthalpy = 0.;
00432 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00433 {
00434 for( ion=0; ion<nelem+2; ++ion )
00435 {
00436 UpstreamIon[nelem][ion] = 0.;
00437 }
00438 }
00439 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00440 {
00441 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00442 {
00443 if( dense.lgElmtOn[nelem] )
00444 {
00445 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
00446 {
00447 UpstreamStatesElem[nelem][nelem-ipISO][level] = 0.;
00448 }
00449 }
00450 }
00451 }
00452
00453 for(mol=0; mol<mole_global.num_calc;mol++)
00454 {
00455 Upstream_molecules[mol] = 0;
00456 }
00457
00458 ipUpstream = -1;
00459 iphUpstream = -1;
00460 ipyUpstream = -1;
00461 return;
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473 upstream = MAX2(Old_depth[0] , radius.depth + Dyn_dr);
00474 hupstream = 0.5*(radius.depth + upstream);
00475
00476
00477
00478 while( (Old_depth[ipUpstream+1] < upstream ) &&
00479 ( ipUpstream < nOld_zone-1 ) )
00480 {
00481 ++ipUpstream;
00482 }
00483 ASSERT( ipUpstream <= nOld_zone-1 );
00484
00485
00486
00487 if( (ipUpstream != -1) && (ipUpstream != nOld_zone-1)&& ((Old_depth[ipUpstream+1] - Old_depth[ipUpstream])> SMALLFLOAT) )
00488 {
00489
00490 frac_next = ( upstream - Old_depth[ipUpstream])/
00491 (Old_depth[ipUpstream+1] - Old_depth[ipUpstream]);
00492 dynamics.Upstream_density = (realnum)(Old_density[ipUpstream] +
00493 (Old_density[ipUpstream+1] - Old_density[ipUpstream])*
00494 frac_next);
00495 dilutionleft = 1./Old_density[ipUpstream];
00496 dilutionright = 1./Old_density[ipUpstream+1];
00497
00498
00499
00500
00501 dilution = 1./dynamics.Upstream_density;
00502
00503
00504 AdvecSpecificEnthalpy = (Old_EnthalpyDensity[ipUpstream]*dilutionleft +
00505 (Old_EnthalpyDensity[ipUpstream+1]*dilutionright - Old_EnthalpyDensity[ipUpstream]*dilutionleft)*
00506 frac_next);
00507
00508 ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] );
00509
00510
00511 realnum lo = (realnum)(Old_EnthalpyDensity[ipUpstream]*dilutionleft);
00512 realnum x = (realnum)AdvecSpecificEnthalpy;
00513 realnum hi = (realnum)(Old_EnthalpyDensity[ipUpstream+1]*dilutionright);
00514 if( ! fp_bound(lo,x,hi) )
00515 {
00516 fprintf(ioQQQ,"PROBLEM interpolated enthalpy density is not within range %.16e\t%.16e\t%.16e\t%e\t%e\n",
00517 lo, x, hi, (hi-x)/(hi-lo), (x-lo)/(hi-lo));
00518 }
00519
00520 ASSERT( fp_bound(lo,x,hi) );
00521
00522 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00523 {
00524 UpstreamElem[nelem] = 0.;
00525 for( ion=0; ion<nelem+2; ++ion )
00526 {
00527
00528
00529
00530
00531
00532
00533
00534 UpstreamIon[nelem][ion] =
00535 Old_xIonDense[ipUpstream][nelem][ion]/Old_density[ipUpstream] +
00536 (Old_xIonDense[ipUpstream+1][nelem][ion]/Old_density[ipUpstream+1] -
00537 Old_xIonDense[ipUpstream][nelem][ion]/Old_density[ipUpstream])*
00538 frac_next;
00539
00540 UpstreamElem[nelem] += UpstreamIon[nelem][ion];
00541 }
00542 }
00543
00544 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00545 {
00546 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00547 {
00548 if( dense.lgElmtOn[nelem] )
00549 {
00550 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
00551 {
00552 UpstreamStatesElem[nelem][nelem-ipISO][level] =
00553 Old_StatesElem[ipUpstream][nelem][nelem-ipISO][level]/Old_density[ipUpstream] +
00554 (Old_StatesElem[ipUpstream+1][nelem][nelem-ipISO][level]/Old_density[ipUpstream+1] -
00555 Old_StatesElem[ipUpstream][nelem][nelem-ipISO][level]/Old_density[ipUpstream])*
00556 frac_next;
00557
00558 ASSERT( !isnan( UpstreamStatesElem[nelem][nelem-ipISO][level] ) );
00559 }
00560 }
00561 }
00562 }
00563
00564 for(mol=0;mol<mole_global.num_calc;mol++)
00565 {
00566 Upstream_molecules[mol] = Old_molecules[ipUpstream][mol]/Old_density[ipUpstream] +
00567 (Old_molecules[ipUpstream+1][mol]/Old_density[ipUpstream+1] -
00568 Old_molecules[ipUpstream][mol]/Old_density[ipUpstream])*
00569 frac_next;
00570
00571 if(mole.species[mol].location == NULL && mole_global.list[mol]->parentLabel.empty())
00572 {
00573 for(molecule::nAtomsMap::iterator atom= mole_global.list[mol]->nAtom.begin();
00574 atom != mole_global.list[mol]->nAtom.end(); ++atom)
00575 {
00576 UpstreamElem[atom->first->el->Z-1] +=
00577 Upstream_molecules[mol]*atom->second;
00578 }
00579 }
00580 }
00581 }
00582 else
00583 {
00584
00585 long ipBound = ipUpstream;
00586 if (ipBound == -1)
00587 ipBound = 0;
00588 dynamics.Upstream_density = Old_density[ipBound];
00589
00590 dilution = 1./dynamics.Upstream_density;
00591
00592 AdvecSpecificEnthalpy = Old_EnthalpyDensity[ipBound]/Old_density[ipBound];
00593 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00594 {
00595 UpstreamElem[nelem] = 0.;
00596 for( ion=0; ion<nelem+2; ++ion )
00597 {
00598
00599 UpstreamIon[nelem][ion] =
00600 Old_xIonDense[ipBound][nelem][ion]/Old_density[ipBound];
00601 UpstreamElem[nelem] += UpstreamIon[nelem][ion];
00602 }
00603 }
00604
00605 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00606 {
00607 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00608 {
00609 if( dense.lgElmtOn[nelem] )
00610 {
00611 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
00612 {
00613 UpstreamStatesElem[nelem][nelem-ipISO][level] =
00614 Old_StatesElem[ipBound][nelem][nelem-ipISO][level]/Old_density[ipBound];
00615
00616 ASSERT( !isnan( UpstreamStatesElem[nelem][nelem-ipISO][level] ) );
00617 }
00618 }
00619 }
00620 }
00621
00622 for(mol=0;mol<mole_global.num_calc;mol++)
00623 {
00624 Upstream_molecules[mol] = Old_molecules[ipBound][mol]/Old_density[ipBound];
00625 if(mole.species[mol].location == NULL && mole_global.list[mol]->parentLabel.empty())
00626 {
00627 for(molecule::nAtomsMap::iterator atom=mole_global.list[mol]->nAtom.begin();
00628 atom != mole_global.list[mol]->nAtom.end(); ++atom)
00629 {
00630 UpstreamElem[atom->first->el->Z-1] +=
00631 Upstream_molecules[mol]*atom->second;
00632 }
00633 }
00634 }
00635 }
00636
00637
00638
00639
00640 while( (Old_depth[iphUpstream+1] < hupstream ) &&
00641 ( iphUpstream < nOld_zone-1 ) )
00642 {
00643 ++iphUpstream;
00644 }
00645 ASSERT( iphUpstream <= nOld_zone-1 );
00646
00647 while( (Old_depth[ipyUpstream+1] < radius.depth ) &&
00648 ( ipyUpstream < nOld_zone-1 ) )
00649 {
00650 ++ipyUpstream;
00651 }
00652 ASSERT( ipyUpstream <= nOld_zone-1 );
00653
00654 dynamics.dRad = BIGFLOAT;
00655
00656 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
00657 hnextfrac = ( hupstream - Old_depth[iphUpstream])/
00658 (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]);
00659 else
00660 hnextfrac = 0.;
00661
00662 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
00663 ynextfrac = ( radius.depth - Old_depth[ipyUpstream])/
00664 (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]);
00665 else
00666 ynextfrac = 0.;
00667
00668
00669 if(ipUpstream != -1 && ipUpstream < nOld_zone-1)
00670 dynamics.dRad = MIN2(dynamics.dRad,
00671 5*(Old_depth[ipUpstream+1] - Old_depth[ipUpstream]));
00672
00673
00674
00675 const double STEP_FACTOR=0.05;
00676
00677 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00678 {
00679 for( ion=0; ion<nelem+2; ++ion )
00680 {
00681 double f1;
00682 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
00683 {
00684 yion =
00685 Old_xIonDense[ipyUpstream][nelem][ion]/Old_density[ipyUpstream] +
00686 (Old_xIonDense[ipyUpstream+1][nelem][ion]/Old_density[ipyUpstream+1] -
00687 Old_xIonDense[ipyUpstream][nelem][ion]/Old_density[ipyUpstream])*
00688 ynextfrac;
00689 }
00690 else
00691 {
00692 long ipBound = ipyUpstream;
00693 if (-1 == ipBound)
00694 ipBound = 0;
00695 yion = Old_xIonDense[ipBound][nelem][ion]/Old_density[ipBound];
00696 }
00697 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
00698 {
00699 hion =
00700 Old_xIonDense[iphUpstream][nelem][ion]/Old_density[iphUpstream] +
00701 (Old_xIonDense[iphUpstream+1][nelem][ion]/Old_density[iphUpstream+1] -
00702 Old_xIonDense[iphUpstream][nelem][ion]/Old_density[iphUpstream])*
00703 hnextfrac;
00704 }
00705 else
00706 {
00707 long ipBound = iphUpstream;
00708 if (-1 == ipBound)
00709 ipBound = 0;
00710 hion = Old_xIonDense[ipBound][nelem][ion]/Old_density[ipBound];
00711 }
00712
00713
00714 f1 = fabs(yion - UpstreamIon[nelem][ion] );
00715 if( f1 > SMALLFLOAT )
00716 {
00717 dynamics.dRad = MIN2(dynamics.dRad,STEP_FACTOR*fabs(Dyn_dr) *
00718
00719 MAX2(yion + UpstreamIon[nelem][ion],1e-6 ) / f1);
00720 }
00721
00722
00723
00724
00725
00726
00727 dynamics.discretization_error += POW2(yion-2.*hion+UpstreamIon[nelem][ion]);
00728 dynamics.error_scale2 += POW2(UpstreamIon[nelem][ion]-yion);
00729 }
00730 }
00731
00732 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00733 {
00734 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00735 {
00736 if( dense.lgElmtOn[nelem] )
00737 {
00738 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
00739 {
00740 double f1;
00741 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 &&
00742 (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
00743 {
00744 yiso =
00745 Old_StatesElem[ipyUpstream][nelem][nelem-ipISO][level]/Old_density[ipyUpstream] +
00746 (Old_StatesElem[ipyUpstream+1][nelem][nelem-ipISO][level]/Old_density[ipyUpstream+1] -
00747 Old_StatesElem[ipyUpstream][nelem][nelem-ipISO][level]/Old_density[ipyUpstream])*
00748 ynextfrac;
00749 }
00750 else
00751 {
00752 long ipBound = ipyUpstream;
00753 if (-1 == ipBound)
00754 ipBound = 0;
00755 yiso = Old_StatesElem[ipBound][nelem][nelem-ipISO][level]/Old_density[ipBound];
00756 }
00757 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 &&
00758 (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
00759 {
00760 hiso =
00761 Old_StatesElem[iphUpstream][nelem][nelem-ipISO][level]/Old_density[iphUpstream] +
00762 (Old_StatesElem[iphUpstream+1][nelem][nelem-ipISO][level]/Old_density[iphUpstream+1] -
00763 Old_StatesElem[iphUpstream][nelem][nelem-ipISO][level]/Old_density[iphUpstream])*
00764 hnextfrac;
00765 }
00766 else
00767 {
00768 long ipBound = iphUpstream;
00769 if (-1 == ipBound)
00770 ipBound = 0;
00771 hiso = Old_StatesElem[ipBound][nelem][nelem-ipISO][level]/Old_density[ipBound];
00772 }
00773
00774
00775 f1 = fabs(yiso - UpstreamStatesElem[nelem][nelem-ipISO][level] );
00776 if( f1 > SMALLFLOAT )
00777 {
00778 dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
00779
00780 MAX2(yiso + UpstreamStatesElem[nelem][nelem-ipISO][level],1e-6 ) / f1);
00781 }
00782
00783
00784
00785
00786
00787 dynamics.discretization_error += POW2(yiso-2.*hiso+UpstreamStatesElem[nelem][nelem-ipISO][level]);
00788 dynamics.error_scale2 += POW2(UpstreamStatesElem[nelem][nelem-ipISO][level]);
00789 }
00790 }
00791 }
00792 }
00793
00794 for(mol=0; mol < mole_global.num_calc; mol++)
00795 {
00796 double f1;
00797 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
00798 {
00799 ymol =
00800 Old_molecules[ipyUpstream][mol]/Old_density[ipyUpstream] +
00801 (Old_molecules[ipyUpstream+1][mol]/Old_density[ipyUpstream+1] -
00802 Old_molecules[ipyUpstream][mol]/Old_density[ipyUpstream])*
00803 ynextfrac;
00804 }
00805 else
00806 {
00807 long ipBound = ipyUpstream;
00808 if (-1 == ipBound)
00809 ipBound = 0;
00810 ymol = Old_molecules[ipBound][mol]/Old_density[ipBound];
00811 }
00812 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
00813 {
00814 hmol =
00815 Old_molecules[iphUpstream][mol]/Old_density[iphUpstream] +
00816 (Old_molecules[iphUpstream+1][mol]/Old_density[iphUpstream+1] -
00817 Old_molecules[iphUpstream][mol]/Old_density[iphUpstream])*
00818 hnextfrac;
00819 }
00820 else
00821 {
00822 long ipBound = iphUpstream;
00823 if (-1 == ipBound)
00824 ipBound = 0;
00825 hmol = Old_molecules[ipBound][mol]/Old_density[ipBound];
00826 }
00827
00828
00829 f1 = fabs(ymol - Upstream_molecules[mol] );
00830 if( f1 > SMALLFLOAT )
00831 {
00832 dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
00833
00834 MAX2(ymol + Upstream_molecules[mol],1e-6 ) / f1 );
00835 }
00836
00837
00838
00839
00840 dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_molecules[mol]);
00841 dynamics.error_scale2 += POW2(Upstream_molecules[mol]-ymol);
00842 }
00843
00844 if( dynamics.lgTracePrint )
00845 {
00846 fprintf(ioQQQ," DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
00847 nzone,dynamics.Rate , dynamics.Source[0][0] , dilution*scalingDensity() );
00848 }
00849 return;
00850 }
00851
00852
00853 void DynaEndZone(void)
00854 {
00855 DEBUG_ENTRY( "DynaEndZone()" );
00856
00857
00858
00859 dynamics.DivergePresInteg += wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad));
00860
00861 if(dynamics.lgTracePrint)
00862 fprintf(ioQQQ,"Check dp: %g %g mom %g %g mas %g\n",
00863 wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad)),
00864 2*wind.windv*DynaFlux(radius.depth)*radius.drad/(1e16-radius.depth),
00865 wind.windv*DynaFlux(radius.depth),
00866 wind.windv*DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth),
00867 DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth));
00868 return;
00869 }
00870
00871
00872
00873
00874 void DynaIterEnd(void)
00875 {
00876
00877
00878
00879
00880 long int i;
00881 static long int nTime_dt_array_element = 0;
00882
00883 DEBUG_ENTRY( "DynaIterEnd()" );
00884
00885
00886
00887
00888
00889
00890
00891 if( iteration == dynamics.n_initial_relax+1)
00892 {
00893 fprintf(ioQQQ,"DYNAMICS DynaIterEnd sets stop radius to %.2e after "
00894 "dynamics.n_initial_relax=%li iterations.\n",
00895 radius.depth,
00896 dynamics.n_initial_relax);
00897 for( i=iteration-1; i<iterations.iter_malloc; ++i )
00898
00899
00900
00901
00902 radius.StopThickness[i] = radius.depth;
00903 }
00904
00905 dynamics.DivergePresInteg = 0.;
00906
00907
00908
00909
00910
00911
00912 if( !dynamics.lgTimeDependentStatic )
00913 {
00914 if(iteration == dynamics.n_initial_relax+1 )
00915 {
00916
00917
00918
00919
00920
00921
00922
00923 if( dynamics.AdvecLengthInit> 0. )
00924 {
00925 Dyn_dr = dynamics.AdvecLengthInit;
00926 }
00927 else
00928 {
00929
00930 Dyn_dr = -dynamics.AdvecLengthInit*radius.depth;
00931 }
00932
00933 if (wind.windv0 > 0)
00934 Dyn_dr = -Dyn_dr;
00935
00936 if( dynamics.lgTracePrint )
00937 {
00938 fprintf(ioQQQ," DynaIterEnd, dr=%.2e \n",
00939 Dyn_dr );
00940 }
00941 }
00942 else if(iteration > dynamics.n_initial_relax+1 )
00943 {
00944
00945 DynaNewStep();
00946 }
00947 }
00948 else
00949 {
00950
00951 static double HeatInitial=-1. , HeatRadiated=-1. ,
00952 DensityInitial=-1;
00953
00954
00955 Dyn_dr = 0.;
00956 fprintf(ioQQQ,
00957 "DEBUG times enter dynamics.timestep %.2e elapsed_time %.2e iteration %li relax %li \n",
00958 dynamics.timestep ,
00959 dynamics.time_elapsed,
00960 iteration , dynamics.n_initial_relax);
00961 if( iteration > dynamics.n_initial_relax )
00962 {
00963
00964 DynaNewStep();
00965
00966
00967
00968 if( dynamics.lg_coronal_time_init )
00969 {
00970 thermal.lgTemperatureConstant = false;
00971 thermal.ConstTemp = 0.;
00972 }
00973
00974
00975
00976
00977 dynamics.time_elapsed += dynamics.timestep;
00978
00979 if( dynamics.timestep_stop > 0 && dynamics.time_elapsed > dynamics.timestep_stop )
00980 {
00981 dynamics.lgStatic_completed = true;
00982 }
00983
00984
00985 if( (phycon.te < StopCalc.TempLoStopIteration) ||
00986 (phycon.te > StopCalc.TempHiStopIteration ) )
00987 dynamics.lgStatic_completed = true;
00988
00989
00990
00991 HeatRadiated += (thermal.ctot-dynamics.Cool()) * dynamics.timestep *
00992 (DensityInitial / scalingDensity());
00993 }
00994 else
00995 {
00996
00997 HeatInitial = 1.5 * pressure.PresGasCurr;
00998 HeatRadiated = 0.;
00999 DensityInitial = scalingDensity();
01000 fprintf(ioQQQ,"DEBUG relaxing times requested %li this is step %li\n",
01001 dynamics.n_initial_relax , iteration);
01002 }
01003 fprintf(ioQQQ,"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n",
01004 HeatInitial , HeatRadiated );
01005
01006
01007
01008 if( dynamics.time_elapsed > time_elapsed_time[nTime_dt_array_element] )
01009 {
01010
01011
01012 ++nTime_dt_array_element;
01013
01014
01015 ASSERT( nTime_dt_array_element < nTime_flux );
01016
01017
01018 if( lgtime_Recom[nTime_dt_array_element] )
01019 {
01020 fprintf(ioQQQ,"DEBUG dynamics turn on recombination logic\n");
01021 dynamics.lgRecom = true;
01022
01023
01024
01025
01026 radius.sdrmax = radius.dr_max_last_iter/3.;
01027 radius.lgSdrmaxRel = false;
01028 }
01029
01030 if( lgtime_dt_specified )
01031 {
01032
01033 fprintf(ioQQQ,"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element,
01034 dynamics.timestep);
01035 dynamics.timestep = time_dt[nTime_dt_array_element];
01036
01037 if( time_dt_scale_factor[nTime_dt_array_element] > 0. )
01038 dynamics.timestep_factor = time_dt_scale_factor[nTime_dt_array_element];
01039 }
01040 }
01041 else if( lgtime_dt_specified )
01042 {
01043
01044 dynamics.timestep *= dynamics.timestep_factor;
01045 fprintf(ioQQQ,"DEBUG lgtimes increment Timeby dynamics.timestep_factor to %li %.2e\n" ,
01046 nTime_dt_array_element,
01047 dynamics.timestep );
01048 if(dynamics.time_elapsed+dynamics.timestep > time_elapsed_time[nTime_dt_array_element] )
01049 {
01050 fprintf(ioQQQ,"DEBUG lgtimes but reset to %.2e\n" ,dynamics.timestep );
01051 dynamics.timestep = 1.0001*(time_elapsed_time[nTime_dt_array_element]-dynamics.time_elapsed);
01052 }
01053 }
01054 else
01055 {
01056
01057 dynamics.timestep = timestep_next();
01058 }
01059
01060 if( cosmology.lgDo )
01061 {
01062 cosmology.redshift_step = -(1.f + cosmology.redshift_current ) *
01063 GetHubbleFactor( cosmology.redshift_current ) * (realnum)dynamics.timestep;
01064 cosmology.redshift_current += cosmology.redshift_step;
01065 }
01066
01067 fprintf(ioQQQ,"DEBUG times exit dynamics.timestep %.2e elapsed_time %.2e scale %.2e ",
01068 dynamics.timestep ,
01069 dynamics.time_elapsed,
01070 rfield.time_continuum_scale);
01071
01072 if( cosmology.lgDo )
01073 {
01074 fprintf(ioQQQ,"redshift %.3e ", cosmology.redshift_current );
01075 }
01076
01077 fprintf(ioQQQ,"\n" );
01078 }
01079
01080
01081 ASSERT( (iteration < dynamics.n_initial_relax+1) ||
01082 Dyn_dr != 0. || (Dyn_dr == 0. && wind.lgStatic()) );
01083
01084
01085 ipUpstream = iphUpstream = ipyUpstream = -1;
01086 dynamics.discretization_error = 0.;
01087 dynamics.error_scale2 = 0.;
01088
01089
01090 DynaSaveLast();
01091 return;
01092 }
01093
01094
01095 STATIC void DynaNewStep(void)
01096 {
01097 long int ilast = 0,
01098 i,
01099 nelem,
01100 ion,
01101 mol;
01102
01103 double frac_next=-BIGFLOAT,
01104 Oldi_density,
01105 Oldi_ion,
01106 Oldi_iso,
01107 Oldi_mol;
01108
01109 DEBUG_ENTRY( "DynaNewStep()" );
01110
01111
01112 dynamics.convergence_error = 0;
01113 dynamics.error_scale1 = 0.;
01114
01115 ASSERT( nzone < struc.nzlim);
01116 for(i=0;i<nzone;++i)
01117 {
01118
01119 while( (Old_depth[ilast] < struc.depth[i] ) &&
01120 ( ilast < nOld_zone-1 ) )
01121 {
01122 ++ilast;
01123 }
01124 ASSERT( ilast <= nOld_zone-1 );
01125
01126 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01127 {
01128 frac_next = ( struc.depth[i] - Old_depth[ilast])/
01129 (Old_depth[ilast+1] - Old_depth[ilast]);
01130 Oldi_density = Old_density[ilast] +
01131 (Old_density[ilast+1] - Old_density[ilast])*
01132 frac_next;
01133 }
01134 else
01135 {
01136 Oldi_density = Old_density[ilast];
01137 }
01138
01139
01140 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01141 {
01142 for( ion=0; ion<nelem+2; ++ion )
01143 {
01144 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01145 {
01146 Oldi_ion = (Old_xIonDense[ilast][nelem][ion] +
01147 (Old_xIonDense[ilast+1][nelem][ion]-Old_xIonDense[ilast][nelem][ion])*
01148 frac_next);
01149 }
01150 else
01151 {
01152 Oldi_ion = Old_xIonDense[ilast][nelem][ion];
01153 }
01154 dynamics.convergence_error += POW2((double)Oldi_ion/Oldi_density-struc.xIonDense[nelem][ion][i]/scalingZoneDensity(i)) ;
01155
01156
01157
01158
01159 dynamics.error_scale1 += POW2((double)struc.xIonDense[nelem][ion][i]/scalingZoneDensity(i));
01160 }
01161 }
01162 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01163 {
01164 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01165 {
01166 if( dense.lgElmtOn[nelem] )
01167 {
01168 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; ++level )
01169 {
01170 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01171 {
01172 Oldi_iso = (Old_StatesElem[ilast][nelem][nelem-ipISO][level] +
01173 (Old_StatesElem[ilast+1][nelem][nelem-ipISO][level]-Old_StatesElem[ilast][nelem][nelem-ipISO][level])*
01174 frac_next);
01175 }
01176 else
01177 {
01178 Oldi_iso = Old_StatesElem[ilast][nelem][nelem-ipISO][level];
01179 }
01180 dynamics.convergence_error += POW2(Oldi_iso/Oldi_density-struc.StatesElem[nelem][nelem-ipISO][level][i]/struc.hden[i]) ;
01181
01182
01183 dynamics.error_scale1 += POW2(struc.StatesElem[nelem][nelem-ipISO][level][i]/struc.hden[i]);
01184 }
01185 }
01186 }
01187 }
01188
01189 for( mol=0; mol < mole_global.num_calc; mol++)
01190 {
01191 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01192 {
01193 Oldi_mol = (Old_molecules[ilast][mol] +
01194 (Old_molecules[ilast+1][mol]-Old_molecules[ilast][mol])*
01195 frac_next);
01196 }
01197 else
01198 {
01199 Oldi_mol = Old_molecules[ilast][mol];
01200 }
01201 dynamics.convergence_error += POW2((double)Oldi_mol/Oldi_density-struc.molecules[mol][i]/scalingZoneDensity(i)) ;
01202
01203
01204
01205 dynamics.error_scale1 += POW2((double)struc.molecules[mol][i]/scalingZoneDensity(i));
01206 }
01207 }
01208
01209
01210
01211
01212
01213
01214
01215 fprintf(ioQQQ,"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
01216 Dyn_dr, dynamics.convergence_error , dynamics.discretization_error ,
01217 dynamics.error_scale1 , dynamics.error_scale2
01218 );
01219
01220
01221 if( dynamics.convergence_error < dynamics.convergence_tolerance*dynamics.discretization_error )
01222 Dyn_dr /= 1.5;
01223 return;
01224 }
01225
01226
01227 STATIC void DynaSaveLast(void)
01228 {
01229 long int i,
01230 ion,
01231 nelem,
01232 mol;
01233
01234 DEBUG_ENTRY( "DynaSaveLast()" );
01235
01236
01237 nOld_zone = nzone;
01238 dynamics.oldFullDepth = struc.depth[nzone-1];
01239 ASSERT( nzone < struc.nzlim );
01240 for( i=0; i<nzone; ++i )
01241 {
01242 Old_histr[i] = struc.histr[i];
01243 Old_depth[i] = struc.depth[i];
01244 Old_xLyman_depth[i] = struc.xLyman_depth[i];
01245
01246 Old_hiistr[i] = struc.hiistr[i];
01247
01248 Old_pressure[i] = struc.pressure[i];
01249
01250 Old_ednstr[i] = struc.ednstr[i];
01251
01252 Old_EnthalpyDensity[i] = EnthalpyDensity[i];
01253 Old_density[i] = scalingZoneDensity(i);
01254 Old_DenMass[i] = struc.DenMass[i];
01255
01256 for(mol=0;mol<mole_global.num_calc;mol++)
01257 {
01258 Old_molecules[i][mol] = struc.molecules[mol][i];
01259 }
01260
01261 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01262 {
01263 Old_gas_phase[i][nelem] = dense.gas_phase[nelem];
01264 for( ion=0; ion<nelem+2; ++ion )
01265 {
01266 Old_xIonDense[i][nelem][ion] = struc.xIonDense[nelem][ion][i];
01267 }
01268 }
01269 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01270 {
01271 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01272 {
01273 if( dense.lgElmtOn[nelem] )
01274 {
01275 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
01276 {
01277 Old_StatesElem[i][nelem][nelem-ipISO][level] = struc.StatesElem[nelem][nelem-ipISO][level][i];
01278 ASSERT( !isnan( Old_StatesElem[i][nelem][nelem-ipISO][level] ) );
01279 }
01280 }
01281 }
01282 }
01283 }
01284 return;
01285 }
01286
01287 realnum DynaFlux(double depth)
01288
01289 {
01290 realnum flux;
01291
01292 DEBUG_ENTRY( "DynaFlux()" );
01293
01294 if(dynamics.FluxIndex == 0)
01295 {
01296 flux = (realnum)dynamics.FluxScale;
01297 }
01298 else
01299 {
01300 flux = (realnum)(dynamics.FluxScale*pow(fabs(depth-dynamics.FluxCenter),dynamics.FluxIndex));
01301 if(depth < dynamics.FluxCenter)
01302 flux = -flux;
01303 }
01304 if(dynamics.lgFluxDScale)
01305 {
01306
01307
01308 flux *= dense.xMassDensity0;
01309 }
01310 return flux;
01311 }
01312
01313
01314
01315
01316 void DynaZero( void )
01317 {
01318 int ipISO;
01319
01320 DEBUG_ENTRY( "DynaZero()" );
01321
01322
01323 nOld_zone = 0;
01324
01325
01326 dynamics.lgAdvection = false;
01327
01328 AdvecSpecificEnthalpy = 0.;
01329 dynamics.Cool_r = 0.;
01330 dynamics.Heat_v = 0.;
01331 dynamics.dHeatdT = 0.;
01332 dynamics.HeatMax = 0.;
01333 dynamics.CoolMax = 0.;
01334 dynamics.Rate = 0.;
01335
01336
01337 dynamics.lgRecom = false;
01338
01339
01340 dynamics.lgEquilibrium = false;
01341
01342
01343 dynamics.lgStatic_completed = false;
01344
01345
01346 dynamics.lgTimeDependentStatic = false;
01347 dynamics.timestep_init = -1.;
01348
01349 dynamics.timestep_factor = 1.2;
01350 dynamics.time_elapsed = 0.;
01351
01352
01353
01354 dynamics.n_initial_relax = 2;
01355
01356
01357
01358 dynamics.AdvecLengthInit = -0.1;
01359
01360
01361 dynamics.convergence_tolerance = 0.1;
01362
01363
01364 dynamics.lgSetPresMode = false;
01365
01366
01367 dynamics.FluxScale = 0.;
01368 dynamics.lgFluxDScale = false;
01369 dynamics.FluxCenter = 0.;
01370 dynamics.FluxIndex = 0.;
01371 dynamics.dRad = BIGFLOAT;
01372
01373 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01374 {
01375
01376
01377
01378 dynamics.lgISO[ipISO] = true;
01379 }
01380
01381 dynamics.lgMETALS = true;
01382
01383 dynamics.lgCoolHeat = true;
01384 dynamics.DivergePresInteg = 0.;
01385
01386 dynamics.discretization_error = 0.;
01387 dynamics.error_scale2 = 0.;
01388 return;
01389 }
01390
01391
01392
01393
01394
01395 void DynaCreateArrays( void )
01396 {
01397 long int nelem,
01398 ns,
01399 i,
01400 ion,
01401 mol;
01402
01403 DEBUG_ENTRY( "DynaCreateArrays()" );
01404
01405 if (mole_global.num_calc != 0)
01406 {
01407 Upstream_molecules = (double*)MALLOC((size_t)mole_global.num_calc*sizeof(double) );
01408
01409 dynamics.molecules = (double*)MALLOC((size_t)mole_global.num_calc*sizeof(double) );
01410 }
01411 else
01412 {
01413 Upstream_molecules = dynamics.molecules = NULL;
01414 }
01415 UpstreamElem = (double*)MALLOC((size_t)LIMELM*sizeof(double) );
01416
01417 dynamics.Source = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
01418 UpstreamIon = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
01419 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01420 {
01421 dynamics.Source[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
01422 UpstreamIon[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
01423 for( ion=0; ion<nelem+2; ++ion )
01424 {
01425 dynamics.Source[nelem][ion] = 0.;
01426 }
01427 }
01428
01429 UpstreamStatesElem = ((double***)MALLOC( (size_t)LIMELM*sizeof(double **) ));
01430 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01431 {
01432 if( dense.lgElmtOn[nelem] )
01433 {
01434 UpstreamStatesElem[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
01435 for( long ion=0; ion<nelem+1; ion++ )
01436 {
01437 long ipISO = nelem-ion;
01438 if( ipISO < NISO )
01439 {
01440 UpstreamStatesElem[nelem][nelem-ipISO] = (double*)MALLOC(sizeof(double)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
01441 }
01442 else
01443 {
01444 fixit();
01445 UpstreamStatesElem[nelem][nelem-ipISO] = NULL;
01446 }
01447 }
01448 }
01449 }
01450
01451
01452 dynamics.StatesElem = ((double***)MALLOC( (size_t)LIMELM*sizeof(double **) ));
01453 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01454 {
01455 if( dense.lgElmtOn[nelem] )
01456 {
01457 dynamics.StatesElem[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
01458 for( long ion=0; ion<nelem+1; ion++ )
01459 {
01460 long ipISO = nelem-ion;
01461 if( ipISO < NISO )
01462 {
01463 dynamics.StatesElem[nelem][nelem-ipISO] = (double*)MALLOC(sizeof(double)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
01464 }
01465 else
01466 {
01467 fixit();
01468 dynamics.StatesElem[nelem][nelem-ipISO] = NULL;
01469 }
01470 }
01471 }
01472 }
01473
01474 dynamics.Rate = 0.;
01475
01476 Old_EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01477
01478 Old_ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01479
01480 EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01481
01482 Old_DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01483
01484 Old_density = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01485
01486 Old_pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01487
01488 Old_histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01489
01490 Old_hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01491
01492 Old_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01493
01494 Old_xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
01495
01496 Old_xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) );
01497
01498 Old_StatesElem = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(struc.nzlim) );
01499
01500 Old_gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
01501
01502 Old_molecules = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
01503
01504
01505 for( ns=0; ns < struc.nzlim; ++ns )
01506 {
01507 Old_xIonDense[ns] =
01508 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM) );
01509
01510 Old_StatesElem[ns] =
01511 (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(LIMELM) );
01512
01513 Old_gas_phase[ns] =
01514 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM) );
01515
01516 if (mole_global.num_calc != 0)
01517 {
01518 Old_molecules[ns] =
01519 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(mole_global.num_calc) );
01520 }
01521 else
01522 {
01523 Old_molecules[ns] = NULL;
01524 }
01525
01526 for( nelem=0; nelem<LIMELM; ++nelem )
01527 {
01528 Old_xIonDense[ns][nelem] =
01529 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+1) );
01530 }
01531
01532 for( nelem=0; nelem< LIMELM; ++nelem )
01533 {
01534 if( dense.lgElmtOn[nelem] )
01535 {
01536 Old_StatesElem[ns][nelem] =
01537 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+1) );
01538 for( ion=0; ion<nelem+1; ion++ )
01539 {
01540 long ipISO = nelem-ion;
01541 if( ipISO < NISO )
01542 {
01543 Old_StatesElem[ns][nelem][ion] =
01544 (realnum*)MALLOC(sizeof(realnum)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
01545 }
01546 else
01547 {
01548 fixit();
01549 Old_StatesElem[ns][nelem][ion] = NULL;
01550 }
01551 }
01552 }
01553 }
01554 }
01555
01556 for( i=0; i < struc.nzlim; i++ )
01557 {
01558
01559 Old_histr[i] = 0.;
01560 Old_xLyman_depth[i] = 0.;
01561 Old_depth[i] = 0.;
01562 dynamics.oldFullDepth = 0.;
01563
01564 Old_hiistr[i] = 0.;
01565
01566 Old_pressure[i] = 0.;
01567
01568 Old_ednstr[i] = 0.;
01569 Old_density[i] = 0.;
01570 Old_DenMass[i] = 0.;
01571 Old_EnthalpyDensity[i] = 0.;
01572 for( nelem=0; nelem<LIMELM; ++nelem )
01573 {
01574 for( ion=0; ion<LIMELM+1; ++ion )
01575 {
01576 Old_xIonDense[i][nelem][ion] = 0.;
01577 }
01578 }
01579 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01580 {
01581 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01582 {
01583 if( dense.lgElmtOn[nelem] )
01584 {
01585 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
01586 {
01587 Old_StatesElem[i][nelem][nelem-ipISO][level] = 0.;
01588 }
01589 }
01590 }
01591 }
01592
01593 for(mol=0;mol<mole_global.num_calc;mol++)
01594 {
01595 Old_molecules[i][mol] = 0.;
01596 }
01597 }
01598 return;
01599 }
01600
01601
01602
01603
01604 STATIC void advection_set_default( bool lgWind )
01605 {
01606
01607 DEBUG_ENTRY( "advection_set_default()" );
01608
01609
01610 dynamics.lgAdvection = true;
01611
01612
01613
01614 thermal.lgPredNextTe = false;
01615
01616
01617
01618
01619
01620
01621 pressure.lgPres_radiation_ON = true;
01622 pressure.lgPres_magnetic_ON = true;
01623 pressure.lgPres_ram_ON = true;
01624
01625
01626 if( lgWind )
01627 {
01628
01629 conv.EdenErrorAllowed = 1e-3;
01630
01631
01632
01633
01634
01635 conv.HeatCoolRelErrorAllowed = 3e-4f;
01636 conv.PressureErrorAllowed = 1e-3f;
01637
01638 if( cosmology.lgDo )
01639 {
01640 conv.EdenErrorAllowed = 1e-5;
01641 conv.PressureErrorAllowed = 1e-5f;
01642 }
01643 }
01644 return;
01645 }
01646
01647
01648
01649 void ParseDynaTime( Parser &p )
01650 {
01651 DEBUG_ENTRY( "ParseDynaTime()" );
01652
01653
01654 dynamics.lgTimeDependentStatic = true;
01655
01656 dynamics.timestep_init = p.getNumberCheckAlwaysLogLim("dynamics.timestep",30.);
01657
01658 dynamics.timestep = dynamics.timestep_init;
01659 if( p.nMatch( "TRAC" ) )
01660 dynamics.lgTracePrint = true;
01661
01662
01663 dynamics.timestep_stop = p.getNumberDefaultAlwaysLog("stop time", -1.);
01664
01665
01666 advection_set_default(false);
01667
01668 wind.windv0 = 0.;
01669 wind.setStatic();
01670 wind.windv = wind.windv0;
01671
01672
01673 time_elapsed_time = (double*)MALLOC((size_t)NTIME*sizeof(double));
01674 time_flux_ratio = (double*)MALLOC((size_t)NTIME*sizeof(double));
01675 time_dt = (double*)MALLOC((size_t)NTIME*sizeof(double));
01676 time_dt_scale_factor = (double*)MALLOC((size_t)NTIME*sizeof(double));
01677 lgtime_Recom = (int*)MALLOC((size_t)NTIME*sizeof(int));
01678
01679
01680 nTime_flux = 0;
01681
01682
01683 p.getline();
01684 if( p.m_lgEOF )
01685 {
01686 fprintf( ioQQQ,
01687 " Hit EOF while reading time-continuum list; use END to end list.\n" );
01688 cdEXIT(EXIT_FAILURE);
01689 }
01690
01691
01692
01693 lgtime_dt_specified = true;
01694
01695 while( p.strcmp("END") != 0 )
01696 {
01697 if( nTime_flux >= NTIME )
01698 {
01699 fprintf( ioQQQ,
01700 " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n",
01701 NTIME );
01702 cdEXIT(EXIT_FAILURE);
01703 }
01704
01705 if( p.nMatch("CYCLE") )
01706 {
01707 double period = p.getNumberCheckAlwaysLog("log time");
01708 ASSERT( period > time_elapsed_time[nTime_flux-1] );
01709 long pointsPerPeriod = nTime_flux;
01710 while( nTime_flux < NTIME - 1 )
01711 {
01712 time_elapsed_time[nTime_flux] = period + time_elapsed_time[nTime_flux-pointsPerPeriod];
01713 time_flux_ratio[nTime_flux] = time_flux_ratio[nTime_flux-pointsPerPeriod];
01714 time_dt[nTime_flux] = time_dt[nTime_flux-pointsPerPeriod];
01715 time_dt_scale_factor[nTime_flux] = time_dt_scale_factor[nTime_flux-pointsPerPeriod];
01716 nTime_flux++;
01717 }
01718
01719 fprintf( ioQQQ, " Adding cycles with period = %e s.\n", period );
01720
01721
01722 p.getline();
01723 if( p.m_lgEOF )
01724 {
01725 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
01726 cdEXIT(EXIT_FAILURE);
01727 }
01728 continue;
01729 }
01730
01731 time_elapsed_time[nTime_flux] = p.getNumberCheckAlwaysLog("log time");
01732 if( nTime_flux >= 1 )
01733 ASSERT( time_elapsed_time[nTime_flux] > time_elapsed_time[nTime_flux-1] );
01734 time_flux_ratio[nTime_flux] = p.getNumberCheckAlwaysLog("log flux ratio");
01735
01736
01737
01738 time_dt[nTime_flux] = p.getNumberDefaultAlwaysLog("log time step",-1.);
01739
01740
01741 if( time_dt[nTime_flux] < 0.0 )
01742 lgtime_dt_specified = false;
01743
01744
01745 time_dt_scale_factor[nTime_flux] = p.getNumberDefaultAlwaysLog(
01746 "scale factor to increase time",-1.);
01747
01748
01749 if( p.nMatch("RECOMBIN") )
01750 {
01751
01752
01753 lgtime_Recom[nTime_flux] = true;
01754 }
01755 else
01756 {
01757 lgtime_Recom[nTime_flux] = false;
01758 }
01759
01760
01761 ++nTime_flux;
01762
01763
01764 p.getline();
01765 if( p.m_lgEOF )
01766 {
01767 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
01768 cdEXIT(EXIT_FAILURE);
01769 }
01770
01771 }
01772
01773 if( nTime_flux < 2 )
01774 {
01775 fprintf( ioQQQ, " At least two instances of time must be specified. There is an implicit instance at t=0.\n" \
01776 " The user must specify at least one additional time. Sorry.\n" );
01777 cdEXIT(EXIT_FAILURE);
01778 }
01779
01780 for( long i=0; i < nTime_flux; i++ )
01781 {
01782 fprintf( ioQQQ, "DEBUG time dep %.2e %.2e %.2e %.2e\n",
01783 time_elapsed_time[i],
01784 time_flux_ratio[i] ,
01785 time_dt[i],
01786 time_dt_scale_factor[i]);
01787 }
01788 fprintf( ioQQQ, "\n" );
01789 return;
01790 }
01791
01792
01793 void ParseDynaWind( Parser &p)
01794 {
01795 int iVelocity_Type;
01796 bool lgModeSet=false;
01797
01798
01799 double dfdr=-BIGDOUBLE;
01800
01801 DEBUG_ENTRY( "ParseDynaWind()" );
01802
01803 if( p.nMatch( "TRAC" ) )
01804 dynamics.lgTracePrint = true;
01805
01806
01807
01808
01809
01810 iVelocity_Type = 0;
01811
01812
01813 if( p.nMatch( "VELO" ) )
01814 {
01815 wind.windv0 = (realnum)(p.getNumberPlain("velocity")*1e5);
01816 wind.windv = wind.windv0;
01817 wind.setDefault();
01818 iVelocity_Type = 1;
01819 }
01820
01821 if( p.nMatch( "BALL" ) )
01822 {
01823 wind.setBallistic();
01824 lgModeSet = true;
01825 }
01826
01827 if( p.nMatch( "STAT" ) )
01828 {
01829 wind.windv0 = 0.;
01830 wind.setStatic();
01831 lgModeSet = true;
01832 iVelocity_Type = 1;
01833 }
01834
01835 if ( 1 == iVelocity_Type && !lgModeSet)
01836 {
01837 if (wind.windv0 > 0.)
01838 {
01839 fprintf(ioQQQ,"Warning, BALListic option needed to switch off pressure gradient terms\n");
01840 }
01841 else if (wind.windv0 == 0.)
01842 {
01843 fprintf(ioQQQ,"Warning, STATic option needed for zero speed solutions\n");
01844 }
01845 }
01846
01847 if( p.nMatch("DFDR") )
01848 {
01849
01850 dfdr = p.getNumberPlain("flux gradient");
01851 iVelocity_Type = 2;
01852 }
01853
01854
01855 if( p.nMatch("CENT") )
01856 {
01857
01858 dynamics.FluxCenter = p.getNumberPlain(
01859 "centre of mass flux distribution");
01860 }
01861
01862
01863 if( p.nMatch("INDE") )
01864 {
01865
01866 dynamics.FluxIndex = p.getNumberPlain(
01867 "power law index of mass flux distribution");
01868 }
01869
01870
01871 if(iVelocity_Type == 1)
01872 {
01873
01874 if(dynamics.FluxIndex == 0)
01875 {
01876 dynamics.FluxScale = wind.windv0;
01877 dynamics.lgFluxDScale = true;
01878
01879
01880
01881 dynamics.FluxCenter = -1.;
01882 }
01883 else
01884 {
01887
01888 dynamics.FluxScale = wind.windv0*
01889 pow(fabs(dynamics.FluxCenter),-dynamics.FluxIndex);
01890
01891 dynamics.lgFluxDScale = true;
01892 if(dynamics.FluxCenter > 0)
01893 {
01894 dynamics.FluxScale = -dynamics.FluxScale;
01895 }
01896 }
01897 }
01898
01899 else if(iVelocity_Type == 2)
01900 {
01901 if(dynamics.FluxIndex == 0)
01902 {
01903 fprintf(ioQQQ,"Can't specify gradient when flux is constant!\n");
01904
01905 cdEXIT(EXIT_FAILURE);
01906 }
01909
01910
01911 dynamics.FluxScale = dfdr/dynamics.FluxIndex*
01912 pow(fabs(dynamics.FluxCenter),1.-dynamics.FluxIndex);
01913 if(dynamics.FluxCenter > 0)
01914 {
01915 dynamics.FluxScale = -dynamics.FluxScale;
01916 }
01917 dynamics.lgFluxDScale = false;
01918
01919
01920
01921 wind.windv0 = -0.01f;
01922 wind.setDefault();
01923 }
01924 else
01925 {
01926
01927
01928
01929 wind.windv0 = (realnum)(p.getNumberCheck("wind velocity")*1e5);
01930 if (wind.windv0 < 0.)
01931 {
01932 wind.setDefault();
01933 }
01934 else if (wind.windv0 > 0.)
01935 {
01936 wind.setBallistic();
01937 }
01938 else
01939 {
01940 wind.setStatic();
01941 }
01942
01943 dynamics.FluxScale = wind.windv0;
01944 dynamics.FluxIndex = 0.;
01945 dynamics.lgFluxDScale = true;
01946
01947
01948
01949 dynamics.FluxCenter = -1.;
01950 }
01951
01952 wind.windv = wind.windv0;
01953
01954 # ifdef FOO
01955 fprintf(ioQQQ,"Scale %g (*%c) Index %g Center %g\n",
01956 dynamics.FluxScale,(dynamics.lgFluxDScale)?'D':'1',
01957 dynamics.FluxIndex,dynamics.FluxCenter);
01958 # endif
01959
01960
01961 if( p.nMatch( "ADVE" ) )
01962 {
01963
01964 advection_set_default(true);
01965 strcpy( dense.chDenseLaw, "DYNA" );
01966 }
01967
01968 else
01969 {
01970
01971 if( wind.windv0 <= 1.e6 )
01972 {
01973
01974 fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
01975 wind.lgWindOK = false;
01976 }
01977
01978
01979 wind.comass = (realnum)p.getNumberDefault("central object mass",1.);
01980
01981
01982
01983 wind.lgDisk = false;
01984 if( p.nMatch( "DISK") )
01985 wind.lgDisk = true;
01986
01987 strcpy( dense.chDenseLaw, "WIND" );
01988 }
01989
01990
01991
01992 if( p.nMatch("NO CO") )
01993 {
01994 pressure.lgContRadPresOn = false;
01995 }
01996 else
01997 {
01998 pressure.lgContRadPresOn = true;
01999 }
02000 return;
02001 }
02002
02003
02004 void DynaPrtZone( void )
02005 {
02006
02007 DEBUG_ENTRY( "DynaPrtZone()" );
02008
02009 ASSERT( nzone>0 && nzone<struc.nzlim );
02010
02011 if( nzone > 0 )
02012 {
02013 fprintf(ioQQQ," DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
02014 timesc.sound_speed_adiabatic/1e5 ,
02015 wind.windv/1e5 ,
02016 dynamics.Cool()/thermal.ctot,
02017 dynamics.Heat()/thermal.ctot);
02018 }
02019
02020 ASSERT( EnthalpyDensity[nzone-1] > 0. );
02021
02022 fprintf(ioQQQ," DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
02023 phycon.EnergyExcitation,
02024 phycon.EnergyIonization,
02025 phycon.EnergyBinding,
02026 0.5*POW2(wind.windv)*dense.xMassDensity,
02027 5./2.*pressure.PresGasCurr ,
02028 EnthalpyDensity[nzone-1]/scalingDensity() , AdvecSpecificEnthalpy
02029 );
02030 return;
02031 }
02032
02033
02034 void DynaPunchTimeDep( FILE* ipPnunit , const char *chJob )
02035 {
02036
02037 DEBUG_ENTRY( "DynaPunchTimeDep()" );
02038
02039 if( strcmp( chJob , "END" ) == 0 )
02040 {
02041 double te_mean,
02042 H2_mean,
02043 H0_mean,
02044 Hp_mean,
02045 Hep_mean;
02046
02047 if( cdTemp(
02048
02049 "HYDR",
02050
02051
02052 2,
02053
02054 &te_mean,
02055
02056 "RADIUS" ) )
02057 {
02058 TotalInsanity();
02059 }
02060 if( cdIonFrac(
02061
02062 "HYDR",
02063
02064
02065 2,
02066
02067 &Hp_mean,
02068
02069 "RADIUS" ,
02070
02071 false ) )
02072 {
02073 TotalInsanity();
02074 }
02075 if( cdIonFrac(
02076
02077 "HYDR",
02078
02079
02080 1,
02081
02082 &H0_mean,
02083
02084 "RADIUS" ,
02085
02086 false ) )
02087 {
02088 TotalInsanity();
02089 }
02090 if( cdIonFrac(
02091
02092 "H2 ",
02093
02094
02095 0,
02096
02097 &H2_mean,
02098
02099 "RADIUS" ,
02100
02101 false ) )
02102 {
02103 TotalInsanity();
02104 }
02105 if( cdIonFrac(
02106
02107 "HELI",
02108
02109
02110 2,
02111
02112 &Hep_mean,
02113
02114 "RADIUS" ,
02115
02116 false ) )
02117 {
02118 TotalInsanity();
02119 }
02120 fprintf( ipPnunit ,
02121 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" ,
02122 dynamics.time_elapsed ,
02123 dynamics.timestep ,
02124 rfield.time_continuum_scale ,
02125 scalingDensity(),
02126 te_mean ,
02127 Hp_mean ,
02128 H0_mean ,
02129 H2_mean ,
02130 Hep_mean ,
02131
02132 findspecieslocal("CO")->column / SDIV( colden.colden[ipCOL_HTOT] ),
02133 cosmology.redshift_current,
02134 dense.eden/dense.gas_phase[ipHYDROGEN]
02135 );
02136 }
02137 else
02138 TotalInsanity();
02139 return;
02140 }
02141
02142
02143 void DynaSave(FILE* ipPnunit , char chJob )
02144 {
02145 DEBUG_ENTRY( "DynaSave()" );
02146
02147 if( chJob=='a' )
02148 {
02149
02150 fprintf( ipPnunit , "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
02151 radius.depth_mid_zone,
02152 thermal.htot ,
02153 dynamics.Cool() ,
02154 dynamics.Heat() ,
02155 dynamics.dCooldT() ,
02156 dynamics.Source[ipHYDROGEN][ipHYDROGEN],
02157 dynamics.Rate,
02158 phycon.EnthalpyDensity/scalingDensity() ,
02159 AdvecSpecificEnthalpy
02160 );
02161 }
02162 else
02163 TotalInsanity();
02164 return;
02165 }
02166
02167 #define MERGE 0
02168 double t_dynamics::Heat()
02169 {
02170 double heat = Heat_v*scalingDensity();
02171 if (MERGE)
02172 {
02173 double cool = Cool_r*phycon.EnthalpyDensity;
02174 if (heat > cool)
02175 return heat-cool;
02176 else
02177 return 0.;
02178 }
02179 return heat;
02180 }
02181
02182 double t_dynamics::Cool()
02183 {
02184 double cool = Cool_r*phycon.EnthalpyDensity;
02185 if (MERGE)
02186 {
02187 double heat = Heat_v*scalingDensity();
02188 if (heat < cool)
02189 return cool-heat;
02190 else
02191 return 0.;
02192 }
02193 return cool;
02194 }
02195 #undef MERGE
02196
02197 double t_dynamics::dCooldT()
02198 {
02199 return Cool_r*5./2.*pressure.PresGasCurr/phycon.te;
02200 }
02201
02202 void DynaIterStart(void)
02203 {
02204 DEBUG_ENTRY( "DynaIterStart()" );
02205
02206 if( 0 == nTime_flux )
02207 {
02208 rfield.time_continuum_scale = 1.;
02209 return;
02210 }
02211 else if( dynamics.time_elapsed <= time_elapsed_time[0] )
02212 {
02213
02214 rfield.time_continuum_scale = (realnum)time_flux_ratio[0];
02215 }
02216 else if( dynamics.time_elapsed > time_elapsed_time[nTime_flux-1] )
02217 {
02218 fprintf( ioQQQ,
02219 " PROBLEM - DynaIterStart - I need the continuum at time %.2e but the table ends at %.2e.\n" ,
02220 dynamics.time_elapsed ,
02221 time_elapsed_time[nTime_flux-1]);
02222 cdEXIT(EXIT_FAILURE);
02223 }
02224 else
02225 {
02226 rfield.time_continuum_scale = (realnum)linint(
02227
02228 time_elapsed_time,
02229
02230 time_flux_ratio,
02231
02232 nTime_flux,
02233
02234 dynamics.time_elapsed);
02235 }
02236
02237 fprintf(ioQQQ,"DEBUG time dep reset continuum iter %ld dynamics.timestep %.2e elapsed time %.2e scale %.2e",
02238 iteration,
02239 dynamics.timestep ,
02240 dynamics.time_elapsed,
02241 rfield.time_continuum_scale);
02242 if( dynamics.lgRecom )
02243 {
02244 fprintf(ioQQQ," recom");
02245 }
02246 fprintf(ioQQQ,"\n");
02247
02248
02249 long int nTimeVary = 0;
02250 for( long int i=0; i < rfield.nShape; i++ )
02251 {
02252
02253
02254 if( rfield.lgTimeVary[i] )
02255 ++nTimeVary;
02256 }
02257
02258 if( hextra.lgTurbHeatVaryTime )
02259 {
02260
02261 hextra.TurbHeat = hextra.TurbHeatSave * rfield.time_continuum_scale;
02262 fprintf(ioQQQ,"DEBUG TurbHeat vary new heat %.2e\n",
02263 hextra.TurbHeat);
02264 }
02265 else if( !nTimeVary )
02266 {
02267 fprintf(ioQQQ," DISASTER - there were no variable continua "
02268 "or heat sources - put TIME option on at least one "
02269 "luminosity or hextra command.\n");
02270 cdEXIT(EXIT_FAILURE);
02271 }
02272 }
02273
02274