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 /*ConvInitSolution drive search for initial temperature, for illuminated face */ 00004 /*lgTempTooHigh returns true if temperature is too high */ 00005 /*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */ 00006 #include "cddefines.h" 00007 #include "physconst.h" 00008 #include "trace.h" 00009 #include "struc.h" 00010 #include "rfield.h" 00011 #include "mole.h" 00012 #include "dense.h" 00013 #include "stopcalc.h" 00014 #include "heavy.h" 00015 #include "wind.h" 00016 #include "geometry.h" 00017 #include "thermal.h" 00018 #include "radius.h" 00019 #include "phycon.h" 00020 #include "pressure.h" 00021 #include "conv.h" 00022 00023 /* derivative of net cooling wrt temperature to check on sign oscillations */ 00024 static double dCoolNetDTOld = 0; 00025 00026 static double OxyInGrains , FracMoleMax; 00027 /* determine whether chemistry is important - what is the largest 00028 * fraction of an atom in molecules? Also is ice formation 00029 * important 00030 * sets above two file static variables */ 00031 00032 /*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */ 00033 STATIC bool lgCoolHeatCheckConverge( 00034 double *CoolNet ) 00035 { 00036 static double HeatOld=-1. , CoolOld=-1.; 00037 DEBUG_ENTRY( "lgCoolHeatCheckConverge()" ); 00038 00039 double HeatChange = thermal.htot - HeatOld; 00040 double CoolChange = thermal.ctot - CoolOld; 00041 /* will use larger of heating or cooling in tests for relative convergence 00042 * because in constant temperature models one may have trivially small value */ 00043 double HeatCoolMax = MAX2( thermal.htot , thermal.ctot ); 00044 00045 /* is the heating / cooling converged? 00046 * get max relative change, determines whether converged heating/cooling */ 00047 double HeatRelChange = fabs(HeatChange)/SDIV( HeatCoolMax ); 00048 double CoolRelChange = fabs(CoolChange)/SDIV( HeatCoolMax ); 00049 bool lgConverged = true; 00050 if( MAX2( HeatRelChange , CoolRelChange ) > conv.HeatCoolRelErrorAllowed ) 00051 lgConverged = false; 00052 00053 *CoolNet = thermal.ctot - thermal.htot; 00054 00055 HeatOld = thermal.htot; 00056 CoolOld = thermal.ctot; 00057 return lgConverged; 00058 } 00059 00060 /* call lgCoolHeatCheckConverge until cooling and heating are converged */ 00061 STATIC bool lgCoolNetConverge( 00062 double *CoolNet , 00063 double *dCoolNetDT ) 00064 { 00065 const long int LOOP_MAX=20; 00066 long int loop = 0; 00067 bool lgConvergeCoolHeat = false; 00068 00069 DEBUG_ENTRY( "lgCoolNetConverge()" ); 00070 00071 while( loop < LOOP_MAX && !lgConvergeCoolHeat && !lgAbort ) 00072 { 00073 if( ConvEdenIoniz() ) 00074 { 00075 /* this is an error return, calculation will immediately stop */ 00076 lgAbort = true; 00077 } 00078 lgConvergeCoolHeat = lgCoolHeatCheckConverge( CoolNet ); 00079 if( trace.lgTrace || trace.nTrConvg>=3 ) 00080 fprintf(ioQQQ," lgCoolNetConverge %li calls to ConvEdenIoniz, converged? %c\n", 00081 loop , TorF( lgConvergeCoolHeat ) ); 00082 ++loop; 00083 } 00084 00085 if( thermal.ConstTemp > 0 ) 00086 { 00087 /* constant temperature model - this is trick so that temperature 00088 * updater uses derivative to go to the set constant temperature */ 00089 *CoolNet = phycon.te - thermal.ConstTemp; 00090 *dCoolNetDT = 1.f; 00091 } 00092 else if( phycon.te < 4000.f ) 00093 { 00094 /* at low temperatures the analytical cooling-heating derivative 00095 * is often of no value - the species densities change when the 00096 * temperature changes. Use simple dCnet/dT=(C-H)/T - this is 00097 * usually too large and results is smaller than optical dT, but 00098 * we do eventually converge */ 00099 *dCoolNetDT = thermal.ctot / phycon.te; 00100 } 00101 else if( thermal.htot / thermal.ctot > 2.f ) 00102 { 00103 /* we are far from the solution and the temperature is much lower than 00104 * equilibrium - analytical derivative from cooling evaluation is probably 00105 * wrong since coolants are not correct */ 00106 *dCoolNetDT = thermal.ctot / phycon.te; 00107 } 00108 else 00109 { 00110 *dCoolNetDT = thermal.dCooldT - thermal.dHeatdT; 00111 if( dCoolNetDTOld * *dCoolNetDT < 0. ) 00112 { 00113 /* derivative is oscillating */ 00114 *dCoolNetDT = *CoolNet / phycon.te; 00115 fprintf(ioQQQ,"PROBLEM CoolNet derivative oscillating in initial solution \n"); 00116 } 00117 } 00118 return lgConvergeCoolHeat; 00119 } 00120 00121 /*ChemImportance find fraction of heavy elements in molecules, O in ices */ 00122 STATIC void ChemImportance( void ) 00123 { 00124 long int i; 00125 00126 DEBUG_ENTRY( "ChemImportance()" ); 00127 00128 FracMoleMax = 0.; 00129 long int ipMax = -1; 00130 for( i=0; i<mole.num_comole_calc; ++i ) 00131 { 00132 /* is element included in the chemistry, and is it a molecule? */ 00133 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] && COmole[i]->n_nuclei>1 ) 00134 { 00135 realnum frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]); 00136 frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol]; 00137 if( frac > FracMoleMax ) 00138 { 00139 /* remember largest fraction */ 00140 FracMoleMax = frac; 00141 ipMax = i; 00142 } 00143 } 00144 } 00145 00146 /* fraction of all oxygen in ices on grains */ 00147 OxyInGrains = 0; 00148 for(i=0;i<mole.num_comole_calc;++i) 00149 { 00150 OxyInGrains += (1 - COmole[i]->lgGas_Phase)*COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN]; 00151 } 00152 /* this is now fraction of O in ices */ 00153 OxyInGrains /= SDIV(dense.gas_phase[ipOXYGEN]); 00154 00155 { 00156 /* option to print out entire matrix */ 00157 enum {DEBUG_LOC=false}; 00158 if( DEBUG_LOC ) 00159 { 00160 fprintf(ioQQQ,"DEBUG fraction molecular %.2e species %li O ices %.2e\n" , 00161 FracMoleMax , ipMax ,OxyInGrains ); 00162 } 00163 } 00164 00165 return; 00166 } 00167 00168 double FindTempChangeFactor(void) 00169 { 00170 double TeChangeFactor; 00171 00172 DEBUG_ENTRY( "FindTempChangeFactor()" ); 00173 00174 /* find fraction of atoms in molecules - need small changes 00175 * in temperature if fully molecular since chemistry 00176 * network is complex and large changes in solution would 00177 * cause linearization to fail */ 00178 /* this evaluates the global variables FracMoleMax and 00179 * OxyInGrains */ 00180 ChemImportance(); 00181 00182 /* Following are series of chemical 00183 * tests that determine the factor by which 00184 * which can change the temperature. must be VERY small when 00185 * ice formation on grains is important due to exponential sublimation 00186 * rate. Also small when molecules are important, to keep chemistry 00187 * within limits of linearized solver 00188 * the complete linearization that is implicit in all these solvers 00189 * will not work when large Te jumps occur when molecules/ices 00190 * are important - solvers can't return to solution if too far away 00191 * keep temp steps small when mole/ice is important */ 00192 if( OxyInGrains > 0.99 ) 00193 { 00194 TeChangeFactor = 0.999; 00195 } 00196 else if( OxyInGrains > 0.9 ) 00197 { 00198 TeChangeFactor = 0.99; 00199 } 00200 /* >>chng 97 mar 3, make TeChangeFactor small when close to lower bound */ 00201 else if( phycon.te < 5. || 00202 /*>>chng 06 jul 30, or if molecules/ices are important */ 00203 OxyInGrains > 0.1 ) 00204 { 00205 TeChangeFactor = 0.97; 00206 } 00207 /*>>chng 07 feb 23, add test of chemistry being very important */ 00208 else if( phycon.te < 20. || FracMoleMax > 0.1 ) 00209 { 00210 /* >>chng 07 feb 24, changed from 0,9 to 0.95 to converge 00211 * pdr_coolbb.in test */ 00212 TeChangeFactor = 0.95; 00213 } 00214 else 00215 { 00216 /* this is the default largest step */ 00217 TeChangeFactor = 0.8; 00218 } 00219 return TeChangeFactor; 00220 } 00221 00222 /*ConvInitSolution drive search for initial temperature, for illuminated face */ 00223 int ConvInitSolution(void) 00224 { 00225 long int i, 00226 ionlup, 00227 nelem , 00228 nflux_old, 00229 nelem_reflux , 00230 ion_reflux; 00231 /* current value of Cooling - Heating */ 00232 bool lgConvergeCoolHeat; 00233 00234 double 00235 TeChangeFactor, 00236 TeBoundLow, 00237 TeBoundHigh; 00238 00239 DEBUG_ENTRY( "ConvInitSolution()" ); 00240 00241 /* set NaN for safety */ 00242 set_NaN( TeBoundLow ); 00243 set_NaN( TeBoundHigh ); 00244 00245 /* this counts number of times ConvBase is called by PressureChange, in current zone */ 00246 conv.nPres2Ioniz = 0; 00247 /* this counts how many times ConvBase is called in this iteration 00248 * and is flag used by various routines to understand they are 00249 * being called the first time*/ 00250 conv.nTotalIoniz = 0; 00251 /* ots rates not oscillating */ 00252 conv.lgOscilOTS = false; 00253 00254 lgAbort = false; 00255 dense.lgEdenBad = false; 00256 dense.nzEdenBad = 0; 00257 /* these are variables to remember the biggest error in the 00258 * electron density, and the zone in which it occurred */ 00259 conv.BigEdenError = 0.; 00260 conv.AverEdenError = 0.; 00261 conv.BigHeatCoolError = 0.; 00262 conv.AverHeatCoolError = 0.; 00263 conv.BigPressError = 0.; 00264 conv.AverPressError = 0.; 00265 00266 /* these are equal if set dr was used, and we know the first dr */ 00267 if( fp_equal( radius.sdrmin, radius.sdrmax ) ) 00268 { 00269 radius.drad = MIN2( radius.sdrmax , radius.drad ); 00270 radius.drad_x_fillfac = radius.drad * geometry.FillFac; 00271 } 00272 00273 /* the H+ density is zero in sims with no H-ionizing radiation and no cosmic rays, 00274 * the code does work in this limit. But it is unphysical for the H0 density 00275 * to be zero. This could only happen in the fully molecular limit, and the 00276 * code does not go to this limit at this stage, due to simple approximations. */ 00277 ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.); 00278 00279 if( trace.nTrConvg ) 00280 { 00281 fprintf( ioQQQ, "\nConvInitSolution entered \n" ); 00282 } 00283 00284 /******************************************************************** 00285 * 00286 * this is second or higher iteration, reestablish original temperature 00287 * 00288 *********************************************************************/ 00289 if( iteration != 1 ) 00290 { 00291 /* this is second or higher iteration on multi-iteration model */ 00292 if( trace.lgTrace || trace.nTrConvg ) 00293 { 00294 fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n", 00295 iteration, struc.testr[0] ); 00296 } 00297 00298 if( trace.lgTrace || trace.nTrConvg ) 00299 { 00300 fprintf( ioQQQ, " search set true\n" ); 00301 } 00302 00303 /* search phase must be turned on so that variables such as the ots rates, 00304 * secondary ionization, and auger yields, can converge more quickly to 00305 * proper values */ 00306 conv.lgSearch = true; 00307 00308 /* this is the temperature and pressure from the previous iteration */ 00309 /*phycon.te = TempPrevIteration;*/ 00310 /*>>chng 06 jun 09, reset to this rather than init value set in this routine */ 00311 TempChange( struc.testr[0] , false); 00312 00313 /* the initial pressure should now be valid */ 00314 /* this sets values of pressure.PresTotlCurr */ 00315 PresTotCurrent(); 00316 00317 /* now get final temperature */ 00318 if( ConvPresTempEdenIoniz() ) 00319 { 00320 /* this is an error return, calculation will immediately stop */ 00321 lgAbort = true; 00322 return(1); 00323 } 00324 00325 if( trace.lgTrace || trace.nTrConvg ) 00326 { 00327 fprintf( ioQQQ, " ========================================================================================================================\n"); 00328 fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te); 00329 fprintf( ioQQQ, " ========================================================================================================================\n"); 00330 } 00331 conv.lgSearch = false; 00332 00333 /* now get final temperature */ 00334 if( ConvTempEdenIoniz() ) 00335 { 00336 /* this is an error return, calculation will immediately stop */ 00337 lgAbort = true; 00338 return(1); 00339 } 00340 00341 /* the initial pressure should now be valid 00342 * this sets values of pressure.PresTotlCurr */ 00343 PresTotCurrent(); 00344 } 00345 00346 else 00347 { 00348 /******************************************************************** 00349 * 00350 * do first te from scratch 00351 * 00352 *********************************************************************/ 00353 00354 /* say that we are in search phase */ 00355 conv.lgSearch = true; 00356 00357 if( trace.lgTrace ) 00358 { 00359 fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" ); 00360 } 00361 00362 /* coming up to here Te is either 4,000 K (usually) or 10^6 */ 00363 TeBoundLow = phycon.TEMP_LIMIT_LOW/1.001; 00364 00365 /* set temperature floor option - StopCalc.TeFloor is usually 00366 * zero but reset this this command - will go over to constant 00367 * temperature calculation if temperature falls below floor */ 00368 TeBoundLow = MAX2( TeBoundLow , StopCalc.TeFloor ); 00369 00370 /* highest possible temperature - must not get this high since 00371 * parts of code will abort if value is reached. 00372 * divide by 1.2 to prevent checking on temp > 1e10 */ 00373 TeBoundHigh = phycon.TEMP_LIMIT_HIGH/1.2; 00374 00375 /* set initial temperature, options are constant temperature, 00376 * or approach equilibrium from either high or low TE */ 00377 double TeFirst; 00378 if( thermal.ConstTemp > 0 ) 00379 { 00380 /* constant temperature , set 4000 K then relax to desired value 00381 * for cold temperatures, to slowly approach fully molecular solution. 00382 * if constant temperature is higher than 4000., go right to it */ 00383 TeFirst = MAX2( 4000. , thermal.ConstTemp ); 00384 /* true allow ionization stage trimming, false blocks it */ 00385 } 00386 else if( thermal.lgTeHigh ) 00387 { 00388 /* approach from high TE */ 00389 TeFirst = 1e6; 00390 } 00391 else 00392 { 00393 /* approach from low TE */ 00394 TeFirst = 4000.; 00395 } 00396 TempChange(TeFirst , false); 00397 if( trace.lgTrace || trace.nTrConvg>=2 ) 00398 fprintf(ioQQQ,"ConvInitSolution doing initial solution with temp=%.2e\n", 00399 phycon.te); 00400 00401 /* this sets values of pressure.PresTotlCurr */ 00402 PresTotCurrent(); 00403 00404 thermal.htot = 1.; 00405 thermal.ctot = 1.; 00406 00407 /* call cooling, heating, opacity, loop to convergence 00408 * this is very first call to it, by default is at 4000K */ 00409 00410 double CoolNet=0, dCoolNetDT=0; 00411 /* do ionization twice to get stable solution, evaluating initial dR each time */ 00412 lgConvergeCoolHeat = false; 00413 for( ionlup=0; ionlup<2; ++ionlup ) 00414 { 00415 if( trace.lgTrace || trace.nTrConvg>=2 ) 00416 fprintf( ioQQQ, "ConvInitSolution calling lgCoolNetConverge " 00417 "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n", 00418 ionlup,phycon.te , thermal.ctot , thermal.htot,dCoolNetDT ); 00419 /* do not flag oscillating d(C-H)/dT until stable */ 00420 dCoolNetDTOld = 0.; 00421 lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , &dCoolNetDT ); 00422 if( lgAbort ) 00423 break; 00424 /* set thickness of first zone, this affects the pumping rates due 00425 * to correction for attenuation across zone, so must be known 00426 * for ConvEdenIoniz to get valid solution - this is why it 00427 * is in this loop */ 00428 radius_first(); 00429 } 00430 if( !lgConvergeCoolHeat ) 00431 fprintf(ioQQQ," PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n"); 00432 00433 if( lgAbort ) 00434 { 00435 /* we hit an abort */ 00436 fprintf( ioQQQ, " DISASTER PROBLEM ConvInitSolution: Search for " 00437 "initial conditions aborted - lgAbort set true.\n" ); 00438 ShowMe(); 00439 /* this is an error return, calculation will immediately stop */ 00440 return(1); 00441 } 00442 00443 /* we now have error in C-H and its derivative - following is better 00444 * derivative for global case where we may be quite far from C==H */ 00445 if( thermal.ConstTemp<=0 ) 00446 dCoolNetDT = thermal.ctot / phycon.te; 00447 bool lgConvergedLoop = false; 00448 long int LoopThermal = 0; 00449 /* initial temperature is usually 4000K, if the dTe scale factor is 0.95 00450 * it will take 140 integrations to lower temperature to 3K, 00451 * and many more if ices are important. */ 00452 const long int LIMIT_THERMAL_LOOP=300; 00453 double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,CoolSave=-1; 00454 while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP ) 00455 { 00456 /* change temperature until we sign of C-H changes */ 00457 CoolMHeatSave = CoolNet; 00458 dCoolNetDTOld = dCoolNetDT; 00459 CoolSave = thermal.ctot; 00460 TempSave = phycon.te; 00461 00462 /* find temperature change factor, a number less than 1*/ 00463 TeChangeFactor = FindTempChangeFactor(); 00464 ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. ); 00465 00466 TeNew = phycon.te - CoolNet / dCoolNetDT; 00467 00468 TeNew = MAX2( phycon.te*TeChangeFactor , TeNew ); 00469 TeNew = MIN2( phycon.te/TeChangeFactor , TeNew ); 00470 /* change temperature */ 00471 TempChange(TeNew , true); 00472 lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , & dCoolNetDT ); 00473 00474 if( trace.lgTrace || trace.nTrConvg>=2 ) 00475 fprintf(ioQQQ,"ConvInitSolution %4li TeNnew=%.3e " 00476 "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e" 00477 " Ograins %.2e FracMoleMax %.2e\n", 00478 LoopThermal , TeNew , TeChangeFactor , 00479 CoolNet/SDIV(thermal.htot) , dCoolNetDT, 00480 dense.eden , OxyInGrains , FracMoleMax ); 00481 if( lgAbort ) 00482 { 00483 /* we hit an abort */ 00484 fprintf( ioQQQ, " Search for initial conditions aborted - " 00485 "lgAbort set true, lgConvergeCoolHeat=%c.\n", 00486 TorF(lgConvergeCoolHeat) ); 00487 /* this is an error return, calculation will immediately stop */ 00488 return(1); 00489 } 00490 00491 /* keep changing temperature until sign of CoolNet changes 00492 * in constant temperature case CoolNet is delta Temperature 00493 * so is zero for last pass in this loop 00494 * this is for constant temperature case */ 00495 if( fabs(CoolNet)< SMALLDOUBLE ) 00496 /* CoolNet is zero */ 00497 lgConvergedLoop = true; 00498 else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0) 00499 /* change in sign of net cooling */ 00500 lgConvergedLoop = true; 00501 else if( fabs(CoolNet)/MAX2( thermal.ctot , thermal.htot ) <conv.HeatCoolRelErrorAllowed*10. ) 00502 lgConvergedLoop = true; 00503 else if( phycon.te <= TeBoundLow || phycon.te>=TeBoundHigh) 00504 lgConvergedLoop = true; 00505 ++LoopThermal; 00506 } 00507 00508 if( !lgConvergedLoop ) 00509 fprintf(ioQQQ,"PROBLEM ConvInitSolution: temperature solution not " 00510 "found in initial search, final Te=%.2e\n", 00511 phycon.te ); 00512 00513 /* interpolate temperature where C=H if not constant temperature sim 00514 * will have set constant temperature mode above if we hit temperature 00515 * floor */ 00516 if( thermal.ConstTemp<=0 && 00517 ! fp_equal( TempSave , phycon.te ) ) 00518 { 00519 if( trace.lgTrace || trace.nTrConvg>=2 ) 00520 fprintf(ioQQQ," Change in sign of C-H found, Te brackets %.3e " 00521 "%.3e Cool brackets %.3e %.3e ", 00522 TempSave , phycon.te , CoolMHeatSave , CoolNet); 00523 /* interpolate new temperature assuming Cool = a T^alpha, 00524 * first find alpha */ 00525 double alpha = log(thermal.ctot/CoolSave) / log( phycon.te/TempSave); 00526 if( fabs(alpha) < SMALLFLOAT ) 00527 /* alpha close to 0 means constant temperature */ 00528 TeNew = phycon.te; 00529 else 00530 { 00531 /* next find log a - work in logs to avoid infinities */ 00532 if( thermal.ctot<0. || thermal.htot<0. ) 00533 TotalInsanity(); 00534 double Alog = log( thermal.ctot ) - alpha * log( phycon.te ); 00535 /* the interpolated temperature where heating equals cooling */ 00536 double TeLn = (log( thermal.htot ) - Alog) / alpha ; 00537 00538 /* a sanity check */ 00539 if( TeLn < log( MIN2(phycon.te , TempSave )) ) 00540 TeNew = MIN2(phycon.te , TempSave ); 00541 else if( TeLn > log( MAX2(phycon.te , TempSave )) ) 00542 TeNew = MAX2(phycon.te , TempSave ); 00543 else 00544 TeNew = pow( EE , TeLn ); 00545 } 00546 00547 ASSERT( TeNew>=MIN2 ( TempSave , phycon.te ) || 00548 TeNew<=MAX2 ( TempSave , phycon.te )); 00549 00550 if( trace.lgTrace || trace.nTrConvg>=2 ) 00551 fprintf(ioQQQ," interpolating to Te %.3e \n\n", 00552 TeNew); 00553 00554 /* change temperature */ 00555 TempChange(TeNew , true); 00556 } 00557 00558 if( ConvTempEdenIoniz() ) 00559 { 00560 /* this is an error return, calculation will immediately stop */ 00561 lgAbort = true; 00562 return(1); 00563 } 00564 00565 /* this sets values of pressure.PresTotlCurr */ 00566 PresTotCurrent(); 00567 00568 if( trace.lgTrace || trace.nTrConvg ) 00569 { 00570 fprintf( ioQQQ, "\nConvInitSolution: end 1st iteration search phase " 00571 "finding Te=%.3e\n\n" , phycon.te); 00572 } 00573 conv.lgSearch = false; 00574 00575 if( trace.lgTrace ) 00576 { 00577 fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n", 00578 phycon.te ); 00579 } 00580 } 00581 00582 /* we now have a fairly good temperature and ionization 00583 * iteration is 1 on first iteration - remember current pressure 00584 * on first iteration so we can hold this constant in constant 00585 * pressure simulation. 00586 * 00587 * flag dense.lgDenseInitConstant false if 00588 * *constant pressure reset* is used - 00589 * default is true, after first iteration initial density is used for 00590 * first zone no matter what pressure results. Small changes occur due 00591 * to radiative transfer convergence, 00592 * when set false with *reset* option the density on later iterations 00593 * can change to keep the pressure constant */ 00594 if( iteration==1 || dense.lgDenseInitConstant ) 00595 { 00596 double PresNew = pressure.PresTotlCurr; 00597 00598 /* option to specify the pressure as option on constant pressure command */ 00599 if( pressure.lgPressureInitialSpecified ) 00600 /* this is log of nT product - if not present then set zero */ 00601 PresNew = pressure.PressureInitialSpecified; 00602 if( trace.lgTrace ) 00603 { 00604 fprintf( ioQQQ, 00605 " PresTotCurrent 1st zn old values of PresTotlInit:%.3e" 00606 " to %.3e \n", 00607 pressure.PresTotlInit, 00608 PresNew); 00609 } 00610 00611 pressure.PresTotlInit = PresNew; 00612 } 00613 00614 /* now bring current pressure to the correct pressure - must do this 00615 * before initial values are saved in iter start/end */ 00616 ConvPresTempEdenIoniz(); 00617 00618 /* this counts number of times ConvBase is called by PressureChange, in 00619 * current zone these are reset here, so that we count from first zone 00620 * not search */ 00621 conv.nPres2Ioniz = 0; 00622 00623 dense.lgEdenBad = false; 00624 dense.nzEdenBad = 0; 00625 00626 /* save counter upon exit so we can see how many iterations were 00627 * needed to do true zones */ 00628 conv.nTotalIoniz_start = conv.nTotalIoniz; 00629 00630 /* these are variables to remember the biggest error in the 00631 * electron density, and the zone in which it occurred */ 00632 conv.BigEdenError = 0.; 00633 conv.AverEdenError = 0.; 00634 conv.BigHeatCoolError = 0.; 00635 conv.AverHeatCoolError = 0.; 00636 conv.BigPressError = 0.; 00637 conv.AverPressError = 0.; 00638 00639 /*>>chng 06 jun 09, only reset on first try - otherwise have stable solution? */ 00640 if( iteration == 1 ) 00641 { 00642 /* now remember some things we may need even in first zone, these 00643 * are normally set towards end of zone calculation in RT_tau_inc */ 00644 struc.testr[0] = (realnum)phycon.te; 00645 /* >>chng 02 May 2001 rjrw: add hden for dilution */ 00646 struc.hden[0] = dense.gas_phase[ipHYDROGEN]; 00647 /* pden is the total number of particles per unit vol */ 00648 struc.DenParticles[0] = dense.pden; 00649 struc.heatstr[0] = thermal.htot; 00650 struc.coolstr[0] = thermal.ctot; 00651 struc.volstr[0] = (realnum)radius.dVeff; 00652 struc.drad_x_fillfac[0] = (realnum)radius.drad_x_fillfac; 00653 struc.histr[0] = dense.xIonDense[ipHYDROGEN][0]; 00654 struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1]; 00655 struc.ednstr[0] = (realnum)dense.eden; 00656 struc.o3str[0] = dense.xIonDense[ipOXYGEN][2]; 00657 struc.DenMass[0] = dense.xMassDensity; 00658 struc.drad[0] = (realnum)radius.drad; 00659 } 00660 00661 /* check that nflux extends above IP of highest ionization species present. 00662 * for collisional case it is possible for species to exist that are higher 00663 * IP than the limit to the continuum. Need continuum to encompass all 00664 * possible emission - to account for diffuse emission 00665 * NB 00666 * on second iteration of multi-iteration model this may result in rfield.nflux increasing 00667 * which can change the solution */ 00668 nflux_old = rfield.nflux; 00669 nelem_reflux = -1; 00670 ion_reflux = -1; 00671 for( nelem=2; nelem < LIMELM; nelem++ ) 00672 { 00673 /* do not include hydrogenic species in following */ 00674 for( i=0; i < nelem; i++ ) 00675 { 00676 if( dense.xIonDense[nelem][i+1] > 0. ) 00677 { 00678 if( Heavy.ipHeavy[nelem][i] > rfield.nflux ) 00679 { 00680 rfield.nflux = Heavy.ipHeavy[nelem][i]; 00681 nelem_reflux = nelem; 00682 ion_reflux = i; 00683 } 00684 } 00685 } 00686 } 00687 00688 /* was the upper limit to the continuum updated? if so need to define 00689 * continuum variables to this new range */ 00690 if( nflux_old != rfield.nflux ) 00691 { 00692 /* zero out parts of rfield arrays that were previously undefined */ 00693 rfield_opac_zero( nflux_old-1 , rfield.nflux ); 00694 00695 /* if continuum reset up, we need to define gaunt factors through high end */ 00696 /*tfidle(false);*/ 00697 /* this calls tfidle, among other things */ 00698 /* this sets values of pressure.PresTotlCurr */ 00699 PresTotCurrent(); 00700 00701 /* redo ionization and update opacities */ 00702 if( ConvBase(1) ) 00703 { 00704 /* this is catastrophic failure */ 00705 lgAbort = true; 00706 return( 1 ); 00707 } 00708 00709 /* need to update continuum opacities */ 00710 if( trace.lgTrace ) 00711 { 00712 fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n", 00713 nflux_old , rfield.nflux , 00714 rfield.anu[nflux_old] , rfield.anu[rfield.nflux] ); 00715 fprintf(ioQQQ," caused by element %li ion %li \n", 00716 nelem_reflux ,ion_reflux ); 00717 } 00718 } 00719 00720 /* wind with negative velocity is flow from PDR - change density slightly 00721 * and call pressure solver to see if it returns to original density */ 00722 if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. ) 00723 { 00724 /* >>chng 04 apr 23, change pressure and let solver come back to correct 00725 * pressure. This trick makes sure we are correctly either super or sub sonic 00726 * since solver will jump from one branch to the other if necessary */ 00727 # define PCHNG 0.98 00728 /* this ignores return condition, assume must be ok if got this far */ 00729 if( ConvPresTempEdenIoniz() ) 00730 { 00731 /* this is an error return, calculation will immediately stop */ 00732 lgAbort = true; 00733 return(1); 00734 } 00735 00736 pressure.PresTotlInit *= PCHNG; 00737 if( ConvPresTempEdenIoniz() ) 00738 { 00739 /* this is an error return, calculation will immediately stop */ 00740 lgAbort = true; 00741 return(1); 00742 } 00743 00744 pressure.PresTotlInit /= PCHNG; 00745 if( ConvPresTempEdenIoniz() ) 00746 { 00747 /* this is an error return, calculation will immediately stop */ 00748 lgAbort = true; 00749 return(1); 00750 } 00751 00752 # undef PCHNG 00753 } 00754 /* this is the only valid return and lgAbort should be false*/ 00755 return( lgAbort ); 00756 }