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