cloudy  trunk
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 
35 STATIC bool lgNucleiConserved(const multi_arr<double,2> &c);
36 
37 void 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  }
60  for(mole_reaction_i p=mole_priv::reactab.begin();
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 
70  for(mole_reaction_i p=mole_priv::reactab.begin();
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  {
142  ASSERT( lgNucleiConserved(c) );
143  }
144 
145  //mole_dominant_rates(findspecies("H+"),ioQQQ);
146  return;
147 }
148 
149 void 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 
266 STATIC bool lgNucleiConserved(const multi_arr<double,2> &c)
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 
335 void 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 
343  for(mole_reaction_i p=mole_priv::reactab.begin();
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 }
UNLIKELY
#define UNLIKELY(x)
Definition: cpu.h:387
mole_reaction::nreactants
int nreactants
Definition: mole_priv.h:52
dense
t_dense dense
Definition: dense.cpp:24
mole_eval_balance
void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
Definition: mole_eval_balance.cpp:37
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
molecule::isMonatomic
bool isMonatomic(void) const
Definition: mole.h:157
realnum
float realnum
Definition: cddefines.h:103
mole_dominant_rates
void mole_dominant_rates(const molecule *debug_species, FILE *ioOut)
Definition: mole_eval_balance.cpp:335
STATIC
#define STATIC
Definition: cddefines.h:97
unresolved_atom_list
ChemAtomList unresolved_atom_list
Definition: mole_species.cpp:70
mole.h
mole_reaction::pvector
molecule * pvector[MAXPRODUCTS]
Definition: mole_priv.h:57
molecule::index
int index
Definition: mole.h:169
multi_arr< double, 2 >
mole_reaction
Definition: mole_priv.h:49
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
trace.h
t_mole_global::list
MoleculeList list
Definition: mole.h:317
molecule::charge
int charge
Definition: mole.h:144
MAXREACTANTS
#define MAXREACTANTS
Definition: mole_priv.h:45
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
atmdat.h
t_mole_local::sink
double ** sink
Definition: mole.h:394
mole_priv::reactab
map< string, count_ptr< mole_reaction > > reactab
Definition: mole_species.cpp:60
dense.h
mole
t_mole_local mole
Definition: mole.cpp:7
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
null_mole
molecule * null_mole
Definition: mole_species.cpp:64
t_mole_local::xMoleChTrRate
realnum *** xMoleChTrRate
Definition: mole.h:396
mole_reaction::index
long index
Definition: mole_priv.h:62
mole_reaction::nproducts
int nproducts
Definition: mole_priv.h:52
chem_atom::ipMl
vector< int > ipMl
Definition: mole.h:45
MAX2
#define MAX2
Definition: cddefines.h:782
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_mole_local::source
double ** source
Definition: mole.h:394
save.h
chem_atom::label
string label(void) const
Definition: mole.h:55
mole_reaction::reactants
molecule * reactants[MAXREACTANTS]
Definition: mole_priv.h:53
mole_reaction::products
molecule * products[MAXPRODUCTS]
Definition: mole_priv.h:56
mole_reaction::label
string label
Definition: mole_priv.h:51
multi_arr::zero
void zero()
Definition: container_classes.h:1051
chem_atom
Definition: mole.h:37
molecule
Definition: mole.h:132
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
molecule::nAtom
nAtomsMap nAtom
Definition: mole.h:143
mole_reaction_i
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
Definition: mole_priv.h:37
lgNucleiConserved
STATIC bool lgNucleiConserved(const multi_arr< double, 2 > &c)
Definition: mole_eval_balance.cpp:266
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
atom_list
ChemAtomList atom_list
Definition: mole_species.cpp:71
t_atmdat::lgCTOn
bool lgCTOn
Definition: atmdat.h:177
molecule::label
string label
Definition: mole.h:142
t_mole_local::reaction_rks
vector< double > reaction_rks
Definition: mole.h:400
mole_priv.h
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
mole_reaction::rvector
molecule * rvector[MAXREACTANTS]
Definition: mole_priv.h:54
mole_eval_sources
void mole_eval_sources(long int num_total)
Definition: mole_eval_balance.cpp:149
Generated by   doxygen 1.8.17