cloudy trunk
Loading...
Searching...
No Matches
mole_eval_balance.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/*CO_step fills in matrix for heavy elements molecular routines */
4#include "cddefines.h"
5#include "mole.h"
6#include "mole_priv.h"
7#include "save.h"
8#include "dense.h"
9#include "atmdat.h"
10#include "trace.h"
11/* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element
12 * molecular network in Cloudy. Before this routine would predict negative abundances if
13 * the fraction of carbon in the form of molecules came close to 100%. A reorganizing of
14 * the reaction network detected several bugs. Treatment of "coupled reactions",
15 * in which both densities in the reaction rate were being predicted by Cloudy, were also
16 * added. Due to these improvements, Cloudy can now perform calculations
17 * where 100% of the carbon is in the form of CO without predicting negative abundances
18 *
19 * Additional changes were made in November of 2003 so that our reaction
20 * network would include all reactions from the TH85 paper. This involved
21 * adding silicon to the chemical network. Also the reaction rates were
22 * labeled to make identification with the reaction easier and the matrix
23 * elements of atomic C, O, and Si are now done in a loop, which makes
24 * the addition of future chemical species (like N or S) easy.
25 * */
26/* Robin Williams in August 2006 onwards reorganized the coding to cut down repetitions.
27 * This isolated several further bugs, and allows a sigificant number of lines of
28 * code to be eliminated. The balance of S2/S2+ amd ClO/ClO+ seems highly sensitive
29 * (with small log scale results varying significantly if the order of arithmetic
30 * operations is changed) -- I suspect this may imply a bug somewhere.
31 * */
32/*lint -e778 constant expression evaluatess to 0 in operation '-' */
33/*=================================================================*/
34
36
37void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c)
38{
39 long int i, j;
40 mole_reaction *rate;
41 double rate_tot, rate_deriv[MAXREACTANTS], rk;
42 molecule *sp;
43
44 DEBUG_ENTRY("mole_eval_balance()");
45 /* zero out array used for formation rates */
46 for( i=0; i < num_total; i++ )
47 {
48 b[i] = 0.;
49 }
50 if (lgJac)
51 c.zero();
52
53 if (trace.lgTrace)
54 {
55 for(long mol=0;mol<mole_global.num_calc;mol++)
56 {
57 fprintf(ioQQQ," TrChemSp %20.13g %s\n",mole.species[mol].den,
58 mole_global.list[mol]->label.c_str());
59 }
61 p != mole_priv::reactab.end(); ++p)
62 {
63 rate = &(*p->second);
64 rk = mole.reaction_rks[ rate->index ];
65 fprintf(ioQQQ," TrChem %20.13g %s\n",rk,rate->label.c_str());
66 }
67 }
68
69
71 p != mole_priv::reactab.end(); ++p)
72 {
73 rate = &(*p->second);
74 rk = mole.reaction_rks[ rate->index ];
75
76 rate_tot = rk;
77 for(i=0;i<rate->nreactants;i++)
78 {
79 rate_tot *= mole.species[ rate->reactants[i]->index ].den;
80 }
81
82 if (trace.lgTrace)
83 {
84 fprintf(ioQQQ," TrChemTot %20.13g %20.13g %s\n",rate_tot,rk,rate->label.c_str());
85 for(i=0;i<rate->nreactants;i++)
86 {
87 fprintf(ioQQQ, " %20.13g", mole.species[ rate->reactants[i]->index ].den);
88 }
89 fprintf(ioQQQ,"\n");
90 }
91 for(i=0;i<rate->nreactants;i++)
92 {
93 sp = rate->reactants[i];
94 if (rate->rvector[i] == NULL)
95 {
96 b[sp->index] -= rate_tot;
97 }
98 }
99
100 for(i=0;i<rate->nproducts;i++)
101 {
102 sp = rate->products[i];
103 if (rate->pvector[i] == NULL)
104 {
105 b[sp->index] += rate_tot;
106 }
107 }
108
109 if (lgJac)
110 {
111 for(i=0;i<rate->nreactants;i++)
112 {
113 rate_deriv[i] = rk;
114 for(j=0;j<rate->nreactants;j++)
115 {
116 if(i!=j)
117 {
118 rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
119 }
120 }
121 }
122 for(j=0;j<rate->nreactants;j++)
123 {
124 sp = rate->reactants[j];
125 const double rated = rate_deriv[j];
126 for(i=0;i<rate->nreactants;i++)
127 {
128 if (rate->rvector[i] == NULL)
129 c[sp->index][rate->reactants[i]->index] -= rated;
130 }
131 for(i=0;i<rate->nproducts;i++)
132 {
133 if (rate->pvector[i] == NULL)
134 c[sp->index][rate->products[i]->index] += rated;
135 }
136 }
137 }
138 }
139
140 if (lgJac)
141 {
143 }
144
145 //mole_dominant_rates(findspecies("H+"),ioQQQ);
146 return;
147}
148
149void mole_eval_sources(long int num_total)
150{
151 long int i, j, nelem, ion, ion2;
152 mole_reaction *rate;
153 double rate_tot, rate_deriv[MAXREACTANTS], rk;
154 molecule *sp;
155
156 DEBUG_ENTRY("mole_eval_sources()");
157 /* zero out array used for formation rates */
158 for( i=0; i < num_total; i++ )
159 {
160 mole.species[i].src = mole.species[i].snk = 0.;
161 }
162
163 for( nelem=0; nelem< LIMELM; ++nelem )
164 {
165 /* these have one more ion than above */
166 for( ion=0; ion<nelem+2; ++ion )
167 {
168 /* zero out the transfer array */
169 for( ion2=0; ion2<nelem+2; ++ion2 )
170 {
171 mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
172 }
173 }
174 }
175
176 for(mole_reaction_i p=mole_priv::reactab.begin();
177 p != mole_priv::reactab.end(); ++p)
178 {
179 rate = &(*p->second);
180 rk = mole.reaction_rks[ rate->index ];
181
182 for(i=0;i<rate->nreactants;i++)
183 {
184 rate_deriv[i] = rk;
185 for(j=0;j<rate->nreactants;j++)
186 {
187 if(i!=j)
188 {
189 rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
190 }
191 }
192 }
193
194 rate_tot = rate_deriv[0] * mole.species[ rate->reactants[0]->index ].den;
195
196 for(i=0;i<rate->nreactants;i++)
197 {
198 sp = rate->reactants[i];
199 if (rate->rvector[i] == NULL)
200 {
201 mole.species[sp->index].snk += rate_deriv[i];
202 }
203 else
204 {
205 // exclude charge exchange with isotopes of same element
206 if( rate->nreactants==2 )
207 {
208 long otherIndex = 1-i;
209 molecule *sp2 = rate->reactants[otherIndex];
210 if( sp->isMonatomic() && sp2->isMonatomic() && sp->nAtom.begin()->first->el == sp2->nAtom.begin()->first->el )
211 continue;
212 }
213
214 if ( atmdat.lgCTOn )
215 {
216 for( ChemAtomList::iterator atom = unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
217 {
218 nelem = (*atom)->el->Z-1;
219 if (sp->nAtom.find(*atom) != sp->nAtom.end() && sp->nAtom[*atom] != 0 && rate->rvector[i]->charge != sp->charge)
220 {
221 mole.xMoleChTrRate[nelem][sp->charge][rate->rvector[i]->charge] +=
222 (realnum) rate_deriv[i];
223 break;
224 }
225 }
226 }
227 }
228 }
229
230 for(i=0;i<rate->nproducts;i++)
231 {
232 sp = rate->products[i];
233 if (rate->pvector[i] == NULL)
234 {
235 mole.species[sp->index].src += rate_tot;
236 }
237 }
238 }
239
240 for (ChemAtomList::iterator atom = unresolved_atom_list.begin();
241 atom != unresolved_atom_list.end(); ++atom)
242 {
243 const long int nelem=(*atom)->el->Z-1;
244 if( !dense.lgElmtOn[nelem] )
245 continue;
246
247 for (long int ion=0;ion<nelem+2;ion++)
248 {
249 if ((*atom)->ipMl[ion] != -1)
250 {
251 mole.source[nelem][ion] = mole.species[(*atom)->ipMl[ion]].src;
252 mole.sink[nelem][ion] = mole.species[(*atom)->ipMl[ion]].snk;
253 }
254 else
255 {
256 mole.source[nelem][ion] = 0.0;
257 mole.sink[nelem][ion] = 0.0;
258 }
259 }
260 }
261
262 //mole_dominant_rates(findspecies("H+"),ioQQQ);
263 return;
264}
265
267{
268 DEBUG_ENTRY ("lgNucleiConserved()");
269 bool checkAllOK = true;
270 multi_arr<double,2> test(atom_list.size(),mole_global.num_calc),
271 tot(atom_list.size(),mole_global.num_calc);
272 vector<double> ccache(mole_global.num_calc);
273 vector<long> ncache(mole_global.num_calc);
274 test.zero();
275 tot.zero();
276
277 map<chem_atom*, long> atom_to_index;
278 for (unsigned long j=0; j<atom_list.size(); ++j )
279 {
280 atom_to_index[atom_list[j].get_ptr()] = j;
281 }
282 for (long j=0;j<mole_global.num_calc;j++)
283 {
284 long nc = 0;
285 for (long i=0;i<mole_global.num_calc;i++)
286 {
287 if (c[i][j] != 0.)
288 {
289 ccache[nc] = c[i][j];
290 ncache[nc] = i;
291 ++nc;
292 }
293 }
294 if (nc > 0)
295 {
296 for (molecule::nAtomsMap::const_iterator el = mole_global.list[j]->nAtom.begin();
297 el != mole_global.list[j]->nAtom.end(); ++el)
298 {
299 long natom = atom_to_index[el->first.get_ptr()];
300 const int nAtomj = el->second;
301 for (long i=0;i<nc;i++)
302 {
303 const double term = ccache[i] * nAtomj;
304 test[natom][ncache[i]] += term;
305 tot[natom][ncache[i]] += fabs(term);
306 }
307 }
308 }
309 }
310
311 for( unsigned long natom=0; natom < atom_list.size(); ++natom)
312 {
313 for (long i=0;i<mole_global.num_calc;i++)
314 {
315 const bool checkOK =
316 ( fabs(test[natom][i]) <= MAX2(1e-10*tot[natom][i], 1e10*DBL_MIN) );
317 if ( UNLIKELY(!checkOK) )
318 {
319 chem_atom *atom = atom_list[natom].get_ptr();
320 fprintf(stdout,"Network conservation error %s %s %g %g %g %g\n",
321 atom->label().c_str(),
322 mole_global.list[i]->label.c_str(),
323 test[natom][i],
324 test[natom][i]/tot[natom][i],
325 mole.species[atom->ipMl[0]].den,
326 mole.species[atom->ipMl[1]].den);
327 //fprintf(stdout,"Problem at %s\n",rate->label);
328 checkAllOK = false;
329 }
330 }
331 }
332 return checkAllOK;
333}
334
335void mole_dominant_rates( const molecule *debug_species, FILE *ioOut )
336{
337 long int i, j;
338 mole_reaction *rate, *ratesnk=NULL, *ratesrc=NULL;
339 double rate_tot, rate_deriv[MAXREACTANTS], rk;
340 molecule *sp;
341 double snkx=0.,srcx=0.;
342
344 p != mole_priv::reactab.end(); ++p)
345 {
346 rate = &(*p->second);
347 rk = mole.reaction_rks[ rate->index ];
348
349 for(i=0;i<rate->nreactants;i++)
350 {
351 rate_deriv[i] = rk;
352 for(j=0;j<rate->nreactants;j++)
353 {
354 if(i!=j)
355 {
356 rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
357 }
358 }
359 }
360
361 rate_tot = rate_deriv[0] * mole.species[ rate->reactants[0]->index ].den;
362
363 if (debug_species != null_mole)
364 {
365 for(i=0;i<rate->nproducts;++i)
366 {
367 sp = rate->products[i];
368 if (sp == debug_species && rate->pvector[i] == NULL)
369 {
370 if (fabs(rate_tot) > srcx)
371 {
372 srcx = rate_tot;
373 ratesrc = rate;
374 }
375 }
376 }
377 for(i=0;i<rate->nreactants;++i)
378 {
379 sp = rate->reactants[i];
380 if (sp == debug_species && rate->rvector[i] == NULL)
381 {
382 if (fabs(rate_deriv[i]) > snkx)
383 {
384 snkx = rate_deriv[i];
385 ratesnk = rate;
386 }
387 }
388 }
389 }
390 }
391
392
393 if (debug_species != null_mole)
394 {
395 if (ratesrc)
396 {
397 fprintf( ioOut, "%20.20s src %13.7g of %13.7g [",
398 ratesrc->label.c_str(),srcx,mole.species[debug_species->index].src);
399 for (j=0;j<ratesrc->nreactants;j++)
400 {
401 if (j)
402 {
403 fprintf( ioOut, "," );
404 }
405 fprintf( ioOut, "%-6.6s %13.7g",
406 ratesrc->reactants[j]->label.c_str(),
407 mole.species[ ratesrc->reactants[j]->index ].den);
408 }
409 fprintf( ioOut, "]" );
410 }
411 if (ratesnk)
412 {
413 fprintf( ioOut, "%20.20s snk %13.7g of %13.7g [",
414 ratesnk->label.c_str(), snkx * mole.species[debug_species->index].den,
415 mole.species[debug_species->index].snk * mole.species[debug_species->index].den);
416 for (j=0;j<ratesnk->nreactants;j++)
417 {
418 if (j)
419 {
420 fprintf( ioOut, "," );
421 }
422 fprintf( ioOut, "%-6.6s %13.7g",
423 ratesnk->reactants[j]->label.c_str(),
424 mole.species[ ratesnk->reactants[j]->index ].den);
425 }
426 fprintf( ioOut, "]" );
427 }
428 }
429 fprintf( ioOut, "\n" );
430
431 return;
432}
t_atmdat atmdat
Definition atmdat.cpp:6
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int LIMELM
Definition cddefines.h:258
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
string label(void) const
Definition mole.h:55
vector< int > ipMl
Definition mole.h:45
molecule * pvector[MAXPRODUCTS]
Definition mole_priv.h:57
molecule * reactants[MAXREACTANTS]
Definition mole_priv.h:53
molecule * rvector[MAXREACTANTS]
Definition mole_priv.h:54
string label
Definition mole_priv.h:51
molecule * products[MAXPRODUCTS]
Definition mole_priv.h:56
int charge
Definition mole.h:144
nAtomsMap nAtom
Definition mole.h:143
string label
Definition mole.h:142
bool isMonatomic(void) const
Definition mole.h:157
int index
Definition mole.h:169
#define UNLIKELY(x)
Definition cpu.h:387
t_dense dense
Definition dense.cpp:24
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
ChemAtomList atom_list
ChemAtomList unresolved_atom_list
molecule * null_mole
void mole_eval_sources(long int num_total)
void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
void mole_dominant_rates(const molecule *debug_species, FILE *ioOut)
STATIC bool lgNucleiConserved(const multi_arr< double, 2 > &c)
#define MAXREACTANTS
Definition mole_priv.h:45
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
Definition mole_priv.h:37
map< string, count_ptr< mole_reaction > > reactab
t_trace trace
Definition trace.cpp:5