cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /* DynaEndIter called at end of iteration when advection is turned on */ 00004 /* DynaStartZone called at start of zone calculation when advection is turned on */ 00005 /* DynaEndZone called at end of zone calculation when advection is turned on */ 00006 /* DynaIonize, called from ionize to evaluate advective terms for current conditions */ 00007 /* DynaPresChngFactor, called from PressureChange to evaluate new density needed for 00008 * current conditions and wind solution, returns ratio of new to old density */ 00009 /* timestep_next - find next time step in time dependent case */ 00010 /* DynaZero zero some dynamics variables, called from zero.c */ 00011 /* DynaCreateArrays allocate some space needed to save the dynamics structure variables, 00012 * called from DynaCreateArrays */ 00013 /* DynaPrtZone - called to print zone results */ 00014 /* DynaPunch punch info related to advection */ 00015 /* DynaPunch, punch output for dynamics solutions */ 00016 /* ParseDynaTime parse the time command, called from ParseCommands */ 00017 /* ParseDynaWind parse the wind command, called from ParseCommands */ 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 "hextra.h" 00026 #include "rfield.h" 00027 #include "iterations.h" 00028 #include "trace.h" 00029 #include "conv.h" 00030 #include "timesc.h" 00031 #include "dense.h" 00032 #include "mole.h" 00033 #include "thermal.h" 00034 #include "pressure.h" 00035 #include "phycon.h" 00036 #include "wind.h" 00037 #include "hmi.h" 00038 #include "dynamics.h" 00039 static int ipUpstream=0,iphUpstream=0,ipyUpstream=0; 00040 00041 /* 00042 * >>chng 01 mar 16, incorporate advection within dynamical solutions 00043 * this file contains the routines that incorporate effects of dynamics and advection upon the 00044 * thermal and ionization solutions. 00045 * 00046 * This code was originally developed in March 2001 during 00047 * a visit to UNAM Morelia, in collaboration with Robin Williams, Jane Arthur, and Will Henney. 00048 * Development was continued by email, and in a meeting in July/Aug 2001 at U Cardiff 00049 * Further development June 2002, UNAM Morelia (Cloudy and the Duendes Verdes) 00050 */ 00051 00052 /* turn this on for dynamics printout */ 00053 #define DIAG_PRINT false 00054 #define MAINPRINT false 00055 00056 /* the interpolated upstream densities of all ionization stages, 00057 * [element][ion] */ 00058 static double **UpstreamIon; 00059 /* total abundance of each element per hydrogen */ 00060 static double *UpstreamElem; 00061 00062 /* hydrogen molecules */ 00063 static double *Upstream_H2_molec; 00064 00065 /* CO molecules */ 00066 static double *Upstream_CO_molec; 00067 00068 /* space to save continuum when time command is used 00069 static realnum *dyna_flux_save;*/ 00070 00071 /* initial timestep (seconds) read off time command, 00072 * each iteration accounts for this much time */ 00073 static double timestep_init, 00074 timestep, 00075 timestep_stop, 00076 timestep_factor; 00077 00078 /* array of times and continuum values we will interpolate upon */ 00079 static double *time_elapsed_time , 00080 *time_flux_ratio , 00081 *time_dt, 00082 *time_dt_scale_factor; 00083 bool lgtime_dt_specified; 00084 int *lgtime_Recom; 00085 #define NTIME 200 00086 00087 /* number of time steps actually read in */ 00088 static long int nTime_flux; 00089 00090 /* routine called at end of iteration to determine new step sizes */ 00091 STATIC void DynaNewStep(void); 00092 00093 /* routine called at end of iteration to save values in previous iteration */ 00094 STATIC void DynaSaveLast(void); 00095 00096 /* routine called to determine mass flux at given distance */ 00097 /* static realnum DynaFlux(double depth); */ 00098 00099 /* lookahead distance, as derived in DynaEndIter */ 00100 static double Dyn_dr; 00101 00102 /* advected work */ 00103 static double AdvecSpecificEnthalpy; 00104 00105 /* the upstream hydrogen atoms per unit vol*/ 00106 static double Upstream_hden; 00107 00108 /* HI ionization structure from previous iteration */ 00109 static realnum *Old_histr/*[NZLIM]*/ , 00110 /* Lyman continuum optical depth from previous iteration */ 00111 *Old_xLyman_depth/*[NZLIM]*/, 00112 /* depth of position in previous iteration */ 00113 *Old_depth/*[NZLIM]*/, 00114 /* old n_p density from previous iteration */ 00115 *Old_hiistr/*[NZLIM]*/, 00116 /* old pressure from previous iteration */ 00117 *Old_pressure/*[NZLIM]*/, 00118 /* H density - particles per unit vol */ 00119 *Old_hden/*[NZLIM]*/ , 00120 /* density - total grams per unit vol */ 00121 *Old_DenMass/*[NZLIM]*/ , 00122 /* sum of enthalpy and kinetic energy per gram */ 00123 *EnthalpyDensity/*[NZLIM]*/, 00124 /* old electron density from previous iteration */ 00125 *Old_ednstr/*[NZLIM]*/, 00126 /* sum of enthalpy and kinetic energy per gram */ 00127 *Old_EnthalpyDensity/*[NZLIM]*/; 00128 00129 static realnum **Old_H2_molec; 00130 static realnum **Old_CO_molec; 00131 00132 /* the ionization fractions from the previous iteration */ 00133 static realnum ***Old_xIonDense; 00134 00135 /* the gas phase abundances from the previous iteration */ 00136 static realnum **Old_gas_phase; 00137 00138 /* the number of zones that were saved in the previous iteration */ 00139 static long int nOld_zone; 00140 00141 /* the number of zones that were saved in the previous iteration */ 00142 static realnum DivergePresInteg; 00143 00144 /*timestep_next - find next time step in time dependent case */ 00145 STATIC double timestep_next( void ) 00146 { 00147 static double te_old=-1; 00148 double timestep_Hp_temp , timestep_return; 00149 00150 DEBUG_ENTRY( "timestep_next()" ); 00151 00152 timestep_return = timestep; 00153 00154 if( dynamics.lgRecom ) 00155 { 00156 double te_new; 00157 if( cdTemp( 00158 /* four char string, null terminzed, giving the element name */ 00159 "hydr", 00160 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */ 00161 2 , 00162 /* will be temperature */ 00163 &te_new, 00164 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 00165 "VOLUME" ) ) 00166 TotalInsanity(); 00167 00168 if( te_old>0 ) 00169 { 00170 double dTdStep = fabs(te_new-te_old)/te_new; 00171 /* this is the change factor to get 0.1 frac change in mean temp */ 00172 double dT_factor = 0.04/SDIV(dTdStep); 00173 dT_factor = MIN2( 2.0 , dT_factor ); 00174 dT_factor = MAX2( 0.01 , dT_factor ); 00175 timestep_Hp_temp = timestep * dT_factor; 00176 } 00177 else 00178 { 00179 timestep_Hp_temp = -1.; 00180 } 00181 te_old = te_new; 00182 if( timestep_Hp_temp > 0. ) 00183 timestep_return = timestep_Hp_temp; 00184 } 00185 else 00186 { 00187 timestep_return = timestep_init; 00188 } 00189 fprintf(ioQQQ,"DEBUG timestep_next returns %.3e, old temp %.2e\n" , timestep_return, te_old ); 00190 00191 return( timestep_return ); 00192 } 00193 00194 /* ============================================================================== */ 00195 /* DynaPresChngFactor, called from PressureChange to evaluate factor needed 00196 * to find new density needed for 00197 * current conditions and wind solution, returns ratio of new to old density, 00198 * called when wind velocity is negative */ 00199 00200 /* object is to get the local ram pressure 00201 * RamNeeded = pressure.PresTotlInit + pressure.PresGasCurr + pressure.PresInteg; 00202 * to equal the sum of the inital pressur at the illuminated face, the local gas pressure, 00203 * and the integrate radiative acceleration from the incident continuum 00204 * 00205 * the local gas pressure is linear in density if the temperature is constant, 00206 * 00207 * the local ram pressure is inversely linear in density because of the relationship 00208 * between velocity and density introduced by the mass flux conservation 00209 */ 00210 #define SUBSONIC 1 00211 #define SUPERSONIC 2 00212 /*#define FREE 3*/ 00213 #define STRONGD 4 00214 #define ORIGINAL 5 00215 #define SHOCK 6 00216 #define ANTISHOCK 7 00217 #define ANTISHOCK2 8 00218 00219 double DynaPresChngFactor(void) 00220 { 00221 double factor, 00222 er, 00223 width; 00224 static double dp = -1., 00225 dpp = -1., 00226 erp = -1., 00227 erpp = -1.; 00230 static int lastzone = -1, 00231 zonepresmode, 00232 globalpresmode; 00233 int ipPRE; 00234 00235 DEBUG_ENTRY( "DynaPresChngFactor()" ); 00236 00237 /* update the current pressure and radiative acceleration */ 00238 PresTotCurrent(); 00239 00240 /* update the current desired pressure */ 00241 pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg*pressure.lgContRadPresOn 00242 + DivergePresInteg; 00243 if( trace.lgTrace ) 00244 { 00245 fprintf( ioQQQ, 00246 " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n", 00247 pressure.PresTotlCorrect , pressure.PresTotlInit , pressure.PresInteg*pressure.lgContRadPresOn, 00248 DivergePresInteg ); 00249 } 00250 00251 if( DIAG_PRINT) 00252 fprintf(stdout,"Pressure: init %g +rad %g +diverge %g = %g cf %g\n", 00253 pressure.PresTotlInit, pressure.PresInteg*pressure.lgContRadPresOn, 00254 DivergePresInteg, pressure.PresTotlCorrect, pressure.PresTotlCurr); 00255 /*fprintf(ioQQQ,"DEBUG\t%.2f\thden\t%.4e\tPtot\t%.4e\tPgas\t%.4e\n", 00256 fnzone, dense.gas_phase[ipHYDROGEN],pressure.PresTotlCurr,pressure.PresGasCurr );*/ 00257 00258 /* dynamics.lgSetPresMode is flag to indicate sane value of dynamics.chPresMode. 00259 * If set true with SET DYNAMICS PRESSURE MODE 00260 * then do not override that choice */ 00261 if( !dynamics.lgSetPresMode ) 00262 { 00263 /* above set true if pressure mode was set with 00264 * SET DYNAMICS PRESSURE MODE - if we got here 00265 * it was not set, and must make a guess */ 00266 if(pressure.PresGasCurr < pressure.PresRamCurr) 00267 strcpy( dynamics.chPresMode , "supersonic" ); 00268 else 00269 strcpy( dynamics.chPresMode , "subsonic" ); 00270 /* clear the flag - pressure mode has been set */ 00271 dynamics.lgSetPresMode = true; 00272 } 00273 00274 /* if globally looking for transonic solution, then locally sometimes 00275 * supersonic, sometimes subsonic - this branch sets global flag, 00276 * which can also be set with SET DYNAMICS PRESSURE MODE. 00277 * Under default conditions, ChPresMode was set in previous branch 00278 * to sub or super sonic depending on current velocity on first time*/ 00279 if( strcmp( dynamics.chPresMode , "original" ) == 0 ) 00280 { 00281 globalpresmode = ORIGINAL; 00282 pressure.lgSonicPointAbortOK = true; 00283 } 00284 else if( strcmp( dynamics.chPresMode , "subsonic" ) == 0 ) 00285 { 00286 globalpresmode = SUBSONIC; 00287 pressure.lgSonicPointAbortOK = true; 00288 } 00289 /* supersonic */ 00290 else if( strcmp( dynamics.chPresMode , "supersonic" ) == 0 ) 00291 { 00292 globalpresmode = SUPERSONIC; 00293 pressure.lgSonicPointAbortOK = true; 00294 } 00295 /* strong d */ 00296 else if( strcmp( dynamics.chPresMode , "strongd" ) == 0 ) 00297 { 00298 globalpresmode = STRONGD; 00299 pressure.lgSonicPointAbortOK = false; 00300 } 00301 else if( strcmp( dynamics.chPresMode , "shock" ) == 0 ) 00302 { 00303 globalpresmode = SHOCK; 00304 pressure.lgSonicPointAbortOK = false; 00305 } 00306 else if( strcmp( dynamics.chPresMode , "antishock-by-mach" ) == 0 ) 00307 { 00308 globalpresmode = ANTISHOCK2; 00309 pressure.lgSonicPointAbortOK = false; 00310 } 00311 else if( strcmp( dynamics.chPresMode , "antishock" ) == 0 ) 00312 { 00313 globalpresmode = ANTISHOCK; 00314 pressure.lgSonicPointAbortOK = false; 00315 } 00316 00317 /* this sets pressure mode for the current zone based on global mode 00318 * and local conditions */ 00319 if(globalpresmode == ORIGINAL) 00320 { 00321 /* in this mode use comparison between ram and gas pressure to determine 00322 * whether sub or super sonic */ 00323 if(pressure.PresGasCurr > pressure.PresRamCurr) 00324 { 00325 zonepresmode = SUBSONIC; 00326 } 00327 else 00328 { 00329 zonepresmode = SUPERSONIC; 00330 } 00331 } 00332 else if(globalpresmode == STRONGD) 00333 { 00334 if(nzone <= 1) 00335 zonepresmode = SUPERSONIC; 00336 } 00337 else if(globalpresmode == SUBSONIC) 00338 { 00339 zonepresmode = SUBSONIC; 00340 } 00341 else if(globalpresmode == SUPERSONIC) 00342 { 00343 zonepresmode = SUPERSONIC; 00344 } 00345 else if(globalpresmode == SHOCK) 00346 { 00347 if(radius.depth < dynamics.ShockDepth) 00348 { 00349 zonepresmode = SUBSONIC; 00350 } 00351 else 00352 { 00353 zonepresmode = SUPERSONIC; 00354 } 00355 } 00356 else if(globalpresmode == ANTISHOCK) 00357 { 00358 if(radius.depth < dynamics.ShockDepth) 00359 { 00360 zonepresmode = SUPERSONIC; 00361 } 00362 else 00363 { 00364 zonepresmode = SUBSONIC; 00365 } 00366 } 00367 else if(globalpresmode == ANTISHOCK2) 00368 { 00369 /* WJH 19 May 2004: This version of the antishock mode will 00370 insert the antishock at the point where the isothermal Mach 00371 number falls below a certain value, dynamics.ShockMach */ 00372 if( fabs(wind.windv) > dynamics.ShockMach*sqrt(pressure.PresGasCurr/dense.xMassDensity)) 00373 { 00374 zonepresmode = SUPERSONIC; 00375 } 00376 else 00377 { 00378 zonepresmode = SUBSONIC; 00379 } 00380 } 00381 else 00382 { 00383 printf("Need to set global pressure mode\n"); 00384 cdEXIT( EXIT_FAILURE ); 00385 } 00386 00387 er = pressure.PresTotlCurr-pressure.PresTotlCorrect; 00388 /* fprintf(ioQQQ,"Ds %ld: %ld; %g %g %g %g %g %g %g %g\n",iteration,nzone,dense.gas_phase[ipHYDROGEN],er,pressure.PresTotlCorrect,pressure.PresTotlCurr,phycon.te,thermal.ctot,thermal.htot,phycon.EnthalpyDensity); */ 00389 /* fprintf(ioQQQ,"Ds %ld: %ld; %g %g %g %g\n",iteration,nzone,dense.gas_phase[ipHYDROGEN],er,pressure.PresTotlCurr,phycon.te); */ 00390 if(globalpresmode == ORIGINAL || lastzone != nzone || fabs(er-erp) < SMALLFLOAT) 00391 { 00392 /* First time in this zone, or when last step made no change, 00393 * take step hopefully in the right direction... 00394 * ...or at least somewhere */ 00395 if( zonepresmode == SUBSONIC ) 00396 { 00397 /* when gas pressure dominated need to increase density to increase pressure */ 00398 factor = pressure.PresTotlCorrect / pressure.PresTotlCurr; 00399 ipPRE = 0; 00400 } 00401 else 00402 { 00403 /* when ram pressure dominated need to decrease density to increase pressure */ 00404 factor = pressure.PresTotlCurr / pressure.PresTotlCorrect; 00405 ipPRE = 1; 00406 } 00407 if(fabs(factor-1.) > 0.01) 00408 { 00409 factor = 1.+sign(0.01,factor-1.); 00410 } 00411 erp = er; 00412 dp = dense.gas_phase[ipHYDROGEN]; 00413 erpp = -1.; 00414 dpp = -1.; 00415 } 00416 else 00417 { 00418 #if 0 00419 printf("Ds: %d; %g %g %g; %g %g %g tot %g\n",nzone,dense.gas_phase[ipHYDROGEN],dp,dpp,er,erp,erpp, 00420 pressure.PresTotlCorrect); 00421 #endif 00422 if(1 || dpp < 0. || fabs((dense.gas_phase[ipHYDROGEN]-dp)*(dp-dpp)*(dpp-dense.gas_phase[ipHYDROGEN])) < SMALLFLOAT) 00423 { 00424 /* second iteration on this zone, do linerar fit to previous Pres - rho curve 00425 * Linear approximation to guess root with two independent samples */ 00426 factor = (dense.gas_phase[ipHYDROGEN]*erp-dp*er)/(erp-er)/dense.gas_phase[ipHYDROGEN]; 00427 /* Limit step length to `reasonable' extrapolation */ 00428 width = fabs(1.-dp/dense.gas_phase[ipHYDROGEN]); 00429 if(width > 1e-2) 00430 width = 1e-2; 00431 00432 /* Subsonic case: pressure ^ with density ^ => increase density further */ 00433 /* Super " case: pressure ^ with density v => decrease density further */ 00434 00435 /* printf("Presmode %d flag %g factor %g\n",zonepresmode,(er-erp)*(dense.gas_phase[ipHYDROGEN]-dp),factor); */ 00436 /* WJH 21 May 2004: I think that this part is to force the solution 00437 * towards the required branch when it appears to have the "wrong" value 00438 * of dP/drho */ 00439 if(zonepresmode == SUBSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) < 0) 00440 { 00441 factor = 1+3*width; 00442 } 00443 else if(zonepresmode == SUPERSONIC && (er-erp)*(dense.gas_phase[ipHYDROGEN]-dp) > 0) 00444 { 00445 factor = 1-3*width; 00446 } 00447 00448 if(fabs(factor-1.) > 3*width) 00449 { 00450 factor = 1.+sign(3*width,factor-1.); 00451 } 00452 ipPRE = 2; 00453 if(fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT) 00454 { 00455 dpp = dp; 00456 erpp = erp; 00457 dp = dense.gas_phase[ipHYDROGEN]; 00458 erp = er; 00459 } 00460 } 00461 else 00462 { 00463 /* 3rd or more solutions in this zone so have extensive data 00464 * on pressure - density relation 00465 * Use quadratic fit to last three errors to estimate optimum */ 00466 double a, b, c, q, dmin, emin, dsol, f1 , f2; 00467 int iroot; 00468 a = er/(dense.gas_phase[ipHYDROGEN]-dp)/(dense.gas_phase[ipHYDROGEN]-dpp) + 00469 erp/(dp-dense.gas_phase[ipHYDROGEN])/(dp-dpp)+ 00470 erpp/(dpp-dense.gas_phase[ipHYDROGEN])/(dpp-dp); 00471 b = (erp-erpp)/(dp-dpp) - a * (dp+dpp); 00472 c = erp - dp*(a*dp+b); 00473 dmin = (-0.5*b/a); 00474 emin = (a*dmin+b)*dmin+c; 00475 if(a < 0) 00476 { 00477 a = -a; 00478 b = -b; 00479 c = -c; 00480 } 00481 #if 0 00482 fprintf(ioQQQ,"Check 1: %g %g\n",er,(a*dense.gas_phase[ipHYDROGEN]+b)*dense.gas_phase[ipHYDROGEN]+c); 00483 fprintf(ioQQQ,"Check 2: %g %g\n",erp,(a*dp+b)*dp+c); 00484 fprintf(ioQQQ,"Check 3: %g %g\n",erpp,(a*dpp+b)*dpp+c); 00485 #endif 00486 q = b*b-4*a*c; 00487 if(q < 0) 00488 { 00489 /* Imaginary root, search for local minimum */ 00490 /* printf("No root at %d (%g cf %g) => %g\n",nzone,q,b*b,dmin); */ 00491 factor = dmin/dense.gas_phase[ipHYDROGEN]; 00492 00496 if(globalpresmode == STRONGD && -q > 1e-3*b*b) 00497 { 00498 zonepresmode = SUBSONIC; 00499 } 00500 } 00501 else 00502 { 00503 /* Look for nearest root */ 00504 /* if(zonepresmode == SUPERSONIC || (zonepresmode != SUBSONIC && (dense.gas_phase[ipHYDROGEN]-dmin) < 0)) */ 00505 00506 if(zonepresmode == SUPERSONIC) 00507 iroot = 1; 00508 else 00509 iroot = 0; 00510 00511 if(emin > 0) 00512 iroot = ! iroot; 00513 00514 if(iroot) 00515 { 00516 /* Look for smaller root */ 00517 if(b > 0) 00518 { 00519 dsol = -(b+sqrt(q))/(2*a); 00520 } 00521 else 00522 { 00523 dsol = 2*c/(-b+sqrt(q)); 00524 } 00525 } 00526 else 00527 { 00528 if(b < 0) 00529 { 00530 dsol = (-b+sqrt(q))/(2*a); 00531 } 00532 else 00533 { 00534 dsol = -2*c/(b+sqrt(q)); 00535 } 00536 } 00537 factor = dsol/dense.gas_phase[ipHYDROGEN]; 00538 #if 0 00539 fprintf(ioQQQ,"Check 4: %g %g %d %g %g\n",dsol,(a*dsol+b)*dsol+c,iroot,dmin,emin); 00540 #endif 00541 } 00542 /* Limit step length */ 00543 f1 = fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]); 00544 f2 = fabs(1.- dp/dense.gas_phase[ipHYDROGEN]); 00545 /*width = MAX2(fabs(1.-dpp/dense.gas_phase[ipHYDROGEN]),fabs(1.-dp/dense.gas_phase[ipHYDROGEN]));*/ 00546 width = MAX2(f1,f2); 00547 /* width = MAX2(width,1e-2); */ 00548 if(fabs(factor-1.) > 3*width) 00549 { 00550 factor = 1.+sign(3*width,factor-1.); 00551 } 00552 ipPRE = 3; 00553 if(fabs(dp-dense.gas_phase[ipHYDROGEN]) > SMALLFLOAT) 00554 { 00555 dpp = dp; 00556 erpp = erp; 00557 dp = dense.gas_phase[ipHYDROGEN]; 00558 erp = er; 00559 } 00560 } 00561 } 00562 00563 #if 0 00564 printf("Out: %d; %g; %g %g; %g %g\n",nzone,factor*dense.gas_phase[ipHYDROGEN],dp,dpp,erp,erpp); 00565 #endif 00566 lastzone = nzone; 00567 00568 if( DIAG_PRINT ) 00569 fprintf(ioQQQ,"windv %li r %g v %g f %g\n", 00570 nzone,radius.depth,wind.windv,DynaFlux(radius.depth)); 00571 00572 /* convergence trace at this level */ 00573 if( trace.nTrConvg>=1 ) 00574 { 00575 fprintf( ioQQQ, 00576 " DynaPresChngFactor mode %s scale factor %.3f vel %.3e\n", 00577 dynamics.chPresMode , factor , wind.windv ); 00578 } 00579 00580 { 00581 /*@-redef@*/ 00582 enum{DEBUG_LOC=false}; 00583 /*@+redef@*/ 00584 if( DEBUG_LOC ) 00585 { 00586 char chPRE[][4] = {"gas" , "ram", "sec", "par" }; 00587 fprintf(ioQQQ, 00588 "pre %s\tfac\t%.5f\n", 00589 chPRE[ipPRE], 00590 factor 00591 ); 00592 } 00593 } 00594 00595 /* identically zero velocities cannot occur */ 00596 ASSERT( wind.windv != 0. || (wind.windv == 0. && dynamics.lgStatic) ); 00597 00598 return factor; 00599 } 00600 00601 /* ============================================================================== */ 00602 /* DynaIonize, called from ConvBase to evaluate advective terms for current conditions, 00603 * calculates terms to add to ionization balance equation */ 00604 void DynaIonize(void) 00605 { 00606 long int nelem, 00607 ion, 00608 mol , 00609 i; 00610 00611 DEBUG_ENTRY( "DynaIonize()" ); 00612 00613 /* the time (s) needed for gas to move Dyn_dr */ 00614 /* >>chng 02 dec 11 rjrw increase timestep when beyond end of previous zone, to allow -> eqm */ 00615 if( !dynamics.lgStatic ) 00616 { 00617 /* in time dependent model timestep only changes once at end of iteration 00618 * and is constant across a model */ 00619 /* usual case - full dynamics */ 00620 timestep = -Dyn_dr/wind.windv; 00621 } 00622 /* printf("%d %g %g %g %g\n",nzone,radius.depth,Dyn_dr,radius.depth-Old_depth[nOld_zone-1],timestep); */ 00623 00624 ASSERT(nzone<struc.nzlim ); 00625 if( nzone > 0 ) 00626 EnthalpyDensity[nzone-1] = (realnum)phycon.EnthalpyDensity; 00627 00628 /* do nothing on first iteration or when looking beyond previous iteration */ 00629 /* >>chng 05 jan 27, from hardwired "iteration < 2" to more general case, 00630 * this is set with SET DYNAMICS RELAX command with the default of 2 */ 00631 /*if( iteration < 2 || radius.depth > dynamics.oldFullDepth)*/ 00632 if( iteration < dynamics.n_initial_relax+1 ) 00633 { 00634 /* first iteration, return zero */ 00635 dynamics.Cool = 0.; 00636 dynamics.Heat = 0.; 00637 dynamics.dCooldT = 0.; 00638 dynamics.dHeatdT = 0.; 00639 00640 dynamics.Rate = 0.; 00641 00642 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 00643 { 00644 for( ion=0; ion<nelem+2; ++ion ) 00645 { 00646 dynamics.Source[nelem][ion] = 0.; 00647 } 00648 } 00649 00650 for(mol=0;mol<N_H_MOLEC;mol++) 00651 { 00652 dynamics.H2_molec[mol] = 0.; 00653 } 00654 for(mol=0;mol<mole.num_comole_calc;mol++) 00655 { 00656 dynamics.CO_molec[mol] = 0.; 00657 } 00658 return; 00659 } 00660 00661 if( MAINPRINT ) 00662 { 00663 fprintf(ioQQQ,"workwork\t%li\t%.3e\t%.3e\t%.3e\n", 00664 nzone, 00665 phycon.EnthalpyDensity, 00666 0.5*POW2(wind.windv)*dense.xMassDensity , 00667 5./2.*pressure.PresGasCurr 00668 ); 00669 } 00670 00671 /* net cooling due to advection */ 00672 /* >>chng 01 aug 01, removed hden from dilution, put back into here */ 00673 /* >>chng 01 sep 15, do not update this variable the very first time this routine 00674 * is called at the new zone. */ 00675 dynamics.Cool = phycon.EnthalpyDensity/timestep*dynamics.lgCoolHeat; 00676 dynamics.Heat = AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep*dynamics.lgCoolHeat; 00677 /* temp deriv of cooling minus heating */ 00678 dynamics.dCooldT = 5./2.*pressure.PresGasCurr/phycon.te/timestep*dynamics.lgCoolHeat; 00679 dynamics.dHeatdT = 0.*dynamics.lgCoolHeat; 00680 00681 #if 0 00682 /* printf("DynaCool %g %g %g\n", dynamics.Cool, phycon.EnthalpyDensity/timestep,AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep); */ 00683 if(dynamics.Cool > 0) { 00684 dynamics.Heat = 0.; 00685 /* temp deriv of cooling minus heating */ 00686 dynamics.dCooldT = 5./2.*pressure.PresGasCurr/phycon.te/timestep*dynamics.lgCoolHeat; 00687 dynamics.dHeatdT = 0.*dynamics.lgCoolHeat; 00688 } else { 00689 dynamics.Heat = -dynamics.Cool; 00690 dynamics.Cool = 0.; 00691 /* temp deriv of cooling minus heating */ 00692 dynamics.dCooldT = 0.*dynamics.lgCoolHeat; 00693 dynamics.dHeatdT = -5./2.*pressure.PresGasCurr/phycon.te/timestep*dynamics.lgCoolHeat; 00694 } 00695 #endif 00696 00697 # if 0 00698 if( MAINPRINT || nzone>17 && iteration == 10) 00699 { 00700 fprintf(ioQQQ, 00701 "dynamics cool-heat\t%li\t%.3e\t%.3e\t%.3e\t%.3e\t %.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 00702 nzone, 00703 phycon.te, 00704 dynamics.CoolHeat, 00705 thermal.htot, 00706 phycon.EnthalpyDensity/timestep, 00707 AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN]/timestep, 00708 phycon.EnthalpyDensity, 00709 AdvecSpecificEnthalpy*dense.gas_phase[ipHYDROGEN], 00710 dense.gas_phase[ipHYDROGEN], 00711 timestep); 00712 } 00713 # endif 00714 00715 /* second or greater iteration, have advective terms */ 00716 /* this will evaluate advective terms for current physical conditions */ 00717 00718 /* the rate (s^-1) atoms drift in from upstream, a source term for the ground */ 00719 00720 /* dynamics.Hatom/timestep is the source (cm^-3 s^-1) of neutrals, 00721 normalized to (s^-1) by the next higher ionization state as for all 00722 other recombination terms */ 00723 00724 /* dynamics.xIonDense[ipHYDROGEN][1]/timestep is the sink (cm^-3 s^-1) of 00725 ions, normalized to (s^-1) by the ionization state as for all other 00726 ionization terms */ 00727 00728 dynamics.Rate = 1./timestep; 00729 00730 for(mol=0;mol<N_H_MOLEC;mol++) 00731 { 00732 dynamics.H2_molec[mol] = Upstream_H2_molec[mol]* dense.gas_phase[ipHYDROGEN]*dynamics.Rate; 00733 } 00734 00735 for(mol=0;mol<mole.num_comole_calc;mol++) 00736 { 00737 dynamics.CO_molec[mol] = Upstream_CO_molec[mol]* dense.gas_phase[ipHYDROGEN]*dynamics.Rate; 00738 } 00739 00740 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 00741 { 00742 if( dense.lgElmtOn[nelem] ) 00743 { 00744 /* check that the total number of each element is conserved along the flow */ 00745 if(fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]>=1e-3) 00746 { 00747 /* sum of all H in standard H chem */ 00748 realnum sumh = 0.; 00749 for(i=0;i<N_H_MOLEC;i++) 00750 { 00751 sumh += hmi.Hmolec[i]*hmi.nProton[i]; 00752 } 00753 fprintf(ioQQQ, 00754 "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", 00755 nzone , 00756 nelem, 00757 UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN], 00758 dense.gas_phase[nelem] , 00759 (UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN]-dense.gas_phase[nelem]) / 00760 (UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN]) , 00761 (dense.gas_phase[ipHYDROGEN]-sumh) / dense.gas_phase[ipHYDROGEN] , 00762 dense.H_sum_in_CO / dense.gas_phase[ipHYDROGEN] ); 00763 } 00764 /* ASSERT( fabs(UpstreamElem[nelem]*dense.gas_phase[ipHYDROGEN] -dense.gas_phase[nelem])/dense.gas_phase[nelem]<1e-3); */ 00765 for( ion=0; ion<dense.IonLow[nelem]; ++ion ) 00766 { 00767 dynamics.Source[nelem][ion] = 0.; 00768 } 00769 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 00770 { 00771 /* Normalize to next higher state in current zone, except at first iteration 00772 where upstream version may be a better estimate (required for 00773 convergence in the small timestep limit) */ 00774 00775 dynamics.Source[nelem][ion] = 00776 /* UpstreamIon is ion number per unit hydrogen because dilution is 1/hden 00777 * this forms the ratio of upstream atom over current ion, per timestep, 00778 * so Recomb has units s-1 */ 00779 UpstreamIon[nelem][ion] * dense.gas_phase[ipHYDROGEN] / timestep; 00780 00781 } 00782 for( ion=dense.IonHigh[nelem]+1;ion<nelem+2; ++ion ) 00783 { 00784 dynamics.Source[nelem][ion] = 0.; 00785 } 00786 } 00787 } 00788 # if 0 00789 fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n", 00790 nzone, 00791 dynamics.Rate, 00792 dynamics.Source[ipHYDROGEN][0], 00793 dynamics.Rate, 00794 dynamics.Source[ipCARBON][3]); 00795 # endif 00796 #if 0 00797 nelem = ipCARBON; 00798 ion = 3; 00799 /*if( nzone > 160 && iteration > 1 )*/ 00800 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", 00801 nzone, 00802 ipUpstream, 00803 radius.depth , 00804 Old_depth[ipUpstream], 00805 dense.xIonDense[nelem][ion], 00806 UpstreamIon[nelem][ion]* dense.gas_phase[ipHYDROGEN], 00807 Old_xIonDense[ipUpstream][nelem][ion] , 00808 dense.xIonDense[nelem][ion+1], 00809 UpstreamIon[nelem][ion+1]* dense.gas_phase[ipHYDROGEN], 00810 Old_xIonDense[ipUpstream][nelem][ion+1] , 00811 timestep, 00812 dynamics.Source[nelem][ion] 00813 ); 00814 #endif 00815 if( MAINPRINT ) 00816 { 00817 fprintf(ioQQQ," DynaIonize, %4li photo=%.2e , H recom= %.2e \n", 00818 nzone,dynamics.Rate , dynamics.Source[0][0] ); 00819 } 00820 return; 00821 } 00822 00823 /* ============================================================================== */ 00824 /* DynaStartZone called at start of zone calculation when advection is turned on */ 00825 void DynaStartZone(void) 00826 { 00827 /* this routine is called at the start of a zone calculation, by ZoneStart: 00828 * 00829 * it sets deduced variables to zero if first iteration, 00830 * 00831 * if upstream depth is is outside the computed structure on previous iteration, 00832 * return value at shielded face 00833 * 00834 * Also calculates discretization_error, an estimate of the accuracy of the source terms. 00835 * 00836 * */ 00837 00838 /* this is index of upstream cell in structure stored from previous iteration */ 00839 double upstream, dilution, dilutionleft, dilutionright, frac_next; 00840 00841 /* Properties for cell half as far upstream, used to converge timestep */ 00842 double hupstream, hnextfrac=-BIGFLOAT, hion, hmol; 00843 00844 /* Properties for cell at present depth, used to converge timestep */ 00845 double ynextfrac=-BIGFLOAT, yion, ymol; 00846 00847 long int nelem , ion, mol; 00848 00849 DEBUG_ENTRY( "DynaStartZone()" ); 00850 00851 /* do nothing on first iteration */ 00852 if( iteration < 2 ) 00853 { 00854 Upstream_hden = 0.; 00855 AdvecSpecificEnthalpy = 0.; 00856 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 00857 { 00858 for( ion=0; ion<nelem+2; ++ion ) 00859 { 00860 UpstreamIon[nelem][ion] = 0.; 00861 } 00862 } 00863 /* hydrogen molecules */ 00864 for(mol=0; mol<N_H_MOLEC;mol++) 00865 { 00866 Upstream_H2_molec[mol] = 0; 00867 } 00868 for(mol=0; mol<mole.num_comole_calc;mol++) 00869 { 00870 Upstream_CO_molec[mol] = 0; 00871 } 00872 00873 ipUpstream = 0; 00874 iphUpstream = 0; 00875 ipyUpstream = 0; 00876 return; 00877 } 00878 00879 /* radius.depth is distance from illuminated face of cloud to outer edge of 00880 * current zone, which has thickness radius.drad */ 00881 00882 /* find where the down stream point is, in previous iteration, 00883 * NB neg sign since velocity is negative, we are looking for point 00884 * where current gas cell was, in previous iteration */ 00885 00886 /* true, both advection and wind solution */ 00887 /* don't interpolate to the illuminated side of the first cell */ 00888 upstream = MAX2(Old_depth[0] , radius.depth + Dyn_dr); 00889 hupstream = 0.5*(radius.depth + upstream); 00890 00891 /* now locate upstream point in previous stored structure, 00892 * will be at the same point or higher than we found previously */ 00893 while( (Old_depth[ipUpstream+1] < upstream ) && 00894 ( ipUpstream < nOld_zone-1 ) ) 00895 { 00896 ++ipUpstream; 00897 } 00898 ASSERT( ipUpstream <= nOld_zone-1 ); 00899 00900 /* above loop will give ipUpstream == nOld_zone-1 if computed structure has been overrun */ 00901 00902 if( (ipUpstream != nOld_zone-1)&& ((Old_depth[ipUpstream+1] - Old_depth[ipUpstream])> SMALLFLOAT) ) 00903 { 00904 /* we have not overrun radius scale of previous iteration */ 00905 frac_next = ( upstream - Old_depth[ipUpstream])/ 00906 (Old_depth[ipUpstream+1] - Old_depth[ipUpstream]); 00907 Upstream_hden = Old_hden[ipUpstream] + 00908 (Old_hden[ipUpstream+1] - Old_hden[ipUpstream])* 00909 frac_next; 00910 dilutionleft = 1./Old_hden[ipUpstream]; 00911 dilutionright = 1./Old_hden[ipUpstream+1]; 00912 00913 /* fractional changes in density from passive advection */ 00914 /* >>chng 01 May 02 rjrw: use hden for dilution */ 00915 /* >>chng 01 aug 01, remove hden here, put back into vars when used in DynaIonize */ 00916 dilution = 1./Upstream_hden; 00917 00918 /* the advected enthalphy */ 00919 AdvecSpecificEnthalpy = (Old_EnthalpyDensity[ipUpstream]*dilutionleft + 00920 (Old_EnthalpyDensity[ipUpstream+1]*dilutionright - Old_EnthalpyDensity[ipUpstream]*dilutionleft)* 00921 frac_next); 00922 00923 ASSERT( Old_depth[ipUpstream] <= upstream && upstream <= Old_depth[ipUpstream+1] ); 00924 if( (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) * 00925 (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) < -SMALLFLOAT ) 00926 { 00927 fprintf(ioQQQ,"PROBLEM advected enthalpy density is not conserved %e\t%e\t%e\t%e\n", 00928 (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) * 00929 (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) , 00930 Old_EnthalpyDensity[ipUpstream]*dilutionleft, 00931 AdvecSpecificEnthalpy, 00932 Old_EnthalpyDensity[ipUpstream+1]*dilutionright); 00933 } 00934 00935 ASSERT( (Old_EnthalpyDensity[ipUpstream]*dilutionleft - AdvecSpecificEnthalpy) * 00936 (AdvecSpecificEnthalpy - Old_EnthalpyDensity[ipUpstream+1]*dilutionright) > -SMALLFLOAT ); 00937 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 00938 { 00939 UpstreamElem[nelem] = 0.; 00940 for( ion=0; ion<nelem+2; ++ion ) 00941 { 00942 /* Robin - I made several changes like the following - it seems easier to 00943 * bring out the division by the old hydrogen density rather than putting in 00944 * dilution and then looking for how diluation is defined. I think the code 00945 * is equivalent */ 00946 /* UpstreamIon is ion number per unit hydrogen, both at the upstream position */ 00947 UpstreamIon[nelem][ion] = 00948 Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream] + 00949 (Old_xIonDense[ipUpstream+1][nelem][ion]/Old_hden[ipUpstream+1] - 00950 Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream])* 00951 frac_next; 00952 00953 UpstreamElem[nelem] += UpstreamIon[nelem][ion]; 00954 } 00955 } 00956 00957 for(mol=0;mol<N_H_MOLEC;mol++) 00958 { 00959 Upstream_H2_molec[mol] = Old_H2_molec[ipUpstream][mol]/Old_hden[ipUpstream] + 00960 (Old_H2_molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] - 00961 Old_H2_molec[ipUpstream][mol]/Old_hden[ipUpstream])* 00962 frac_next; 00963 /* test on mol > 1, first two elements are H0 and H+, which are alread 00964 * counted in upstreamElem */ 00965 if(mol > 1) 00966 UpstreamElem[ipHYDROGEN] += Upstream_H2_molec[mol]*hmi.nProton[mol]; 00967 } 00968 /* loop is only up to mole.num_heavy_molec, not mole.num_comole_calc, since do not want 00969 * to consider atoms and ions that have also been done */ 00970 for(mol=0;mol<mole.num_comole_calc;mol++) 00971 { 00972 Upstream_CO_molec[mol] = Old_CO_molec[ipUpstream][mol]/Old_hden[ipUpstream] + 00973 (Old_CO_molec[ipUpstream+1][mol]/Old_hden[ipUpstream+1] - 00974 Old_CO_molec[ipUpstream][mol]/Old_hden[ipUpstream])* 00975 frac_next; 00976 00977 if(COmole[mol]->n_nuclei > 1) 00978 { 00979 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00980 { 00981 if( mole.lgElem_in_chemistry[nelem] ) 00982 { 00983 UpstreamElem[nelem] += 00984 COmole[mol]->nElem[nelem]*Upstream_CO_molec[mol]; 00985 } 00986 } 00987 } 00988 } 00989 } 00990 else 00991 { 00992 /* SPECIAL CASE - we have overrun the previous iteration's radius */ 00993 Upstream_hden = Old_hden[ipUpstream]; 00994 /* fractional changes in density from passive advection */ 00995 dilution = 1./Upstream_hden; 00996 /* AdvecSpecificEnthalpy enters as heat term */ 00997 AdvecSpecificEnthalpy = Old_EnthalpyDensity[ipUpstream]/Old_hden[ipUpstream]; 00998 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 00999 { 01000 UpstreamElem[nelem] = 0.; 01001 for( ion=0; ion<nelem+2; ++ion ) 01002 { 01003 /* UpstreamIon is ion number per unit hydrogen */ 01004 UpstreamIon[nelem][ion] = 01005 Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream]; 01006 UpstreamElem[nelem] += UpstreamIon[nelem][ion]; 01007 } 01008 } 01009 01010 for(mol=0;mol<N_H_MOLEC;mol++) 01011 { 01012 Upstream_H2_molec[mol] = Old_H2_molec[ipUpstream][mol]/Old_hden[ipUpstream]; 01013 if(mol > 1) 01014 UpstreamElem[ipHYDROGEN] += Upstream_H2_molec[mol]*hmi.nProton[mol]; 01015 } 01016 for(mol=0;mol<mole.num_comole_calc;mol++) 01017 { 01018 Upstream_CO_molec[mol] = Old_CO_molec[ipUpstream][mol]/Old_hden[ipUpstream]; 01019 if(COmole[mol]->n_nuclei > 1) 01020 { 01021 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 01022 { 01023 if( mole.lgElem_in_chemistry[nelem] ) 01024 { 01025 UpstreamElem[nelem] += 01026 Upstream_CO_molec[mol]*COmole[mol]->nElem[nelem]; 01027 } 01028 } 01029 } 01030 } 01031 } 01032 01033 /* Repeat enough of the above for half-step and no-step to judge convergence: 01034 the result of this code is the increment of discretization_error */ 01035 01036 while( (Old_depth[iphUpstream+1] < hupstream ) && 01037 ( iphUpstream < nOld_zone-1 ) ) 01038 { 01039 ++iphUpstream; 01040 } 01041 ASSERT( iphUpstream <= nOld_zone-1 ); 01042 01043 while( (Old_depth[ipyUpstream+1] < radius.depth ) && 01044 ( ipyUpstream < nOld_zone-1 ) ) 01045 { 01046 ++ipyUpstream; 01047 } 01048 ASSERT( ipyUpstream <= nOld_zone-1 ); 01049 01050 dynamics.dRad = BIGFLOAT; 01051 01052 if(iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT)) 01053 { 01054 hnextfrac = ( hupstream - Old_depth[iphUpstream])/ 01055 (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]); 01056 /* >>chng 02 nov 11, hhden never appears on rhs of = */ 01057 /*hhden = Old_hden[iphUpstream] + 01058 (Old_hden[iphUpstream+1] - Old_hden[iphUpstream])* 01059 hnextfrac;*/ 01060 } 01061 else 01062 { 01063 /* >>chng 06 mar 17, set this to zero */ 01064 hnextfrac = 0.; 01065 /* >>chng 02 nov 11, hhden never appears on rhs of = */ 01066 /*hhden = Old_hden[iphUpstream];*/ 01067 } 01068 if(ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT)) 01069 { 01070 ynextfrac = ( radius.depth - Old_depth[ipyUpstream])/ 01071 (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]); 01072 /* >>chng 02 nov 11, yhden never appears on rhsof = */ 01073 /*yhden = Old_hden[ipyUpstream] + 01074 (Old_hden[ipyUpstream+1] - Old_hden[ipyUpstream])* 01075 ynextfrac;*/ 01076 } 01077 else 01078 { 01079 /* >>chng 06 mar 17, set this to zero */ 01080 ynextfrac = 0.; 01081 /* >>chng 02 nov 11, yhden never appears on rhsof = */ 01082 /*yhden = Old_hden[ipyUpstream];*/ 01083 } 01084 01085 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 01086 { 01087 for( ion=0; ion<nelem+2; ++ion ) 01088 { 01089 double f1; 01090 if(ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT)) 01091 { 01092 yion = 01093 Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream] + 01094 (Old_xIonDense[ipyUpstream+1][nelem][ion]/Old_hden[ipyUpstream+1] - 01095 Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream])* 01096 ynextfrac; 01097 } 01098 else 01099 { 01100 yion = Old_xIonDense[ipyUpstream][nelem][ion]/Old_hden[ipyUpstream]; 01101 } 01102 if(iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT)) 01103 { 01104 hion = 01105 Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream] + 01106 (Old_xIonDense[iphUpstream+1][nelem][ion]/Old_hden[iphUpstream+1] - 01107 Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream])* 01108 hnextfrac; 01109 } 01110 else 01111 { 01112 hion = Old_xIonDense[iphUpstream][nelem][ion]/Old_hden[iphUpstream]; 01113 } 01114 01115 /* the proposed thickness of the next zone, there will be a scale factor, something 01116 * like 1/500, that will be applied when this is used in nextdr */ 01117 f1 = fabs(yion - UpstreamIon[nelem][ion] ); 01118 f1 = SDIV( f1 ); 01119 dynamics.dRad = MIN2(dynamics.dRad,Dyn_dr * 01120 /* don't pay attention to species with abundance relative to H less tghan 1e-6 */ 01121 MAX2(yion + UpstreamIon[nelem][ion],1e-6 ) / f1); 01122 /*MAX2(fabs(yion - UpstreamIon[nelem][ion] ),SMALLFLOAT));*/ 01123 /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */ 01124 01125 /* Must be consistent with convergence_error below */ 01126 /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */ 01127 /* >>chngf 02 aug 01, multiply by cell width */ 01128 /* this error is error due to the advection length not being zero - a finite 01129 * advection length. no need to bring convergence error to below 01130 * discretization error. when convergece error is lower than a fraction of this 01131 * errror we reduce the advection length. */ 01132 dynamics.discretization_error += POW2(yion-2.*hion+UpstreamIon[nelem][ion]); 01133 dynamics.error_scale2 += POW2(UpstreamIon[nelem][ion]-yion); 01134 } 01135 } 01136 01137 for(mol=0; mol < N_H_MOLEC; mol++) 01138 { 01139 double f1; 01140 if(ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT)) 01141 { 01142 ymol = 01143 Old_H2_molec[ipyUpstream][mol]/Old_hden[ipyUpstream] + 01144 (Old_H2_molec[ipyUpstream+1][mol]/Old_hden[ipyUpstream+1] - 01145 Old_H2_molec[ipyUpstream][mol]/Old_hden[ipyUpstream])* 01146 ynextfrac; 01147 } 01148 else 01149 { 01150 ymol = Old_H2_molec[ipyUpstream][mol]/Old_hden[ipyUpstream]; 01151 } 01152 if(iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT)) 01153 { 01154 hmol = 01155 Old_H2_molec[iphUpstream][mol]/Old_hden[iphUpstream] + 01156 (Old_H2_molec[iphUpstream+1][mol]/Old_hden[iphUpstream+1] - 01157 Old_H2_molec[iphUpstream][mol]/Old_hden[iphUpstream])* 01158 hnextfrac; 01159 } 01160 else 01161 { 01162 hmol = Old_H2_molec[iphUpstream][mol]/Old_hden[iphUpstream]; 01163 } 01164 01165 /* the proposed thickness of the next zone, there will be a scale factor, something 01166 * like 1/500, that will be applied when this is used in nextdr */ 01167 f1 = fabs(ymol - Upstream_H2_molec[mol] ); 01168 f1 = SDIV( f1 ); 01169 dynamics.dRad = MIN2(dynamics.dRad,Dyn_dr * 01170 /* don't pay attention to species with abundance relative to H less tghan 1e-6 */ 01171 MAX2(ymol + Upstream_H2_molec[mol],1e-6 ) / f1 ); 01172 /*MAX2(fabs(ymol - Upstream_H2_molec[mol] ),SMALLFLOAT));*/ 01173 /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */ 01174 01175 /* Must be consistent with convergence_error below */ 01176 /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */ 01177 /* >>chngf 02 aug 01, multiply by cell width */ 01178 dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_H2_molec[mol]); 01179 dynamics.error_scale2 += POW2(Upstream_H2_molec[mol]-ymol); 01180 } 01181 01182 for(mol=0; mol < mole.num_comole_calc; mol++) 01183 { 01184 double f1; 01185 if(ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT)) 01186 { 01187 ymol = 01188 Old_CO_molec[ipyUpstream][mol]/Old_hden[ipyUpstream] + 01189 (Old_CO_molec[ipyUpstream+1][mol]/Old_hden[ipyUpstream+1] - 01190 Old_CO_molec[ipyUpstream][mol]/Old_hden[ipyUpstream])* 01191 ynextfrac; 01192 } 01193 else 01194 { 01195 ymol = Old_CO_molec[ipyUpstream][mol]/Old_hden[ipyUpstream]; 01196 } 01197 if(iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT)) 01198 { 01199 hmol = 01200 Old_CO_molec[iphUpstream][mol]/Old_hden[iphUpstream] + 01201 (Old_CO_molec[iphUpstream+1][mol]/Old_hden[iphUpstream+1] - 01202 Old_CO_molec[iphUpstream][mol]/Old_hden[iphUpstream])* 01203 hnextfrac; 01204 } 01205 else 01206 { 01207 hmol = Old_CO_molec[iphUpstream][mol]/Old_hden[iphUpstream]; 01208 } 01209 01210 /* the proposed thickness of the next zone, there will be a scale factor, something 01211 * like 1/500, that will be applied when this is used in nextdr */ 01212 f1 = fabs(ymol - Upstream_CO_molec[mol] ); 01213 f1 = SDIV( f1 ); 01214 dynamics.dRad = MIN2(dynamics.dRad,Dyn_dr * 01215 /* don't pay attention to species with abundance relative to H less tghan 1e-6 */ 01216 MAX2(ymol + Upstream_CO_molec[mol],1e-6 ) / f1 ); 01217 /*MAX2(fabs(ymol - Upstream_H2_molec[mol] ),SMALLFLOAT));*/ 01218 /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */ 01219 01220 /* Must be consistent with convergence_error below */ 01221 /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */ 01222 /* >>chngf 02 aug 01, multiply by cell width */ 01223 dynamics.discretization_error += POW2(ymol-2.*hmol+Upstream_CO_molec[mol]); 01224 dynamics.error_scale2 += POW2(Upstream_CO_molec[mol]-ymol); 01225 } 01226 01227 if( MAINPRINT ) 01228 { 01229 fprintf(ioQQQ," DynaStartZone, %4li photo=%.2e , H recom= %.2e dil %.2e \n", 01230 nzone,dynamics.Rate , dynamics.Source[0][0] , dilution*dense.gas_phase[ipHYDROGEN] ); 01231 } 01232 return; 01233 } 01234 01235 /* DynaEndZone called at end of zone calculation when advection is turned on */ 01236 void DynaEndZone(void) 01237 { 01238 DEBUG_ENTRY( "DynaEndZone()" ); 01239 01240 /* this routine is called at the end of a zone calculation, by ZoneEnd */ 01241 01242 DivergePresInteg += wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad)); 01243 if(DIAG_PRINT) 01244 fprintf(ioQQQ,"Check dp: %g %g mom %g %g mas %g\n", 01245 wind.windv*(DynaFlux(radius.depth)-DynaFlux(radius.depth-radius.drad)), 01246 2*wind.windv*DynaFlux(radius.depth)*radius.drad/(1e16-radius.depth), 01247 wind.windv*DynaFlux(radius.depth), 01248 wind.windv*DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth), 01249 DynaFlux(radius.depth)*(1e16-radius.depth)*(1e16-radius.depth)); 01250 return; 01251 } 01252 01253 01254 /* ============================================================================== */ 01255 /* DynaEndIter called at end of iteration when advection is turned on */ 01256 void DynaEndIter(void) 01257 { 01258 /* this routine is called by IterRestart at the end of an iteration 01259 * when advection is turned on. currently it only derives a 01260 * timestep by looking at the spatial derivative of 01261 * some stored quantities */ 01262 long int i; 01263 static long int nTime_dt_array_element = 0; 01264 01265 DEBUG_ENTRY( "DynaEndIter()" ); 01266 01267 /* set stopping outer radius to current outer radius once we have let 01268 * solution relax dynamics.n_initial_relax times 01269 * note the off by one confusion - relax is 2 by default, 01270 * want to do two static iterations then start dynamics 01271 * iteration was incremented before this call so iteration == 2 at 01272 * end of first iteration */ 01273 if( iteration == dynamics.n_initial_relax+1) 01274 { 01275 fprintf(ioQQQ,"DYNAMICS DynaEndIter sets stop radius to %.2e after" 01276 "dynamics.n_initial_relax=%li iterations.\n", 01277 radius.depth, 01278 dynamics.n_initial_relax); 01279 for( i=iteration-1; i<iterations.iter_malloc; ++i ) 01280 /* set stopping radius to current radius, this stops 01281 * dynamical solutions from overrunning the upstream scale 01282 * and extending to very large radius due to unphysical heat 01283 * appearing to advect into region */ 01284 radius.router[i] = radius.depth; 01285 } 01286 01287 DivergePresInteg = 0.; 01288 01289 /* This routine is only called if advection is turned on at the end of 01290 * each iteration. The test 01291 * is to check whether wind velocity is also set by dynamics code */ 01292 01293 /* !dynamics.lgStatic true - a true dynamical model */ 01294 if( !dynamics.lgStatic ) 01295 { 01296 if(iteration == dynamics.n_initial_relax+1 ) 01297 { 01298 /* let model settle down for n_initial_relax iterations before we begin 01299 * dynamics.n_initial_relax set with "set dynamics relax" 01300 * command - it gives the first iteration where we do dynamics 01301 * note the off by one confusion - relax is 2 by default, 01302 * want to do two static iterations then start dynamics 01303 * iteration was incremented before this call so iteration == 2 01304 * at end of first iteration */ 01305 if( dynamics.AdvecLengthInit> 0. ) 01306 { 01307 Dyn_dr = dynamics.AdvecLengthInit; 01308 } 01309 else 01310 { 01311 /* -ve dynamics.adveclengthlimit sets length as fraction of first iter total depth */ 01312 Dyn_dr = -dynamics.AdvecLengthInit*radius.depth; 01313 } 01314 01315 if( MAINPRINT ) 01316 { 01317 fprintf(ioQQQ," DynaEndIter, dr=%.2e \n", 01318 Dyn_dr ); 01319 } 01320 } 01321 else if(iteration > dynamics.n_initial_relax+1 ) 01322 { 01323 /* evaluate errors and update Dyn_dr */ 01324 DynaNewStep(); 01325 } 01326 } 01327 else 01328 { 01329 /* this is time-dependent static model */ 01330 static double HeatInitial=-1. , HeatRadiated=-1. , 01331 DensityInitial=-1; 01332 /* n_initial_relax is number of time-steady models before we start 01333 * to evolve, set with "set dynamics relax" command */ 01334 Dyn_dr = 0.; 01335 fprintf(ioQQQ, 01336 "DEBUG times enter timestep %.2e elapsed_time %.2e iteration %li relax %li \n", 01337 timestep , 01338 dynamics.time_elapsed, 01339 iteration , dynamics.n_initial_relax); 01340 if( iteration > dynamics.n_initial_relax ) 01341 { 01342 /* this is into the changing part of the simulation */ 01343 long int nTimeVary; 01344 01345 /* evaluate errors */ 01346 DynaNewStep(); 01347 01348 /* this is set true on CORONAL XXX TIME INIT command, says to use constant 01349 * temperature for first n_initial_relax iterations, then let run free */ 01350 if( dynamics.lg_coronal_time_init ) 01351 { 01352 thermal.lgTemperatureConstant = false; 01353 thermal.ConstTemp = 0.; 01354 } 01355 01356 /* time variable branch, now without dynamics */ 01357 /* elapsed time - don't increase dynamics.time_elapsed during first two 01358 * two iterations since this sets static model */ 01359 dynamics.time_elapsed += timestep; 01360 /* timestep_stop is -1 if not set with explicit stop time */ 01361 if( timestep_stop > 0 && dynamics.time_elapsed > timestep_stop ) 01362 { 01363 dynamics.lgStatic_completed = true; 01364 } 01365 01366 if( dynamics.time_elapsed < time_elapsed_time[0] ) 01367 { 01368 /* if very early times not specified assume no flux variation yet */ 01369 rfield.time_continuum_scale = 1.; 01370 } 01371 else if( dynamics.time_elapsed > time_elapsed_time[nTime_flux-1] ) 01372 { 01373 fprintf( ioQQQ, 01374 " PROBLEM - DynaEnditer - I need the continuum at time %.2e but the table ends at %.2e.\n" , 01375 dynamics.time_elapsed , 01376 time_elapsed_time[nTime_flux-1]); 01377 cdEXIT(EXIT_FAILURE); 01378 } 01379 else 01380 { 01381 rfield.time_continuum_scale = (realnum)linint( 01382 /* the times in seconds */ 01383 time_elapsed_time, 01384 /* the rfield.time_continuum_scale factors */ 01385 time_flux_ratio, 01386 /* the number of rfield.time_continuum_scale factors */ 01387 nTime_flux, 01388 /* the desired time */ 01389 dynamics.time_elapsed); 01390 } 01391 01392 fprintf(ioQQQ,"DEBUG time dep reset continuum zero timestep %.2e elapsed time %.2e scale %.2e", 01393 timestep , 01394 dynamics.time_elapsed, 01395 rfield.time_continuum_scale); 01396 if( dynamics.lgRecom ) 01397 { 01398 fprintf(ioQQQ," recom"); 01399 } 01400 fprintf(ioQQQ,"\n"); 01401 01402 /* make sure that at least one continuum source is variable */ 01403 nTimeVary = 0; 01404 for( i=0; i < rfield.nspec; i++ ) 01405 { 01406 /* this is set true if particular continuum source can vary with time 01407 * set true if TIME appears on intensity / luminosity command line */ 01408 if( rfield.lgTimeVary[i] ) 01409 ++nTimeVary; 01410 } 01411 01412 if( hextra.lgTurbHeatVaryTime ) 01413 { 01414 /* vary extra heating */ 01415 hextra.TurbHeat = hextra.TurbHeatSave * rfield.time_continuum_scale; 01416 fprintf(ioQQQ,"DEBUG TurbHeat vary new heat %.2e\n", 01417 hextra.TurbHeat); 01418 } 01419 # if 0 01420 /* >>chng 07 may 23, rm this time dependent logic from here 01421 * time dependent continua implemented when continuum is reset 01422 * in inter_startend */ 01423 else if( nTimeVary ) 01424 { 01425 for( i=0; i<rfield.nupper; ++i ) 01426 { 01427 /* >>chng 06 mar 22, break into constant and time dependent parts */ 01428 rfield. flux[i] = rfield.flux_beam_const[i] + 01429 rfield.flux_beam_time[i]*rfield.time_continuum_scale + 01430 rfield.flux_isotropic[i]; 01431 rfield.flux_total_incident[i] = rfield.flux[i]; 01432 } 01433 } 01434 # endif 01435 else if( !nTimeVary ) 01436 { 01437 fprintf(ioQQQ," DISASTER - there were no variable continua " 01438 "or heat sources - put TIME option on at least one " 01439 "luminosity or hextra command.\n"); 01440 cdEXIT(EXIT_FAILURE); 01441 } 01442 /* this is heat radiated, after correction for change of H density in constant 01443 * pressure cloud */ 01444 HeatRadiated += (thermal.ctot-dynamics.Cool) * timestep * 01445 (DensityInitial / dense.gas_phase[ipHYDROGEN]); 01446 } 01447 else 01448 { 01449 /* this branch, during initial relaxation of solution */ 01450 HeatInitial = 1.5 * pressure.PresGasCurr; 01451 HeatRadiated = 0.; 01452 DensityInitial = dense.gas_phase[ipHYDROGEN]; 01453 fprintf(ioQQQ,"DEBUG relaxing times requested %li this is step %li\n", 01454 dynamics.n_initial_relax , iteration); 01455 } 01456 fprintf(ioQQQ,"DEBUG heat conser HeatInitial=%.2e HeatRadiated=%.2e\n", 01457 HeatInitial , HeatRadiated ); 01458 01459 /* at this point dynamics.time_elapsed is the time at the end of the 01460 * previous iteration. We need dt for the next iteration */ 01461 if( dynamics.time_elapsed > time_elapsed_time[nTime_dt_array_element] ) 01462 { 01463 /* time has advanced to next table point, 01464 * set timestep to specified value */ 01465 ++nTime_dt_array_element; 01466 /* this is an assert since we already qualified the array 01467 * element above */ 01468 ASSERT( nTime_dt_array_element < nTime_flux ); 01469 01470 /* option to set flag for recombination logic */ 01471 if( lgtime_Recom[nTime_dt_array_element] ) 01472 { 01473 fprintf(ioQQQ,"DEBUG dynamics turn on recombination logic\n"); 01474 dynamics.lgRecom = true; 01475 /* set largest possible zone thickness to value on previous 01476 * iteration when light was on - during recombination conditions 01477 * become much more homogeneous and dr can get too large, 01478 * crashing into H i-front */ 01479 radius.sdrmax = radius.dr_max_last_iter/3.; 01480 } 01481 01482 if( lgtime_dt_specified ) 01483 { 01484 /* this is the new timestep */ 01485 fprintf(ioQQQ,"DEBUG lgtimes increment Time to %li %.2e\n" ,nTime_dt_array_element, 01486 timestep); 01487 timestep = time_dt[nTime_dt_array_element]; 01488 /* option to change time step factor - default is 1.2 set in DynaZero */ 01489 if( time_dt_scale_factor[nTime_dt_array_element] > 0. ) 01490 timestep_factor = time_dt_scale_factor[nTime_dt_array_element]; 01491 } 01492 } 01493 else if( lgtime_dt_specified ) 01494 { 01495 /* we are between two points in the table, increase timestep */ 01496 timestep *= timestep_factor; 01497 fprintf(ioQQQ,"DEBUG lgtimes increment Timeby timestep_factor to %li %.2e\n" , 01498 nTime_dt_array_element, 01499 timestep ); 01500 if(dynamics.time_elapsed+timestep > time_elapsed_time[nTime_dt_array_element] ) 01501 { 01502 fprintf(ioQQQ,"DEBUG lgtimes but reset to %.2e\n" ,timestep ); 01503 timestep = 1.0001*(time_elapsed_time[nTime_dt_array_element]-dynamics.time_elapsed); 01504 } 01505 } 01506 else 01507 { 01508 /* time not specified in third column, so use initial */ 01509 timestep = timestep_next(); 01510 } 01511 fprintf(ioQQQ,"DEBUG times exit timestep %.2e elapsed_time %.2e scale %.2e \n", 01512 timestep , 01513 dynamics.time_elapsed, 01514 rfield.time_continuum_scale); 01515 } 01516 01517 /* Dyn_dr == 0 is for static time dependent cloud */ 01518 ASSERT( (iteration < dynamics.n_initial_relax+1) || 01519 Dyn_dr > 0. || 01520 (Dyn_dr == 0. && wind.windv==0.) ); 01521 01522 /* reset the upstream counters */ 01523 ipUpstream = iphUpstream = ipyUpstream = 0; 01524 dynamics.discretization_error = 0.; 01525 dynamics.error_scale2 = 0.; 01526 01527 /* save results from previous iteration */ 01528 DynaSaveLast(); 01529 return; 01530 } 01531 01532 /*DynaNewStep work out convergence errors */ 01533 STATIC void DynaNewStep(void) 01534 { 01535 long int ilast = 0, 01536 i, 01537 nelem, 01538 ion, 01539 mol; 01540 01541 double frac_next=-BIGFLOAT, 01542 Oldi_hden, 01543 Oldi_ion, 01544 Oldi_mol; 01545 01546 DEBUG_ENTRY( "DynaNewStep()" ); 01547 01548 /*n = MIN2(nzone, NZLIM-1);*/ 01549 dynamics.convergence_error = 0; 01550 dynamics.error_scale1 = 0.; 01551 01552 ASSERT( nzone < struc.nzlim); 01553 for(i=0;i<nzone;++i) 01554 { 01555 /* Interpolate for present position in previous solution */ 01556 while( (Old_depth[ilast] < struc.depth[i] ) && 01557 ( ilast < nOld_zone-1 ) ) 01558 { 01559 ++ilast; 01560 } 01561 ASSERT( ilast <= nOld_zone-1 ); 01562 01563 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) ) 01564 { 01565 frac_next = ( struc.depth[i] - Old_depth[ilast])/ 01566 (Old_depth[ilast+1] - Old_depth[ilast]); 01567 Oldi_hden = Old_hden[ilast] + 01568 (Old_hden[ilast+1] - Old_hden[ilast])* 01569 frac_next; 01570 } 01571 else 01572 { 01573 Oldi_hden = Old_hden[ilast]; 01574 } 01575 /* Must be consistent with discretization_error above */ 01576 /* >>chngf 02 aug 01, multiply by cell width */ 01577 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 01578 { 01579 for( ion=0; ion<nelem+2; ++ion ) 01580 { 01581 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) ) 01582 { 01583 Oldi_ion = (Old_xIonDense[ilast][nelem][ion] + 01584 (Old_xIonDense[ilast+1][nelem][ion]-Old_xIonDense[ilast][nelem][ion])* 01585 frac_next); 01586 } 01587 else 01588 { 01589 Oldi_ion = Old_xIonDense[ilast][nelem][ion]; 01590 } 01591 dynamics.convergence_error += POW2(Oldi_ion/Oldi_hden-struc.xIonDense[nelem][ion][i]/struc.hden[i]) /* *struc.dr[i] */; 01592 01593 /* >>chng 02 nov 11, add first error scale estimate from Robin */ 01594 dynamics.error_scale1 += POW2(struc.xIonDense[nelem][ion][i]/struc.hden[i]); 01595 } 01596 } 01597 for( mol=0; mol < N_H_MOLEC; mol++) 01598 { 01599 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) ) 01600 { 01601 Oldi_mol = (Old_H2_molec[ilast][mol] + 01602 (Old_H2_molec[ilast+1][mol]-Old_H2_molec[ilast][mol])* 01603 frac_next); 01604 } 01605 else 01606 { 01607 Oldi_mol = Old_H2_molec[ilast][mol]; 01608 } 01609 dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.H2_molec[mol][i]/struc.hden[i]) /* *struc.dr[i] */; 01610 01611 /* >>chng 02 nov 11, add first error scale estimate from Robin 01612 * used to normalize the above convergence_error */ 01613 dynamics.error_scale1 += POW2(struc.H2_molec[mol][i]/struc.hden[i]); 01614 } 01615 for( mol=0; mol < mole.num_comole_calc; mol++) 01616 { 01617 if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) ) 01618 { 01619 Oldi_mol = (Old_CO_molec[ilast][mol] + 01620 (Old_CO_molec[ilast+1][mol]-Old_CO_molec[ilast][mol])* 01621 frac_next); 01622 } 01623 else 01624 { 01625 Oldi_mol = Old_CO_molec[ilast][mol]; 01626 } 01627 dynamics.convergence_error += POW2(Oldi_mol/Oldi_hden-struc.CO_molec[mol][i]/struc.hden[i]); 01628 01629 /* >>chng 02 nov 11, add first error scale estimate from Robin 01630 * used to normalize the above convergence_error */ 01631 dynamics.error_scale1 += POW2(struc.CO_molec[mol][i]/struc.hden[i]); 01632 } 01633 } 01634 01635 /* convergence_error is an estimate of the convergence of the solution from its change during the last iteration, 01636 discretization_error is an estimate of the accuracy of the advective terms, calculated in DynaStartZone above: 01637 if dominant error is from the advective terms, need to make them more accurate. 01638 */ 01639 01640 /* report properties of previous iteration */ 01641 fprintf(ioQQQ,"DYNAMICS DynaNewStep: Dyn_dr %.2e convergence_error %.2e discretization_error %.2e error_scale1 %.2e error_scale2 %.2e\n", 01642 Dyn_dr, dynamics.convergence_error , dynamics.discretization_error , 01643 dynamics.error_scale1 , dynamics.error_scale2 01644 ); 01645 01646 /* >>chng 02 nov 29, dynamics.convergence_tolerance is now set to 0.1 in init routine */ 01647 if( dynamics.convergence_error < dynamics.convergence_tolerance*dynamics.discretization_error ) 01648 Dyn_dr /= 1.5; 01649 return; 01650 } 01651 01652 /*DynaSaveLast save results from previous iteration */ 01653 STATIC void DynaSaveLast(void) 01654 { 01655 long int i, 01656 ion, 01657 nelem, 01658 mol; 01659 01660 DEBUG_ENTRY( "DynaSaveLast()" ); 01661 01662 /* Save results from previous iteration */ 01663 nOld_zone = nzone; 01664 dynamics.oldFullDepth = struc.depth[nzone-1]; 01665 ASSERT( nzone < struc.nzlim ); 01666 for( i=0; i<nzone; ++i ) 01667 { 01668 Old_histr[i] = struc.histr[i]; 01669 Old_depth[i] = struc.depth[i]; 01670 Old_xLyman_depth[i] = struc.xLyman_depth[i]; 01671 /* old n_p density from previous iteration */ 01672 Old_hiistr[i] = struc.hiistr[i]; 01673 /* old pressure from previous iteration */ 01674 Old_pressure[i] = struc.pressure[i]; 01675 /* old electron density from previous iteration */ 01676 Old_ednstr[i] = struc.ednstr[i]; 01677 /* energy term */ 01678 Old_EnthalpyDensity[i] = EnthalpyDensity[i]; 01679 /* >>chng 02 May 2001 rjrw: add hden for dilution */ 01680 Old_hden[i] = struc.hden[i]; 01681 Old_DenMass[i] = struc.DenMass[i]; 01682 01683 for(mol=0;mol<N_H_MOLEC;mol++) 01684 { 01685 Old_H2_molec[i][mol] = struc.H2_molec[mol][i]; 01686 } 01687 for(mol=0;mol<mole.num_comole_calc;mol++) 01688 { 01689 Old_CO_molec[i][mol] = struc.CO_molec[mol][i]; 01690 } 01691 01692 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 01693 { 01694 Old_gas_phase[i][nelem] = dense.gas_phase[nelem]; 01695 for( ion=0; ion<nelem+2; ++ion ) 01696 { 01697 Old_xIonDense[i][nelem][ion] = struc.xIonDense[nelem][ion][i]; 01698 } 01699 } 01700 } 01701 return; 01702 } 01703 01704 realnum DynaFlux(double depth) 01705 01706 { 01707 realnum flux; 01708 01709 DEBUG_ENTRY( "DynaFlux()" ); 01710 01711 if(dynamics.FluxIndex == 0) 01712 { 01713 flux = (realnum)dynamics.FluxScale; 01714 } 01715 else 01716 { 01717 flux = (realnum)(dynamics.FluxScale*pow(fabs(depth-dynamics.FluxCenter),dynamics.FluxIndex)); 01718 if(depth < dynamics.FluxCenter) 01719 flux = -flux; 01720 } 01721 if(dynamics.lgFluxDScale) 01722 { 01723 /*flux *= struc.DenMass[0]; */ 01724 /* WJH 21 may 04, changed to use dense.xMassDensity0, which should be strictly constant */ 01725 flux *= dense.xMassDensity0; 01726 } 01727 return flux; 01728 } 01729 01730 /* ============================================================================== */ 01731 /*DynaZero zero some dynamics variables, called from zero.c, 01732 * before parsing commands */ 01733 void DynaZero( void ) 01734 { 01735 int ipISO; 01736 01737 DEBUG_ENTRY( "DynaZero()" ); 01738 01739 /* the number of zones in the previous iteration */ 01740 nOld_zone = 0; 01741 01742 /* by default advection is turned off */ 01743 dynamics.lgAdvection = false; 01744 /*dynamics.Velocity = 0.;*/ 01745 AdvecSpecificEnthalpy = 0.; 01746 dynamics.Cool = 0.; 01747 dynamics.Heat = 0.; 01748 dynamics.dCooldT = 0.; 01749 dynamics.dHeatdT = 0.; 01750 dynamics.HeatMax = 0.; 01751 dynamics.CoolMax = 0.; 01752 dynamics.Rate = 0.; 01753 01754 /* sets recombination logic, keyword RECOMBINATION on a time step line */ 01755 dynamics.lgRecom = false; 01756 01757 /* set true if time dependent calculation is finished */ 01758 dynamics.lgStatic_completed = false; 01759 01760 /* vars that determine whether time dependent soln only - set with time command */ 01761 dynamics.lgStatic = false; 01762 timestep_init = -1.; 01763 /* this factor multiplies the time step */ 01764 timestep_factor = 1.2; 01765 dynamics.time_elapsed = 0.; 01766 01767 /* set the first iteration to include dynamics rather than constant pressure */ 01768 /* iteration number, initial iteration is 1, default is 2 - changed with SET DYNAMICS FIRST command */ 01769 dynamics.n_initial_relax = 2; 01770 01771 /* set initial value of the advection length, 01772 * neg => fraction of depth of init model, + length cm */ 01773 dynamics.AdvecLengthInit = -0.1; 01774 01775 /* this is a tolerance for determining whether dynamics has converged */ 01776 dynamics.convergence_tolerance = 0.1; 01777 01778 /* this says that set dynamics pressure mode was set */ 01779 dynamics.lgSetPresMode = false; 01780 01781 /* set default values for uniform mass flux */ 01782 dynamics.FluxScale = 0.; 01783 dynamics.lgFluxDScale = false; 01784 dynamics.FluxCenter = 0.; 01785 dynamics.FluxIndex = 0.; 01786 dynamics.dRad = BIGFLOAT; 01787 01788 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 01789 { 01790 /* factor to allow turning off advection for one of the iso seq, 01791 * this is done with command "no advection h-like" or "he-like" 01792 * only for testing */ 01793 dynamics.lgISO[ipISO] = true; 01794 } 01795 /* turn off advection for rest of ions, command "no advection metals" */ 01796 dynamics.lgMETALS = true; 01797 /* turn off thermal effects of advection, command "no advection cooling" */ 01798 dynamics.lgCoolHeat = true; 01799 DivergePresInteg = 0.; 01800 01801 dynamics.discretization_error = 0.; 01802 dynamics.error_scale2 = 0.; 01803 return; 01804 } 01805 01806 01807 /* ============================================================================== */ 01808 /* DynaCreateArrays allocate some space needed to save the dynamics structure variables, 01809 * called from DynaCreateArrays */ 01810 void DynaCreateArrays( void ) 01811 { 01812 long int nelem, 01813 ns, 01814 i, 01815 ion, 01816 mol; 01817 01818 DEBUG_ENTRY( "DynaCreateArrays()" ); 01819 01820 Upstream_H2_molec = (double*)MALLOC((size_t)N_H_MOLEC*sizeof(double) ); 01821 Upstream_CO_molec = (double*)MALLOC((size_t)mole.num_comole_calc*sizeof(double) ); 01822 01823 dynamics.H2_molec = (double*)MALLOC((size_t)N_H_MOLEC*sizeof(double) ); 01824 dynamics.CO_molec = (double*)MALLOC((size_t)mole.num_comole_calc*sizeof(double) ); 01825 01826 UpstreamElem = (double*)MALLOC((size_t)LIMELM*sizeof(double) ); 01827 01828 dynamics.Source = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) )); 01829 UpstreamIon = ((double**)MALLOC( (size_t)LIMELM*sizeof(double *) )); 01830 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 01831 { 01832 dynamics.Source[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) )); 01833 UpstreamIon[nelem] = ((double*)MALLOC( (size_t)(nelem+2)*sizeof(double) )); 01834 for( ion=0; ion<nelem+2; ++ion ) 01835 { 01836 dynamics.Source[nelem][ion] = 0.; 01837 } 01838 } 01839 dynamics.Rate = 0.; 01840 01841 Old_EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01842 01843 Old_ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01844 01845 EnthalpyDensity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01846 01847 Old_DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01848 01849 Old_hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01850 01851 Old_pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01852 01853 Old_histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01854 01855 Old_hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01856 01857 Old_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01858 01859 Old_xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 01860 01861 Old_xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) ); 01862 01863 Old_gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) ); 01864 01865 Old_H2_molec = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) ); 01866 01867 Old_CO_molec = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) ); 01868 01869 /* now create diagonal of space for ionization arrays */ 01870 for( ns=0; ns < struc.nzlim; ++ns ) 01871 { 01872 Old_xIonDense[ns] = 01873 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+3) ); 01874 01875 Old_gas_phase[ns] = 01876 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+3) ); 01877 01878 Old_H2_molec[ns] = 01879 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(N_H_MOLEC) ); 01880 01881 Old_CO_molec[ns] = 01882 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(mole.num_comole_calc) ); 01883 01884 for( nelem=0; nelem< (LIMELM+3);++nelem ) 01885 { 01886 Old_xIonDense[ns][nelem] = 01887 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+1) ); 01888 } 01889 } 01890 01891 for( i=0; i < struc.nzlim; i++ ) 01892 { 01893 /* these are values if H0 and tau_912 from previous iteration */ 01894 Old_histr[i] = 0.; 01895 Old_xLyman_depth[i] = 0.; 01896 Old_depth[i] = 0.; 01897 dynamics.oldFullDepth = 0.; 01898 /* old n_p density from previous iteration */ 01899 Old_hiistr[i] = 0.; 01900 /* old pressure from previous iteration */ 01901 Old_pressure[i] = 0.; 01902 /* old electron density from previous iteration */ 01903 Old_ednstr[i] = 0.; 01904 Old_hden[i] = 0.; 01905 Old_DenMass[i] = 0.; 01906 Old_EnthalpyDensity[i] = 0.; 01907 for( nelem=0; nelem< (LIMELM+3);++nelem ) 01908 { 01909 for( ion=0; ion<LIMELM+1; ++ion ) 01910 { 01911 Old_xIonDense[i][nelem][ion] = 0.; 01912 } 01913 } 01914 for(mol=0;mol<N_H_MOLEC;mol++) 01915 { 01916 Old_H2_molec[i][mol] = 0.; 01917 } 01918 for(mol=0;mol<mole.num_comole_calc;mol++) 01919 { 01920 Old_CO_molec[i][mol] = 0.; 01921 } 01922 } 01923 return; 01924 } 01925 01926 /*advection_set_detault - called to set default conditions 01927 * when time and wind commands are parsed, 01928 * lgWind is true if dynamics, false if time dependent */ 01929 STATIC void advection_set_detault( bool lgWind ) 01930 { 01931 01932 DEBUG_ENTRY( "advection_set_detault()" ); 01933 01934 /* turn on advection */ 01935 dynamics.lgAdvection = true; 01936 01937 /* turn off prediction of next zone's temperature, as guessed in ZoneStart, 01938 * also set with no tepredictor */ 01939 thermal.lgPredNextTe = false; 01940 01941 /* turn off both CO and H2 networks since advection not included, also say physical 01942 * conditions are not ok*/ 01943 /* the heavy element molecules are turned off for the moment, 01944 * until I am confident that the H-only models work robustly */ 01945 /* co.lgNoH2Mole = true; */ 01947 /* >>chng 06 jun 29, add advective terms to CO solver */ 01948 co.lgNoCOMole = true; 01949 # if 0 01950 co.lgNoCOMole = true; /* >>chng 04 apr 23, tried to rm this line - problems */ 01951 phycon.lgPhysOK = false;/* >>chng 04 apr 23, rm this line */ 01952 # endif 01953 /* >>chng 06 nov 29, there is a conservation problem in the ionization - 01954 * molecular solvers that is demonstrated by the */ 01957 /* use the new temperature solver 01958 strcpy( conv.chSolverEden , "new" ); */ 01959 01960 /* constant total pressure, gas+rad+incident continuum 01961 * turn on radiation pressure */ 01962 pressure.lgPres_radiation_ON = true; 01963 pressure.lgPres_magnetic_ON = true; 01964 pressure.lgPres_ram_ON = true; 01965 01966 /* we need to make the solvers much more exact when advection is in place */ 01967 if( lgWind ) 01968 { 01969 /* increse precision of solution */ 01970 conv.EdenErrorAllowed = 1e-3; 01971 /* the actual relative error is relative to the total heating and cooling, 01972 * which include the dynamics.heat and .cool, which are the advected heating/cooling. 01973 * the two terms can be large and nearly cancel, what is written to the .heat and cool files 01974 * by punch files has had the smaller of the two subtracted, leaving only the net advected 01975 * heating and cooling */ 01976 conv.HeatCoolRelErrorAllowed = 3e-4f; 01977 conv.PressureErrorAllowed = 1e-3f; 01978 } 01979 return; 01980 } 01981 01982 /* ============================================================================== */ 01983 /* ParseDynaTime parse the time command, called from ParseCommands */ 01984 void ParseDynaTime( char *chCard ) 01985 { 01986 long int i; 01987 bool lgEOL; 01988 char chCap[INPUT_LINE_LENGTH]; 01989 01990 DEBUG_ENTRY( "ParseDynaTime()" ); 01991 01992 /*flag set true when time dependent only */ 01993 dynamics.lgStatic = true; 01994 01995 i = 5; 01996 timestep_init = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01997 /* this is log of timestep */ 01998 if( timestep_init > 30. ) 01999 { 02000 /* this number is large and will almost certainly crash */ 02001 fprintf(ioQQQ,"WARNING - the log of the initial time step is too " 02002 "large, I shall probably crash. The value was log t = %.2e\n", 02003 timestep_init ); 02004 fflush(ioQQQ); 02005 } 02006 timestep_init = pow( 10. , timestep_init ); 02007 timestep = timestep_init; 02008 if( lgEOL ) 02009 { 02010 NoNumb(chCard); 02011 } 02012 02013 /* this is the stop time and is optional */ 02014 timestep_stop = pow( 10. , FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) ); 02015 if( lgEOL ) 02016 { 02017 timestep_stop = -1.; 02018 } 02019 02020 /* set default flags - false says that time dependent, not dynamical solution */ 02021 advection_set_detault(false); 02022 02023 wind.windv0 = 0.; 02024 wind.windv = wind.windv0; 02025 02026 /* this is used in convpres to say wind solution - both cases use this*/ 02027 strcpy( dense.chDenseLaw, "WIND" ); 02028 02029 /* create time step and flux arrays */ 02030 time_elapsed_time = (double*)MALLOC((size_t)NTIME*sizeof(double)); 02031 time_flux_ratio = (double*)MALLOC((size_t)NTIME*sizeof(double)); 02032 time_dt = (double*)MALLOC((size_t)NTIME*sizeof(double)); 02033 time_dt_scale_factor = (double*)MALLOC((size_t)NTIME*sizeof(double)); 02034 lgtime_Recom = (int*)MALLOC((size_t)NTIME*sizeof(int)); 02035 02036 /* number of lines we will save */ 02037 nTime_flux = 0; 02038 02039 /* get the next line, and check for eof */ 02040 input_readarray(chCard,&lgEOL); 02041 if( lgEOL ) 02042 { 02043 fprintf( ioQQQ, 02044 " Hit EOF while reading time-continuum list; use END to end list.\n" ); 02045 cdEXIT(EXIT_FAILURE); 02046 } 02047 02048 /* convert line to caps */ 02049 strcpy( chCap, chCard ); 02050 caps(chCap); 02051 02052 /* third column might set dt - if any third column is missing then 02053 * this is set false and only time on command line is used */ 02054 lgtime_dt_specified = true; 02055 02056 while( strncmp(chCap, "END" ,3 ) != 0 ) 02057 { 02058 if( nTime_flux >= NTIME ) 02059 { 02060 fprintf( ioQQQ, 02061 " Too many time points have been entered; the limit is %d. Increase variable NTIME in dynamics.c.\n", 02062 NTIME ); 02063 cdEXIT(EXIT_FAILURE); 02064 } 02065 02066 i = 1; 02067 time_elapsed_time[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 02068 time_flux_ratio[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 02069 if( lgEOL ) 02070 { 02071 fprintf( ioQQQ, 02072 " each line must have two numbers, log of time (s0 and flux ratio.\n"); 02073 cdEXIT(EXIT_FAILURE); 02074 } 02075 /* this is optional dt to set time step - if not given then initial 02076 * time step is always used */ 02077 time_dt[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 02078 02079 /* if any of these are not specified then do not use times array */ 02080 if( lgEOL ) 02081 lgtime_dt_specified = false; 02082 02083 /* this is optional scale factor to increase time */ 02084 time_dt_scale_factor[nTime_flux] = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 02085 if( lgEOL ) 02086 time_dt_scale_factor[nTime_flux] = -1.; 02087 02088 /* turn on recombination front logic */ 02089 if( nMatch("RECO" , chCap ) && nMatch("MBIN" , chCap ) ) 02090 { 02091 /* this sets flag dynamics.lgRecom true so that all of code knows recombination 02092 * is taking place */ 02093 lgtime_Recom[nTime_flux] = true; 02094 } 02095 else 02096 { 02097 lgtime_Recom[nTime_flux] = false; 02098 } 02099 02100 /* this is total number stored so far */ 02101 ++nTime_flux; 02102 02103 /* get next line and check for eof */ 02104 input_readarray(chCard,&lgEOL); 02105 if( lgEOL ) 02106 { 02107 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" ); 02108 cdEXIT(EXIT_FAILURE); 02109 } 02110 02111 /* convert it to caps */ 02112 strcpy( chCap, chCard ); 02113 caps(chCap); 02114 } 02115 02116 for( i=0; i < nTime_flux; i++ ) 02117 { 02118 fprintf( ioQQQ, "DEBUG time dep %.2e %.2e %.2e %.2e\n", 02119 time_elapsed_time[i], 02120 time_flux_ratio[i] , 02121 time_dt[i], 02122 time_dt_scale_factor[i]); 02123 } 02124 fprintf( ioQQQ, "\n" ); 02125 return; 02126 } 02127 /* ============================================================================== */ 02128 /* ParseDynaWind parse the wind command, called from ParseCommands */ 02129 void ParseDynaWind( char *chCard ) 02130 { 02131 long int i; 02132 int iVelocity_Type; 02133 bool lgEOL; 02134 /* compiler flagged possible paths where dfdr used but not set - 02135 * this is for safety/keep it happy */ 02136 double dfdr=-BIGDOUBLE; 02137 02138 DEBUG_ENTRY( "ParseDynaWind()" ); 02139 02140 /* Flag for type of velocity law: 02141 * 1 is original, give initial velocity at illuminated face 02142 * 2 is face flux gradient (useful if face velocity is zero), 02143 * set to zero, but will be reset if velocity specified */ 02144 iVelocity_Type = 0; 02145 /* wind structure, parameters are initial velocity and optional mass 02146 * v read in in km s-1 and convert to cm s-1, mass in solar masses */ 02147 if( (i = nMatch("VELO",chCard))>0 ) 02148 { 02149 i += 5; 02150 wind.windv0 = (realnum)(FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)*1e5); 02151 wind.windv = wind.windv0; 02152 iVelocity_Type = 1; 02153 } 02154 02155 if( (i = nMatch("DFDR",chCard))>0 ) 02156 { 02157 /* velocity not specified, rather mass flux gradient */ 02158 i += 5; 02159 dfdr = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02160 iVelocity_Type = 2; 02161 } 02162 02163 /* center option, gives xxx */ 02164 if( (i = nMatch("CENT",chCard))>0 ) 02165 { 02166 /* physical length in cm, can be either sign */ 02167 i += 5; 02168 dynamics.FluxCenter = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02169 } 02170 02171 /* flux index */ 02172 if( (i = nMatch("INDE",chCard))>0 ) 02173 { 02174 /* power law index */ 02175 i += 5; 02176 dynamics.FluxIndex = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02177 } 02178 02179 /* the case where velocity was set */ 02180 if(iVelocity_Type == 1) 02181 { 02182 /* was flux index also set? */ 02183 if(dynamics.FluxIndex == 0) 02184 { 02185 dynamics.FluxScale = wind.windv0; 02186 dynamics.lgFluxDScale = true; 02187 /* Center doesn't mean much in this case -- make sure it's 02188 * in front of grid so DynaFlux doesn't swap signs where 02189 * it shouldn't */ 02190 dynamics.FluxCenter = -1.; 02191 } 02192 else 02193 { 02196 /* velocity was set but flux index was not set - estimate it */ 02197 dynamics.FluxScale = wind.windv0* 02198 pow(fabs(dynamics.FluxCenter),-dynamics.FluxIndex); 02199 02200 dynamics.lgFluxDScale = true; 02201 if(dynamics.FluxCenter > 0) 02202 { 02203 dynamics.FluxScale = -dynamics.FluxScale; 02204 } 02205 } 02206 } 02207 /* the case where flux gradient is set */ 02208 else if(iVelocity_Type == 2) 02209 { 02210 if(dynamics.FluxIndex == 0) 02211 { 02212 fprintf(ioQQQ,"Can't specify gradient when flux is constant!\n"); 02213 /* use this exit handler, which closes out MPI when multiprocessing */ 02214 cdEXIT(EXIT_FAILURE); 02215 } 02218 /* Can't specify FluxScale from dvdr rather than dfdr, as 02219 * d(rho)/dr != 0 */ 02220 dynamics.FluxScale = dfdr/dynamics.FluxIndex* 02221 pow(fabs(dynamics.FluxCenter),1.-dynamics.FluxIndex); 02222 if(dynamics.FluxCenter > 0) 02223 { 02224 dynamics.FluxScale = -dynamics.FluxScale; 02225 } 02226 dynamics.lgFluxDScale = false; 02227 02228 /* put in bogus value simply as flag -- assume that surface velocity 02229 * is small or we wouldn't be using this to specify. */ 02230 wind.windv0 = -0.01f; 02231 } 02232 else 02233 { 02234 /* assume old-style velocity-only specification */ 02235 /* wind structure, parameters are initial velocity and optional mass 02236 * v in km/sec, mass in solar masses */ 02237 i = 5; 02238 wind.windv0 = (realnum)(FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)*1e5); 02239 dynamics.FluxScale = wind.windv0; 02240 dynamics.FluxIndex = 0.; 02241 dynamics.lgFluxDScale = true; 02242 /* Center doesn't mean much in this case -- make sure it's 02243 * in front of grid so DynaFlux doesn't swap signs where 02244 * it shouldn't */ 02245 dynamics.FluxCenter = -1.; 02246 if( lgEOL ) 02247 { 02248 NoNumb(chCard); 02249 } 02250 } 02251 02252 wind.windv = wind.windv0; 02253 02254 # ifdef FOO 02255 fprintf(ioQQQ,"Scale %g (*%c) Index %g Center %g\n", 02256 dynamics.FluxScale,(dynamics.lgFluxDScale)?'D':'1', 02257 dynamics.FluxIndex,dynamics.FluxCenter); 02258 # endif 02259 02260 /* option to include advection */ 02261 if( nMatch( "ADVE" , chCard ) ) 02262 { 02263 /* set default flags - true says dynamical solution */ 02264 advection_set_detault(true); 02265 } 02266 02267 else 02268 { 02269 /* this is usual hypersonic outflow */ 02270 if( wind.windv0 <= 1.e6 ) 02271 { 02272 /* speed of sound roughly 10 km/s */ 02273 fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" ); 02274 wind.lgWindOK = false; 02275 } 02276 02277 /* set the central object mass, in solar masses */ 02278 wind.comass = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02279 /* default is one solar mass */ 02280 if( lgEOL ) 02281 wind.comass = 1.; 02282 02283 /* option for rotating disk, keyword is disk */ 02284 wind.lgDisk = false; 02285 if( nMatch( "DISK" , chCard ) ) 02286 wind.lgDisk = true; 02287 02288 } 02289 02290 /* this is used in convpres to say wind solution - both cases use this*/ 02291 strcpy( dense.chDenseLaw, "WIND" ); 02292 02293 /* option to turn off continuum radiative acceleration */ 02294 if( nMatch("O CO",chCard) ) 02295 { 02296 pressure.lgContRadPresOn = false; 02297 } 02298 else 02299 { 02300 pressure.lgContRadPresOn = true; 02301 } 02302 return; 02303 } 02304 02305 /*DynaPrtZone called to print zone results */ 02306 void DynaPrtZone( void ) 02307 { 02308 02309 DEBUG_ENTRY( "DynaPrtZone()" ); 02310 02311 ASSERT( nzone>0 && nzone<struc.nzlim ); 02312 02313 if( nzone > 0 ) 02314 { 02315 fprintf(ioQQQ," DYNAMICS Advection: Uad %.2f Uwd%.2e FRCcool: %4.2f Heat %4.2f\n", 02316 timesc.sound_speed_adiabatic/1e5 , 02317 wind.windv/1e5 , 02318 dynamics.Cool/thermal.ctot, 02319 dynamics.Heat/thermal.ctot); 02320 } 02321 02322 ASSERT( EnthalpyDensity[nzone-1] > 0. ); 02323 02324 fprintf(ioQQQ," DYNAMICS Eexcit:%.4e Eion:%.4e Ebin:%.4e Ekin:%.4e ET+pdv %.4e EnthalpyDensity/rho%.4e AdvSpWork%.4e\n", 02325 phycon.EnergyExcitation, 02326 phycon.EnergyIonization, 02327 phycon.EnergyBinding, 02328 0.5*POW2(wind.windv)*dense.xMassDensity, 02329 5./2.*pressure.PresGasCurr , 02330 EnthalpyDensity[nzone-1]/dense.gas_phase[ipHYDROGEN] , AdvecSpecificEnthalpy 02331 ); 02332 return; 02333 } 02334 02335 /*DynaPunchTimeDep - punch info about time dependent solution */ 02336 void DynaPunchTimeDep( FILE* ipPnunit , const char *chJob ) 02337 { 02338 02339 DEBUG_ENTRY( "DynaPunchTimeDep()" ); 02340 02341 if( strcmp( chJob , "END" ) == 0 ) 02342 { 02343 double te_mean, 02344 H2_mean, 02345 H0_mean, 02346 Hp_mean, 02347 Hep_mean; 02348 /* punch info at end */ 02349 if( cdTemp( 02350 /* four char string, null terminated, giving the element name */ 02351 "HYDR", 02352 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 02353 * 0 means that chLabel is a special case */ 02354 2, 02355 /* will be temperature */ 02356 &te_mean, 02357 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02358 "RADIUS" ) ) 02359 { 02360 TotalInsanity(); 02361 } 02362 if( cdIonFrac( 02363 /* four char string, null terminated, giving the element name */ 02364 "HYDR", 02365 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 02366 * 0 says special case */ 02367 2, 02368 /* will be fractional ionization */ 02369 &Hp_mean, 02370 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02371 "RADIUS" , 02372 /* if true then weighting also has electron density, if false then only volume or radius */ 02373 false ) ) 02374 { 02375 TotalInsanity(); 02376 } 02377 if( cdIonFrac( 02378 /* four char string, null terminated, giving the element name */ 02379 "HYDR", 02380 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 02381 * 0 says special case */ 02382 1, 02383 /* will be fractional ionization */ 02384 &H0_mean, 02385 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02386 "RADIUS" , 02387 /* if true then weighting also has electron density, if false then only volume or radius */ 02388 false ) ) 02389 { 02390 TotalInsanity(); 02391 } 02392 if( cdIonFrac( 02393 /* four char string, null terminated, giving the element name */ 02394 "H2 ", 02395 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 02396 * 0 says special case */ 02397 0, 02398 /* will be fractional ionization */ 02399 &H2_mean, 02400 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02401 "RADIUS" , 02402 /* if true then weighting also has electron density, if false then only volume or radius */ 02403 false ) ) 02404 { 02405 TotalInsanity(); 02406 } 02407 if( cdIonFrac( 02408 /* four char string, null terminated, giving the element name */ 02409 "HELI", 02410 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 02411 * 0 says special case */ 02412 2, 02413 /* will be fractional ionization */ 02414 &Hep_mean, 02415 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02416 "RADIUS" , 02417 /* if true then weighting also has electron density, if false then only volume or radius */ 02418 false ) ) 02419 { 02420 TotalInsanity(); 02421 } 02422 fprintf( ipPnunit , 02423 "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n" , 02424 dynamics.time_elapsed , 02425 timestep , 02426 rfield.time_continuum_scale , 02427 dense.gas_phase[ipHYDROGEN], 02428 te_mean , 02429 Hp_mean , 02430 H0_mean , 02431 H2_mean , 02432 Hep_mean , 02433 /* ratio of CO to total H column densities */ 02434 findspecies("CO")->hevcol / SDIV( colden.colden[ipCOL_HTOT] )); 02435 } 02436 else 02437 TotalInsanity(); 02438 return; 02439 } 02440 02441 /*DynaPunch punch dynamics - info related to advection */ 02442 void DynaPunch(FILE* ipPnunit , char chJob ) 02443 { 02444 DEBUG_ENTRY( "DynaPunch()" ); 02445 02446 if( chJob=='a' ) 02447 { 02448 /* this is punch dynamics advection, the only punch dynamics */ 02449 fprintf( ipPnunit , "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 02450 radius.depth_mid_zone, 02451 thermal.htot , 02452 dynamics.Cool , 02453 dynamics.Heat , 02454 dynamics.dCooldT , 02455 dynamics.Source[ipHYDROGEN][ipHYDROGEN], 02456 dynamics.Rate, 02457 phycon.EnthalpyDensity/dense.gas_phase[ipHYDROGEN] , 02458 AdvecSpecificEnthalpy 02459 ); 02460 } 02461 else 02462 TotalInsanity(); 02463 return; 02464 }