cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
dense.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cddefines.h"
4 
5 #include "dense.h"
6 
7 #include "abund.h"
8 #include "colden.h"
9 #include "conv.h"
10 #include "dynamics.h"
11 #include "elementnames.h"
12 #include "deuterium.h"
13 #include "hmi.h"
14 #include "phycon.h"
15 #include "radius.h"
16 #include "struc.h"
17 #include "thermal.h"
18 #include "trace.h"
19 #include "iso.h"
20 #include "h2.h"
21 #include "mole.h"
22 #include "save.h"
23 
25 {
27 }
28 
30 {
31  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
32  {
33  dense.SetGasPhaseDensity( nelem, 0. );
34  for( long ion=0; ion < LIMELM+1; ion++ )
35  {
36  dense.xIonDense[nelem][ion] = 0.;
37  }
38  }
39  for (long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem)
40  m_xMolecules[nelem] = 0.f;
41  // needed for TempChange to work but arrays needed for EdenChange to
42  // work are not yet defined
43  dense.eden = 1.;
44 
45  dense.xMassTotal = 0.;
46  dense.xNucleiTotal = 1.;
47  /* WJH */
48  dense.xMassDensity0 = -1.0f;
49 
50  /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the
51  * electron density to account for collisions with atomic H. it is an order of
52  * magnitude guess, so this command provides ability to see whether it affects results */
53  dense.HCorrFac = 1.f;
54 }
55 
57 {
58  double edensave = dense.eden;
59 
60  for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
61  {
62  if (dense.lgElmtOn[nelem])
63  {
64  ScaleIonDensities( nelem, factor );
65  dense.SetGasPhaseDensity( nelem, dense.gas_phase[nelem] * factor);
66  }
67  }
68 
69  EdenChange( dense.eden * factor );
70 
71  if( trace.lgTrace && trace.lgNeBug )
72  {
73  fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n",
74  edensave, dense.eden, edensave/dense.eden );
75  }
76 
77  hmi.H2_total *= factor;
78  h2.ortho_density *= factor;
79  h2.para_density *= factor;
80  for( long mol=0; mol < mole_global.num_total; mol++ )
81  {
82  mole.species[mol].den *= factor;
83  }
86 
88 }
89 
90 void ScaleIonDensities( const long nelem, const realnum factor )
91 {
92  double renorm;
93  for( long ion=0; ion < (nelem + 2); ion++ )
94  {
95  dense.xIonDense[nelem][ion] *= factor;
96  if (nelem-ion >= 0 && nelem-ion < NISO)
97  iso_renorm(nelem, nelem-ion, renorm);
98  }
99 
100  if( nelem==ipHYDROGEN && deut.lgElmtOn )
101  ScaleDensitiesDeuterium( factor );
102 
103  return;
104 }
105 
106 void t_dense::SetGasPhaseDensity( const long nelem, const realnum density )
107 {
108  gas_phase[nelem] = density;
109 
110  if( nelem==ipHYDROGEN && deut.lgElmtOn )
111  {
112  // fractionation is taken care of inside called routine
113  SetGasPhaseDeuterium( density );
114  }
115 
116  return;
117 }
118 
119 bool lgElemsConserved (void)
120 {
121  bool lgOK=true;
122 
123  /* this checks all molecules and ions add up to abundance */
124  for( ChemNuclideList::iterator nuc = nuclide_list.begin(); nuc != nuclide_list.end(); ++nuc)
125  {
126  long nelem = (*nuc)->el->Z-1;
127 
128  double sum_monatomic = 0.;
129  realnum sum_molecules = 0.f;
130  realnum sum_gas_phase = 0.f;
131 
132  if( (*nuc)->label() == "D" )
133  {
134  sum_monatomic = deut.xIonDense[0] + deut.xIonDense[1];
135  sum_molecules = deut.xMolecules();
136  sum_gas_phase = deut.gas_phase;
137 
138  }
139  /* check that species add up if element is turned on */
140  else if( dense.lgElmtOn[nelem] )
141  {
142  /* this sum of over the molecules did not include the atom and first
143  * ion that is calculated in the molecular solver */
144  for( long ion=0; ion<nelem+2; ++ion )
145  {
146  sum_monatomic += dense.xIonDense[nelem][ion] * (*nuc)->frac;
147  }
148  sum_molecules = dense.xMolecules(nelem) * (*nuc)->frac;
149  sum_gas_phase = dense.gas_phase[nelem] * (*nuc)->frac;
150  }
151  else
152  continue;
153 
154 
155  if( sum_monatomic + sum_molecules <= SMALLFLOAT &&
156  sum_gas_phase > SMALLFLOAT)
157  {
158  fprintf(ioQQQ,"PROBLEM non-conservation of nuclei %s\tions %g moles %g error %g of %g\n",
159  (*nuc)->label().c_str(),
160  sum_monatomic,
161  sum_molecules,
162  sum_monatomic + sum_molecules-sum_gas_phase,
163  sum_gas_phase );
164  lgOK=false;
165  }
166 
167  /* check that sum agrees with total gas phase abundance */
171  if( fabs( sum_monatomic + sum_molecules - sum_gas_phase ) >
172  conv.GasPhaseAbundErrorAllowed * sum_gas_phase )
173  {
174  fprintf(ioQQQ,"PROBLEM non-conservation of nuclei %s\t nzone %li atoms %.12e moles %.12e sum %.12e tot gas %.12e rel err %.3e\n",
175  (*nuc)->label().c_str(),
176  nzone,
177  sum_monatomic,
178  sum_molecules,
179  sum_monatomic + sum_molecules,
180  sum_gas_phase,
181  (sum_monatomic + sum_molecules - sum_gas_phase)/sum_gas_phase );
182  lgOK=false;
183  }
184 
185  if( 0 )
186  {
187  // print reactions involving the neutral
189  }
190  }
191 
192  return lgOK;
193 }
194 
195 void lgStatesConserved ( long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion )
196 {
197  if( 0 && conv.lgConvIoniz() &&
198  dense.lgElmtOn[nelem] &&
199  dense.xIonDense[nelem][ionStage]/dense.eden > conv.EdenErrorAllowed/10. &&
200  dense.xIonDense[nelem][ionStage] > SMALLFLOAT )
201  {
202  double abund = 0.;
203  for (long ipLev=0; ipLev<numStates; ipLev++)
204  {
205  abund += states[ipLev].Pop();
206  }
207 
208  double rel_err = fabs(abund-dense.xIonDense[nelem][ionStage])/
209  SDIV(dense.xIonDense[nelem][ionStage]);
210 
211  //fprintf(ioQQQ,"Iso consistency error %ld %ld %g\n",nelem,ionStage,rel_err);
212 
213  if( rel_err > err_tol )
214  {
215  /* we'll set the flag for non-convergence, but don't bother complaining in printout for first sweeps */
216  if( 0 && nzone > 0 && loop_ion > 0 )
217  {
218  fprintf(ioQQQ,"PROBLEM Inconsistent states/stage pops nzone %3ld loop_ion %2ld nelem %2ld ion %2ld states = %e stage = %e error = %e\n",
219  nzone,
220  loop_ion,
221  nelem,
222  ionStage,
223  abund,
224  dense.xIonDense[nelem][ionStage],
225  rel_err );
226  }
227  char chConvIoniz[INPUT_LINE_LENGTH];
228  sprintf( chConvIoniz , "States!=stage pops nelem %ld ion %ld ", nelem, ionStage );
229  conv.setConvIonizFail( chConvIoniz, abund,
230  dense.xIonDense[nelem][ionStage] );
231  }
232  }
233 }
234 
235 void SumDensities( void )
236 {
237  realnum den_mole,
238  DenseAtomsIons;
239 
240  /* find total number of atoms and ions
241  * at this point does not include molecules */
242  DenseAtomsIons = 0.;
243  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
244  {
245  if( dense.lgElmtOn[nelem] )
246  {
247  for( long ion=0; ion<=nelem+1; ++ion )
248  DenseAtomsIons += dense.xIonDense[nelem][ion];
249  }
250  }
251 
252  /* find number of molecules of the heavy elements in the gas phase.
253  * Atoms and ions were counted above. Do not include ices, solids
254  * on grains */
255  den_mole = total_molecules_gasphase();
256 
257  /* total number of atoms, ions, and molecules in gas phase per unit vol */
258  dense.xNucleiTotal = den_mole + DenseAtomsIons;
259  if( dense.xNucleiTotal > BIGFLOAT )
260  {
261  fprintf(ioQQQ, "PROBLEM DISASTER SumDensities has found "
262  "dense.xNucleiTotal with an insane density.\n");
263  fprintf(ioQQQ, "The density was %.2e\n",
265  TotalInsanity();
266  }
267  ASSERT( dense.xNucleiTotal > 0. );
268 
269  /* particle density that enters into the pressure includes electrons
270  * cm-3 */
272 
273  /* dense.wmole is molecular weight, AtomicMassUnit per particle */
274  dense.wmole = 0.;
275  for( long i=0; i < LIMELM; i++ )
276  {
278  }
279  /* dense.wmole is now molecular weight, AtomicMassUnit per particle */
280  dense.wmole /= dense.pden;
281  ASSERT( dense.wmole > 0. && dense.pden > 0. );
282 
283  /* xMassDensity is density in grams per cc */
286  dense.xMassDensity = (realnum)(dense.wmole*ATOMIC_MASS_UNIT*dense.pden);
287 
288  /*>>chng 04 may 25, introduce following test on xMassDensity0
289  * WJH 21 may 04, this is the mass density that corresponds to the hden
290  * specified in the init file. It is used to calculate the mass flux in the
291  * dynamics models. It may not necessarily be the same as struc.DenMass[0],
292  * since the pressure solver may have jumped to a different density at the
293  * illuminated face from that specified.*/
294  if( dense.xMassDensity0 < 0.0 )
296 
297  return;
298 }
299 
300 bool AbundChange( )
301 {
302  DEBUG_ENTRY( "AbundChange()" );
303 
304  fixit("AbundChange breaks conservation if molecular abundance is finite");
305 
306  // Suggested improved procedure:
307  // 1. Scale all molecules by the scaling for their most restrictive constituent
308  // 2. dense.updateXMolecules();
309  // 3. Scale ionic species so as to maintain conservation
310  // This will probably lead to material being dumped into the ionic species,
311  // but this should slosh back at the next molecular network solve. This
312  // procedure is at least a) conservative; b) consistent with the pure-ion case.
313 
314  /* this will be set true if density or abundances change in this zone */
315  bool lgChange = false;
316 
317  /* >>chng 97 jun 03, added variable abundances for element table command
318  * option to get abundances off a table with element read command */
319  if( abund.lgAbTaON )
320  {
321  lgChange = true;
322  for( long nelem=1; nelem < LIMELM; nelem++ )
323  {
324  if( abund.lgAbunTabl[nelem] )
325  {
326  double abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
328 
329  double hold = abun/dense.gas_phase[nelem];
330  dense.SetGasPhaseDensity( nelem, (realnum)abun );
331 
332  for( long ion=0; ion < (nelem + 2); ion++ )
333  {
334  dense.xIonDense[nelem][ion] *= (realnum)hold;
335  }
336  }
337  }
338  }
339 
340  /* this is set false if fluctuations abundances command entered,
341  * when density variations are in place this is true,
342  * and dense.chDenseLaw is "SINE" */
343  double FacAbun = 1.;
344  if( !dense.lgDenFlucOn )
345  {
346  static double FacAbunSav;
347  double OldAbun = 0.0;
348  /* abundance variations are in place */
349  lgChange = true;
350  if( nzone > 1 )
351  {
352  OldAbun = FacAbunSav;
353  }
354 
355  /* rapid abundances fluctuation */
356  if( dense.lgDenFlucRadius )
357  {
358  /* cycle over radius */
359  FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
361  }
362  else
363  {
364  /* cycle over column density */
365  FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
367  }
368 
369  if( nzone > 1 )
370  {
371  FacAbun = FacAbunSav/OldAbun;
372  }
373  }
374 
375  if( FacAbun != 1. )
376  {
377  ASSERT(!dynamics.lgAdvection); // Local abundance change doesn't make sense with advective transport
378 
379  /* Li on up is affect by abundance variations */
380  for( long nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
381  {
382  if (dense.lgElmtOn[nelem])
383  {
384  dense.SetGasPhaseDensity(nelem, dense.gas_phase[nelem] * (realnum) FacAbun);
385  ScaleIonDensities( nelem, (realnum)(FacAbun) );
386  }
387  }
388 
389  for( long mol=0; mol < mole_global.num_total; mol++ )
390  {
391  mole.species[mol].den *= (realnum)(FacAbun);
392  }
393  }
394 
395  if (lgChange)
396  {
397  /* must call TempChange since ionization has changed, there are some
398  * terms that affect collision rates (H0 term in electron collision) */
399  TempChange(phycon.te , false);
400  }
401 
402  return lgChange;
403 }
404 
405 namespace {
406  const int SCALE_NEW = 1;
407 }
408 
410 {
411  if (SCALE_NEW)
412  return dense.xMassDensity/(realnum)ATOMIC_MASS_UNIT;
413  else
414  return dense.gas_phase[ipHYDROGEN];
415 }
417 {
418  if (SCALE_NEW)
419  return struc.DenMass[i]/(realnum)ATOMIC_MASS_UNIT;
420  else
421  return struc.hden[i];
422 }
realnum * hden
Definition: struc.h:25
void ScaleIonDensities(const long nelem, const realnum factor)
Definition: dense.cpp:90
t_mole_global mole_global
Definition: mole.cpp:7
double Radius
Definition: radius.h:31
double depth
Definition: radius.h:31
void SetGasPhaseDensity(const long nelem, const realnum density)
Definition: dense.cpp:106
t_colden colden
Definition: colden.cpp:5
realnum xMolecules(long nelem)
Definition: dense.h:88
bool lgElmtOn
Definition: deuterium.h:21
realnum flcPhase
Definition: dense.h:265
double EdenErrorAllowed
Definition: conv.h:263
realnum csecnd
Definition: dense.h:264
int num_total
Definition: mole.h:362
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
t_struc struc
Definition: struc.cpp:6
const realnum SMALLFLOAT
Definition: cpu.h:246
realnum xNucleiTotal
Definition: dense.h:114
const int NISO
Definition: cddefines.h:311
realnum HCorrFac
Definition: dense.h:121
bool lgDenFlucRadius
Definition: dense.h:259
bool lgAdvection
Definition: dynamics.h:66
t_conv conv
Definition: conv.cpp:5
void updateXMolecules()
Definition: dense.cpp:24
t_phycon phycon
Definition: phycon.cpp:6
realnum GasPhaseAbundErrorAllowed
Definition: conv.h:281
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
bool lgConvIoniz() const
Definition: conv.h:108
ChemNuclideList nuclide_list
t_dynamics dynamics
Definition: dynamics.cpp:42
void lgStatesConserved(long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion)
Definition: dense.cpp:195
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:32
t_dense dense
Definition: global.cpp:15
bool lgNeBug
Definition: trace.h:112
realnum xMassDensity0
Definition: dense.h:105
t_elementnames elementnames
Definition: elementnames.cpp:5
realnum * DenMass
Definition: struc.h:25
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_trace trace
Definition: trace.cpp:5
double ortho_density
Definition: h2_priv.h:326
t_abund abund
Definition: abund.cpp:5
void ScaleAllDensities(realnum factor)
Definition: dense.cpp:56
realnum pden
Definition: dense.h:108
double para_density
Definition: h2_priv.h:326
bool lgTrace
Definition: trace.h:12
t_mole_local mole
Definition: mole.cpp:8
molecule * findspecies(const char buf[])
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
const realnum BIGFLOAT
Definition: cpu.h:244
realnum m_xMolecules[LIMELM]
Definition: dense.h:85
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
realnum flong
Definition: dense.h:262
realnum total_molecules_gasphase(void)
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
double AbundancesTable(double r0, double depth, long int iel)
Definition: abundances.cpp:464
bool AbundChange()
Definition: dense.cpp:300
void SetGasPhaseDeuterium(const realnum &Hdensity)
Definition: deuterium.cpp:65
realnum gas_phase
Definition: deuterium.h:22
realnum scalingZoneDensity(long i)
Definition: dense.cpp:416
t_radius radius
Definition: radius.cpp:5
double xIonDense[2]
Definition: deuterium.h:23
bool lgElmtOn[LIMELM]
Definition: dense.h:160
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
realnum gas_phase[LIMELM]
Definition: dense.h:76
realnum wmole
Definition: dense.h:111
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
void ScaleDensitiesDeuterium(const realnum &factor)
Definition: deuterium.cpp:21
#define ASSERT(exp)
Definition: cddefines.h:613
double density(const genericState &gs)
bool lgAbTaON
Definition: abund.h:218
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
realnum scalingDensity(void)
Definition: dense.cpp:409
realnum xMassDensity
Definition: dense.h:101
double H2_total
Definition: hmi.h:25
t_deuterium deut
Definition: deuterium.cpp:7
double eden
Definition: dense.h:201
void iso_renorm(long nelem, long ipISO, double &renorm)
Definition: iso_solve.cpp:310
realnum xMassTotal
Definition: dense.h:117
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
bool lgDenFlucOn
Definition: dense.h:255
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
void setConvIonizFail(const char *reason, double oldval, double newval)
Definition: conv.h:100
void zero()
Definition: dense.cpp:29
void SumDensities(void)
Definition: dense.cpp:235
realnum xMolecules(void)
Definition: deuterium.h:30
void EdenChange(double EdenNew)
Definition: eden_change.cpp:11
#define fixit(a)
Definition: cddefines.h:417
bool lgElemsConserved(void)
Definition: dense.cpp:119
t_hmi hmi
Definition: hmi.cpp:5
double te
Definition: phycon.h:21
realnum cfirst
Definition: dense.h:263
const int ipHYDROGEN
Definition: cddefines.h:349
const int ipLITHIUM
Definition: cddefines.h:351
void total_molecule_elems(realnum total[LIMELM])
realnum colden[NCOLD]
Definition: colden.h:32
bool lgAbunTabl[LIMELM]
Definition: abund.h:218
void mole_print_species_reactions(molecule *speciesToPrint)
void updateXMolecules()
Definition: deuterium.cpp:16