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 /*PressureChange called by ConvPresTempEdenIoniz 00004 * evaluate the current pressure, change needed to get it to converge, 00005 * the global static variable pressure_change_factor 00006 * applies this correction factor to all gas constituents, 00007 * sets conv.lgConvPres true if good pressure, false if pressure change capped */ 00008 /*lgConvPres finds amount by which density should change to move towards pressure equilibrium 00009 * returns true if pressure is converged */ 00010 #include "cddefines.h" 00011 #include "abund.h" 00012 #include "hmi.h" 00013 #include "struc.h" 00014 #include "trace.h" 00015 #include "wind.h" 00016 #include "phycon.h" 00017 #include "thermal.h" 00018 #include "dense.h" 00019 #include "geometry.h" 00020 #include "radius.h" 00021 #include "mole.h" 00022 #include "dynamics.h" 00023 #include "pressure.h" 00024 #include "colden.h" 00025 #include "conv.h" 00026 00027 /* this is the pressure change pressure_change_factor, needed to move current pressure to correct pressure */ 00028 static double pressure_change_factor; 00029 00030 /*PressureChange evaluate the current pressure, and change needed to 00031 * get it to PresTotlInit, 00032 * return value is true is density was changed, false if no changes were necessary */ 00033 int PressureChange( 00034 /* this is change factor, 1 at first, becomes smaller as oscillations occur */ 00035 double dP_chng_factor ) 00036 { 00037 long int ion, 00038 nelem, 00039 lgChange, 00040 mol; 00041 00042 double abun, 00043 edensave, 00044 hold; 00045 00046 /* biggest multiplicative change in the pressure that we will allow */ 00047 /* allowed error in the pressure is conv.PressureErrorAllowed*/ 00048 double pdelta; 00049 00050 static double FacAbun, 00051 FacAbunSav, 00052 OldAbun; 00053 00054 DEBUG_ENTRY( "PressureChange()" ); 00055 00056 edensave = dense.eden; 00057 00058 /* first evaluate total pressure for this location, and current conditions 00059 * CurrentPressure is just sum of gas plus local line radiation pressure */ 00060 /* this sets values of pressure.PresTotlCurr */ 00061 PresTotCurrent(); 00062 00063 /* this is equal to zero upon first call in this iteration, 00064 * init some vars if this is first call */ 00065 if( !conv.nTotalIoniz ) 00066 { 00067 /* one time initialization of variables */ 00068 conv.hist_pres_npres = -1; 00069 conv.hist_pres_nzone = -1; 00070 conv.hist_pres_limit = 0; 00071 } 00072 00073 /* this will save the history of density - pressure relationship 00074 * for the current zone */ 00075 if( nzone!=conv.hist_pres_nzone ) 00076 { 00077 /* first time in this zone - reset counters */ 00078 conv.hist_pres_nzone = nzone; 00079 /* this counter will be updated after vars are saved so will be 00080 * total number of saved values */ 00081 conv.hist_pres_npres = 0; 00082 } 00083 00084 /* do we need to allocate, or reallocate, memory for the vectors */ 00085 if( conv.hist_pres_limit==0 ) 00086 { 00087 conv.hist_pres_limit = 200; 00088 conv.hist_pres_density = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) ); 00089 conv.hist_pres_current = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) ); 00090 conv.hist_pres_correct = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) ); 00091 } 00092 00093 /* ran out of space - allocate more */ 00094 if( (conv.hist_pres_npres+1) >= conv.hist_pres_limit ) 00095 { 00096 conv.hist_pres_limit *= 3; 00097 conv.hist_pres_density = (double *)REALLOC( conv.hist_pres_density, (unsigned)(conv.hist_pres_limit)*sizeof(double)); 00098 conv.hist_pres_current = (double *)REALLOC( conv.hist_pres_current, (unsigned)(conv.hist_pres_limit)*sizeof(double)); 00099 conv.hist_pres_correct = (double *)REALLOC( conv.hist_pres_correct, (unsigned)(conv.hist_pres_limit)*sizeof(double)); 00100 } 00101 00102 /* >>chng 04 feb 11, add option to remember current density and pressure */ 00103 conv.hist_pres_density[conv.hist_pres_npres] = dense.gas_phase[ipHYDROGEN]; 00104 conv.hist_pres_current[conv.hist_pres_npres] = pressure.PresTotlCurr; 00105 conv.hist_pres_correct[conv.hist_pres_npres] = pressure.PresTotlCorrect; 00106 ++conv.hist_pres_npres; 00107 00108 /* remember old hydrogen density */ 00109 hold = dense.gas_phase[ipHYDROGEN]; 00110 00111 /* this will be set true if density or abundances change in this zone */ 00112 lgChange = false; 00113 00114 /* this evaluates current pressure, sets pressure_change_factor and 00115 * pressure.PresTotlCorrect, updates velocity, 00116 * and returns true if pressure has converged 00117 * sets pressure.PresTotlCorrect */ 00118 conv.lgConvPres = lgConvPres(); 00119 /*fprintf(ioQQQ,"DEBUG pressure nH %.3e init %.3e correct %.3e currnt %.3e \n", 00120 dense.gas_phase[ipHYDROGEN] , 00121 pressure.PresTotlInit, 00122 pressure.PresTotlCorrect , 00123 pressure.PresTotlCurr);*/ 00124 00125 /* if convergence is OK at present state, so no change reqd, simply return 00126 * >>chng 05 feb 04, cannot do this test here since variable abundances, test has 00127 * not been done, so variable abundances did not work, 00128 * caught by Marcelo Castellanos and David Valls-Gabaud 00129 if( conv.lgConvPres ) 00130 return false;*/ 00131 00132 /* >> chng 02 dec 13 rjrw: short-circuit if nothing changes */ 00133 if( pressure_change_factor != 1. ) 00134 { 00135 lgChange = true; 00136 } 00137 00138 /* allow 3 percent changes, 00139 * dP_chng_factor is initially 1, becomes smaller if sign of pressure change, changes */ 00140 pdelta = 0.03 * dP_chng_factor; 00141 00142 /* make sure that change is not too extreme */ 00143 pressure_change_factor = MIN2(pressure_change_factor,1.+pdelta); 00144 pressure_change_factor = MAX2(pressure_change_factor,1.-pdelta); 00145 00146 { 00147 /*@-redef@*/ 00148 enum{DEBUG_LOC=false}; 00149 static long int nsave=-1; 00150 /*@+redef@*/ 00151 if( DEBUG_LOC /*&& nzone > 150 && iteration > 1*/ ) 00152 { 00153 if( nsave-nzone ) fprintf(ioQQQ,"\n"); 00154 nsave = nzone; 00155 fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\n", 00156 nzone, 00157 pressure_change_factor, 00158 /* when this is negative we need to raise the density */ 00159 (pressure.PresTotlCorrect-pressure.PresTotlCurr)*100./pressure.PresTotlCorrect, 00160 pressure.PresTotlCorrect, 00161 pressure.PresTotlCurr, 00162 pressure.PresGasCurr, 00163 pressure.pres_radiation_lines_curr ); 00164 } 00165 } 00166 00167 /* >>chng 97 jun 03, added variable abundances for element table command 00168 * option to get abundances off a table with element read command */ 00169 if( abund.lgAbTaON ) 00170 { 00171 lgChange = true; 00172 for( nelem=1; nelem < LIMELM; nelem++ ) 00173 { 00174 if( abund.lgAbunTabl[nelem] ) 00175 { 00176 abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)* 00177 dense.gas_phase[ipHYDROGEN]; 00178 00179 hold = abun/dense.gas_phase[nelem]; 00180 dense.gas_phase[nelem] = (realnum)abun; 00181 00182 for( ion=0; ion < (nelem + 2); ion++ ) 00183 { 00184 dense.xIonDense[nelem][ion] *= (realnum)hold; 00185 } 00186 } 00187 } 00188 } 00189 00190 /* this is set false if fluctuations abundances command entered, 00191 * when density variations are in place this is true, 00192 * and dense.chDenseLaw is "SINE" */ 00193 if( !dense.lgDenFlucOn ) 00194 { 00195 /* abundance variations are in place */ 00196 lgChange = true; 00197 if( nzone <= 1 ) 00198 { 00199 OldAbun = 1.; 00200 FacAbun = 1.; 00201 if( dense.lgDenFlucRadius ) 00202 { 00203 /* cycle over radius */ 00204 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+ 00205 dense.flcPhase) + dense.csecnd; 00206 } 00207 else 00208 { 00209 /* cycle over column density */ 00210 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+ 00211 dense.flcPhase) + dense.csecnd; 00212 } 00213 } 00214 else 00215 { 00216 OldAbun = FacAbunSav; 00217 /* rapid abundances fluctuation */ 00218 if( dense.lgDenFlucRadius ) 00219 { 00220 /* cycle over radius */ 00221 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+ 00222 dense.flcPhase) + dense.csecnd; 00223 } 00224 else 00225 { 00226 /* cycle over column density */ 00227 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+ 00228 dense.flcPhase) + dense.csecnd; 00229 } 00230 FacAbun = FacAbunSav/OldAbun; 00231 } 00232 } 00233 else 00234 { 00235 /* abundance variations are NOT in place */ 00236 FacAbun = 1.; 00237 } 00238 00239 /* chng 02 dec 11 rjrw -- remove test, saves marginal time could generate nasty intermittent bug when 00240 * ( pressure_change_factor*FacAbun == 1 ) && (pressure_change_factor != 1) */ 00241 if( lgChange ) 00242 { 00243 /* H, He not affected by abundance fluctuations, so only change these 00244 * by the pressure change factor */ 00245 for( nelem=ipHYDROGEN; nelem <= ipHELIUM; ++nelem ) 00246 { 00247 dense.xMolecules[nelem] *= (realnum)pressure_change_factor; 00248 dense.gas_phase[nelem] *= (realnum)pressure_change_factor; 00249 for( ion=0; ion < (nelem + 2); ion++ ) 00250 { 00251 /* >>chng 96 jul 12 had only multiplied total abun, not ions */ 00252 dense.xIonDense[nelem][ion] *= (realnum)pressure_change_factor; 00253 } 00254 } 00255 00256 /* Li on up is affect by both pressure and abundance variations, 00257 * so multiply by both factors */ 00258 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ ) 00259 { 00260 dense.xMolecules[nelem] *= (realnum)(pressure_change_factor*FacAbun); 00261 dense.gas_phase[nelem] *= (realnum)(pressure_change_factor*FacAbun); 00262 for( ion=0; ion < (nelem + 2); ion++ ) 00263 { 00264 /* >>chng 96 jul 12 had only multiplied total abun, not ions */ 00265 dense.xIonDense[nelem][ion] *= (realnum)(pressure_change_factor*FacAbun); 00266 } 00267 } 00268 00269 dense.eden *= pressure_change_factor; 00270 00271 /* must call TempChange since ionization has changed, there are some 00272 * terms that affect collision rates (H0 term in electron collision) */ 00273 TempChange(phycon.te , false); 00274 00275 /* molecules done in hmole, only change pressure, not abundances */ 00276 for(mol = 0; mol < N_H_MOLEC; mol++) 00277 { 00278 hmi.Hmolec[mol] *= (realnum)pressure_change_factor; 00279 } 00280 hmi.H2_total *= (realnum)pressure_change_factor; 00281 00282 /* molecules done in comole */ 00283 /* >>chng 03 sep 15, upper limit had not included the C, O atoms/ions */ 00284 /*for( ion=0; ion < NUM_HEAVY_MOLEC; ion++ )*/ 00285 for( ion=0; ion < mole.num_comole_calc; ion++ ) 00286 { 00287 COmole[ion]->hevmol *= (realnum)(pressure_change_factor*FacAbun); 00288 00289 /* check for NaN */ 00290 ASSERT( !isnan( COmole[ion]->hevmol ) ); 00291 } 00292 } 00293 00294 if( trace.lgTrace ) 00295 { 00296 fprintf( ioQQQ, 00297 " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n", 00298 hold, dense.gas_phase[ipHYDROGEN], geometry.FillFac ); 00299 00300 if( trace.lgNeBug ) 00301 { 00302 fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n", 00303 edensave, dense.eden, edensave/dense.eden ); 00304 } 00305 } 00306 00307 { 00308 /*@-redef@*/ 00309 enum {DEBUG_LOC=false}; 00310 /*@+redef@*/ 00311 if( DEBUG_LOC && (nzone>215) ) 00312 { 00313 fprintf( ioQQQ, 00314 "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n", 00315 radius.depth, 00316 pressure.PresTotlCurr, 00317 pressure.PresTotlInit + pressure.PresInteg, 00318 pressure.PresTotlInit, 00319 pressure.PresGasCurr, 00320 pressure.PresRamCurr, 00321 pressure.pres_radiation_lines_curr, 00322 /* subtract continuum rad pres which has already been added on */ 00323 pressure.PresInteg - pressure.pinzon, 00324 wind.windv/1e5, 00325 sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5, 00326 TorF(conv.lgConvPres) ); 00327 } 00328 } 00329 00330 return lgChange; 00331 } 00332 00333 /*lgConvPres finds amount by which density should change to move towards pressure equilibrium 00334 * sets pressure.PresTotlCorrect 00335 * returns true if pressure is converged */ 00336 bool lgConvPres(void) 00337 { 00338 double dnew, 00339 term; 00340 bool lgRet; 00341 00342 DEBUG_ENTRY( "lgConvPres()" ); 00343 00344 /* make sure this is set by one of the following branches - set to zero here, 00345 * then assert that it is greater than zero at end */ 00346 pressure.PresTotlCorrect = 0.; 00347 00348 /* evaluate a series of possible pressure options, and set the file static variable 00349 * pressure_change_factor */ 00350 /* inside out globule */ 00351 if( strcmp(dense.chDenseLaw,"GLOB") == 0 ) 00352 { 00353 /* GLBDST is distance from globule, or glbrad-DEPTH */ 00354 if( radius.glbdst < 0. ) 00355 { 00356 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" ); 00357 fprintf( ioQQQ, " This is routine lgConvPres, GLBDST is%10.2e\n", 00358 radius.glbdst ); 00359 cdEXIT(EXIT_FAILURE); 00360 } 00361 pressure_change_factor = (radius.glbden*pow(radius.glbrad/(radius.glbdst),radius.glbpow))/ 00362 dense.gas_phase[ipHYDROGEN]; 00363 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00364 } 00365 00366 /* this is impossible - wind with identically zero velocity */ 00367 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && dynamics.lgStatic ) 00368 { 00369 /* this is default case, keep initial H den constant at first zone, 00370 * from iteration to iteration */ 00371 if( dense.lgDenseInitConstant ) 00372 { 00373 pressure_change_factor = 1.; 00374 } 00375 else 00376 { 00377 /* this option is if constant pressure set and flag also 00378 * set */ 00379 pressure.PresTotlCorrect = pressure.PresTotlInit; 00380 /* ratio of correct to current pressures */ 00381 pressure_change_factor = pressure.PresTotlCorrect / pressure.PresTotlCurr; 00382 } 00383 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00384 } 00385 00386 /* this is positive wind velocity the outflowing wind beyond sonic point */ 00387 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv > 0. ) 00388 { 00389 00390 /* this is logic for supersonic outflowing wind solution, 00391 * which assumes positive velocity, well above sonic point. 00392 * following makes sure wind v only updated once per zone */ 00393 if( nzone > 1 ) 00394 { 00395 /* Wind model */ 00396 00397 if( trace.lgTrace && trace.lgWind ) 00398 { 00399 fprintf(ioQQQ," lgConvPres sets AccelGravity %.3e lgDisk?%c\n", 00400 wind.AccelGravity , 00401 TorF(wind.lgDisk) ); 00402 } 00403 00404 # if 0 00405 /* acceleration due to grad P(rad), xMassDensity is gm per cc */ 00406 wind.AccelPres = (realnum)(-(pressure.pres_radiation_lines_curr - pold)/radius.drad/ 00407 dense.xMassDensity); 00408 /* is this the first zone? */ 00409 if( nzone > 2 ) 00410 { 00411 /* acceleration due to grad P(rad), xMassDensity is gm per cc */ 00412 wind.AccelPres = (realnum)(-(pressure.pres_radiation_lines_curr - pold)/radius.drad/ 00413 dense.xMassDensity); 00414 } 00415 else 00416 { 00417 wind.AccelPres = 0.; 00418 } 00419 # endif 00420 00421 # if 0 00422 /* this is evaluated in pressure_total*/ 00423 /* this is numerically unstable */ 00424 wind.AccelPres = 0.; 00425 00426 /* AccelXXX is computed in radinc, is continuum and line acceleration 00427 * AccelCont is continuum acceleration 00428 * AccelLine is line acceleration 00429 * AccelPres is gradient of local gas pressure */ 00430 wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres; 00431 00432 if( nzone==NzoneEval+1) 00433 { 00434 double da = AccelNet - AccelNetLast; 00435 /* we have previous acceleration in this iteration, use it 00436 * a = a_i + (a_i - a_i-1)/(r_i - r_i-1) * (r_i+1 - r_i)/2 */ 00437 AccelNet += da* (radius.drNext/radius.drad) /2.; 00438 00439 } 00440 # endif 00441 00442 /* following is form of energy equation 00443 * struc.windv[nzone-2] is velocity of previous zone 00444 * this increments that velocity to form square of new wind 00445 * velocity for outer edge of this zone */ 00446 term = POW2(struc.windv[nzone-2]) + 2.*(wind.AccelTot - wind.AccelGravity)* radius.drad; 00447 00448 /* increment velocity if it is substantially positive */ 00449 if( term <= 1e3 ) 00450 { 00451 /* wind velocity is well below sonic point, give up, 00452 * do not change velocity */ 00453 wind.lgVelPos = false; 00454 } 00455 else 00456 { 00457 /* struc.windv[nzone-2] is velocity of previous zone 00458 * this increments that velocity to form square of new wind 00459 * velocity for outer edge of this zone 00460 double windnw = (double)POW2(struc.windv[nzone-2]) + 00461 (double)(2.*(wind.AccelTot-wind.AccelGravity))*radius.drad;*/ 00462 00463 /* wind.windv is velocity at OUTER edge of this zone */ 00464 wind.windv = (realnum)sqrt(term); 00465 wind.lgVelPos = true; 00466 } 00467 00468 /* following needed for expansion cooling */ 00469 dynamics.dDensityDT = (realnum)(wind.AccelTot/wind.windv + 2.*wind.windv/ 00470 radius.Radius); 00471 00472 if( trace.lgTrace && trace.lgWind ) 00473 { 00474 fprintf(ioQQQ," lgConvPres new wind V zn%li %.3e AccelTot %.3e AccelGravity %.3e\n", 00475 nzone,wind.windv, wind.AccelTot, wind.AccelGravity ); 00476 } 00477 } 00478 else 00479 { 00480 pressure_change_factor = 1.; 00481 } 00482 00483 /* conservation of mass sets density here */ 00484 pressure_change_factor = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq; 00485 pressure.PresTotlCorrect = pressure.PresTotlCurr * pressure_change_factor; 00486 } 00487 00488 /* this is negative wind velocity the new dynamics */ 00489 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. ) 00490 { 00491 /* sets pressure.PresTotlCorrect */ 00492 pressure_change_factor = DynaPresChngFactor(); 00493 } 00494 00495 /* this is impossible - wind with identically zero velocity */ 00496 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv == 0. ) 00497 { 00498 /* declare insanity */ 00499 fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" ); 00500 /* TotalInsanity announces fatal problem, ShowMe, then cdEXIT with failure */ 00501 TotalInsanity(); 00502 } 00503 00504 else if( strcmp(dense.chDenseLaw,"SINE") == 0 ) 00505 { 00506 /* rapid density fluctuation */ 00507 if( dense.lgDenFlucRadius ) 00508 { 00509 pressure_change_factor = (dense.cfirst*cos(radius.depth*dense.flong+dense.flcPhase) + 00510 dense.csecnd)/dense.gas_phase[ipHYDROGEN]; 00511 } 00512 else 00513 { 00514 pressure_change_factor = (dense.cfirst*cos(colden.colden[ipCOL_HTOT]*dense.flong+dense.flcPhase) + 00515 dense.csecnd)/dense.gas_phase[ipHYDROGEN]; 00516 } 00517 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00518 } 00519 00520 else if( strcmp(dense.chDenseLaw,"POWR") == 0 ) 00521 { 00522 /* power law function of radius */ 00523 dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower); 00524 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN]; 00525 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00526 } 00527 00528 else if( strcmp(dense.chDenseLaw,"POWD") == 0 ) 00529 { 00530 /* power law function of depth */ 00531 dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower); 00532 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN]; 00533 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00534 } 00535 00536 else if( strcmp(dense.chDenseLaw,"POWC") == 0 ) 00537 { 00538 /* power law function of column density */ 00539 dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/ 00540 dense.rscale,dense.DensityPower); 00541 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN]; 00542 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00543 } 00544 00545 else if( strcmp(dense.chDenseLaw,"CPRE") == 0 ) 00546 { 00547 /* constant pressure */ 00548 if( pressure.lgContRadPresOn ) 00549 { 00550 /* >>chng 01 oct 31, replace pneed with CorretPressure */ 00551 /* this has pressure due to incident continuum */ 00552 pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg; 00553 } 00554 else 00555 { 00556 /* this does not have pressure due to incident continuum*/ 00557 pressure.PresTotlCorrect = pressure.PresTotlInit* 00558 /* following term normally unity, power law set with option par on cmmnd*/ 00559 pow(radius.Radius/radius.rinner,(double)pressure.PresPowerlaw); 00560 } 00561 00562 /* ratio of correct to current pressures */ 00563 pressure_change_factor = pressure.PresTotlCorrect / pressure.PresTotlCurr; 00564 } 00565 00566 else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 ) 00567 { 00568 if( strcmp(dense.chDenseLaw,"DLW1") == 0 ) 00569 { 00570 /* call ACF sub */ 00571 pressure_change_factor = dense_fabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN]; 00572 } 00573 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 ) 00574 { 00575 /* call table interpolation subroutine 00576 * >>chng 96 nov 29, added dense_tabden */ 00577 pressure_change_factor = dense_tabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN]; 00578 } 00579 else 00580 { 00581 fprintf( ioQQQ, " Insanity, lgConvPres gets chCPres=%4.4s\n", 00582 dense.chDenseLaw ); 00583 cdEXIT(EXIT_FAILURE); 00584 } 00585 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00586 /* >>chng 06 feb 21, from above to below. we want to look at change, 00587 * not the value. Bug found by Valentina Luridiana */ 00588 /*if( pressure.PresTotlCorrect > 3. || pressure.PresTotlCorrect< 1./3 )*/ 00589 if( pressure_change_factor > 3. || pressure_change_factor< 1./3 ) 00590 { 00591 static bool lgWARN2BIG=false; 00592 if( !lgWARN2BIG ) 00593 { 00594 lgWARN2BIG = true; 00595 fprintf(ioQQQ,"\n\n >========== Warning! The tabulated or functional change in density as a function of depth was VERY large. This is zone %li.\n",nzone); 00596 fprintf(ioQQQ," >========== Warning! This will cause convergence problems. \n"); 00597 fprintf(ioQQQ," >========== Warning! The current radius is %.3e. \n",radius.Radius); 00598 fprintf(ioQQQ," >========== Warning! The current depth is %.3e. \n",radius.depth); 00599 fprintf(ioQQQ," >========== Warning! Consider using a more modest change in density vs radius. \n\n\n"); 00600 } 00601 } 00602 } 00603 00604 else if( strcmp(dense.chDenseLaw,"CDEN") == 0 ) 00605 { 00606 /* this is the default, constant density */ 00607 pressure_change_factor = 1.; 00608 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor; 00609 } 00610 00611 else 00612 { 00613 fprintf( ioQQQ, " Unknown pressure law=%s= This is a major internal error.\n", 00614 dense.chDenseLaw ); 00615 ShowMe(); 00616 cdEXIT(EXIT_FAILURE); 00617 } 00618 00619 /* one of the branches above must have reset this variable, 00620 * and it was init to 0 at start. Confirm that non-zero */ 00621 ASSERT( pressure.PresTotlCorrect > FLT_MIN ); 00622 00623 /* now see whether current pressure is within error bounds */ 00624 if( pressure_change_factor > 1. + conv.PressureErrorAllowed || pressure_change_factor < 1. - conv.PressureErrorAllowed ) 00625 { 00626 lgRet = false; 00627 conv.lgConvPres = false; 00628 } 00629 else 00630 { 00631 lgRet = true; 00632 conv.lgConvPres = true; 00633 } 00634 00635 return lgRet; 00636 }