cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pressure_change.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*PressureChange called by ConvPresTempEdenIoniz
4  * evaluate the current pressure, change needed to get it to converge,
5  * the global static variable pressure_change_factor
6  * applies this correction factor to all gas constituents,
7  * sets conv.lgConvPres true if good pressure, false if pressure change capped */
8 /*lgConvPres finds amount by which density should change to move towards pressure equilibrium
9  * returns true if pressure is converged */
10 #include "cddefines.h"
11 #include "abund.h"
12 #include "hmi.h"
13 #include "struc.h"
14 #include "trace.h"
15 #include "wind.h"
16 #include "phycon.h"
17 #include "thermal.h"
18 #include "dense.h"
19 #include "geometry.h"
20 #include "radius.h"
21 #include "mole.h"
22 #include "dynamics.h"
23 #include "pressure.h"
24 #include "colden.h"
25 #include "conv.h"
26 
27 /* this is the pressure change pressure_change_factor, needed to move current pressure to correct pressure */
28 static double pressure_change_factor;
29 
30 /*PressureChange evaluate the current pressure, and change needed to
31  * get it to PresTotlInit,
32  * return value is true is density was changed, false if no changes were necessary */
34  /* this is change factor, 1 at first, becomes smaller as oscillations occur */
35  double dP_chng_factor )
36 {
37  long int ion,
38  nelem,
39  lgChange,
40  mol;
41 
42  double abun,
43  edensave,
44  hold;
45 
46  /* biggest multiplicative change in the pressure that we will allow */
47  /* allowed error in the pressure is conv.PressureErrorAllowed*/
48  double pdelta;
49 
50  static double FacAbun,
51  FacAbunSav,
52  OldAbun;
53 
54  DEBUG_ENTRY( "PressureChange()" );
55 
56  edensave = dense.eden;
57 
58  /* first evaluate total pressure for this location, and current conditions
59  * CurrentPressure is just sum of gas plus local line radiation pressure */
60  /* this sets values of pressure.PresTotlCurr */
62 
63  /* this is equal to zero upon first call in this iteration,
64  * init some vars if this is first call */
65  if( !conv.nTotalIoniz )
66  {
67  /* one time initialization of variables */
68  conv.hist_pres_npres = -1;
69  conv.hist_pres_nzone = -1;
71  }
72 
73  /* this will save the history of density - pressure relationship
74  * for the current zone */
76  {
77  /* first time in this zone - reset counters */
79  /* this counter will be updated after vars are saved so will be
80  * total number of saved values */
82  }
83 
84  /* do we need to allocate, or reallocate, memory for the vectors */
85  if( conv.hist_pres_limit==0 )
86  {
87  conv.hist_pres_limit = 200;
88  conv.hist_pres_density = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
89  conv.hist_pres_current = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
90  conv.hist_pres_correct = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
91  }
92 
93  /* ran out of space - allocate more */
95  {
96  conv.hist_pres_limit *= 3;
97  conv.hist_pres_density = (double *)REALLOC( conv.hist_pres_density, (unsigned)(conv.hist_pres_limit)*sizeof(double));
98  conv.hist_pres_current = (double *)REALLOC( conv.hist_pres_current, (unsigned)(conv.hist_pres_limit)*sizeof(double));
99  conv.hist_pres_correct = (double *)REALLOC( conv.hist_pres_correct, (unsigned)(conv.hist_pres_limit)*sizeof(double));
100  }
101 
102  /* >>chng 04 feb 11, add option to remember current density and pressure */
107 
108  /* remember old hydrogen density */
109  hold = dense.gas_phase[ipHYDROGEN];
110 
111  /* this will be set true if density or abundances change in this zone */
112  lgChange = false;
113 
114  /* this evaluates current pressure, sets pressure_change_factor and
115  * pressure.PresTotlCorrect, updates velocity,
116  * and returns true if pressure has converged
117  * sets pressure.PresTotlCorrect */
119  /*fprintf(ioQQQ,"DEBUG pressure nH %.3e init %.3e correct %.3e currnt %.3e \n",
120  dense.gas_phase[ipHYDROGEN] ,
121  pressure.PresTotlInit,
122  pressure.PresTotlCorrect ,
123  pressure.PresTotlCurr);*/
124 
125  /* if convergence is OK at present state, so no change reqd, simply return
126  * >>chng 05 feb 04, cannot do this test here since variable abundances, test has
127  * not been done, so variable abundances did not work,
128  * caught by Marcelo Castellanos and David Valls-Gabaud
129  if( conv.lgConvPres )
130  return false;*/
131 
132  /* >> chng 02 dec 13 rjrw: short-circuit if nothing changes */
133  if( pressure_change_factor != 1. )
134  {
135  lgChange = true;
136  }
137 
138  /* allow 3 percent changes,
139  * dP_chng_factor is initially 1, becomes smaller if sign of pressure change, changes */
140  pdelta = 0.03 * dP_chng_factor;
141 
142  /* make sure that change is not too extreme */
145 
146  {
147  /*@-redef@*/
148  enum{DEBUG_LOC=false};
149  static long int nsave=-1;
150  /*@+redef@*/
151  if( DEBUG_LOC /*&& nzone > 150 && iteration > 1*/ )
152  {
153  if( nsave-nzone ) fprintf(ioQQQ,"\n");
154  nsave = nzone;
155  fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\n",
156  nzone,
158  /* when this is negative we need to raise the density */
164  }
165  }
166 
167  /* >>chng 97 jun 03, added variable abundances for element table command
168  * option to get abundances off a table with element read command */
169  if( abund.lgAbTaON )
170  {
171  lgChange = true;
172  for( nelem=1; nelem < LIMELM; nelem++ )
173  {
174  if( abund.lgAbunTabl[nelem] )
175  {
176  abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
178 
179  hold = abun/dense.gas_phase[nelem];
180  dense.gas_phase[nelem] = (realnum)abun;
181 
182  for( ion=0; ion < (nelem + 2); ion++ )
183  {
184  dense.xIonDense[nelem][ion] *= (realnum)hold;
185  }
186  }
187  }
188  }
189 
190  /* this is set false if fluctuations abundances command entered,
191  * when density variations are in place this is true,
192  * and dense.chDenseLaw is "SINE" */
193  if( !dense.lgDenFlucOn )
194  {
195  /* abundance variations are in place */
196  lgChange = true;
197  if( nzone <= 1 )
198  {
199  OldAbun = 1.;
200  FacAbun = 1.;
201  if( dense.lgDenFlucRadius )
202  {
203  /* cycle over radius */
204  FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
206  }
207  else
208  {
209  /* cycle over column density */
210  FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
212  }
213  }
214  else
215  {
216  OldAbun = FacAbunSav;
217  /* rapid abundances fluctuation */
218  if( dense.lgDenFlucRadius )
219  {
220  /* cycle over radius */
221  FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
223  }
224  else
225  {
226  /* cycle over column density */
227  FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
229  }
230  FacAbun = FacAbunSav/OldAbun;
231  }
232  }
233  else
234  {
235  /* abundance variations are NOT in place */
236  FacAbun = 1.;
237  }
238 
239  /* chng 02 dec 11 rjrw -- remove test, saves marginal time could generate nasty intermittent bug when
240  * ( pressure_change_factor*FacAbun == 1 ) && (pressure_change_factor != 1) */
241  if( lgChange )
242  {
243  /* H, He not affected by abundance fluctuations, so only change these
244  * by the pressure change factor */
245  for( nelem=ipHYDROGEN; nelem <= ipHELIUM; ++nelem )
246  {
249  for( ion=0; ion < (nelem + 2); ion++ )
250  {
251  /* >>chng 96 jul 12 had only multiplied total abun, not ions */
253  }
254  }
255 
256  /* Li on up is affect by both pressure and abundance variations,
257  * so multiply by both factors */
258  for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
259  {
260  dense.xMolecules[nelem] *= (realnum)(pressure_change_factor*FacAbun);
261  dense.gas_phase[nelem] *= (realnum)(pressure_change_factor*FacAbun);
262  for( ion=0; ion < (nelem + 2); ion++ )
263  {
264  /* >>chng 96 jul 12 had only multiplied total abun, not ions */
265  dense.xIonDense[nelem][ion] *= (realnum)(pressure_change_factor*FacAbun);
266  }
267  }
268 
270 
271  /* must call TempChange since ionization has changed, there are some
272  * terms that affect collision rates (H0 term in electron collision) */
273  TempChange(phycon.te , false);
274 
275  /* molecules done in hmole, only change pressure, not abundances */
276  for(mol = 0; mol < N_H_MOLEC; mol++)
277  {
279  }
281 
282  /* molecules done in comole */
283  /* >>chng 03 sep 15, upper limit had not included the C, O atoms/ions */
284  /*for( ion=0; ion < NUM_HEAVY_MOLEC; ion++ )*/
285  for( ion=0; ion < mole.num_comole_calc; ion++ )
286  {
287  COmole[ion]->hevmol *= (realnum)(pressure_change_factor*FacAbun);
288 
289  /* check for NaN */
290  ASSERT( !isnan( COmole[ion]->hevmol ) );
291  }
292  }
293 
294  if( trace.lgTrace )
295  {
296  fprintf( ioQQQ,
297  " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
299 
300  if( trace.lgNeBug )
301  {
302  fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n",
303  edensave, dense.eden, edensave/dense.eden );
304  }
305  }
306 
307  {
308  /*@-redef@*/
309  enum {DEBUG_LOC=false};
310  /*@+redef@*/
311  if( DEBUG_LOC && (nzone>215) )
312  {
313  fprintf( ioQQQ,
314  "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n",
315  radius.depth,
322  /* subtract continuum rad pres which has already been added on */
324  wind.windv/1e5,
325  sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5,
326  TorF(conv.lgConvPres) );
327  }
328  }
329 
330  return lgChange;
331 }
332 
333 /*lgConvPres finds amount by which density should change to move towards pressure equilibrium
334  * sets pressure.PresTotlCorrect
335  * returns true if pressure is converged */
336 bool lgConvPres(void)
337 {
338  double dnew,
339  term;
340  bool lgRet;
341 
342  DEBUG_ENTRY( "lgConvPres()" );
343 
344  /* make sure this is set by one of the following branches - set to zero here,
345  * then assert that it is greater than zero at end */
347 
348  /* evaluate a series of possible pressure options, and set the file static variable
349  * pressure_change_factor */
350  /* inside out globule */
351  if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
352  {
353  /* GLBDST is distance from globule, or glbrad-DEPTH */
354  if( radius.glbdst < 0. )
355  {
356  fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
357  fprintf( ioQQQ, " This is routine lgConvPres, GLBDST is%10.2e\n",
358  radius.glbdst );
359  cdEXIT(EXIT_FAILURE);
360  }
364  }
365 
366  /* this is impossible - wind with identically zero velocity */
367  else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && dynamics.lgStatic )
368  {
369  /* this is default case, keep initial H den constant at first zone,
370  * from iteration to iteration */
372  {
374  }
375  else
376  {
377  /* this option is if constant pressure set and flag also
378  * set */
380  /* ratio of correct to current pressures */
382  }
384  }
385 
386  /* this is positive wind velocity the outflowing wind beyond sonic point */
387  else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv > 0. )
388  {
389 
390  /* this is logic for supersonic outflowing wind solution,
391  * which assumes positive velocity, well above sonic point.
392  * following makes sure wind v only updated once per zone */
393  if( nzone > 1 )
394  {
395  /* Wind model */
396 
397  if( trace.lgTrace && trace.lgWind )
398  {
399  fprintf(ioQQQ," lgConvPres sets AccelGravity %.3e lgDisk?%c\n",
400  wind.AccelGravity ,
401  TorF(wind.lgDisk) );
402  }
403 
404 # if 0
405  /* acceleration due to grad P(rad), xMassDensity is gm per cc */
408  /* is this the first zone? */
409  if( nzone > 2 )
410  {
411  /* acceleration due to grad P(rad), xMassDensity is gm per cc */
414  }
415  else
416  {
417  wind.AccelPres = 0.;
418  }
419 # endif
420 
421 # if 0
422  /* this is evaluated in pressure_total*/
423  /* this is numerically unstable */
424  wind.AccelPres = 0.;
425 
426  /* AccelXXX is computed in radinc, is continuum and line acceleration
427  * AccelCont is continuum acceleration
428  * AccelLine is line acceleration
429  * AccelPres is gradient of local gas pressure */
431 
432  if( nzone==NzoneEval+1)
433  {
434  double da = AccelNet - AccelNetLast;
435  /* we have previous acceleration in this iteration, use it
436  * a = a_i + (a_i - a_i-1)/(r_i - r_i-1) * (r_i+1 - r_i)/2 */
437  AccelNet += da* (radius.drNext/radius.drad) /2.;
438 
439  }
440 # endif
441 
442  /* following is form of energy equation
443  * struc.windv[nzone-2] is velocity of previous zone
444  * this increments that velocity to form square of new wind
445  * velocity for outer edge of this zone */
446  term = POW2(struc.windv[nzone-2]) + 2.*(wind.AccelTot - wind.AccelGravity)* radius.drad;
447 
448  /* increment velocity if it is substantially positive */
449  if( term <= 1e3 )
450  {
451  /* wind velocity is well below sonic point, give up,
452  * do not change velocity */
453  wind.lgVelPos = false;
454  }
455  else
456  {
457  /* struc.windv[nzone-2] is velocity of previous zone
458  * this increments that velocity to form square of new wind
459  * velocity for outer edge of this zone
460  double windnw = (double)POW2(struc.windv[nzone-2]) +
461  (double)(2.*(wind.AccelTot-wind.AccelGravity))*radius.drad;*/
462 
463  /* wind.windv is velocity at OUTER edge of this zone */
464  wind.windv = (realnum)sqrt(term);
465  wind.lgVelPos = true;
466  }
467 
468  /* following needed for expansion cooling */
470  radius.Radius);
471 
472  if( trace.lgTrace && trace.lgWind )
473  {
474  fprintf(ioQQQ," lgConvPres new wind V zn%li %.3e AccelTot %.3e AccelGravity %.3e\n",
476  }
477  }
478  else
479  {
481  }
482 
483  /* conservation of mass sets density here */
486  }
487 
488  /* this is negative wind velocity the new dynamics */
489  else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
490  {
491  /* sets pressure.PresTotlCorrect */
493  }
494 
495  /* this is impossible - wind with identically zero velocity */
496  else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv == 0. )
497  {
498  /* declare insanity */
499  fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
500  /* TotalInsanity announces fatal problem, ShowMe, then cdEXIT with failure */
501  TotalInsanity();
502  }
503 
504  else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
505  {
506  /* rapid density fluctuation */
507  if( dense.lgDenFlucRadius )
508  {
511  }
512  else
513  {
516  }
518  }
519 
520  else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
521  {
522  /* power law function of radius */
523  dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower);
526  }
527 
528  else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
529  {
530  /* power law function of depth */
531  dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower);
534  }
535 
536  else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
537  {
538  /* power law function of column density */
539  dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/
543  }
544 
545  else if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
546  {
547  /* constant pressure */
549  {
550  /* >>chng 01 oct 31, replace pneed with CorretPressure */
551  /* this has pressure due to incident continuum */
553  }
554  else
555  {
556  /* this does not have pressure due to incident continuum*/
558  /* following term normally unity, power law set with option par on cmmnd*/
560  }
561 
562  /* ratio of correct to current pressures */
564  }
565 
566  else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 )
567  {
568  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
569  {
570  /* call ACF sub */
572  }
573  else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
574  {
575  /* call table interpolation subroutine
576  * >>chng 96 nov 29, added dense_tabden */
578  }
579  else
580  {
581  fprintf( ioQQQ, " Insanity, lgConvPres gets chCPres=%4.4s\n",
582  dense.chDenseLaw );
583  cdEXIT(EXIT_FAILURE);
584  }
586  /* >>chng 06 feb 21, from above to below. we want to look at change,
587  * not the value. Bug found by Valentina Luridiana */
588  /*if( pressure.PresTotlCorrect > 3. || pressure.PresTotlCorrect< 1./3 )*/
590  {
591  static bool lgWARN2BIG=false;
592  if( !lgWARN2BIG )
593  {
594  lgWARN2BIG = true;
595  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);
596  fprintf(ioQQQ," >========== Warning! This will cause convergence problems. \n");
597  fprintf(ioQQQ," >========== Warning! The current radius is %.3e. \n",radius.Radius);
598  fprintf(ioQQQ," >========== Warning! The current depth is %.3e. \n",radius.depth);
599  fprintf(ioQQQ," >========== Warning! Consider using a more modest change in density vs radius. \n\n\n");
600  }
601  }
602  }
603 
604  else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
605  {
606  /* this is the default, constant density */
609  }
610 
611  else
612  {
613  fprintf( ioQQQ, " Unknown pressure law=%s= This is a major internal error.\n",
614  dense.chDenseLaw );
615  ShowMe();
616  cdEXIT(EXIT_FAILURE);
617  }
618 
619  /* one of the branches above must have reset this variable,
620  * and it was init to 0 at start. Confirm that non-zero */
621  ASSERT( pressure.PresTotlCorrect > FLT_MIN );
622 
623  /* now see whether current pressure is within error bounds */
625  {
626  lgRet = false;
627  conv.lgConvPres = false;
628  }
629  else
630  {
631  lgRet = true;
632  conv.lgConvPres = true;
633  }
634 
635  return lgRet;
636 }

Generated for cloudy by doxygen 1.8.1.1