00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "physconst.h"
00012 #include "abund.h"
00013 #include "hmi.h"
00014 #include "trace.h"
00015 #include "wind.h"
00016 #include "thermal.h"
00017 #include "dense.h"
00018 #include "geometry.h"
00019 #include "radius.h"
00020 #include "mole.h"
00021 #include "dynamics.h"
00022 #include "pressure.h"
00023 #include "colden.h"
00024 #include "rt.h"
00025 #include "conv.h"
00026
00027
00028 static double pressure_change_factor;
00029
00030
00031
00032
00033 int PressureChange(
00034
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
00047
00048 double pdelta;
00049
00050 static double FacAbun,
00051 FacAbunSav,
00052 OldAbun;
00053
00054 DEBUG_ENTRY( "PressureChange()" );
00055
00056 edensave = dense.eden;
00057
00058
00059
00060
00061 PresTotCurrent();
00062
00063 # if 0
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 if( nzone <= 1 && (iteration==1 || dense.lgDenseInitConstant) )
00074 {
00075
00076
00077
00078 pressure.PresTotlInit = pressure.PresTotlCurr;
00079 pressure.PresRamInit = pressure.PresRamCurr;
00080 if( trace.lgTrace )
00081 {
00082 fprintf( ioQQQ,
00083 " PressureChange called, this is first zone, so reseting pressure to %e\n",
00084 pressure.PresTotlInit );
00085 }
00086 }
00087 # endif
00088
00089
00090
00091 if( !conv.nTotalIoniz )
00092 {
00093
00094 conv.hist_pres_npres = -1;
00095 conv.hist_pres_nzone = -1;
00096 conv.hist_pres_limit = 0;
00097 }
00098
00099
00100
00101 if( nzone!=conv.hist_pres_nzone )
00102 {
00103
00104 conv.hist_pres_nzone = nzone;
00105
00106
00107 conv.hist_pres_npres = 0;
00108 }
00109
00110
00111 if( conv.hist_pres_limit==0 )
00112 {
00113 conv.hist_pres_limit = 200;
00114 conv.hist_pres_density = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
00115 conv.hist_pres_current = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
00116 conv.hist_pres_correct = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
00117 }
00118
00119
00120 if( (conv.hist_pres_npres+1) >= conv.hist_pres_limit )
00121 {
00122 conv.hist_pres_limit *= 3;
00123 conv.hist_pres_density = (double *)REALLOC( conv.hist_pres_density, (unsigned)(conv.hist_pres_limit)*sizeof(double));
00124 conv.hist_pres_current = (double *)REALLOC( conv.hist_pres_current, (unsigned)(conv.hist_pres_limit)*sizeof(double));
00125 conv.hist_pres_correct = (double *)REALLOC( conv.hist_pres_correct, (unsigned)(conv.hist_pres_limit)*sizeof(double));
00126 }
00127
00128
00129 conv.hist_pres_density[conv.hist_pres_npres] = dense.gas_phase[ipHYDROGEN];
00130 conv.hist_pres_current[conv.hist_pres_npres] = pressure.PresTotlCurr;
00131 conv.hist_pres_correct[conv.hist_pres_npres] = pressure.PresTotlCorrect;
00132 ++conv.hist_pres_npres;
00133
00134
00135 hold = dense.gas_phase[ipHYDROGEN];
00136
00137
00138 lgChange = false;
00139
00140
00141
00142
00143
00144 conv.lgConvPres = lgConvPres();
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 if( pressure_change_factor != 1. )
00160 {
00161 lgChange = true;
00162 }
00163
00164
00165
00166 pdelta = 0.03 * dP_chng_factor;
00167
00168
00169 pressure_change_factor = MIN2(pressure_change_factor,1.+pdelta);
00170 pressure_change_factor = MAX2(pressure_change_factor,1.-pdelta);
00171
00172 {
00173
00174 enum{DEBUG_LOC=false};
00175 static long int nsave=-1;
00176
00177 if( DEBUG_LOC )
00178 {
00179 if( nsave-nzone ) fprintf(ioQQQ,"\n");
00180 nsave = nzone;
00181 fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\n",
00182 nzone,
00183 pressure_change_factor,
00184
00185 (pressure.PresTotlCorrect-pressure.PresTotlCurr)*100./pressure.PresTotlCorrect,
00186 pressure.PresTotlCorrect,
00187 pressure.PresTotlCurr,
00188 pressure.PresGasCurr,
00189 pressure.pres_radiation_lines_curr );
00190 }
00191 }
00192
00193
00194
00195 if( abund.lgAbTaON )
00196 {
00197 lgChange = true;
00198 for( nelem=1; nelem < LIMELM; nelem++ )
00199 {
00200 if( abund.lgAbunTabl[nelem] )
00201 {
00202 abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
00203 dense.gas_phase[ipHYDROGEN];
00204
00205 hold = abun/dense.gas_phase[nelem];
00206 dense.gas_phase[nelem] = (float)abun;
00207
00208 for( ion=0; ion < (nelem + 2); ion++ )
00209 {
00210 dense.xIonDense[nelem][ion] *= (float)hold;
00211 }
00212 }
00213 }
00214 }
00215
00216
00217
00218
00219 if( !dense.lgDenFlucOn )
00220 {
00221
00222 lgChange = true;
00223 if( nzone <= 1 )
00224 {
00225 OldAbun = 1.;
00226 FacAbun = 1.;
00227 if( dense.lgDenFlucRadius )
00228 {
00229
00230 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
00231 dense.flcPhase) + dense.csecnd;
00232 }
00233 else
00234 {
00235
00236 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
00237 dense.flcPhase) + dense.csecnd;
00238 }
00239 }
00240 else
00241 {
00242 OldAbun = FacAbunSav;
00243
00244 if( dense.lgDenFlucRadius )
00245 {
00246
00247 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
00248 dense.flcPhase) + dense.csecnd;
00249 }
00250 else
00251 {
00252
00253 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
00254 dense.flcPhase) + dense.csecnd;
00255 }
00256 FacAbun = FacAbunSav/OldAbun;
00257 }
00258 }
00259 else
00260 {
00261
00262 FacAbun = 1.;
00263 }
00264
00265
00266
00267 if( lgChange )
00268 {
00269
00270
00271 for( nelem=ipHYDROGEN; nelem <= ipHELIUM; ++nelem )
00272 {
00273 dense.xMolecules[nelem] *= (float)pressure_change_factor;
00274 dense.gas_phase[nelem] *= (float)pressure_change_factor;
00275 for( ion=0; ion < (nelem + 2); ion++ )
00276 {
00277
00278 dense.xIonDense[nelem][ion] *= (float)pressure_change_factor;
00279 }
00280 }
00281
00282
00283
00284 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
00285 {
00286 dense.xMolecules[nelem] *= (float)(pressure_change_factor*FacAbun);
00287 dense.gas_phase[nelem] *= (float)(pressure_change_factor*FacAbun);
00288 for( ion=0; ion < (nelem + 2); ion++ )
00289 {
00290
00291 dense.xIonDense[nelem][ion] *= (float)(pressure_change_factor*FacAbun);
00292 }
00293 }
00294
00295 dense.eden *= pressure_change_factor;
00296
00297
00298 tfidle(false);
00299
00300
00301 for(mol = 0; mol < N_H_MOLEC; mol++)
00302 {
00303 hmi.Hmolec[mol] *= (float)pressure_change_factor;
00304 }
00305 hmi.H2_total *= (float)pressure_change_factor;
00306
00307
00308
00309
00310 for( ion=0; ion < mole.num_comole_calc; ion++ )
00311 {
00312 COmole[ion]->hevmol *= (float)(pressure_change_factor*FacAbun);
00313
00314
00315 ASSERT( !isnan( COmole[ion]->hevmol ) );
00316 }
00317 }
00318
00319 if( trace.lgTrace )
00320 {
00321 fprintf( ioQQQ,
00322 " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
00323 hold, dense.gas_phase[ipHYDROGEN], geometry.FillFac );
00324
00325 if( trace.lgNeBug )
00326 {
00327 fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n",
00328 edensave, dense.eden, edensave/dense.eden );
00329 }
00330 }
00331
00332 {
00333
00334 enum {DEBUG_LOC=false};
00335
00336 if( DEBUG_LOC && (nzone>215) )
00337 {
00338 fprintf( ioQQQ,
00339 "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n",
00340 radius.depth,
00341 pressure.PresTotlCurr,
00342 pressure.PresTotlInit + pressure.PresInteg,
00343 pressure.PresTotlInit,
00344 pressure.PresGasCurr,
00345 pressure.PresRamCurr,
00346 pressure.pres_radiation_lines_curr,
00347
00348 pressure.PresInteg - pressure.pinzon,
00349 wind.windv/1e5,
00350 sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5,
00351 TorF(conv.lgConvPres) );
00352 }
00353 }
00354
00355 DEBUG_EXIT( "PressureChange()" );
00356
00357 return lgChange;
00358 }
00359
00360
00361
00362
00363 bool lgConvPres(void)
00364 {
00365
00366 static double pold = 0.;
00367
00368 double windnw,
00369 dnew,
00370 term;
00371 bool lgRet;
00372
00373
00374
00375 pressure.PresTotlCorrect = 0.;
00376
00377
00378
00379
00380 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00381 {
00382
00383 if( radius.glbdst < 0. )
00384 {
00385 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
00386 fprintf( ioQQQ, " This is routine PressureChange, GLBDST is%10.2e\n",
00387 radius.glbdst );
00388 puts( "[Stop in PressureChange]" );
00389 cdEXIT(EXIT_FAILURE);
00390 }
00391 pressure_change_factor = (radius.glbden*pow(radius.glbrad/(radius.glbdst),radius.glbpow))/
00392 dense.gas_phase[ipHYDROGEN];
00393 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00394 }
00395
00396
00397 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && dynamics.lgStatic )
00398 {
00399
00400
00401 if( dense.lgDenseInitConstant )
00402 {
00403 pressure_change_factor = 1.;
00404 }
00405 else
00406 {
00407
00408
00409 pressure.PresTotlCorrect = pressure.PresTotlInit;
00410
00411 pressure_change_factor = pressure.PresTotlCorrect / pressure.PresTotlCurr;
00412 }
00413 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00414 }
00415
00416
00417 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv > 0. )
00418 {
00419
00420 static double rsave = 0.;
00421
00422
00423
00424
00425 if( radius.depth != rsave )
00426
00427 {
00428 rsave = radius.depth;
00429
00430 RT_radiative_acceleration();
00431
00432
00433 wind.agrav = (float)((6.67e-8*wind.comass)*(SOLAR_MASS/radius.Radius)/
00434 radius.Radius);
00435
00436
00437 if( nzone > 2 )
00438 {
00439
00440 wind.AccelPres = (float)(-(pressure.pres_radiation_lines_curr - pold)/radius.drad/
00441 dense.xMassDensity);
00442 }
00443 else
00444 {
00445 wind.AccelPres = 0.;
00446 }
00447
00448
00449 wind.AccelPres = 0.;
00450
00451
00452
00453
00454
00455 wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres;
00456
00457
00458 wind.AccelMax = (float)MAX2(wind.AccelMax,wind.AccelTot);
00459
00460
00461 wind.AccelAver += wind.AccelTot*(float)radius.drad_x_fillfac;
00462 wind.acldr += (float)radius.drad_x_fillfac;
00463
00464
00465 term = POW2(wind.windv) + 2.*(wind.AccelTot - wind.agrav)*
00466 radius.drad;
00467
00468
00469
00470 if( term <= 1e3 )
00471 {
00472
00473 wind.lgVelPos = false;
00474 }
00475 else
00476 {
00477
00478 windnw = (double)(wind.windv*wind.windv) + (double)(2.*
00479 (wind.AccelTot-wind.agrav))*radius.drad;
00480
00481
00482 wind.windv = (float)sqrt(windnw);
00483 wind.lgVelPos = true;
00484 }
00485
00486
00487 dynamics.dDensityDT = (float)(wind.AccelTot/wind.windv + 2.*wind.windv/
00488 radius.Radius);
00489
00490
00491
00492
00493 }
00494 else
00495 {
00496 pressure_change_factor = 1.;
00497 }
00498
00499 pressure_change_factor = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq;
00500 pold = pressure.pres_radiation_lines_curr;
00501 pressure.PresTotlCorrect = pressure.PresTotlCurr * pressure_change_factor;
00502 }
00503
00504
00505 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
00506 {
00507
00508 pressure_change_factor = DynaPresChngFactor();
00509 }
00510
00511
00512 else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv == 0. )
00513 {
00514
00515 fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
00516
00517 TotalInsanity();
00518 }
00519
00520 else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
00521 {
00522
00523 if( dense.lgDenFlucRadius )
00524 {
00525 pressure_change_factor = (dense.cfirst*cos(radius.depth*dense.flong+dense.flcPhase) +
00526 dense.csecnd)/dense.gas_phase[ipHYDROGEN];
00527 }
00528 else
00529 {
00530 pressure_change_factor = (dense.cfirst*cos(colden.colden[ipCOL_HTOT]*dense.flong+dense.flcPhase) +
00531 dense.csecnd)/dense.gas_phase[ipHYDROGEN];
00532 }
00533 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00534 }
00535
00536 else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
00537 {
00538
00539 dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower);
00540 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00541 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00542 }
00543
00544 else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
00545 {
00546
00547 dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower);
00548 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00549 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00550 }
00551
00552 else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
00553 {
00554
00555 dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/
00556 dense.rscale,dense.DensityPower);
00557 pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
00558 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00559 }
00560
00561 else if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
00562 {
00563
00564 if( pressure.lgContRadPresOn )
00565 {
00566
00567
00568 pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg;
00569 }
00570 else
00571 {
00572
00573 pressure.PresTotlCorrect = pressure.PresTotlInit*
00574
00575 pow(radius.Radius/radius.rinner,(double)pressure.PresPowerlaw);
00576 }
00577
00578
00579 pressure_change_factor = pressure.PresTotlCorrect / pressure.PresTotlCurr;
00580 }
00581
00582 else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 )
00583 {
00584 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00585 {
00586
00587 pressure_change_factor = dense_fabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN];
00588 }
00589 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00590 {
00591
00592
00593 pressure_change_factor = dense_tabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN];
00594 }
00595 else
00596 {
00597 fprintf( ioQQQ, " Insanity, PressureChange gets chCPres=%4.4s\n",
00598 dense.chDenseLaw );
00599 puts( "[Stop in PressureChange]" );
00600 cdEXIT(EXIT_FAILURE);
00601 }
00602 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00603
00604
00605
00606 if( pressure_change_factor > 3. || pressure_change_factor< 1./3 )
00607 {
00608 static bool lgWARN2BIG=false;
00609 if( !lgWARN2BIG )
00610 {
00611 lgWARN2BIG = true;
00612 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);
00613 fprintf(ioQQQ," >========== Warning! This will cause convergence problems. \n");
00614 fprintf(ioQQQ," >========== Warning! The current radius is %.3e. \n",radius.Radius);
00615 fprintf(ioQQQ," >========== Warning! The current depth is %.3e. \n",radius.depth);
00616 fprintf(ioQQQ," >========== Warning! Consider using a more modest change in density vs radius. \n\n\n");
00617 }
00618 }
00619 }
00620
00621 else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
00622 {
00623
00624 pressure_change_factor = 1.;
00625 pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
00626 }
00627
00628 else
00629 {
00630 fprintf( ioQQQ, " Unknown pressure law=%s= This is a major internal error.\n",
00631 dense.chDenseLaw );
00632 puts( "[Stop in PressureChange]" );
00633 ShowMe();
00634 cdEXIT(EXIT_FAILURE);
00635 }
00636
00637
00638
00639 ASSERT( pressure.PresTotlCorrect > FLT_MIN );
00640
00641
00642 if( pressure_change_factor > 1. + conv.PressureErrorAllowed || pressure_change_factor < 1. - conv.PressureErrorAllowed )
00643 {
00644 lgRet = false;
00645 conv.lgConvPres = false;
00646 }
00647 else
00648 {
00649 lgRet = true;
00650 conv.lgConvPres = true;
00651 }
00652
00653 return lgRet;
00654 }