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 static int ipUpstream=-1,iphUpstream=-1,ipyUpstream=-1;
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 static double **UpstreamIon;
00060 static double ***UpstreamStatesElemNEW;
00061
00062 static double *UpstreamElem;
00063
00064
00065 static double *Upstream_H2_molec;
00066
00067
00068 static double *Upstream_CO_molec;
00069
00070
00071
00072
00073
00074 static double *time_elapsed_time ,
00075 *time_flux_ratio ,
00076 *time_dt,
00077 *time_dt_scale_factor;
00078 bool lgtime_dt_specified;
00079 int *lgtime_Recom;
00080 #define NTIME 200
00081
00082
00083 static long int nTime_flux=0;
00084
00085
00086 STATIC void DynaNewStep(void);
00087
00088
00089 STATIC void DynaSaveLast(void);
00090
00091
00092
00093
00094
00095 static double Dyn_dr;
00096
00097
00098 static double AdvecSpecificEnthalpy;
00099
00100
00101 static realnum *Old_histr ,
00102
00103 *Old_xLyman_depth,
00104
00105 *Old_depth,
00106
00107 *Old_hiistr,
00108
00109 *Old_pressure,
00110
00111 *Old_hden ,
00112
00113 *Old_DenMass ,
00114
00115 *EnthalpyDensity,
00116
00117 *Old_ednstr,
00118
00119 *Old_EnthalpyDensity;
00120
00121 static realnum **Old_H2_molec;
00122 static realnum **Old_CO_molec;
00123
00124
00125 static realnum ***Old_xIonDense;
00126
00127
00128 static realnum **Old_gas_phase;
00129
00130
00131 static realnum ****Old_StatesElemNEW;
00132
00133
00134 static long int nOld_zone;
00135
00136
00137 static realnum DivergePresInteg;
00138
00139
00140 STATIC double timestep_next( void )
00141 {
00142 static double te_old=-1;
00143 double timestep_Hp_temp , timestep_return;
00144
00145 DEBUG_ENTRY( "timestep_next()" );
00146
00147 timestep_return = dynamics.timestep;
00148
00149 if( dynamics.lgRecom )
00150 {
00151 double te_new;
00152 if( cdTemp(
00153
00154 "hydr",
00155
00156 2 ,
00157
00158 &te_new,
00159
00160 "VOLUME" ) )
00161 TotalInsanity();
00162
00163 if( te_old>0 )
00164 {
00165 double dTdStep = fabs(te_new-te_old)/te_new;
00166
00167 double dT_factor = 0.04/SDIV(dTdStep);
00168 dT_factor = MIN2( 2.0 , dT_factor );
00169 dT_factor = MAX2( 0.01 , dT_factor );
00170 timestep_Hp_temp = dynamics.timestep * dT_factor;
00171 }
00172 else
00173 {
00174 timestep_Hp_temp = -1.;
00175 }
00176 te_old = te_new;
00177 if( timestep_Hp_temp > 0. )
00178 timestep_return = timestep_Hp_temp;
00179 }
00180 else
00181 {
00182 timestep_return = dynamics.timestep_init;
00183 }
00184 fprintf(ioQQQ,"DEBUG timestep_next returns %.3e, old temp %.2e\n" , timestep_return, te_old );
00185
00186 return( timestep_return );
00187 }
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 #define SUBSONIC 1
00206 #define SUPERSONIC 2
00207
00208 #define STRONGD 4
00209 #define ORIGINAL 5
00210 #define SHOCK 6
00211 #define ANTISHOCK 7
00212 #define ANTISHOCK2 8
00213
00214 double DynaPresChngFactor(void)
00215 {
00216 double factor,
00217 er,
00218 width;
00219 static double dp = -1.,
00220 dpp = -1.,
00221 erp = -1.,
00222 erpp = -1.;
00225 static int lastzone = -1,
00226 zonepresmode,
00227 globalpresmode;
00228 int ipPRE;
00229
00230 DEBUG_ENTRY( "DynaPresChngFactor()" );
00231
00232
00233 PresTotCurrent();
00234
00235
00236 pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg*pressure.lgContRadPresOn
00237 + DivergePresInteg;
00238 if( trace.lgTrace )
00239 {
00240 fprintf( ioQQQ,
00241 " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n",
00242 pressure.PresTotlCorrect , pressure.PresTotlInit , pressure.PresInteg*pressure.lgContRadPresOn,
00243 DivergePresInteg );
00244 }
00245
00246 if( dynamics.lgTracePrint)
00247 fprintf(ioQQQ,"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
00248 pressure.PresTotlInit, pressure.PresInteg*pressure.lgContRadPresOn,
00249 DivergePresInteg, pressure.PresTotlCorrect, pressure.PresTotlCurr);
00250
00251
00252
00253
00254
00255
00256 if( !dynamics.lgSetPresMode )
00257 {
00258
00259
00260
00261 if(pressure.PresGasCurr < pressure.PresRamCurr)
00262 strcpy( dynamics.chPresMode , "supersonic" );
00263 else
00264 strcpy( dynamics.chPresMode , "subsonic" );
00265
00266 dynamics.lgSetPresMode = true;
00267 }
00268
00269
00270
00271
00272
00273
00274 if( strcmp( dynamics.chPresMode , "original" ) == 0 )
00275 {
00276 globalpresmode = ORIGINAL;
00277 pressure.lgSonicPointAbortOK = true;
00278 }
00279 else if( strcmp( dynamics.chPresMode , "subsonic" ) == 0 )
00280 {
00281 globalpresmode = SUBSONIC;
00282 pressure.lgSonicPointAbortOK = true;
00283 }
00284
00285 else if( strcmp( dynamics.chPresMode , "supersonic" ) == 0 )
00286 {
00287 globalpresmode = SUPERSONIC;
00288 pressure.lgSonicPointAbortOK = true;
00289 }
00290
00291 else if( strcmp( dynamics.chPresMode , "strongd" ) == 0 )
00292 {
00293 globalpresmode = STRONGD;
00294 pressure.lgSonicPointAbortOK = false;
00295 }
00296 else if( strcmp( dynamics.chPresMode , "shock" ) == 0 )
00297 {
00298 globalpresmode = SHOCK;
00299 pressure.lgSonicPointAbortOK = false;
00300 }
00301 else if( strcmp( dynamics.chPresMode , "antishock-by-mach" ) == 0 )
00302 {
00303 globalpresmode = ANTISHOCK2;
00304 pressure.lgSonicPointAbortOK = false;
00305 }
00306 else if( strcmp( dynamics.chPresMode , "antishock" ) == 0 )
00307 {
00308 globalpresmode = ANTISHOCK;
00309 pressure.lgSonicPointAbortOK = false;
00310 }
00311
00312
00313
00314 if(globalpresmode == ORIGINAL)
00315 {
00316
00317
00318 if(pressure.PresGasCurr > pressure.PresRamCurr)
00319 {
00320 zonepresmode = SUBSONIC;
00321 }
00322 else
00323 {
00324 zonepresmode = SUPERSONIC;
00325 }
00326 }
00327 else if(globalpresmode == STRONGD)
00328 {
00329 if(nzone <= 1)
00330 zonepresmode = SUPERSONIC;
00331 }
00332 else if(globalpresmode == SUBSONIC)
00333 {
00334 zonepresmode = SUBSONIC;
00335 }
00336 else if(globalpresmode == SUPERSONIC)
00337 {
00338 zonepresmode = SUPERSONIC;
00339 }
00340 else if(globalpresmode == SHOCK)
00341 {
00342 if(radius.depth < dynamics.ShockDepth)
00343 {
00344 zonepresmode = SUBSONIC;
00345 }
00346 else
00347 {
00348 zonepresmode = SUPERSONIC;
00349 }
00350 }
00351 else if(globalpresmode == ANTISHOCK)
00352 {
00353 if(radius.depth < dynamics.ShockDepth)
00354 {
00355 zonepresmode = SUPERSONIC;
00356 }
00357 else
00358 {
00359 zonepresmode = SUBSONIC;
00360 }
00361 }
00362 else if(globalpresmode == ANTISHOCK2)
00363 {
00364
00365
00366
00367 if( fabs(wind.windv) > dynamics.ShockMach*sqrt(pressure.PresGasCurr/dense.xMassDensity))
00368 {
00369 zonepresmode = SUPERSONIC;
00370 }
00371 else
00372 {
00373 zonepresmode = SUBSONIC;
00374 }
00375 }
00376 else
00377 {
00378 printf("Need to set global pressure mode\n");
00379 cdEXIT( EXIT_FAILURE );
00380 }
00381
00382 er = pressure.PresTotlCurr-pressure.PresTotlCorrect;
00383
00384
00385 if(globalpresmode == ORIGINAL || lastzone != nzone || fabs(er-erp) < SMALLFLOAT)
00386 {
00387
00388
00389
00390 if( zonepresmode == SUBSONIC )
00391 {
00392
00393 factor = pressure.PresTotlCorrect / pressure.PresTotlCurr;
00394 ipPRE = 0;
00395 }
00396 else
00397 {
00398
00399 factor = pressure.PresTotlCurr / pressure.PresTotlCorrect;
00400 ipPRE = 1;
00401 }
00402 if(fabs(factor-1.) > 0.01)
00403 {
00404 factor = 1.+sign(0.01,factor-1.);
00405 }
00406 erp = er;
00407 dp = dense.gas_phase[ipHYDROGEN];
00408 erpp = -1.;
00409 dpp = -1.;
00410 }
00411 else
00412 {
00413 #if 0
00414 printf("Ds: %d; %g %g %g; %g %g %g tot %g\n",nzone,dense.gas_phase[ipHYDROGEN],dp,dpp,er,erp,erpp,
00415 pressure.PresTotlCorrect);
00416 #endif
00417 if(1 || dpp < 0. || fabs((dense.gas_phase[ipHYDROGEN]-dp)*(dp-dpp)*(dpp-dense.gas_phase[ipHYDROGEN])) < SMALLFLOAT)
00418 {
00419
00420
00421 factor = (dense.gas_phase[ipHYDROGEN]*erp-dp*er)/(erp-er)/dense.gas_phase[ipHYDROGEN];
00422
00423 width = fabs(1.-dp/dense.gas_phase[ipHYDROGEN]);
00424 if(width > 1e-2)
00425 width = 1e-2;
00426
00427
00428
00429
00430
00431
00432
00433
00434 if(zonepresmode == SUBSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) < 0)
00435 {
00436 factor = 1+3*width;
00437 }
00438 else if(zonepresmode == SUPERSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) > 0)
00439 {
00440 factor = 1-3*width;
00441 }
00442
00443 if(fabs(factor-1.) > 3*width)
00444 {
00445 factor = 1.+sign(3*width,factor-1.);
00446 }
00447 ipPRE = 2;
00448 if(fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT)
00449 {
00450 dpp = dp;
00451 erpp = erp;
00452 dp = dense.gas_phase[ipHYDROGEN];
00453 erp = er;
00454 }
00455 }
00456 else
00457 {
00458
00459
00460
00461 double a, b, c, q, dmin, emin, dsol, f1 , f2;
00462 int iroot;
00463 a = er/(dense.gas_phase[ipHYDROGEN]-dp)/(dense.gas_phase[ipHYDROGEN]-dpp) +
00464 erp/(dp-dense.gas_phase[ipHYDROGEN])/(dp-dpp)+
00465 erpp/(dpp-dense.gas_phase[ipHYDROGEN])/(dpp-dp);
00466 b = (erp-erpp)/(dp-dpp) - a * (dp+dpp);
00467 c = erp - dp*(a*dp+b);
00468 dmin = (-0.5*b/a);
00469 emin = (a*dmin+b)*dmin+c;
00470 if(a < 0)
00471 {
00472 a = -a;
00473 b = -b;
00474 c = -c;
00475 }
00476 #if 0
00477 fprintf(ioQQQ,"Check 1: %g %g\n",er,(a*dense.gas_phase[ipHYDROGEN]+b)*dense.gas_phase[ipHYDROGEN]+c);
00478 fprintf(ioQQQ,"Check 2: %g %g\n",erp,(a*dp+b)*dp+c);
00479 fprintf(ioQQQ,"Check 3: %g %g\n",erpp,(a*dpp+b)*dpp+c);
00480 #endif
00481 q = b*b-4*a*c;
00482 if(q < 0)
00483 {
00484
00485
00486 factor = dmin/dense.gas_phase[ipHYDROGEN];
00487
00491 if(globalpresmode == STRONGD && -q > 1e-3*b*b)
00492 {
00493 zonepresmode = SUBSONIC;
00494 }
00495 }
00496 else
00497 {
00498
00499
00500
00501 if(zonepresmode == SUPERSONIC)
00502 iroot = 1;
00503 else
00504 iroot = 0;
00505
00506 if(emin > 0)
00507 iroot = ! iroot;
00508
00509 if(iroot)
00510 {
00511
00512 if(b > 0)
00513 {
00514 dsol = -(b+sqrt(q))/(2*a);
00515 }
00516 else
00517 {
00518 dsol = 2*c/(-b+sqrt(q));
00519 }
00520 }
00521 else
00522 {
00523 if(b < 0)
00524 {
00525 dsol = (-b+sqrt(q))/(2*a);
00526 }
00527 else
00528 {
00529 dsol = -2*c/(b+sqrt(q));
00530 }
00531 }
00532 factor = dsol/dense.gas_phase[ipHYDROGEN];
00533 #if 0
00534 fprintf(ioQQQ,"Check 4: %g %g %d %g %g\n",dsol,(a*dsol+b)*dsol+c,iroot,dmin,emin);
00535 #endif
00536 }
00537
00538 f1 = fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]);
00539 f2 = fabs(1.- dp/dense.gas_phase[ipHYDROGEN]);
00540
00541 width = MAX2(f1,f2);
00542
00543 if(fabs(factor-1.) > 3*width)
00544 {
00545 factor = 1.+sign(3*width,factor-1.);
00546 }
00547 ipPRE = 3;
00548 if(fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT)
00549 {
00550 dpp = dp;
00551 erpp = erp;
00552 dp = dense.gas_phase[ipHYDROGEN];
00553 erp = er;
00554 }
00555 }
00556 }
00557
00558 #if 0
00559 printf("Out: %d; %g; %g %g; %g %g\n",nzone,factor*dense.gas_phase[ipHYDROGEN],dp,dpp,erp,erpp);
00560 #endif
00561 lastzone = nzone;
00562
00563 if( dynamics.lgTracePrint )
00564 fprintf(ioQQQ,"windv %li r %g v %g f %g\n",
00565 nzone,radius.depth,wind.windv,DynaFlux(radius.depth));
00566
00567
00568 if( trace.nTrConvg>=1 )
00569 {
00570 fprintf( ioQQQ,
00571 " DynaPresChngFactor mode %s scale factor %.3f vel %.3e\n",
00572 dynamics.chPresMode , factor , wind.windv );
00573 }
00574
00575 {
00576
00577 enum{DEBUG_LOC=false};
00578
00579 if( DEBUG_LOC )
00580 {
00581 char chPRE[][4] = {"gas" , "ram", "sec", "par" };
00582 fprintf(ioQQQ,
00583 "pre %s\tfac\t%.5f\n",
00584 chPRE[ipPRE],
00585 factor
00586 );
00587 }
00588 }
00589
00590
00591 ASSERT( wind.windv != 0. || (wind.windv == 0. && dynamics.lgTimeDependentStatic) );
00592
00593 return factor;
00594 }
00595
00596
00597
00598
00599 void DynaIonize(void)
00600 {
00601 long int nelem,
00602 ion,
00603 mol ,
00604 i;
00605
00606 DEBUG_ENTRY( "DynaIonize()" );
00607
00608
00609
00610 if( !dynamics.lgTimeDependentStatic )
00611 {
00612
00613
00614
00615 dynamics.timestep = -Dyn_dr/wind.windv;
00616 }
00617
00618
00619 ASSERT(nzone<struc.nzlim );
00620 if( nzone > 0 )
00621 EnthalpyDensity[nzone-1] = (realnum)phycon.EnthalpyDensity;
00622
00623
00624
00625
00626
00627
00628
00629 double depth = radius.depth;
00630 if( iteration < dynamics.n_initial_relax+1 ||
00631 ( ! dynamics.lgTimeDependentStatic &&
00632 ( depth < 0 || depth > dynamics.oldFullDepth ) ) )
00633 {
00634
00635 dynamics.Cool_r = 0.;
00636 dynamics.Heat_v = 0.;
00637 dynamics.dHeatdT = 0.;
00638
00639 dynamics.Rate = 0.;
00640
00641 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00642 {
00643 for( ion=0; ion<nelem+2; ++ion )
00644 {
00645 dynamics.Source[nelem][ion] = 0.;
00646 }
00647 }
00648 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00649 {
00650 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00651 {
00652 if( dense.lgElmtOn[nelem] )
00653 {
00654 for( long level=0; level < iso.numLevels_local[ipISO][nelem]; ++level )
00655 {
00656 dynamics.StatesElemNEW[nelem][nelem-ipISO][level] = 0.;
00657 }
00658 }
00659 }
00660 }
00661 for(mol=0;mol<N_H_MOLEC;mol++)
00662 {
00663 dynamics.H2_molec[mol] = 0.;
00664 }
00665 for(mol=0;mol<mole.num_comole_calc;mol++)
00666 {
00667 dynamics.CO_molec[mol] = 0.;
00668 }
00669 return;
00670 }
00671
00672 if( dynamics.lgTracePrint )
00673 {
00674 fprintf(ioQQQ,"workwork\t%li\t%.3e\t%.3e\t%.3e\n",
00675 nzone,
00676 phycon.EnthalpyDensity,
00677 0.5*POW2(wind.windv)*dense.xMassDensity ,
00678 5./2.*pressure.PresGasCurr
00679 );
00680 }
00681
00682
00683
00684
00685
00686 dynamics.Cool_r = 1./dynamics.timestep*dynamics.lgCoolHeat;
00687 dynamics.Heat_v = AdvecSpecificEnthalpy/dynamics.timestep*dynamics.lgCoolHeat;
00688 dynamics.dHeatdT = 0.*dynamics.lgCoolHeat;
00689
00690 #if 0
00691
00692 if(dynamics.Cool > 0) {
00693 dynamics.Heat = 0.;
00694
00695 dynamics.dCooldT = 5./2.*pressure.PresGasCurr/phycon.te/dynamics.timestep*dynamics.lgCoolHeat;
00696 dynamics.dHeatdT = 0.*dynamics.lgCoolHeat;
00697 } else {
00698 dynamics.Heat = -dynamics.Cool;
00699 dynamics.Cool = 0.;
00700
00701 dynamics.dCooldT = 0.*dynamics.lgCoolHeat;
00702 dynamics.dHeatdT = -5./2.*pressure.PresGasCurr/phycon.te/dynamics.timestep*dynamics.lgCoolHeat;
00703 }
00704 #endif
00705
00706 # if 0
00707 if( dynamics.lgTracePrint || nzone>17 && iteration == 10)
00708 {
00709 fprintf(ioQQQ,
00710 "dynamics cool-heat\t%li\t%.3e\t%i\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
00711 nzone,
00712 phycon.te,
00713 dynamics.lgCoolHeat,
00714 thermal.htot,
00715 phycon.EnthalpyDensity/dynamics.timestep,
00716 AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/dynamics.timestep,
00717 phycon.EnthalpyDensity,
00718 AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN],
00719 dense.gas_phase[ipHYDROGEN],
00720 dynamics.timestep);
00721 }
00722 # endif
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737 dynamics.Rate = 1./dynamics.timestep;
00738
00739 for(mol=0;mol<N_H_MOLEC;mol++)
00740 {
00741 dynamics.H2_molec[mol] = Upstream_H2_molec[mol]* dense.gas_phase[ipHYDROGEN]*dynamics.Rate;
00742 }
00743
00744 for(mol=0;mol<mole.num_comole_calc;mol++)
00745 {
00746 dynamics.CO_molec[mol] = Upstream_CO_molec[mol]* dense.gas_phase[ipHYDROGEN]*dynamics.Rate;
00747 }
00748
00749 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00750 {
00751 if( dense.lgElmtOn[nelem] )
00752 {
00753
00754 if(fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]>=1e-3)
00755 {
00756
00757 realnum sumh = 0.;
00758 for(i=0;i<N_H_MOLEC;i++)
00759 {
00760 sumh += hmi.Hmolec[i]*hmi.nProton[i];
00761 }
00762 fprintf(ioQQQ,
00763 "PROBLEM conservation error: zn %li elem %li upstream %.8e abund %.8e (up-ab)/up %.2e f1(H n CO) %.2e f2(H n CO) %.2e\n",
00764 nzone ,
00765 nelem,
00766 UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN],
00767 dense.gas_phase[nelem] ,
00768 (UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN]-dense.gas_phase[nelem]) /
00769 (UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN]) ,
00770 (dense.gas_phase[ipHYDROGEN]-sumh) / dense.gas_phase[ipHYDROGEN] ,
00771 dense.H_sum_in_CO / dense.gas_phase[ipHYDROGEN] );
00772 }
00773
00774 for( ion=0; ion<dense.IonLow[nelem]; ++ion )
00775 {
00776 dynamics.Source[nelem][ion] = 0.;
00777 }
00778 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00779 {
00780
00781
00782
00783
00784 dynamics.Source[nelem][ion] =
00785
00786
00787
00788 UpstreamIon[nelem][ion] * dense.gas_phase[ipHYDROGEN] / dynamics.timestep;
00789
00790 }
00791 for( ion=dense.IonHigh[nelem]+1;ion<nelem+2; ++ion )
00792 {
00793 dynamics.Source[nelem][ion] = 0.;
00794 }
00795 }
00796 }
00797
00798 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00799 {
00800 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00801 {
00802 if( dense.lgElmtOn[nelem] )
00803 {
00804 for( long level=0; level < iso.numLevels_local[ipISO][nelem]; ++level )
00805 {
00806 dynamics.StatesElemNEW[nelem][nelem-ipISO][level] =
00807 UpstreamStatesElemNEW[nelem][nelem-ipISO][level] * dense.gas_phase[ipHYDROGEN] / dynamics.timestep;
00808 }
00809 }
00810 }
00811 }
00812
00813 # if 0
00814 fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
00815 nzone,
00816 dynamics.Rate,
00817 dynamics.Source[ipHYDROGEN][0],
00818 dynamics.Rate,
00819 dynamics.Source[ipCARBON][3]);
00820 # endif
00821 #if 0
00822 nelem = ipCARBON;
00823 ion = 3;
00824
00825 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",
00826 nzone,
00827 ipUpstream,
00828 radius.depth ,
00829 Old_depth[ipUpstream],
00830 dense.xIonDense[nelem][ion],
00831 UpstreamIon[nelem][ion]* dense.gas_phase[ipHYDROGEN],
00832 Old_xIonDense[ipUpstream][nelem][ion] ,
00833 dense.xIonDense[nelem][ion+1],
00834 UpstreamIon[nelem][ion+1]* dense.gas_phase[ipHYDROGEN],
00835 Old_xIonDense[ipUpstream][nelem][ion+1] ,
00836 dynamics.timestep,
00837 dynamics.Source[nelem][ion]
00838 );
00839 #endif
00840 if( dynamics.lgTracePrint )
00841 {
00842 fprintf(ioQQQ," DynaIonize, %4li photo=%.2e , H recom= %.2e \n",
00843 nzone,dynamics.Rate , dynamics.Source[0][0] );
00844 }
00845 return;
00846 }
00847
00848
00849
00850 void DynaStartZone(void)
00851 {
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 double upstream, dilution, dilutionleft, dilutionright, frac_next;
00865
00866
00867 double hupstream, hnextfrac=-BIGFLOAT, hion, hmol, hiso;
00868
00869
00870 double ynextfrac=-BIGFLOAT, yion, ymol, yiso;
00871
00872 long int nelem , ion, mol;
00873
00874 DEBUG_ENTRY( "DynaStartZone()" );
00875
00876
00877 if( iteration < 2 )
00878 {
00879 dynamics.Upstream_hden = 0.;
00880 AdvecSpecificEnthalpy = 0.;
00881 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00882 {
00883 for( ion=0; ion<nelem+2; ++ion )
00884 {
00885 UpstreamIon[nelem][ion] = 0.;
00886 }
00887 }
00888 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00889 {
00890 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00891 {
00892 if( dense.lgElmtOn[nelem] )
00893 {
00894 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
00895 {
00896 UpstreamStatesElemNEW[nelem][nelem-ipISO][level] = 0.;
00897 }
00898 }
00899 }
00900 }
00901
00902 for(mol=0; mol<N_H_MOLEC;mol++)
00903 {
00904 Upstream_H2_molec[mol] = 0;
00905 }
00906 for(mol=0; mol<mole.num_comole_calc;mol++)
00907 {
00908 Upstream_CO_molec[mol] = 0;
00909 }
00910
00911 ipUpstream = -1;
00912 iphUpstream = -1;
00913 ipyUpstream = -1;
00914 return;
00915 }
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926 upstream = MAX2(Old_depth[0] , radius.depth + Dyn_dr);
00927 hupstream = 0.5*(radius.depth + upstream);
00928
00929
00930
00931 while( (Old_depth[ipUpstream+1] < upstream ) &&
00932 ( ipUpstream < nOld_zone-1 ) )
00933 {
00934 ++ipUpstream;
00935 }
00936 ASSERT( ipUpstream <= nOld_zone-1 );
00937
00938
00939
00940 if( (ipUpstream != -1) && (ipUpstream != nOld_zone-1)&& ((Old_depth[ipUpstream+1] - Old_depth[ipUpstream])> SMALLFLOAT) )
00941 {
00942
00943 frac_next = ( upstream - Old_depth[ipUpstream])/
00944 (Old_depth[ipUpstream+1] - Old_depth[ipUpstream]);
00945 dynamics.Upstream_hden = (realnum)(Old_hden[ipUpstream] +
00946 (Old_hden[ipUpstream+1] - Old_hden[ipUpstream])*
00947 frac_next);
00948 dilutionleft = 1./Old_hden[ipUpstream];
00949 dilutionright = 1./Old_hden[ipUpstream+1];
00950
00951
00952
00953
00954 dilution = 1./dynamics.Upstream_hden;
00955
00956
00957 AdvecSpecificEnthalpy = (Old_EnthalpyDensity[ipUpstream]*dilutionleft +
00958 (Old_EnthalpyDensity[ipUpstream+1]*dilutionright - Old_EnthalpyDensity[ipUpstream]*dilutionleft)*
00959 frac_next);
00960
00961 ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] );
00962
00963
00964 realnum lo = (realnum)(Old_EnthalpyDensity[ipUpstream]*dilutionleft);
00965 realnum x = (realnum)AdvecSpecificEnthalpy;
00966 realnum hi = (realnum)(Old_EnthalpyDensity[ipUpstream+1]*dilutionright);
00967 if( ! fp_bound(lo,x,hi) )
00968 {
00969 fprintf(ioQQQ,"PROBLEM interpolated enthalpy density is not within range %.16e\t%.16e\t%.16e\t%e\t%e\n",
00970 lo, x, hi, (hi-x)/(hi-lo), (x-lo)/(hi-lo));
00971 }
00972
00973 ASSERT( fp_bound(lo,x,hi) );
00974
00975 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00976 {
00977 UpstreamElem[nelem] = 0.;
00978 for( ion=0; ion<nelem+2; ++ion )
00979 {
00980
00981
00982
00983
00984
00985 UpstreamIon[nelem][ion] =
00986 Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream] +
00987 (Old_xIonDense[ipUpstream+1][nelem][ion]/Old_hden[ipUpstream+1] -
00988 Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream])*
00989 frac_next;
00990
00991 UpstreamElem[nelem] += UpstreamIon[nelem][ion];
00992 }
00993 }
00994
00995 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00996 {
00997 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00998 {
00999 if( dense.lgElmtOn[nelem] )
01000 {
01001 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
01002 {
01003 UpstreamStatesElemNEW[nelem][nelem-ipISO][level] =
01004 Old_StatesElemNEW[ipUpstream][nelem][nelem-ipISO][level]/Old_hden[ipUpstream] +
01005 (Old_StatesElemNEW[ipUpstream+1][nelem][nelem-ipISO][level]/Old_hden[ipUpstream+1] -
01006 Old_StatesElemNEW[ipUpstream][nelem][nelem-ipISO][level]/Old_hden[ipUpstream])*
01007 frac_next;
01008
01009 ASSERT( !isnan( UpstreamStatesElemNEW[nelem][nelem-ipISO][level] ) );
01010 }
01011 }
01012 }
01013 }
01014
01015 for(mol=0;mol<N_H_MOLEC;mol++)
01016 {
01017 Upstream_H2_molec[mol] = Old_H2_molec[ipUpstream][mol]/Old_hden[ipUpstream] +
01018 (Old_H2_molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] -
01019 Old_H2_molec[ipUpstream][mol]/Old_hden[ipUpstream])*
01020 frac_next;
01021
01022
01023 if(mol > 1)
01024 UpstreamElem[ipHYDROGEN] += Upstream_H2_molec[mol]*hmi.nProton[mol];
01025 }
01026
01027
01028 for(mol=0;mol<mole.num_comole_calc;mol++)
01029 {
01030 Upstream_CO_molec[mol] = Old_CO_molec[ipUpstream][mol]/Old_hden[ipUpstream] +
01031 (Old_CO_molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] -
01032 Old_CO_molec[ipUpstream][mol]/Old_hden[ipUpstream])*
01033 frac_next;
01034
01035 if(COmole[mol]->n_nuclei > 1)
01036 {
01037 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01038 {
01039 if( mole.lgElem_in_chemistry[nelem] )
01040 {
01041 UpstreamElem[nelem] +=
01042 COmole[mol]->nElem[nelem]*Upstream_CO_molec[mol];
01043 }
01044 }
01045 }
01046 }
01047 }
01048 else
01049 {
01050
01051 long ipBound = ipUpstream;
01052 if (ipBound == -1)
01053 ipBound = 0;
01054 dynamics.Upstream_hden = Old_hden[ipBound];
01055
01056 dilution = 1./dynamics.Upstream_hden;
01057
01058 AdvecSpecificEnthalpy = Old_EnthalpyDensity[ipBound]/Old_hden[ipBound];
01059 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01060 {
01061 UpstreamElem[nelem] = 0.;
01062 for( ion=0; ion<nelem+2; ++ion )
01063 {
01064
01065 UpstreamIon[nelem][ion] =
01066 Old_xIonDense[ipBound][nelem][ion]/Old_hden[ipBound];
01067 UpstreamElem[nelem] += UpstreamIon[nelem][ion];
01068 }
01069 }
01070
01071 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01072 {
01073 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01074 {
01075 if( dense.lgElmtOn[nelem] )
01076 {
01077 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
01078 {
01079 UpstreamStatesElemNEW[nelem][nelem-ipISO][level] =
01080 Old_StatesElemNEW[ipBound][nelem][nelem-ipISO][level]/Old_hden[ipBound];
01081
01082 ASSERT( !isnan( UpstreamStatesElemNEW[nelem][nelem-ipISO][level] ) );
01083 }
01084 }
01085 }
01086 }
01087
01088 for(mol=0;mol<N_H_MOLEC;mol++)
01089 {
01090 Upstream_H2_molec[mol] = Old_H2_molec[ipBound][mol]/Old_hden[ipBound];
01091 if(mol > 1)
01092 UpstreamElem[ipHYDROGEN] += Upstream_H2_molec[mol]*hmi.nProton[mol];
01093 }
01094 for(mol=0;mol<mole.num_comole_calc;mol++)
01095 {
01096 Upstream_CO_molec[mol] = Old_CO_molec[ipBound][mol]/Old_hden[ipBound];
01097 if(COmole[mol]->n_nuclei > 1)
01098 {
01099 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01100 {
01101 if( mole.lgElem_in_chemistry[nelem] )
01102 {
01103 UpstreamElem[nelem] +=
01104 Upstream_CO_molec[mol]*COmole[mol]->nElem[nelem];
01105 }
01106 }
01107 }
01108 }
01109 }
01110
01111
01112
01113
01114 while( (Old_depth[iphUpstream+1] < hupstream ) &&
01115 ( iphUpstream < nOld_zone-1 ) )
01116 {
01117 ++iphUpstream;
01118 }
01119 ASSERT( iphUpstream <= nOld_zone-1 );
01120
01121 while( (Old_depth[ipyUpstream+1] < radius.depth ) &&
01122 ( ipyUpstream < nOld_zone-1 ) )
01123 {
01124 ++ipyUpstream;
01125 }
01126 ASSERT( ipyUpstream <= nOld_zone-1 );
01127
01128 dynamics.dRad = BIGFLOAT;
01129
01130 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
01131 hnextfrac = ( hupstream - Old_depth[iphUpstream])/
01132 (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]);
01133 else
01134 hnextfrac = 0.;
01135
01136 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
01137 ynextfrac = ( radius.depth - Old_depth[ipyUpstream])/
01138 (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]);
01139 else
01140 ynextfrac = 0.;
01141
01142
01143 if(ipUpstream != -1 && ipUpstream < nOld_zone-1)
01144 dynamics.dRad = MIN2(dynamics.dRad,
01145 5*(Old_depth[ipUpstream+1] - Old_depth[ipUpstream]));
01146
01147
01148
01149 const double STEP_FACTOR=0.05;
01150
01151 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01152 {
01153 for( ion=0; ion<nelem+2; ++ion )
01154 {
01155 double f1;
01156 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
01157 {
01158 yion =
01159 Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream] +
01160 (Old_xIonDense[ipyUpstream+1][nelem][ion]/Old_hden[ipyUpstream+1] -
01161 Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream])*
01162 ynextfrac;
01163 }
01164 else
01165 {
01166 long ipBound = ipyUpstream;
01167 if (-1 == ipBound)
01168 ipBound = 0;
01169 yion = Old_xIonDense[ipBound][nelem][ion]/Old_hden[ipBound];
01170 }
01171 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
01172 {
01173 hion =
01174 Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream] +
01175 (Old_xIonDense[iphUpstream+1][nelem][ion]/Old_hden[iphUpstream+1] -
01176 Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream])*
01177 hnextfrac;
01178 }
01179 else
01180 {
01181 long ipBound = iphUpstream;
01182 if (-1 == ipBound)
01183 ipBound = 0;
01184 hion = Old_xIonDense[ipBound][nelem][ion]/Old_hden[ipBound];
01185 }
01186
01187
01188 f1 = fabs(yion - UpstreamIon[nelem][ion] );
01189 if( f1 > SMALLFLOAT )
01190 {
01191 dynamics.dRad = MIN2(dynamics.dRad,STEP_FACTOR*fabs(Dyn_dr) *
01192
01193 MAX2(yion + UpstreamIon[nelem][ion],1e-6 ) / f1);
01194 }
01195
01196
01197
01198
01199
01200
01201 dynamics.discretization_error += POW2(yion-2.*hion+UpstreamIon[nelem][ion]);
01202 dynamics.error_scale2 += POW2(UpstreamIon[nelem][ion]-yion);
01203 }
01204 }
01205
01206 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01207 {
01208 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01209 {
01210 if( dense.lgElmtOn[nelem] )
01211 {
01212 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
01213 {
01214 double f1;
01215 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 &&
01216 (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
01217 {
01218 yiso =
01219 Old_StatesElemNEW[ipyUpstream][nelem][nelem-ipISO][level]/Old_hden[ipyUpstream] +
01220 (Old_StatesElemNEW[ipyUpstream+1][nelem][nelem-ipISO][level]/Old_hden[ipyUpstream+1] -
01221 Old_StatesElemNEW[ipyUpstream][nelem][nelem-ipISO][level]/Old_hden[ipyUpstream])*
01222 ynextfrac;
01223 }
01224 else
01225 {
01226 long ipBound = ipyUpstream;
01227 if (-1 == ipBound)
01228 ipBound = 0;
01229 yiso = Old_StatesElemNEW[ipBound][nelem][nelem-ipISO][level]/Old_hden[ipBound];
01230 }
01231 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 &&
01232 (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
01233 {
01234 hiso =
01235 Old_StatesElemNEW[iphUpstream][nelem][nelem-ipISO][level]/Old_hden[iphUpstream] +
01236 (Old_StatesElemNEW[iphUpstream+1][nelem][nelem-ipISO][level]/Old_hden[iphUpstream+1] -
01237 Old_StatesElemNEW[iphUpstream][nelem][nelem-ipISO][level]/Old_hden[iphUpstream])*
01238 hnextfrac;
01239 }
01240 else
01241 {
01242 long ipBound = iphUpstream;
01243 if (-1 == ipBound)
01244 ipBound = 0;
01245 hiso = Old_StatesElemNEW[ipBound][nelem][nelem-ipISO][level]/Old_hden[ipBound];
01246 }
01247
01248
01249 f1 = fabs(yiso - UpstreamStatesElemNEW[nelem][nelem-ipISO][level] );
01250 if( f1 > SMALLFLOAT )
01251 {
01252 dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
01253
01254 MAX2(yiso + UpstreamStatesElemNEW[nelem][nelem-ipISO][level],1e-6 ) / f1);
01255 }
01256
01257
01258
01259
01260
01261
01262 dynamics.discretization_error += POW2(yiso-2.*hiso+UpstreamStatesElemNEW[nelem][nelem-ipISO][level]);
01263 dynamics.error_scale2 += POW2(UpstreamStatesElemNEW[nelem][nelem-ipISO][level]);
01264 }
01265 }
01266 }
01267 }
01268
01269
01270 for(mol=0; mol < N_H_MOLEC; mol++)
01271 {
01272 double f1;
01273 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
01274 {
01275 ymol =
01276 Old_H2_molec[ipyUpstream][mol]/Old_hden[ipyUpstream] +
01277 (Old_H2_molec[ipyUpstream+1][mol]/Old_hden[ipyUpstream+1] -
01278 Old_H2_molec[ipyUpstream][mol]/Old_hden[ipyUpstream])*
01279 ynextfrac;
01280 }
01281 else
01282 {
01283 long ipBound = ipyUpstream;
01284 if (-1 == ipBound)
01285 ipBound = 0;
01286 ymol = Old_H2_molec[ipBound][mol]/Old_hden[ipBound];
01287 }
01288 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
01289 {
01290 hmol =
01291 Old_H2_molec[iphUpstream][mol]/Old_hden[iphUpstream] +
01292 (Old_H2_molec[iphUpstream+1][mol]/Old_hden[iphUpstream+1] -
01293 Old_H2_molec[iphUpstream][mol]/Old_hden[iphUpstream])*
01294 hnextfrac;
01295 }
01296 else
01297 {
01298 long ipBound = iphUpstream;
01299 if (-1 == ipBound)
01300 ipBound = 0;
01301 hmol = Old_H2_molec[ipBound][mol]/Old_hden[ipBound];
01302 }
01303
01304
01305 f1 = fabs(ymol - Upstream_H2_molec[mol] );
01306 if( f1 > SMALLFLOAT )
01307 {
01308 dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
01309
01310 MAX2(ymol + Upstream_H2_molec[mol],1e-6 ) / f1 );
01311 }
01312
01313
01314 dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_H2_molec[mol]);
01315 dynamics.error_scale2 += POW2(Upstream_H2_molec[mol]-ymol);
01316 }
01317
01318 for(mol=0; mol < mole.num_comole_calc; mol++)
01319 {
01320 double f1;
01321 if(ipyUpstream != -1 && ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
01322 {
01323 ymol =
01324 Old_CO_molec[ipyUpstream][mol]/Old_hden[ipyUpstream] +
01325 (Old_CO_molec[ipyUpstream+1][mol]/Old_hden[ipyUpstream+1] -
01326 Old_CO_molec[ipyUpstream][mol]/Old_hden[ipyUpstream])*
01327 ynextfrac;
01328 }
01329 else
01330 {
01331 long ipBound = ipyUpstream;
01332 if (-1 == ipBound)
01333 ipBound = 0;
01334 ymol = Old_CO_molec[ipBound][mol]/Old_hden[ipBound];
01335 }
01336 if(iphUpstream != -1 && iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
01337 {
01338 hmol =
01339 Old_CO_molec[iphUpstream][mol]/Old_hden[iphUpstream] +
01340 (Old_CO_molec[iphUpstream+1][mol]/Old_hden[iphUpstream+1] -
01341 Old_CO_molec[iphUpstream][mol]/Old_hden[iphUpstream])*
01342 hnextfrac;
01343 }
01344 else
01345 {
01346 long ipBound = iphUpstream;
01347 if (-1 == ipBound)
01348 ipBound = 0;
01349 hmol = Old_CO_molec[ipBound][mol]/Old_hden[ipBound];
01350 }
01351
01352
01353 f1 = fabs(ymol - Upstream_CO_molec[mol] );
01354 if( f1 > SMALLFLOAT )
01355 {
01356 dynamics.dRad = MIN2(dynamics.dRad,fabs(Dyn_dr*STEP_FACTOR) *
01357
01358 MAX2(ymol + Upstream_CO_molec[mol],1e-6 ) / f1 );
01359 }
01360
01361
01362 dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_CO_molec[mol]);
01363 dynamics.error_scale2 += POW2(Upstream_CO_molec[mol]-ymol);
01364 }
01365
01366 if( dynamics.lgTracePrint )
01367 {
01368 fprintf(ioQQQ," DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n",
01369 nzone,dynamics.Rate , dynamics.Source[0][0] , dilution*dense.gas_phase[ipHYDROGEN] );
01370 }
01371 return;
01372 }
01373
01374
01375 void DynaEndZone(void)
01376 {
01377 DEBUG_ENTRY( "DynaEndZone()" );
01378
01379
01380
01381 DivergePresInteg += wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad));
01382
01383 if(dynamics.lgTracePrint)
01384 fprintf(ioQQQ,"Check dp: %g %g mom %g %g mas %g\n",
01385 wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad)),
01386 2*wind.windv*DynaFlux(radius.depth)*radius.drad/(1e16-radius.depth),
01387 wind.windv*DynaFlux(radius.depth),
01388 wind.windv*DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth),
01389 DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth));
01390 return;
01391 }
01392
01393
01394
01395
01396 void DynaIterEnd(void)
01397 {
01398
01399
01400
01401
01402 long int i;
01403 static long int nTime_dt_array_element = 0;
01404
01405 DEBUG_ENTRY( "DynaIterEnd()" );
01406
01407
01408
01409
01410
01411
01412
01413 if( iteration == dynamics.n_initial_relax+1)
01414 {
01415 fprintf(ioQQQ,"DYNAMICS DynaIterEnd sets stop radius to %.2e after "
01416 "dynamics.n_initial_relax=%li iterations.\n",
01417 radius.depth,
01418 dynamics.n_initial_relax);
01419 for( i=iteration-1; i<iterations.iter_malloc; ++i )
01420
01421
01422
01423
01424 radius.StopThickness[i] = radius.depth;
01425 }
01426
01427 DivergePresInteg = 0.;
01428
01429
01430
01431
01432
01433
01434 if( !dynamics.lgTimeDependentStatic )
01435 {
01436 if(iteration == dynamics.n_initial_relax+1 )
01437 {
01438
01439
01440
01441
01442
01443
01444
01445 if( dynamics.AdvecLengthInit> 0. )
01446 {
01447 Dyn_dr = dynamics.AdvecLengthInit;
01448 }
01449 else
01450 {
01451
01452 Dyn_dr = -dynamics.AdvecLengthInit*radius.depth;
01453 }
01454
01455 if (wind.windv0 > 0)
01456 Dyn_dr = -Dyn_dr;
01457
01458 if( dynamics.lgTracePrint )
01459 {
01460 fprintf(ioQQQ," DynaIterEnd, dr=%.2e \n",
01461 Dyn_dr );
01462 }
01463 }
01464 else if(iteration > dynamics.n_initial_relax+1 )
01465 {
01466
01467 DynaNewStep();
01468 }
01469 }
01470 else
01471 {
01472
01473 static double HeatInitial=-1. , HeatRadiated=-1. ,
01474 DensityInitial=-1;
01475
01476
01477 Dyn_dr = 0.;
01478 fprintf(ioQQQ,
01479 "DEBUG times enter dynamics.timestep %.2e elapsed_time %.2e iteration %li relax %li \n",
01480 dynamics.timestep ,
01481 dynamics.time_elapsed,
01482 iteration , dynamics.n_initial_relax);
01483 if( iteration > dynamics.n_initial_relax )
01484 {
01485
01486 DynaNewStep();
01487
01488
01489
01490 if( dynamics.lg_coronal_time_init )
01491 {
01492 thermal.lgTemperatureConstant = false;
01493 thermal.ConstTemp = 0.;
01494 }
01495
01496
01497
01498
01499 dynamics.time_elapsed += dynamics.timestep;
01500
01501 if( dynamics.timestep_stop > 0 && dynamics.time_elapsed > dynamics.timestep_stop )
01502 {
01503 dynamics.lgStatic_completed = true;
01504 }
01505
01506
01507 if( (phycon.te < StopCalc.TempLoStopIteration) ||
01508 (phycon.te > StopCalc.TempHiStopIteration ) )
01509 dynamics.lgStatic_completed = true;
01510
01511
01512
01513 HeatRadiated += (thermal.ctot-dynamics.Cool()) * dynamics.timestep *
01514 (DensityInitial / dense.gas_phase[ipHYDROGEN]);
01515 }
01516 else
01517 {
01518
01519 HeatInitial = 1.5 * pressure.PresGasCurr;
01520 HeatRadiated = 0.;
01521 DensityInitial = dense.gas_phase[ipHYDROGEN];
01522 fprintf(ioQQQ,"DEBUG relaxing times requested %li this is step %li\n",
01523 dynamics.n_initial_relax , iteration);
01524 }
01525 fprintf(ioQQQ,"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n",
01526 HeatInitial , HeatRadiated );
01527
01528
01529
01530 if( dynamics.time_elapsed > time_elapsed_time[nTime_dt_array_element] )
01531 {
01532
01533
01534 ++nTime_dt_array_element;
01535
01536
01537 ASSERT( nTime_dt_array_element < nTime_flux );
01538
01539
01540 if( lgtime_Recom[nTime_dt_array_element] )
01541 {
01542 fprintf(ioQQQ,"DEBUG dynamics turn on recombination logic\n");
01543 dynamics.lgRecom = true;
01544
01545
01546
01547
01548 radius.sdrmax = radius.dr_max_last_iter/3.;
01549 radius.lgSdrmaxRel = false;
01550 }
01551
01552 if( lgtime_dt_specified )
01553 {
01554
01555 fprintf(ioQQQ,"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element,
01556 dynamics.timestep);
01557 dynamics.timestep = time_dt[nTime_dt_array_element];
01558
01559 if( time_dt_scale_factor[nTime_dt_array_element] > 0. )
01560 dynamics.timestep_factor = time_dt_scale_factor[nTime_dt_array_element];
01561 }
01562 }
01563 else if( lgtime_dt_specified )
01564 {
01565
01566 dynamics.timestep *= dynamics.timestep_factor;
01567 fprintf(ioQQQ,"DEBUG lgtimes increment Timeby dynamics.timestep_factor to %li %.2e\n" ,
01568 nTime_dt_array_element,
01569 dynamics.timestep );
01570 if(dynamics.time_elapsed+dynamics.timestep > time_elapsed_time[nTime_dt_array_element] )
01571 {
01572 fprintf(ioQQQ,"DEBUG lgtimes but reset to %.2e\n" ,dynamics.timestep );
01573 dynamics.timestep = 1.0001*(time_elapsed_time[nTime_dt_array_element]-dynamics.time_elapsed);
01574 }
01575 }
01576 else
01577 {
01578
01579 dynamics.timestep = timestep_next();
01580 }
01581
01582 if( cosmology.lgDo )
01583 {
01584 cosmology.redshift_step = -(1.f + cosmology.redshift_current ) *
01585 GetHubbleFactor( cosmology.redshift_current ) * (realnum)dynamics.timestep;
01586 cosmology.redshift_current += cosmology.redshift_step;
01587 }
01588
01589 fprintf(ioQQQ,"DEBUG times exit dynamics.timestep %.2e elapsed_time %.2e scale %.2e ",
01590 dynamics.timestep ,
01591 dynamics.time_elapsed,
01592 rfield.time_continuum_scale);
01593
01594 if( cosmology.lgDo )
01595 {
01596 fprintf(ioQQQ,"redshift %.3e ", cosmology.redshift_current );
01597 }
01598
01599 fprintf(ioQQQ,"\n" );
01600 }
01601
01602
01603 ASSERT( (iteration < dynamics.n_initial_relax+1) ||
01604 Dyn_dr != 0. || (Dyn_dr == 0. && wind.lgStatic()) );
01605
01606
01607 ipUpstream = iphUpstream = ipyUpstream = -1;
01608 dynamics.discretization_error = 0.;
01609 dynamics.error_scale2 = 0.;
01610
01611
01612 DynaSaveLast();
01613 return;
01614 }
01615
01616
01617 STATIC void DynaNewStep(void)
01618 {
01619 long int ilast = 0,
01620 i,
01621 nelem,
01622 ion,
01623 mol;
01624
01625 double frac_next=-BIGFLOAT,
01626 Oldi_hden,
01627 Oldi_ion,
01628 Oldi_iso,
01629 Oldi_mol;
01630
01631 DEBUG_ENTRY( "DynaNewStep()" );
01632
01633
01634 dynamics.convergence_error = 0;
01635 dynamics.error_scale1 = 0.;
01636
01637 ASSERT( nzone < struc.nzlim);
01638 for(i=0;i<nzone;++i)
01639 {
01640
01641 while( (Old_depth[ilast] < struc.depth[i] ) &&
01642 ( ilast < nOld_zone-1 ) )
01643 {
01644 ++ilast;
01645 }
01646 ASSERT( ilast <= nOld_zone-1 );
01647
01648 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01649 {
01650 frac_next = ( struc.depth[i] - Old_depth[ilast])/
01651 (Old_depth[ilast+1] - Old_depth[ilast]);
01652 Oldi_hden = Old_hden[ilast] +
01653 (Old_hden[ilast+1] - Old_hden[ilast])*
01654 frac_next;
01655 }
01656 else
01657 {
01658 Oldi_hden = Old_hden[ilast];
01659 }
01660
01661
01662 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01663 {
01664 for( ion=0; ion<nelem+2; ++ion )
01665 {
01666 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01667 {
01668 Oldi_ion = (Old_xIonDense[ilast][nelem][ion] +
01669 (Old_xIonDense[ilast+1][nelem][ion]-Old_xIonDense[ilast][nelem][ion])*
01670 frac_next);
01671 }
01672 else
01673 {
01674 Oldi_ion = Old_xIonDense[ilast][nelem][ion];
01675 }
01676 dynamics.convergence_error += POW2(Oldi_ion/Oldi_hden-struc.xIonDense[nelem][ion][i]/struc.hden[i]) ;
01677
01678
01679 dynamics.error_scale1 += POW2(struc.xIonDense[nelem][ion][i]/struc.hden[i]);
01680 }
01681 }
01682 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01683 {
01684 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01685 {
01686 if( dense.lgElmtOn[nelem] )
01687 {
01688 for( long level=0; level < iso.numLevels_local[ipISO][nelem]; ++level )
01689 {
01690 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01691 {
01692 Oldi_iso = (Old_StatesElemNEW[ilast][nelem][nelem-ipISO][level] +
01693 (Old_StatesElemNEW[ilast+1][nelem][nelem-ipISO][level]-Old_StatesElemNEW[ilast][nelem][nelem-ipISO][level])*
01694 frac_next);
01695 }
01696 else
01697 {
01698 Oldi_iso = Old_StatesElemNEW[ilast][nelem][nelem-ipISO][level];
01699 }
01700 dynamics.convergence_error += POW2(Oldi_iso/Oldi_hden-struc.StatesElemNEW[nelem][nelem-ipISO][level][i]/struc.hden[i]) ;
01701
01702
01703 dynamics.error_scale1 += POW2(struc.StatesElemNEW[nelem][nelem-ipISO][level][i]/struc.hden[i]);
01704 }
01705 }
01706 }
01707 }
01708 for( mol=0; mol < N_H_MOLEC; mol++)
01709 {
01710 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01711 {
01712 Oldi_mol = (Old_H2_molec[ilast][mol] +
01713 (Old_H2_molec[ilast+1][mol]-Old_H2_molec[ilast][mol])*
01714 frac_next);
01715 }
01716 else
01717 {
01718 Oldi_mol = Old_H2_molec[ilast][mol];
01719 }
01720 dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.H2_molec[mol][i]/struc.hden[i]) ;
01721
01722
01723
01724 dynamics.error_scale1 += POW2(struc.H2_molec[mol][i]/struc.hden[i]);
01725 }
01726 for( mol=0; mol < mole.num_comole_calc; mol++)
01727 {
01728 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
01729 {
01730 Oldi_mol = (Old_CO_molec[ilast][mol] +
01731 (Old_CO_molec[ilast+1][mol]-Old_CO_molec[ilast][mol])*
01732 frac_next);
01733 }
01734 else
01735 {
01736 Oldi_mol = Old_CO_molec[ilast][mol];
01737 }
01738 dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.CO_molec[mol][i]/struc.hden[i]);
01739
01740
01741
01742 dynamics.error_scale1 += POW2(struc.CO_molec[mol][i]/struc.hden[i]);
01743 }
01744 }
01745
01746
01747
01748
01749
01750
01751
01752 fprintf(ioQQQ,"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n",
01753 Dyn_dr, dynamics.convergence_error , dynamics.discretization_error ,
01754 dynamics.error_scale1 , dynamics.error_scale2
01755 );
01756
01757
01758 if( dynamics.convergence_error < dynamics.convergence_tolerance*dynamics.discretization_error )
01759 Dyn_dr /= 1.5;
01760 return;
01761 }
01762
01763
01764 STATIC void DynaSaveLast(void)
01765 {
01766 long int i,
01767 ion,
01768 nelem,
01769 mol;
01770
01771 DEBUG_ENTRY( "DynaSaveLast()" );
01772
01773
01774 nOld_zone = nzone;
01775 dynamics.oldFullDepth = struc.depth[nzone-1];
01776 ASSERT( nzone < struc.nzlim );
01777 for( i=0; i<nzone; ++i )
01778 {
01779 Old_histr[i] = struc.histr[i];
01780 Old_depth[i] = struc.depth[i];
01781 Old_xLyman_depth[i] = struc.xLyman_depth[i];
01782
01783 Old_hiistr[i] = struc.hiistr[i];
01784
01785 Old_pressure[i] = struc.pressure[i];
01786
01787 Old_ednstr[i] = struc.ednstr[i];
01788
01789 Old_EnthalpyDensity[i] = EnthalpyDensity[i];
01790
01791 Old_hden[i] = struc.hden[i];
01792 Old_DenMass[i] = struc.DenMass[i];
01793
01794 for(mol=0;mol<N_H_MOLEC;mol++)
01795 {
01796 Old_H2_molec[i][mol] = struc.H2_molec[mol][i];
01797 }
01798 for(mol=0;mol<mole.num_comole_calc;mol++)
01799 {
01800 Old_CO_molec[i][mol] = struc.CO_molec[mol][i];
01801 }
01802
01803 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01804 {
01805 Old_gas_phase[i][nelem] = dense.gas_phase[nelem];
01806 for( ion=0; ion<nelem+2; ++ion )
01807 {
01808 Old_xIonDense[i][nelem][ion] = struc.xIonDense[nelem][ion][i];
01809 }
01810 }
01811 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01812 {
01813 for( nelem=ipISO; nelem<LIMELM; ++nelem)
01814 {
01815 if( dense.lgElmtOn[nelem] )
01816 {
01817 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
01818 {
01819 Old_StatesElemNEW[i][nelem][nelem-ipISO][level] = struc.StatesElemNEW[nelem][nelem-ipISO][level][i];
01820 ASSERT( !isnan( Old_StatesElemNEW[i][nelem][nelem-ipISO][level] ) );
01821 }
01822 }
01823 }
01824 }
01825 }
01826 return;
01827 }
01828
01829 realnum DynaFlux(double depth)
01830
01831 {
01832 realnum flux;
01833
01834 DEBUG_ENTRY( "DynaFlux()" );
01835
01836 if(dynamics.FluxIndex == 0)
01837 {
01838 flux = (realnum)dynamics.FluxScale;
01839 }
01840 else
01841 {
01842 flux = (realnum)(dynamics.FluxScale*pow(fabs(depth-dynamics.FluxCenter),dynamics.FluxIndex));
01843 if(depth < dynamics.FluxCenter)
01844 flux = -flux;
01845 }
01846 if(dynamics.lgFluxDScale)
01847 {
01848
01849
01850 flux *= dense.xMassDensity0;
01851 }
01852 return flux;
01853 }
01854
01855
01856
01857
01858 void DynaZero( void )
01859 {
01860 int ipISO;
01861
01862 DEBUG_ENTRY( "DynaZero()" );
01863
01864
01865 nOld_zone = 0;
01866
01867
01868 dynamics.lgAdvection = false;
01869
01870 AdvecSpecificEnthalpy = 0.;
01871 dynamics.Cool_r = 0.;
01872 dynamics.Heat_v = 0.;
01873 dynamics.dHeatdT = 0.;
01874 dynamics.HeatMax = 0.;
01875 dynamics.CoolMax = 0.;
01876 dynamics.Rate = 0.;
01877
01878
01879 dynamics.lgRecom = false;
01880
01881
01882 dynamics.lgEquilibrium = false;
01883
01884
01885 dynamics.lgStatic_completed = false;
01886
01887
01888 dynamics.lgTimeDependentStatic = false;
01889 dynamics.timestep_init = -1.;
01890
01891 dynamics.timestep_factor = 1.2;
01892 dynamics.time_elapsed = 0.;
01893
01894
01895
01896 dynamics.n_initial_relax = 2;
01897
01898
01899
01900 dynamics.AdvecLengthInit = -0.1;
01901
01902
01903 dynamics.convergence_tolerance = 0.1;
01904
01905
01906 dynamics.lgSetPresMode = false;
01907
01908
01909 dynamics.FluxScale = 0.;
01910 dynamics.lgFluxDScale = false;
01911 dynamics.FluxCenter = 0.;
01912 dynamics.FluxIndex = 0.;
01913 dynamics.dRad = BIGFLOAT;
01914
01915 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01916 {
01917
01918
01919
01920 dynamics.lgISO[ipISO] = true;
01921 }
01922
01923 dynamics.lgMETALS = true;
01924
01925 dynamics.lgCoolHeat = true;
01926 DivergePresInteg = 0.;
01927
01928 dynamics.discretization_error = 0.;
01929 dynamics.error_scale2 = 0.;
01930 return;
01931 }
01932
01933
01934
01935
01936
01937 void DynaCreateArrays( void )
01938 {
01939 long int nelem,
01940 ns,
01941 i,
01942 ion,
01943 mol;
01944
01945 DEBUG_ENTRY( "DynaCreateArrays()" );
01946
01947 Upstream_H2_molec = (double*)MALLOC((size_t)N_H_MOLEC*sizeof(double) );
01948 Upstream_CO_molec = (double*)MALLOC((size_t)mole.num_comole_calc*sizeof(double) );
01949
01950 dynamics.H2_molec = (double*)MALLOC((size_t)N_H_MOLEC*sizeof(double) );
01951 dynamics.CO_molec = (double*)MALLOC((size_t)mole.num_comole_calc*sizeof(double) );
01952
01953 UpstreamElem = (double*)MALLOC((size_t)LIMELM*sizeof(double) );
01954
01955 dynamics.Source = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
01956 UpstreamIon = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) ));
01957 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
01958 {
01959 dynamics.Source[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
01960 UpstreamIon[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) ));
01961 for( ion=0; ion<nelem+2; ++ion )
01962 {
01963 dynamics.Source[nelem][ion] = 0.;
01964 }
01965 }
01966
01967 UpstreamStatesElemNEW = ((double***)MALLOC( (size_t)LIMELM*sizeof(double **) ));
01968 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01969 {
01970 if( dense.lgElmtOn[nelem] )
01971 {
01972 UpstreamStatesElemNEW[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
01973 for( long ion=0; ion<nelem+1; ion++ )
01974 {
01975 long ipISO = nelem-ion;
01976 if( ipISO < NISO )
01977 {
01978 UpstreamStatesElemNEW[nelem][nelem-ipISO] = (double*)MALLOC(sizeof(double)*(unsigned)iso.numLevels_max[ipISO][nelem]);
01979 }
01980 else
01981 {
01982 fixit();
01983 UpstreamStatesElemNEW[nelem][nelem-ipISO] = NULL;
01984 }
01985 }
01986 }
01987 }
01988
01989
01990 dynamics.StatesElemNEW = ((double***)MALLOC( (size_t)LIMELM*sizeof(double **) ));
01991 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
01992 {
01993 if( dense.lgElmtOn[nelem] )
01994 {
01995 dynamics.StatesElemNEW[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
01996 for( long ion=0; ion<nelem+1; ion++ )
01997 {
01998 long ipISO = nelem-ion;
01999 if( ipISO < NISO )
02000 {
02001 dynamics.StatesElemNEW[nelem][nelem-ipISO] = (double*)MALLOC(sizeof(double)*(unsigned)iso.numLevels_max[ipISO][nelem]);
02002 }
02003 else
02004 {
02005 fixit();
02006 dynamics.StatesElemNEW[nelem][nelem-ipISO] = NULL;
02007 }
02008 }
02009 }
02010 }
02011
02012 dynamics.Rate = 0.;
02013
02014 Old_EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02015
02016 Old_ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02017
02018 EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02019
02020 Old_DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02021
02022 Old_hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02023
02024 Old_pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02025
02026 Old_histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02027
02028 Old_hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02029
02030 Old_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02031
02032 Old_xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
02033
02034 Old_xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) );
02035
02036 Old_StatesElemNEW = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(struc.nzlim) );
02037
02038 Old_gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
02039
02040 Old_H2_molec = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
02041
02042 Old_CO_molec = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
02043
02044
02045 for( ns=0; ns < struc.nzlim; ++ns )
02046 {
02047 Old_xIonDense[ns] =
02048 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+3) );
02049
02050 Old_StatesElemNEW[ns] =
02051 (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(LIMELM) );
02052
02053 Old_gas_phase[ns] =
02054 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+3) );
02055
02056 Old_H2_molec[ns] =
02057 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(N_H_MOLEC) );
02058
02059 Old_CO_molec[ns] =
02060 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(mole.num_comole_calc) );
02061
02062 for( nelem=0; nelem< (LIMELM+3);++nelem )
02063 {
02064 Old_xIonDense[ns][nelem] =
02065 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+1) );
02066 }
02067
02068 for( nelem=0; nelem< LIMELM; ++nelem )
02069 {
02070 if( dense.lgElmtOn[nelem] )
02071 {
02072 Old_StatesElemNEW[ns][nelem] =
02073 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+1) );
02074 for( ion=0; ion<nelem+1; ion++ )
02075 {
02076 long ipISO = nelem-ion;
02077 if( ipISO < NISO )
02078 {
02079 Old_StatesElemNEW[ns][nelem][ion] =
02080 (realnum*)MALLOC(sizeof(realnum)*(unsigned)iso.numLevels_max[ipISO][nelem]);
02081 }
02082 else
02083 {
02084 fixit();
02085 Old_StatesElemNEW[ns][nelem][ion] = NULL;
02086 }
02087 }
02088 }
02089 }
02090 }
02091
02092 for( i=0; i < struc.nzlim; i++ )
02093 {
02094
02095 Old_histr[i] = 0.;
02096 Old_xLyman_depth[i] = 0.;
02097 Old_depth[i] = 0.;
02098 dynamics.oldFullDepth = 0.;
02099
02100 Old_hiistr[i] = 0.;
02101
02102 Old_pressure[i] = 0.;
02103
02104 Old_ednstr[i] = 0.;
02105 Old_hden[i] = 0.;
02106 Old_DenMass[i] = 0.;
02107 Old_EnthalpyDensity[i] = 0.;
02108 for( nelem=0; nelem< (LIMELM+3);++nelem )
02109 {
02110 for( ion=0; ion<LIMELM+1; ++ion )
02111 {
02112 Old_xIonDense[i][nelem][ion] = 0.;
02113 }
02114 }
02115 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
02116 {
02117 for( nelem=ipISO; nelem<LIMELM; ++nelem)
02118 {
02119 if( dense.lgElmtOn[nelem] )
02120 {
02121 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
02122 {
02123 Old_StatesElemNEW[i][nelem][nelem-ipISO][level] = 0.;
02124 }
02125 }
02126 }
02127 }
02128 for(mol=0;mol<N_H_MOLEC;mol++)
02129 {
02130 Old_H2_molec[i][mol] = 0.;
02131 }
02132 for(mol=0;mol<mole.num_comole_calc;mol++)
02133 {
02134 Old_CO_molec[i][mol] = 0.;
02135 }
02136 }
02137 return;
02138 }
02139
02140
02141
02142
02143 STATIC void advection_set_default( bool lgWind )
02144 {
02145
02146 DEBUG_ENTRY( "advection_set_default()" );
02147
02148
02149 dynamics.lgAdvection = true;
02150
02151
02152
02153 thermal.lgPredNextTe = false;
02154
02155
02156
02157
02158
02159
02161
02162 co.lgNoCOMole = true;
02163 # if 0
02164 co.lgNoCOMole = true;
02165 phycon.lgPhysOK = false;
02166 # endif
02167
02168
02171
02172
02173
02174
02175
02176 pressure.lgPres_radiation_ON = true;
02177 pressure.lgPres_magnetic_ON = true;
02178 pressure.lgPres_ram_ON = true;
02179
02180
02181 if( lgWind )
02182 {
02183
02184 conv.EdenErrorAllowed = 1e-3;
02185
02186
02187
02188
02189
02190 conv.HeatCoolRelErrorAllowed = 3e-4f;
02191 conv.PressureErrorAllowed = 1e-3f;
02192
02193 if( cosmology.lgDo )
02194 {
02195 conv.EdenErrorAllowed = 1e-5;
02196 conv.PressureErrorAllowed = 1e-5f;
02197 }
02198 }
02199 return;
02200 }
02201
02202
02203
02204 void ParseDynaTime( Parser &p )
02205 {
02206 DEBUG_ENTRY( "ParseDynaTime()" );
02207
02208
02209 dynamics.lgTimeDependentStatic = true;
02210
02211 dynamics.timestep_init = p.getNumberCheckAlwaysLogLim("dynamics.timestep",30.);
02212
02213 dynamics.timestep = dynamics.timestep_init;
02214 if( p.nMatch( "TRAC" ) )
02215 dynamics.lgTracePrint = true;
02216
02217
02218 dynamics.timestep_stop = p.getNumberDefaultAlwaysLog("stop time", -1.);
02219
02220
02221 advection_set_default(false);
02222
02223 wind.windv0 = 0.;
02224 wind.setStatic();
02225 wind.windv = wind.windv0;
02226
02227
02228
02229
02230
02231 time_elapsed_time = (double*)MALLOC((size_t)NTIME*sizeof(double));
02232 time_flux_ratio = (double*)MALLOC((size_t)NTIME*sizeof(double));
02233 time_dt = (double*)MALLOC((size_t)NTIME*sizeof(double));
02234 time_dt_scale_factor = (double*)MALLOC((size_t)NTIME*sizeof(double));
02235 lgtime_Recom = (int*)MALLOC((size_t)NTIME*sizeof(int));
02236
02237
02238 nTime_flux = 0;
02239
02240
02241 p.getline();
02242 if( p.m_lgEOF )
02243 {
02244 fprintf( ioQQQ,
02245 " Hit EOF while reading time-continuum list; use END to end list.\n" );
02246 cdEXIT(EXIT_FAILURE);
02247 }
02248
02249
02250
02251 lgtime_dt_specified = true;
02252
02253 while( p.strcmp("END") != 0 )
02254 {
02255 if( nTime_flux >= NTIME )
02256 {
02257 fprintf( ioQQQ,
02258 " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n",
02259 NTIME );
02260 cdEXIT(EXIT_FAILURE);
02261 }
02262
02263 if( p.nMatch("CYCLE") )
02264 {
02265 double period = p.getNumberCheckAlwaysLog("log time");
02266 ASSERT( period > time_elapsed_time[nTime_flux-1] );
02267 long pointsPerPeriod = nTime_flux;
02268 while( nTime_flux < NTIME - 1 )
02269 {
02270 time_elapsed_time[nTime_flux] = period + time_elapsed_time[nTime_flux-pointsPerPeriod];
02271 time_flux_ratio[nTime_flux] = time_flux_ratio[nTime_flux-pointsPerPeriod];
02272 time_dt[nTime_flux] = time_dt[nTime_flux-pointsPerPeriod];
02273 time_dt_scale_factor[nTime_flux] = time_dt_scale_factor[nTime_flux-pointsPerPeriod];
02274 nTime_flux++;
02275 }
02276
02277 fprintf( ioQQQ, " Adding cycles with period = %e s.\n", period );
02278
02279
02280 p.getline();
02281 if( p.m_lgEOF )
02282 {
02283 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
02284 cdEXIT(EXIT_FAILURE);
02285 }
02286 continue;
02287 }
02288
02289 time_elapsed_time[nTime_flux] = p.getNumberCheckAlwaysLog("log time");
02290 time_flux_ratio[nTime_flux] = p.getNumberCheckAlwaysLog("log flux ratio");
02291
02292
02293
02294 time_dt[nTime_flux] = p.getNumberDefaultAlwaysLog("log time step",-1.);
02295
02296
02297 if( time_dt[nTime_flux] < 0.0 )
02298 lgtime_dt_specified = false;
02299
02300
02301 time_dt_scale_factor[nTime_flux] = p.getNumberDefaultAlwaysLog(
02302 "scale factor to increase time",-1.);
02303
02304
02305 if( p.nMatch("RECOMBIN") )
02306 {
02307
02308
02309 lgtime_Recom[nTime_flux] = true;
02310 }
02311 else
02312 {
02313 lgtime_Recom[nTime_flux] = false;
02314 }
02315
02316
02317 ++nTime_flux;
02318
02319
02320 p.getline();
02321 if( p.m_lgEOF )
02322 {
02323 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
02324 cdEXIT(EXIT_FAILURE);
02325 }
02326
02327 }
02328
02329 if( nTime_flux < 2 )
02330 {
02331 fprintf( ioQQQ, " At least two instances of time must be specified. There is an implicit instance at t=0.\n" \
02332 " The user must specify at least one additional time. Sorry.\n" );
02333 cdEXIT(EXIT_FAILURE);
02334 }
02335
02336 for( long i=0; i < nTime_flux; i++ )
02337 {
02338 fprintf( ioQQQ, "DEBUG time dep %.2e %.2e %.2e %.2e\n",
02339 time_elapsed_time[i],
02340 time_flux_ratio[i] ,
02341 time_dt[i],
02342 time_dt_scale_factor[i]);
02343 }
02344 fprintf( ioQQQ, "\n" );
02345 return;
02346 }
02347
02348
02349 void ParseDynaWind( Parser &p)
02350 {
02351 int iVelocity_Type;
02352 bool lgModeSet=false;
02353
02354
02355 double dfdr=-BIGDOUBLE;
02356
02357 DEBUG_ENTRY( "ParseDynaWind()" );
02358
02359 if( p.nMatch( "TRAC" ) )
02360 dynamics.lgTracePrint = true;
02361
02362
02363
02364
02365
02366 iVelocity_Type = 0;
02367
02368
02369 if( p.nMatch( "VELO" ) )
02370 {
02371 wind.windv0 = (realnum)(p.getNumberPlain("velocity")*1e5);
02372 wind.windv = wind.windv0;
02373 wind.setDefault();
02374 iVelocity_Type = 1;
02375 }
02376
02377 if( p.nMatch( "BALL" ) )
02378 {
02379 wind.setBallistic();
02380 lgModeSet = true;
02381 }
02382
02383 if( p.nMatch( "STAT" ) )
02384 {
02385 wind.windv0 = 0.;
02386 wind.setStatic();
02387 lgModeSet = true;
02388 iVelocity_Type = 1;
02389 }
02390
02391 if ( 1 == iVelocity_Type && !lgModeSet)
02392 {
02393 if (wind.windv0 > 0.)
02394 {
02395 fprintf(ioQQQ,"Warning, BALListic option needed to switch off pressure gradient terms\n");
02396 }
02397 else if (wind.windv0 == 0.)
02398 {
02399 fprintf(ioQQQ,"Warning, STATic option needed for zero speed solutions\n");
02400 }
02401 }
02402
02403 if( p.nMatch("DFDR") )
02404 {
02405
02406 dfdr = p.getNumberPlain("flux gradient");
02407 iVelocity_Type = 2;
02408 }
02409
02410
02411 if( p.nMatch("CENT") )
02412 {
02413
02414 dynamics.FluxCenter = p.getNumberPlain(
02415 "centre of mass flux distribution");
02416 }
02417
02418
02419 if( p.nMatch("INDE") )
02420 {
02421
02422 dynamics.FluxIndex = p.getNumberPlain(
02423 "power law index of mass flux distribution");
02424 }
02425
02426
02427 if(iVelocity_Type == 1)
02428 {
02429
02430 if(dynamics.FluxIndex == 0)
02431 {
02432 dynamics.FluxScale = wind.windv0;
02433 dynamics.lgFluxDScale = true;
02434
02435
02436
02437 dynamics.FluxCenter = -1.;
02438 }
02439 else
02440 {
02443
02444 dynamics.FluxScale = wind.windv0*
02445 pow(fabs(dynamics.FluxCenter),-dynamics.FluxIndex);
02446
02447 dynamics.lgFluxDScale = true;
02448 if(dynamics.FluxCenter > 0)
02449 {
02450 dynamics.FluxScale = -dynamics.FluxScale;
02451 }
02452 }
02453 }
02454
02455 else if(iVelocity_Type == 2)
02456 {
02457 if(dynamics.FluxIndex == 0)
02458 {
02459 fprintf(ioQQQ,"Can't specify gradient when flux is constant!\n");
02460
02461 cdEXIT(EXIT_FAILURE);
02462 }
02465
02466
02467 dynamics.FluxScale = dfdr/dynamics.FluxIndex*
02468 pow(fabs(dynamics.FluxCenter),1.-dynamics.FluxIndex);
02469 if(dynamics.FluxCenter > 0)
02470 {
02471 dynamics.FluxScale = -dynamics.FluxScale;
02472 }
02473 dynamics.lgFluxDScale = false;
02474
02475
02476
02477 wind.windv0 = -0.01f;
02478 wind.setDefault();
02479 }
02480 else
02481 {
02482
02483
02484
02485 wind.windv0 = (realnum)(p.getNumberCheck("wind velocity")*1e5);
02486 if (wind.windv0 < 0.)
02487 {
02488 wind.setDefault();
02489 }
02490 else if (wind.windv0 > 0.)
02491 {
02492 wind.setBallistic();
02493 }
02494 else
02495 {
02496 wind.setStatic();
02497 }
02498
02499 dynamics.FluxScale = wind.windv0;
02500 dynamics.FluxIndex = 0.;
02501 dynamics.lgFluxDScale = true;
02502
02503
02504
02505 dynamics.FluxCenter = -1.;
02506 }
02507
02508 wind.windv = wind.windv0;
02509
02510 # ifdef FOO
02511 fprintf(ioQQQ,"Scale %g (*%c) Index %g Center %g\n",
02512 dynamics.FluxScale,(dynamics.lgFluxDScale)?'D':'1',
02513 dynamics.FluxIndex,dynamics.FluxCenter);
02514 # endif
02515
02516
02517 if( p.nMatch( "ADVE" ) )
02518 {
02519
02520 advection_set_default(true);
02521 }
02522
02523 else
02524 {
02525
02526 if( wind.windv0 <= 1.e6 )
02527 {
02528
02529 fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" );
02530 wind.lgWindOK = false;
02531 }
02532
02533
02534 wind.comass = (realnum)p.getNumberDefault("central object mass",1.);
02535
02536
02537
02538 wind.lgDisk = false;
02539 if( p.nMatch( "DISK") )
02540 wind.lgDisk = true;
02541
02542 }
02543
02544
02545 strcpy( dense.chDenseLaw, "WIND" );
02546
02547
02548 if( p.nMatch("NO CO") )
02549 {
02550 pressure.lgContRadPresOn = false;
02551 }
02552 else
02553 {
02554 pressure.lgContRadPresOn = true;
02555 }
02556 return;
02557 }
02558
02559
02560 void DynaPrtZone( void )
02561 {
02562
02563 DEBUG_ENTRY( "DynaPrtZone()" );
02564
02565 ASSERT( nzone>0 && nzone<struc.nzlim );
02566
02567 if( nzone > 0 )
02568 {
02569 fprintf(ioQQQ," DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n",
02570 timesc.sound_speed_adiabatic/1e5 ,
02571 wind.windv/1e5 ,
02572 dynamics.Cool()/thermal.ctot,
02573 dynamics.Heat()/thermal.ctot);
02574 }
02575
02576 ASSERT( EnthalpyDensity[nzone-1] > 0. );
02577
02578 fprintf(ioQQQ," DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n",
02579 phycon.EnergyExcitation,
02580 phycon.EnergyIonization,
02581 phycon.EnergyBinding,
02582 0.5*POW2(wind.windv)*dense.xMassDensity,
02583 5./2.*pressure.PresGasCurr ,
02584 EnthalpyDensity[nzone-1]/dense.gas_phase[ipHYDROGEN] , AdvecSpecificEnthalpy
02585 );
02586 return;
02587 }
02588
02589
02590 void DynaPunchTimeDep( FILE* ipPnunit , const char *chJob )
02591 {
02592
02593 DEBUG_ENTRY( "DynaPunchTimeDep()" );
02594
02595 if( strcmp( chJob , "END" ) == 0 )
02596 {
02597 double te_mean,
02598 H2_mean,
02599 H0_mean,
02600 Hp_mean,
02601 Hep_mean;
02602
02603 if( cdTemp(
02604
02605 "HYDR",
02606
02607
02608 2,
02609
02610 &te_mean,
02611
02612 "RADIUS" ) )
02613 {
02614 TotalInsanity();
02615 }
02616 if( cdIonFrac(
02617
02618 "HYDR",
02619
02620
02621 2,
02622
02623 &Hp_mean,
02624
02625 "RADIUS" ,
02626
02627 false ) )
02628 {
02629 TotalInsanity();
02630 }
02631 if( cdIonFrac(
02632
02633 "HYDR",
02634
02635
02636 1,
02637
02638 &H0_mean,
02639
02640 "RADIUS" ,
02641
02642 false ) )
02643 {
02644 TotalInsanity();
02645 }
02646 if( cdIonFrac(
02647
02648 "H2 ",
02649
02650
02651 0,
02652
02653 &H2_mean,
02654
02655 "RADIUS" ,
02656
02657 false ) )
02658 {
02659 TotalInsanity();
02660 }
02661 if( cdIonFrac(
02662
02663 "HELI",
02664
02665
02666 2,
02667
02668 &Hep_mean,
02669
02670 "RADIUS" ,
02671
02672 false ) )
02673 {
02674 TotalInsanity();
02675 }
02676 fprintf( ipPnunit ,
02677 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" ,
02678 dynamics.time_elapsed ,
02679 dynamics.timestep ,
02680 rfield.time_continuum_scale ,
02681 dense.gas_phase[ipHYDROGEN],
02682 te_mean ,
02683 Hp_mean ,
02684 H0_mean ,
02685 H2_mean ,
02686 Hep_mean ,
02687
02688 findspecies("CO")->hevcol / SDIV( colden.colden[ipCOL_HTOT] ),
02689 cosmology.redshift_current,
02690 dense.eden/dense.gas_phase[ipHYDROGEN]
02691 );
02692 }
02693 else
02694 TotalInsanity();
02695 return;
02696 }
02697
02698
02699 void DynaSave(FILE* ipPnunit , char chJob )
02700 {
02701 DEBUG_ENTRY( "DynaSave()" );
02702
02703 if( chJob=='a' )
02704 {
02705
02706 fprintf( ipPnunit , "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
02707 radius.depth_mid_zone,
02708 thermal.htot ,
02709 dynamics.Cool() ,
02710 dynamics.Heat() ,
02711 dynamics.dCooldT() ,
02712 dynamics.Source[ipHYDROGEN][ipHYDROGEN],
02713 dynamics.Rate,
02714 phycon.EnthalpyDensity/dense.gas_phase[ipHYDROGEN] ,
02715 AdvecSpecificEnthalpy
02716 );
02717 }
02718 else
02719 TotalInsanity();
02720 return;
02721 }
02722
02723 double t_dynamics::Heat()
02724 {
02725 return Heat_v*dense.gas_phase[ipHYDROGEN];
02726 }
02727
02728 double t_dynamics::Cool()
02729 {
02730 return Cool_r*phycon.EnthalpyDensity;
02731 }
02732
02733 double t_dynamics::dCooldT()
02734 {
02735 return Cool_r*5./2.*pressure.PresGasCurr/phycon.te;
02736 }
02737
02738 void DynaIterStart(void)
02739 {
02740 DEBUG_ENTRY( "DynaIterStart()" );
02741
02742 if( 0 == nTime_flux )
02743 {
02744 rfield.time_continuum_scale = 1.;
02745 return;
02746 }
02747 else if( dynamics.time_elapsed <= time_elapsed_time[0] )
02748 {
02749
02750 rfield.time_continuum_scale = (realnum)time_flux_ratio[0];
02751 }
02752 else if( dynamics.time_elapsed > time_elapsed_time[nTime_flux-1] )
02753 {
02754 fprintf( ioQQQ,
02755 " PROBLEM - DynaIterStart - I need the continuum at time %.2e but the table ends at %.2e.\n" ,
02756 dynamics.time_elapsed ,
02757 time_elapsed_time[nTime_flux-1]);
02758 cdEXIT(EXIT_FAILURE);
02759 }
02760 else
02761 {
02762 rfield.time_continuum_scale = (realnum)linint(
02763
02764 time_elapsed_time,
02765
02766 time_flux_ratio,
02767
02768 nTime_flux,
02769
02770 dynamics.time_elapsed);
02771 }
02772
02773 fprintf(ioQQQ,"DEBUG time dep reset continuum iter %ld dynamics.timestep %.2e elapsed time %.2e scale %.2e",
02774 iteration,
02775 dynamics.timestep ,
02776 dynamics.time_elapsed,
02777 rfield.time_continuum_scale);
02778 if( dynamics.lgRecom )
02779 {
02780 fprintf(ioQQQ," recom");
02781 }
02782 fprintf(ioQQQ,"\n");
02783
02784
02785 long int nTimeVary = 0;
02786 for( long int i=0; i < rfield.nShape; i++ )
02787 {
02788
02789
02790 if( rfield.lgTimeVary[i] )
02791 ++nTimeVary;
02792 }
02793
02794 if( hextra.lgTurbHeatVaryTime )
02795 {
02796
02797 hextra.TurbHeat = hextra.TurbHeatSave * rfield.time_continuum_scale;
02798 fprintf(ioQQQ,"DEBUG TurbHeat vary new heat %.2e\n",
02799 hextra.TurbHeat);
02800 }
02801 else if( !nTimeVary )
02802 {
02803 fprintf(ioQQQ," DISASTER - there were no variable continua "
02804 "or heat sources - put TIME option on at least one "
02805 "luminosity or hextra command.\n");
02806 cdEXIT(EXIT_FAILURE);
02807 }
02808 }
02809