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