cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pressure_total.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 /*PresTotCurrent determine the gas and line radiation pressures for current conditions,
4  * this sets values of pressure.PresTotlCurr, also calls tfidle */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "taulines.h"
8 #include "opacity.h"
9 #include "hextra.h"
10 #include "elementnames.h"
11 #include "hydrogenic.h"
12 #include "conv.h"
13 #include "iso.h"
14 #include "wind.h"
15 #include "magnetic.h"
16 #include "doppvel.h"
17 #include "rfield.h"
18 #include "phycon.h"
19 #include "thermal.h"
20 #include "hmi.h"
21 #include "h2.h"
22 #include "dense.h"
23 #include "atomfeii.h"
24 #include "mole.h"
25 #include "dynamics.h"
26 #include "trace.h"
27 #include "rt.h"
28 #include "atmdat.h"
29 #include "lines_service.h"
30 #include "pressure.h"
31 #include "radius.h"
32 
33 /* this sets values of pressure.PresTotlCurr, also calls tfidle */
34 void PresTotCurrent(void)
35 {
36  long int i,
37  ion,
38  ipISO ,
39  ipHi,
40  ipLo,
41  nelem;
42 
43  static long int
44  /* used in debug print statement to see which line adds most pressure */
45  ipLineTypePradMax=-1 ,
46  /* points to line causing radiation pressure */
47  ipLinePradMax=-LONG_MAX,
48  ip2=-LONG_MAX,
49  ip3=-LONG_MAX,
50  ip4=-LONG_MAX;
51 
52  /* the line radiation pressure variables that must be preserved since
53  * a particular line radiation pressure may not be evaluated if it is
54  * not important */
55  static double Piso_seq[NISO]={0.,0.},
56  PLevel1=0.,
57  PLevel2=0.,
58  PHfs=0.,
59  PCO=0.,
60  P_H2=0.,
61  P_FeII=0.;
62 
63  double
64  RadPres1,
65  TotalPressure_v,
66  pmx;
67  realnum den_hmole ,
68  den_comole ,
69  DenseAtomsIons,
70  DensElements;
71 
72  DEBUG_ENTRY( "PresTotCurrent()" );
73 
74  if( !conv.nTotalIoniz )
75  {
76  /* resetting ipLinePradMax, necessary for
77  * multiple cloudy calls in single coreload. */
78  ipLinePradMax=-LONG_MAX;
79  //pressure.PresTotlCurr = 0.;
80  }
81 
82  /* PresTotCurrent - evaluate total pressure, dyne cm^-2
83  * and radiative acceleration */
84 
85  /* several loops over the emission lines for radiation pressure and
86  * radiative acceleration are expensive due to the number of lines and
87  * the memory they occupy. Many equations of state do not include
88  * radiation pressure or radiative acceleration. Only update these
89  * when included in EOS - when not included only evaluate once per zone.
90  * do it once per zone since we will still report these quantities.
91  * this flag says whether we must update all terms */
92  bool lgMustEvaluate = false;
93 
94  /* this says we already have an evaluation in this zone so do not
95  * evaluate terms known to be trivial, even when reevaluating major terms */
96  bool lgMustEvaluateTrivial = false;
97  /* if an individual agent is larger than this fraction of the total current
98  * radiation pressure then it is not trivial
99  * conv.PressureErrorAllowed is relative error allowed in pressure */
100  double TrivialLineRadiationPressure = conv.PressureErrorAllowed/10. *
102 
103  /* reevaluate during search phase when conditions are changing dramatically */
104  if( conv.lgSearch )
105  {
106  lgMustEvaluate = true;
107  lgMustEvaluateTrivial = true;
108  }
109 
110  /* reevaluate if zone or iteration has changed */
111  static long int nzoneEvaluated=0, iterationEvaluated=0;
112  if( nzone!=nzoneEvaluated || iteration!=iterationEvaluated )
113  {
114  lgMustEvaluate = true;
115  lgMustEvaluateTrivial = true;
116  /* this flag says which source of radiation pressure was strongest */
117  ipLineTypePradMax = -1;
118  }
119 
120  /* constant pressure or dynamical sim - we must reevaluate terms
121  * but do not need to reevaluate trivial contributors */
122  if( (strcmp(dense.chDenseLaw,"WIND") == 0 ) ||
123  (strcmp(dense.chDenseLaw,"CPRE") == 0 ) )
124  lgMustEvaluate = true;
125 
126  if( lgMustEvaluate )
127  {
128  /* we are reevaluating quantities in this zone and iteration,
129  * remember that we did it */
130  nzoneEvaluated = nzone;
131  iterationEvaluated = iteration;
132  }
133 
134  /* update density - temperature variables which may have changed */
135  /* must call TempChange since ionization has changed, there are some
136  * terms that affect collision rates (H0 term in electron collision) */
137  TempChange(phycon.te , false);
138 
139  /* find total number of atoms and ions */
140  DenseAtomsIons = 0.;
141  DensElements = 0.;
142  realnum factor = 1.1f;
143  if( conv.lgSearch )
144  factor = 2.f;
145  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
146  {
147  /* only check element solution if it is on, and ionization
148  * has not been set with TABLE ELEMENT command */
149  if( dense.lgElmtOn[nelem] )
150  {
151  /* gas_phase is density of all atoms and ions, but not molecules */
152  DensElements += dense.gas_phase[nelem];
153  realnum DenseOne = 0;
154  for( ion=0; ion<=nelem+1; ++ion )
155  DenseOne += dense.xIonDense[nelem][ion];
156 
157  /* during search phase sums may not add up correctly, and during
158  * early parts of search phase may be badly off. This test
159  * was introduced as result off fpe due to another
160  * reason - Te had changed too much during initial search
161  * for sim in which chemistry was important - fix was to not
162  * let Te change. During resulting insanity caused by large
163  * change, linearization did not work, CO and ionization solvers
164  * could diverge leading to molecule or ionization population
165  * larger than 1e38 limit to a realnum, fpe followed. this is to
166  * catch that in its early stages */
167  if( conv.nTotalIoniz>5 && !dense.lgSetIoniz[nelem] &&
168  DenseOne > dense.gas_phase[nelem]*factor )
169  {
170  /* species is not conserved */
171  fprintf(ioQQQ,"\n\n PROBLEM DISASTER PressureTotal: the chemical species "
172  "are not conserved.\n");
173  fprintf(ioQQQ," The sum of the densities of the atoms and ions for "
174  "%s is %.3e which is greater than the gas-phase abundance "
175  "of the element, %.3e.\n",
176  elementnames.chElementName[nelem],DenseOne,
177  dense.gas_phase[nelem]);
178  /* not including cosmic rays is the most likely cause of this */
179  if( hextra.cryden == 0. )
180  {
181  fprintf(ioQQQ," The chemistry network is known to fail this way"
182  " in cold molecular environments when cosmic rays are not"
183  " included. They were not - try including them with the"
184  " COSMIC RAY BACKBROUND command.\n");
185  fprintf(ioQQQ," The calculation is aborting. conv.nTotalIoniz=%li\n "
186  "Sorry.\n\n", conv.nTotalIoniz );
187  lgAbort = true;
188  }
189  /* if they were included then something seriously wrong has happened*/
190  else
191  TotalInsanity();
192  }
193  DenseAtomsIons += DenseOne;
194  }
195  }
196  /* DensElements is sum of all gas-phase densities of all elements,
197  * DenseAtomsIons is sum of density of atoms and ions, does not
198  * include molecules */
199  ASSERT( conv.nTotalIoniz<5 || DenseAtomsIons <= DensElements*factor );
200  ASSERT( DenseAtomsIons > 0. );
201 
202  /* evaluate the total ionization energy on a scale where neutral atoms
203  * correspond to an energy of zero, so the ions have a positive energy */
205 #if 0
206  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
207  {
208  for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
209  {
210  /* lgMETALS is true, set false with "no advection metals" command */
211  int kadvec = dynamics.lgMETALS;
212  /* this gives the iso sequence for this ion - should it be included in the
213  * advected energy equation? lgISO is true, set false with
214  * "no advection H-like" and "he-like" - for testing*/
215  ipISO = nelem-ion;
216  fixit(); /* should this be kadvec = kadvec && dynamics.lgISO[ipISO]; ? */
217  if( ipISO >= 0 && ipISO<NISO )
218  kadvec = dynamics.lgISO[ipISO];
219  for( i=1; i<=ion; ++i )
220  {
221  long int nelec = nelem+1 - i;
222  /* this is the sum of all the energy needed to bring the atom up
223  * to the ion+1 level of ionization - at this point a positive number */
224  phycon.EnergyIonization += dense.xIonDense[nelem][ion] *
225  t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD* 0.9998787*kadvec;
226  }
227  }
228  }
229  /* convert to ergs/cm^3 */
231 #endif
232 
236  phycon.EnergyBinding = 0.;
237 
238  /* find number of molecules of the heavy elements in the gas phase.
239  * Atoms and ions were counted above. Do not include ices, solids
240  * on grains */
241  den_comole = 0.;
242  for( i=0; i < mole.num_comole_calc; i++ )
243  {
244  /* term COmole[i]->lgGas_Phase is 0 or 1 depending on whether molecule
245  * is on grain or in gas phase
246  * n_nuclei is number of nuclei in molecule, this tests whether
247  * actually an atom or molecule - atoms were already counted above */
248  if(COmole[i]->n_nuclei > 1)
249  den_comole += COmole[i]->hevmol * COmole[i]->lgGas_Phase;
250  }
251 
252  /* number of hydrogen molecules, loop over all H molecular species,
253  * do not include H0, H+ */
254  den_hmole = 0.;
255  for( i=0; i<N_H_MOLEC; ++i )
256  {
257  if( i!=ipMH && i!=ipMHp )
258  den_hmole += hmi.Hmolec[i];
259  }
260 
261  /* total number of atoms, ions, and molecules in gas phase per unit vol */
262  dense.xNucleiTotal = den_hmole + den_comole + DenseAtomsIons;
263  if( dense.xNucleiTotal > BIGFLOAT )
264  {
265  fprintf(ioQQQ, "PROBLEM DISASTER pressure_total has found "
266  "dense.xNucleiTotal with an insane density.\n");
267  fprintf(ioQQQ, "The density was %.2e\n",
269  TotalInsanity();
270  }
271  ASSERT( dense.xNucleiTotal > 0. );
272 
273  /* particle density that enters into the pressure includes electrons
274  * cm-3 */
276 
277  /* the current gas pressure */
279  /*fprintf(ioQQQ,"DEBUG gassss %.2f %.4e %.4e %.4e \n",
280  fnzone, pressure.PresGasCurr , dense.pden , pressure.PresInteg );*/
281 
282  /* dense.wmole is molecular weight, AtomicMassUnit per particle */
283  dense.wmole = 0.;
284  for( i=0; i < LIMELM; i++ )
285  {
287  }
288 
289  /* dense.wmole is now molecular weight, AtomicMassUnit per particle */
290  dense.wmole /= dense.pden;
291  ASSERT( dense.wmole > 0. && dense.pden > 0. );
292 
293  /* xMassDensity is density in grams per cc */
297 
298  /*>>chng 04 may 25, introduce following test on xMassDensity0
299  * WJH 21 may 04, this is the mass density that corresponds to the hden
300  * specified in the init file. It is used to calculate the mass flux in the
301  * dynamics models. It may not necessarily be the same as struc.DenMass[0],
302  * since the pressure solver may have jumped to a different density at the
303  * illuminated face from that specified.*/
304  if( dense.xMassDensity0 < 0.0 )
306 
307  /* >>chng 01 nov 02, move to here from dynamics routines */
308  /* >>chng 02 jun 18 (rjrw), fix to use local values */
309  /* WJH 21 may 04, now recalculate wind v for the first zone too.
310  * This is necessary when we are forcing the sub or supersonic branch */
311  if(wind.windv < 0)
313 
314  /* this is the current ram pressure, at this location */
316 
317  /* this is the current turbulent pressure, not now added to equation of state
318  * >>chng 06 mar 29, add Heiles_Troland_F / 6. term*/
321 
324  /* radiative acceleration, lines and continuum */
325  if( lgMustEvaluate )
326  {
327  /* continuous radiative acceleration */
328  double rforce = 0.;
329  for( i=0; i < rfield.nflux; i++ )
330  {
331  rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
332  rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]);
333  }
334 
335  /* line acceleration; xMassDensity is gm per cc */
338  /* this is numerically unstable */
339  wind.AccelPres = 0.;
340  /* total acceleration */
342  /* G, COMASS = mass of star in solar masses */
343  double reff = radius.Radius-radius.drad/2.;
345  (1.-wind.DiskRadius/reff) );
346 # if 0
347  if( fudge(-1) )
348  fprintf(ioQQQ,"DEBUG pressure_total updates AccelTot to %.2e grav %.2e\n",
350 # endif
351  }
352 
353  /* must always evaluate H La width since used elsewhere */
354  if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc > 0. )
355  {
357 # if 0
358  fprintf(ioQQQ,"DEBUG widLya\t%li\t%.2f\t%.3e",
359  iteration,
360  fnzone,
361  hydro.HLineWidth);
363  RT_LyaWidth(
364  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn,
365  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot,
366  4.72e-2/phycon.sqrte,
368  fprintf(ioQQQ,"\t%.3e\n",
369  hydro.HLineWidth);
370 # endif
371  }
372  else
373  hydro.HLineWidth = 0.;
374 
375 
376  /* find max rad pressure even if capped
377  * lgLineRadPresOn is turned off by NO RADIATION PRESSURE command */
378  if( trace.lgTrace )
379  {
380  fprintf(ioQQQ,
381  " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
382  "lgMustEvaluateTrivial:%c "
383  "pressure.lgLineRadPresOn:%c "
384  "rfield.lgDoLineTrans:%c \n",
385  nzone , iteration , TorF(lgMustEvaluate) , TorF(lgMustEvaluateTrivial),
387  }
388 
389  if( lgMustEvaluate && pressure.lgLineRadPresOn && rfield.lgDoLineTrans )
390  {
391  /* RadPres is pressure due to lines, lgPres_radiation_ON turns off or on */
393  /* used to remember largest radiation pressure source */
394  pmx = 0.;
395  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
396  {
397  if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
398  {
399  Piso_seq[ipISO] = 0.;
400  for( nelem=ipISO; nelem < LIMELM; nelem++ )
401  {
402  /* does this ion stage exist? */
403  if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
404  {
405  /* do not include highest levels since maser can occur due to topoff,
406  * and pops are set to small number in this case */
407  for( ipHi=1; ipHi <iso.numLevels_local[ipISO][nelem] - iso.nCollapsed_local[ipISO][nelem]; ipHi++ )
408  {
409  for( ipLo=0; ipLo < ipHi; ipLo++ )
410  {
411  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
412  continue;
413 
414  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > iso.SmallA );
415 
416  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc > SMALLFLOAT &&
417  /* test that have not overrun optical depth scale */
418  ( (Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot - Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn) > SMALLFLOAT ) &&
419  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc > FLT_EPSILON*100. )
420  {
421  RadPres1 = PressureRadiationLine( &Transitions[ipISO][nelem][ipHi][ipLo], dense.xIonDense[nelem][nelem+1-ipISO] );
422 
423  if( RadPres1 > pmx )
424  {
425  ipLineTypePradMax = 2;
426  pmx = RadPres1;
427  ip4 = ipISO;
428  ip3 = nelem;
429  ip2 = ipHi;
430  ipLinePradMax = ipLo;
431  }
432  Piso_seq[ipISO] += RadPres1;
433  {
434  /* option to print particulars of some line when called */
435  enum {DEBUG_LOC=false};
436  if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 )
437  {
438  fprintf(ioQQQ,
439  "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
440  fnzone,
441  RadPres1,
442  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc,
443  StatesElem[ipISO][nelem][ipLo].Pop,
444  StatesElem[ipISO][nelem][ipHi].Pop,
445  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc);
446  }
447  }
448  }
449  }
450  }
451  }
452  }
453  ASSERT( Piso_seq[ipISO] >= 0. );
454  }
455  pressure.pres_radiation_lines_curr += Piso_seq[ipISO];
456  }
457  {
458  /* option to print particulars of some line when called */
459  enum {DEBUG_LOC=false};
460 # if 0
461  if( DEBUG_LOC /*&& iteration > 1*/ && nzone > 260 )
462  {
463  fprintf(ioQQQ,
464  "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
465  nzone,
466  pmx,
467  Transitions[ipISO][ip3][ipLinePradMax][ip2].Emis->PopOpc,
468  StatesElem[ipISO][ip3][0].Pop,
469  StatesElem[ipISO][ip3][2].Pop,
470  StatesElem[ipISO][ip3][3].Pop,
471  StatesElem[ipISO][ip3][4].Pop,
472  StatesElem[ipISO][ip3][5].Pop,
473  StatesElem[ipISO][ip3][6].Pop);
474  }
475 # endif
476  if( DEBUG_LOC /*&& iteration > 1 && nzone > 150 */)
477  {
478  fprintf(ioQQQ,
479  "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
480  pmx,
481  ip4,
482  ip3,
483  ip2,
484  ipLinePradMax,
485  Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->PopOpc,
486  StatesElem[ip4][ip3][ip2].Pop,
487  1.-Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc );
488  }
489  }
490 
491  if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure )
492  {
493  PLevel1 = 0.;
494  /* line radiation pressure from large set of level 1 lines */
495  for( i=1; i <= nLevel1; i++ )
496  {
497  if( TauLines[i].Hi->Pop > 1e-30 )
498  {
499  RadPres1 = PressureRadiationLine( &TauLines[i], 1. );
500 
501  if( RadPres1 > pmx )
502  {
503  ipLineTypePradMax = 3;
504  pmx = RadPres1;
505  ipLinePradMax = i;
506  }
507  PLevel1 += RadPres1;
508  }
509  }
510  ASSERT( PLevel1 >= 0. );
511  }
513 
514  if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
515  {
516  /* level 2 lines */
517  PLevel2 = 0.;
518  for( i=0; i < nWindLine; i++ )
519  {
520  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
521  {
522  if( TauLine2[i].Hi->Pop > 1e-30 )
523  {
524  RadPres1 = PressureRadiationLine( &TauLine2[i], 1. );
525 
526  PLevel2 += RadPres1;
527  if( RadPres1 > pmx )
528  {
529  ipLineTypePradMax = 4;
530  pmx = RadPres1;
531  ipLinePradMax = i;
532  }
533  }
534  }
535  }
536  ASSERT( PLevel2 >= 0. );
537  }
539 
540  /* fine structure lines */
541  if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
542  {
543  PHfs = 0.;
544  for( i=0; i < nHFLines; i++ )
545  {
546  if( HFLines[i].Hi->Pop > 1e-30 )
547  {
548  RadPres1 = PressureRadiationLine( &HFLines[i], 1. );
549 
550  PHfs += RadPres1;
551  if( RadPres1 > pmx )
552  {
553  ipLineTypePradMax = 8;
554  pmx = RadPres1;
555  ipLinePradMax = i;
556  }
557  }
558  }
559  ASSERT( PHfs >= 0. );
560  }
562 
563  /* radiation pressure due to H2 lines */
564  if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
565  {
566  P_H2 = H2_RadPress();
567  /* flag to remember H2 radiation pressure */
568  if( P_H2 > pmx )
569  {
570  pmx = P_H2;
571  ipLineTypePradMax = 9;
572  }
573  ASSERT( P_H2 >= 0. );
574  }
576 
577  /* radiation pressure due to FeII lines and large atom */
578  if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure )
579  {
580  P_FeII = FeIIRadPress();
581  if( P_FeII > pmx )
582  {
583  pmx = P_FeII;
584  ipLineTypePradMax = 7;
585  }
586  ASSERT( P_FeII >= 0. );
587  }
589 
590  /* co carbon monoxide lines */
591  if( lgMustEvaluateTrivial || PCO > TrivialLineRadiationPressure )
592  {
593  PCO = 0.;
594  for( i=0; i < nCORotate; i++ )
595  {
596  if( C12O16Rotate[i].Hi->Pop > 1e-30 )
597  {
598  RadPres1 = PressureRadiationLine( &C12O16Rotate[i], 1. );
599 
600  PCO += RadPres1;
601  if( RadPres1 > pmx )
602  {
603  ipLineTypePradMax = 5;
604  pmx = RadPres1;
605  ipLinePradMax = i;
606  }
607  }
608  if( C13O16Rotate[i].Hi->Pop > 1e-30 )
609  {
610  RadPres1 = PressureRadiationLine( &C13O16Rotate[i], 1. );
611 
612  PCO += RadPres1;
613  if( RadPres1 > pmx )
614  {
615  ipLineTypePradMax = 6;
616  pmx = RadPres1;
617  ipLinePradMax = i;
618  }
619  }
620  }
621  ASSERT( PCO >= 0. );
622  }
624  }
626  {
627  /* case where radiation pressure is not turned on */
628  ipLinePradMax = -1;
629  ipLineTypePradMax = 0;
630  }
631 
632  /* the ratio of radiation to total (gas + continuum + etc) pressure */
634 
635  /* this would be a major logical error */
637  {
638  fprintf(ioQQQ,
639  "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
640  Piso_seq[ipH_LIKE],
641  Piso_seq[ipHE_LIKE],
642  PLevel1,
643  PLevel2,
644  PCO);
645  ShowMe();
646  cdEXIT(EXIT_FAILURE);
647  }
648 
649  /* following test will never succeed, here to trick lint, ipLineTypePradMax is only used
650  * when needed for debug */
651  if( trace.lgTrace && ipLineTypePradMax <0 )
652  {
653  fprintf(ioQQQ,
654  " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
655  pressure.pbeta,ipLineTypePradMax,ipLinePradMax );
656  }
657 
658  /* this is the total internal energy of the gas, the amount of energy needed
659  * to produce the current level populations, relative to ground,
660  * only used for energy terms in advection */
662 #if 0
663  fixit(); /* the quantities phycon.EnergyExcitation, phycon.EnergyIonization, phycon.EnergyBinding
664  * are not used anywhere, except in print statements... */
665  broken(); /* the code below contains serious bugs! It is supposed to implement loops
666  * over quantum states in order to evaluate the potential energy stored in
667  * excited states of all atoms, ions, and molecules. The code below however
668  * often implements loops over all combinations of upper and lower levels!
669  * This code needs to be rewritten once quantumstates are fully implemented. */
670  ipLo = 0;
671  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
672  {
673  for( nelem=ipISO; nelem < LIMELM; nelem++ )
674  {
675  if( dense.IonHigh[nelem] == nelem + 1 - ipISO )
676  {
677  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
678  for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
679  {
681  Transitions[ipISO][nelem][ipHi][ipLo].Hi->Pop *
682  Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg*
683  /* last term is option to turn off energy term, no advection hlike, he-like */
684  dense.xIonDense[nelem][nelem+1-ipISO]*dynamics.lgISO[ipISO];
685  }
686  }
687  }
688  }
689 
690  if( dynamics.lgISO[ipH_LIKE] )
691  /* internal energy of H2 */
693 
694  /* this is option to turn off energy effects of advection, no advection metals */
695  if( dynamics.lgMETALS )
696  {
697  for( i=1; i <= nLevel1; i++ )
698  {
700  TauLines[i].Hi->Pop* TauLines[i].EnergyErg;
701  }
702  for( i=0; i < nWindLine; i++ )
703  {
704  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
705  {
707  TauLine2[i].Hi->Pop* TauLine2[i].EnergyErg;
708  }
709  }
710  for( i=0; i < nHFLines; i++ )
711  {
713  HFLines[i].Hi->Pop* HFLines[i].EnergyErg;
714  }
715  double Energy12 = 0.;
716  double Energy13 = 0.;
717  for( i=0; i < nCORotate; i++ )
718  {
719  Energy12 += C12O16Rotate[i].EnergyErg;
721  C12O16Rotate[i].Hi->Pop* Energy12;
722 
723  Energy13 += C13O16Rotate[i].EnergyErg;
725  C13O16Rotate[i].Hi->Pop* Energy13;
726  }
727 
728  /* internal energy of large FeII atom */
730  }
731 #endif
732 
733  /* ==================================================================
734  * end internal energy of atoms and molecules */
735 
736  /* evaluate some parameters to do with magnetic field */
738 
739  /*lint -e644 Piso_seq not init */
741  {
742  fprintf(ioQQQ,
743  " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e "
744  "gas parts; H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n",
745  fnzone,
751  Piso_seq[ipH_LIKE]/pressure.PresTotlCurr,
752  Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr,
753  PLevel1/pressure.PresTotlCurr,
754  PLevel2/pressure.PresTotlCurr,
756  PHfs/pressure.PresTotlCurr,
757  P_H2/pressure.PresTotlCurr,
759  /*lint +e644 Piso_seq not initialized */
760  }
761 
762  /* Conserved quantity in steady-state energy equation for dynamics:
763  * thermal energy, since recombination is treated as cooling
764  * (i.e. is loss of electron k.e. to emitted photon), so don't
765  * include
766  * ...phycon.EnergyExcitation + phycon.EnergyIonization + phycon.EnergyBinding
767  * */
768 
769  /* constant density case is also hypersonic case */
771  {
772  /* this branch is time dependent AND constant density */
773  /*fprintf(ioQQQ,"DEBUG enthalpy HIT1\n");*/
774  /* this is the time-varying case where density is constant */
775  /* \todo 1 this has 3/2 on the PresGasCurr while true dynamics case below
776  * has 5/2 - so this is not really the enthalpy density - better
777  * would be to always use this term and add the extra PresGasCurr
778  * when the enthalpy is actually needed */
780  0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */
781  3./2.*pressure.PresGasCurr + /* thermal plus PdV work */
782  magnetic.EnthalpyDensity; /* magnetic terms */
783  }
784  else
785  {
786  /* this branch either advective or constant pressure */
787  /*fprintf(ioQQQ,"DEBUG enthalpy HIT2\n");*/
788  /* this is usual dynamics case, or time-varying case where pressure
789  * is kept constant and PdV work is added to the cell */
791  0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */
792  5./2.*pressure.PresGasCurr + /* thermal plus PdV work */
793  magnetic.EnthalpyDensity; /* magnetic terms */
794  }
795 
796  /* stop model from running away on first iteration, when line optical
797  * depths are not defined correctly anyway.
798  * if( iter.le.itermx .and. RadPres.ge.GasPres ) then
799  * >>chng 97 jul 23, only cap radiation pressure on first iteration */
800  if( iteration <= 1 && pressure.pres_radiation_lines_curr >= pressure.PresGasCurr )
801  {
802  /* stop RadPres from exceeding the gas pressure on first iteration */
805  pressure.lgPradCap = true;
806  }
807 
808  /* remember the globally most important line, in the entire model
809  * test on nzone so we only do test after solution is stable */
811  {
813  pressure.ipPradMax_line = ipLinePradMax;
815  if( ipLineTypePradMax == 2 )
816  {
817  /* hydrogenic */
818  strcpy( pressure.chLineRadPres, "ISO ");
819  ASSERT( ip4 < NISO && ip3<LIMELM );
820  ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
821  strcat( pressure.chLineRadPres, chLineLbl(&Transitions[ip4][ip3][ip2][ipLinePradMax]) );
822  {
823  /* option to print particulars of some line when called */
824  enum {DEBUG_LOC=false};
825  /*lint -e644 Piso_seq not initialized */
826  /* trace serious constituents in radiation pressure */
827  if( DEBUG_LOC )
828  {
829  fprintf(ioQQQ,"DEBUG iso prad\t%li\t%li\tISO,nelem=\t%li\t%li\tnu, no=\t%li\t%li\t%.4e\t%.4e\t%e\t%e\t%e\n",
830  iteration, nzone,
831  ip4,ip3,ip2,ipLinePradMax,
832  Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauIn,
833  Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauTot,
834  Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc,
835  Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pelec_esc,
836  Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pdest);
837  if( ip2==5 && ipLinePradMax==4 )
838  {
839  double width;
840  fprintf(ioQQQ,"hit it\n");
841  width = RT_LineWidth(&Transitions[ip4][ip3][ip2][ipLinePradMax]);
842  fprintf(ioQQQ,"DEBUG width %e\n", width);
843  }
844  }
845  }
846  }
847  else if( ipLineTypePradMax == 3 )
848  {
849  /* level 1 lines */
850  ASSERT( ipLinePradMax>=0 );
851  strcpy( pressure.chLineRadPres, "Level1 ");
852  strcat( pressure.chLineRadPres, chLineLbl(&TauLines[ipLinePradMax]) );
853  }
854  else if( ipLineTypePradMax == 4 )
855  {
856  /* level 2 lines */
857  ASSERT( ipLinePradMax>=0 );
858  strcpy( pressure.chLineRadPres, "Level2 ");
859  strcat( pressure.chLineRadPres, chLineLbl(&TauLine2[ipLinePradMax]) );
860  }
861  else if( ipLineTypePradMax == 5 )
862  {
863  /* c12o16 carbon monoxide lines */
864  ASSERT( ipLinePradMax>=0 );
865  strcpy( pressure.chLineRadPres, "12CO ");
866  strcat( pressure.chLineRadPres, chLineLbl(&C12O16Rotate[ipLinePradMax]) );
867  }
868  else if( ipLineTypePradMax == 6 )
869  {
870  /* c13o16 carbon monoxide lines */
871  ASSERT( ipLinePradMax>=0 );
872  strcpy( pressure.chLineRadPres, "13CO ");
873  strcat( pressure.chLineRadPres, chLineLbl(&C13O16Rotate[ipLinePradMax]) );
874  }
875  else if( ipLineTypePradMax == 7 )
876  {
877  /* FeII lines */
878  strcpy( pressure.chLineRadPres, "Fe II ");
879  }
880  else if( ipLineTypePradMax == 8 )
881  {
882  /* hyperfine struct lines */
883  strcpy( pressure.chLineRadPres, "hyperf ");
884  ASSERT( ipLinePradMax>=0 );
885  strcat( pressure.chLineRadPres, chLineLbl(&HFLines[ipLinePradMax]) );
886  }
887  else if( ipLineTypePradMax == 9 )
888  {
889  /* large H2 lines */
890  strcpy( pressure.chLineRadPres, "large H2 ");
891  }
892  else
893  {
894  fprintf(ioQQQ," PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
895  ShowMe();
896  cdEXIT(EXIT_FAILURE);
897  }
898  }
899 
900  if( trace.lgTrace && pressure.pbeta > 0.5 )
901  {
902  fprintf(ioQQQ,
903  " PresTotCurrent Pbeta:%.2f due to %s\n",
904  pressure.pbeta ,
906  );
907  }
908 
909  /* >>chng 02 jun 27 - add in magnetic pressure */
910  /* start to bring total pressure together */
911  TotalPressure_v = pressure.PresGasCurr;
912 
913  /* these flags are set false by default since constant density is default,
914  * set true for constant pressure or dynamics */
915  TotalPressure_v += pressure.PresRamCurr * pressure.lgPres_ram_ON;
916 
917  /* magnetic pressure, evaluated in magnetic.c - this can be negative for an ordered field!
918  * option is on by default, turned off with constant density, or constant gas pressure, cases */
921  TotalPressure_v += magnetic.pressure * pressure.lgPres_magnetic_ON;
922 
923  /* turbulent pressure
924  * >>chng 06 mar 24, include this by default */
925  TotalPressure_v += pressure.PresTurbCurr * DoppVel.lgTurb_pressure;
926 
927  /* radiation pressure only included in total eqn of state when this flag is
928  * true, set with constant pressure command */
929  /* option to add in internal line radiation pressure */
931 
932  {
933  enum{DEBUG_LOC=false};
934  if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT /*&& iteration > 1*/ )
935  {
936  fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
937  nzone,
940  TotalPressure_v ,
944  );
945  }
946  }
947 
948  if( TotalPressure_v< 0. && magnetic.pressure < 0. )
949  {
950  /* negative pressure due to ordered field overwhelms total pressure - punt */
951  fprintf(ioQQQ," The negative pressure due to ordered magnetic field overwhelms the total pressure - please reconsider the geometry & field.\n");
952  cdEXIT(EXIT_FAILURE);
953  }
954 
955  ASSERT( TotalPressure_v > 0. );
956 
957  /* remember highest pressure anywhere */
958  pressure.PresMax = MAX2(pressure.PresMax,(realnum)TotalPressure_v);
959 
960  /* this is what we came for - set the current pressure */
961  pressure.PresTotlCurr = TotalPressure_v;
962 
963 # if 0
964  /* this is special case where we are working on first zone, at
965  * illuminated face of the cloud. Here we want to remember the
966  * initial pressure in case constant pressure is needed */
967  /* >>chng 05 jan 10, chng from nzone==1 to nzone<=1 so pressure not changed
968  * during search phase */
969  /*>>chng 06 jun 20, add test on first iteration, or we are holding
970  * density constant - flag dense.lgDenseInitConstant false if
971  * constant pressure reset is used - default is true, after first iteration
972  * initial density is kept constant, when set false with reset option on
973  * constant pressure density on first iteration is allowed to change to keep
974  * pressure constant */
975  if( nzone <= 1 && (iteration==1 || dense.lgDenseInitConstant) )
976  {
977  double PresTotlInitSave;
978  double PresRamInitSave;
979  /* this is first zone, lock onto pressure */
980  if( conv.nTotalIoniz )
981  {
982  PresTotlInitSave = pressure.PresTotlInit;
983  PresRamInitSave = pressure.PresRamInit;
984  }
985  else
986  {
987  PresTotlInitSave = 0.;
988  PresRamInitSave = 0.;
989  }
991  pressure.PresRamInit = pressure.PresRamCurr;
992  if( trace.lgTrace )
993  {
994  fprintf( ioQQQ,
995  " PresTotCurrent 1st zn reset PresTotlInit to PresTotlCurr=%.3e "
996  "PresRamInit to PresRamCurr=%.3e old tot=%.3e old ram %.3e hden=%.3e\n",
998  pressure.PresRamInit,
999  PresTotlInitSave,PresRamInitSave,
1001  }
1002  }
1003 # endif
1004  return;
1005 }

Generated for cloudy by doxygen 1.8.1.1