cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_co_solve.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 /*CO_step fills in matrix for heavy elements molecular routines */
4 #include "cddefines.h"
5 #include "taulines.h"
6 #include "dense.h"
7 #include "ionbal.h"
8 #include "thermal.h"
9 #include "phycon.h"
10 #include "hmi.h"
11 #include "dynamics.h"
12 #include "conv.h"
13 #include "trace.h"
14 #include "coolheavy.h"
15 #include "timesc.h"
16 #include "thirdparty.h"
17 #include "mole.h"
18 #include "mole_co_priv.h"
19 #include "mole_co_atom.h"
20 /* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element
21  * molecular network in Cloudy. Before this routine would predict negative abundances if
22  * the fraction of carbon in the form of molecules came close to 100%. A reorganizing of
23  * the reaction network detected several bugs. Treatment of "coupled reactions",
24  * in which both densities in the reaction rate were being predicted by Cloudy, were also
25  * added. Due to these improvements, Cloudy can now perform calculations
26  * where 100% of the carbon is in the form of CO without predicting negative abundances
27  *
28  * Additional changes were made in November of 2003 so that our reaction
29  * network would include all reactions from the TH85 paper. This involved
30  * adding silicon to the chemical network. Also the reaction rates were
31  * labeled to make identification with the reaction easier and the matrix
32  * elements of atomic C, O, and Si are now done in a loop, which makes
33  * the addition of future chemical species (like N or S) easy.
34  * */
35 
36 void CO_solve(
37  /* set true if we found neg pops */
38  bool *lgNegPop,
39  /* set true if we tried to compute the pops, but some were zero */
40  bool *lgZerPop )
41 {
42 
43  int32 merror;
44  long int i, j, k, n,
45  nelem , ion , ion2;
46  double
47  co_denominator;
48  realnum cartot_mol, oxytot_mol;
49  realnum abundan;
50  struct chem_element_s *element;
51  double sum;
52 
53  DEBUG_ENTRY( "CO_solve()" );
54 
55  CO_step(); /* Calculate the matrix elements */
56 
57  /* Ugly hack to deal with species which have become uncoupled */
58  for(i=0;i<mole.num_comole_calc;i++)
59  {
60  sum = 0.;
61  for(j=0;j<mole.num_comole_calc;j++)
62  {
63  sum = sum+fabs(mole.c[i][j]);
64  }
65  if(sum < SMALLFLOAT && fabs(mole.b[i]) < SMALLFLOAT)
66  {
67  mole.c[i][i] = 1.;
68  }
69  }
70 
71  cartot_mol = dense.xMolecules[ipCARBON] + findspecies("C")->hevmol + findspecies("C+")->hevmol;
72  oxytot_mol = dense.xMolecules[ipOXYGEN] + findspecies("O")->hevmol + findspecies("O+")->hevmol;
73 
74  ASSERT( cartot_mol >= 0. && oxytot_mol >= 0.);
75 
76  *lgNegPop = false;
77  *lgZerPop = false;
78 
79  /* zero out molecular charge transfer rates */
80  for(nelem=ipLITHIUM; nelem < LIMELM; ++nelem)
81  {
82  for(ion=0; ion < nelem+2; ++ion)
83  {
84  /*zero out the arrays */
85  mole.sink[nelem][ion] = 0.;
86  mole.source[nelem][ion] = 0.;
87 
88  for(ion2=0; ion2< nelem+2; ++ion2)
89  {
90  mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
91  }
92  }
93  }
94 
95 
96  /* >>06 jun 29, add advective terms here */
97  /* Add rate terms for dynamics to equilibrium, makes c[][] non-singular
98  * do before cross talk with heavy elements is done to not double count charge
99  * transfer */
100  if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection /*&& radius.depth < dynamics.oldFullDepth*/ )
101  {
102  for(i=0;i<mole.num_comole_calc;i++)
103  {
104  if(COmole[i]->n_nuclei > 1) {
105  mole.c[i][i] -= dynamics.Rate;
106  /* this is the net gain of the species due to advection -
107  * in this form we have the net gain as the difference between
108  * current species flowing out of this region, downstream, and
109  * new advected material coming into this region from upstream */
110  mole.b[i] -= (COmole[i]->hevmol*dynamics.Rate-dynamics.CO_molec[i]);
111  }
112  }
113  }
114 
115  /* >>chng Oct. 21, 2004 -- Terms that contribute to the ionization balance of C, O, S, Si, and N
116  will now be inserted directly into the ion solver. Before, we just corrected the recombination rate
117  by scaling saying that IONIZATION_RATE = RECOMBINATION_RATE * [n(X+)/n(X0)]. This code follows the logic
118  of hmole_step, written by Robin Williams. */
119 
120 
121  /* sink terms should include only terms that either form or destroy an atomic or singly ionized
122  element in the network, but not both. This means that terms that destroy X at the expense of
123  forming X+ can't be included. The rate of destruction of X or X+ is equal to the formation of
124  the inverse. We therefore add this rate to sink, which gets rid of this dependence */
125 
126  /* The following code is all generalized. After all the molecules, the total number being mole.num_heavy_molec,
127  there are 2*mole.num_elements remaining. The first set, equal to mole.num_elements, are the ionized elemental
128  species. The second set, also equal to mole.num_elements, are the atomic species. They are in order such that
129 
130  array number for X = array number for (X+) + mole.num_elements
131 
132  In the future, if one wants to add another element to the network, such as Na the macros must be changed.
133  Also, the array number for Na and Na+ must obey the equation above. If this is done, then the sources,
134  sinks, and molecular recombination terms are automatically added!!! */
135 
136  /* COmole[i]->nelem_hevmol[j] is just the dominant element for a given molecule. For an element this
137  string is just the element itself!! */
138 
139  /* >> chng 2006 Sep 02 rjrw: Change to using element properties rather than properties of elements as molecules, for ease of
140  generalization -- ordering of atoms and ions no longer hardwired */
141 
142  for(i=0;i<mole.num_elements; ++i)
143  {
144  element = chem_element[i];
145  mole.sink[element->ipCl][0] -= (mole.c[element->ipMl][element->ipMl] + mole.c[element->ipMl][element->ipMlP])*dense.lgElmtOn[element->ipCl];
146  mole.sink[element->ipCl][1] -= (mole.c[element->ipMlP][element->ipMlP] + mole.c[element->ipMlP][element->ipMl] )*dense.lgElmtOn[element->ipCl];
147  }
148 
149  /* source terms must be multiplied by the density of the corresponding matrix element */
150 
151  for(j=0;j < mole.num_elements; ++j)
152  {
153  element = chem_element[j];
154  for(i=0; i < mole.num_comole_calc; ++i)
155  {
156  mole.source[element->ipCl][1] += COmole[i]->hevmol*mole.c[i][element->ipMlP]*dense.lgElmtOn[element->ipCl];
157  mole.source[element->ipCl][0] += COmole[i]->hevmol*mole.c[i][element->ipMl]*dense.lgElmtOn[element->ipCl];
158  }
159  }
160 
161  /* subtract out diagonal terms, they are already in sink. Also subtract recombination terms,
162  as these are done by mole.xMoleChTrRate */
163 
164  for(i=0;i < mole.num_elements; ++i)
165  {
166  element = chem_element[i];
167  mole.source[element->ipCl][1] -= ( mole.c[element->ipMlP][element->ipMlP]*COmole[element->ipMlP]->hevmol +
168  mole.c[element->ipMl][element->ipMlP]*COmole[element->ipMl]->hevmol)*dense.lgElmtOn[element->ipCl];
169 
170  mole.source[element->ipCl][0] -= ( mole.c[element->ipMl][element->ipMl]*COmole[element->ipMl]->hevmol +
171  mole.c[element->ipMlP][element->ipMl]*COmole[element->ipMlP]->hevmol)*dense.lgElmtOn[element->ipCl];
172 
173  /* The source terms, as they are right now, include negative contributions. This is because the
174  linearization "trick" creates source terms. Take for example the reaction:
175  C+ H2O => HCO+ H
176  This reaction destroys C+, but in the act of linearizing, there will be the following term:
177  mole.c[ipH2O][ipCP]
178  This is only a numerical "trick" and not a real matrix term. All terms like this
179  have to be removed to get the "true" source terms. Fortunately, the sum of all these terms
180  is just mole.b[i]. To remove these terms, we have to add mole.b[i] if the overall contribution from these terms
181  was negative, subtract if their contribution was positive. This is done by subtracting mole.b[i]
182  from the current value of source. */
183 
184  mole.source[element->ipCl][1] = mole.source[element->ipCl][1] - mole.b[element->ipMlP];
185  mole.source[element->ipCl][0] = mole.source[element->ipCl][0] - mole.b[element->ipMl];
186 
187  }
188 
189  /* negative source terms are actually destruction mechanisms, so should go into the sink vector.
190  source and sinks have different units, so if source is negative divide by the density of
191  the corresponding element to get same units. For example, if:
192 
193  mole.source[ipCARBON][0]
194 
195  is negative, we divide by dense.xIonDense[ipCARBON][0] to get rate in same units as mole.sink */
196 
197  for(i=2; i < LIMELM; ++i)
198  {
199  for(j=0; j < 2; ++j)
200  {
201  /*if source term is negative, make it a sink and then set to zero*/
202  if(mole.source[i][j] < 0)
203  {
204  mole.sink[i][j] -= (mole.source[i][j]/SDIV(dense.xIonDense[i][j]));
205  mole.source[i][j] = 0;
206  }
207  }
208  }
209  /* These are rates of recombination for atomic and singly ionized species that are due to molecular processes
210  This will be added to ion_solver to get correct ionization balance */
211 
212  for(i=0;i < mole.num_elements; ++i)
213  {
214  element = chem_element[i];
215  mole.xMoleChTrRate[element->ipCl][1][0] = (realnum)(mole.c[element->ipMlP][element->ipMl] -
216  ionbal.RateRecomTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl];
217 
218  mole.xMoleChTrRate[element->ipCl][0][1] = (realnum)(mole.c[element->ipMl][element->ipMlP] -
219  ionbal.RateIonizTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl];
220  }
221 
222  /* If rate for going from 1-0 or 0-1 is negative then it should be added to the inverse process */
223  for(i=2; i < LIMELM; ++i)
224  {
225  for(j=0; j < 2; ++j)
226  {
227  for(k=0; k< 2; ++k)
228  {
229  /*only possible charge transfers are 0-1 or 1-0 */
230  if(j != k)
231  {
232  if( mole.xMoleChTrRate[i][j][k] < 0. )
233  {
234  mole.xMoleChTrRate[i][k][j] -= mole.xMoleChTrRate[i][j][k];
235  mole.xMoleChTrRate[i][j][k] = 0.;
236  }
237  }
238  }
239  }
240  }
241 
242  for(n=0;n<mole.num_elements;n++)
243  {
245  for( i=2; i < chem_element[n]->ipCl+2; i++ )
246  {
247  tot_ion[n] -= dense.xIonDense[chem_element[n]->ipCl][i];
248  }
249  }
250 
251 
252  /* at this point the cartot_mol, sum of molecules and atom/first ion,
253  * should be equal to the gas_phase minus double and higher ions */
254  /*fprintf(ioQQQ," dbuggggas\t %.2f\t%f\t%f\n",fnzone,
255  cartot_mol/cartot_ion,
256  oxytot_mol/oxytot_ion);*/
257 
258  /* these will be used in population equation in case of homogeneous equation */
259 
260 
261  /* rjrw 2006 Aug 08: messing with the matrix like this will likely break advection --
262  was done correctly in hmole */
263 
265  fixit();
266 
267  for(n=0;n<mole.num_elements;n++)
268  {
269  mole.b[chem_element[n]->ipMl] = tot_ion[n];
270  }
271 
272  /* <<chng 03 Nov 11--Nick Abel, Set up the non-zero matrix elements that go into the conservation equations
273  for atomic C, O, and Si, this is now set up by looping over the atomic species in
274  co.h and setting the number of atoms of C, O, and Si equal to the variable findspecies("CARBON")->nElem,
275  findspecies("OXYGEN")->nElem, and findspecies("SILICON")->nElem respectively. For every element (excluding Hydrogen) not in
276  the network an if statement will be needed */
277 
278  /* loop over last mole.num_elements elements in co vector - these are atoms of the mole.num_elements */
279  for( j=0; j < mole.num_comole_calc; j++ )
280  {
281  for(i=0;i<mole.num_elements;i++)
282  {
283  element = chem_element[i];
284  mole.c[j][element->ipMl] = COmole[j]->nElem[element->ipCl];
285  }
286  }
287 
288  /*--------------------------------------------------------------------
289  * */
290  /*printf( " Here are all matrix elements\n ");
291 
292  for(i=0; i < mole.num_comole_calc; i++)
293 
294  {
295 
296 
297  for(j=0; j < mole.num_comole_calc; j++)
298  {
299  printf( "%4.4s", COmole[i]->label );
300  printf( "%4.4s\t", COmole[j]->label );
301  printf( "%.8e,\n", mole.c[j][i] );
302  }
303  }
304  printf( "b's are:\n" );
305  for(i=0; i < mole.num_comole_calc; i++)
306  {
307  printf( "%.8e,\n", mole.b[i] );
308  }*/
309 
310 
312  {
313  fprintf( ioQQQ, " COMOLE matrix\n " );
314 
315 # if 0
316  for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
317  {
318  fprintf( ioQQQ, "%4.4s", COmole[i]->label );
319  }
320  fprintf( ioQQQ, " \n" );
321 
322  for( j=0; j < mole.num_comole_calc; j++ )
323  {
324  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
325  fprintf( ioQQQ, " " );
326  for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
327  {
328  fprintf( ioQQQ, "%9.1e", mole.c[i][j] );
329  }
330  fprintf( ioQQQ, "\n" );
331  }
332 # endif
333 
334  fprintf( ioQQQ, " COMOLE matrix\n " );
335 
336  for( i=0; i < mole.num_comole_calc; i++ )
337  {
338  fprintf( ioQQQ, "%4.4s", COmole[i]->label );
339  }
340  fprintf( ioQQQ, " \n" );
341 
342  for( j=0; j < mole.num_comole_calc; j++ )
343  {
344  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
345  fprintf( ioQQQ, " " );
346  for( i=0; i < mole.num_comole_calc; i++ )
347  {
348  /*fprintf( ioQQQ, "%9.1e", mole.c[i][j] );*/
349  fprintf( ioQQQ, ", %.15e", mole.c[i][j] );
350  }
351  fprintf( ioQQQ, "\n" );
352  }
353 
354  fprintf( ioQQQ, " COMOLE b\n " );
355  for( j=0; j < mole.num_comole_calc; j++ )
356  {
357  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
358  fprintf( ioQQQ, ", %.15e\n",mole.b[j] );
359  }
360 
361 #if 0
362  fprintf( ioQQQ,
363  " COMOLE H2 den:%10.3e, H2,3+=%10.2e%10.2e Carb sum=%10.3e Oxy sum=%10.3e\n",
364  hmi.H2_total,
365  hmi.Hmolec[ipMH2p],
366  hmi.Hmolec[ipMH3p],
367  mole.c[mole.num_heavy_molec][findspecies("C")->index],
368  mole.c[mole.num_heavy_molec][findspecies("O")->index] );
369 
370 #endif
371  }
372 
373  /* remember longest molecule destruction timescales */
374  for( j=0; j < mole.num_comole_calc; j++ )
375  {
376  if( -mole.c[j][j] > SMALLFLOAT )
377  {
378  /* this is rate CO is destroyed, equal to formation rate in equilibrium */
379  timesc.AgeCOMoleDest[j] = -1./mole.c[j][j];
380  /* moved to radinc */
381  /*timesc.BigCOMoleForm = (realnum)MAX2(timesc.BigCOMoleForm,timesc.AgeCOMoleForm);*/
382  }
383  else
384  {
385  timesc.AgeCOMoleDest[j] = 0.;
386  }
387  }
388 
389  /* copy to amat, saving c for later use */
390  for( j=0; j < mole.num_comole_calc; j++ )
391  {
392  for( i=0; i < mole.num_comole_calc; i++ )
393  {
394  mole.amat[i][j] = mole.c[i][j];
395  }
396  }
397 
398  merror = 0;
399 
401 
402  if( merror != 0 )
403  {
404  fprintf( ioQQQ, " CO_solve getrf_wrapper error\n" );
405  cdEXIT(EXIT_FAILURE);
406  }
407 
409 
410  if( merror != 0 )
411  {
412  fprintf( ioQQQ, " CO_solve: dgetrs finds singular or ill-conditioned matrix\n" );
413  cdEXIT(EXIT_FAILURE);
414  }
415 
416  /* check for negative populations, which happens when 100% co */
417  *lgNegPop = false;
418  *lgZerPop = false;
419  co_denominator = 0;
420 # if 0
421  if(fnzone > 1.02)
422  {
423  for( i=0; i < mole.num_comole_calc; i++)
424  {
425 
426  for( j=61; j < 62; j++ )
427  {
428 
429  printf("ZONE\t%.5e\tSPECIES_1\t%s\tSPECIES_2\t%s\tDEST_RATE\t%.3e\tbvec\t%.3e\n",
430  fnzone, COmole[i]->label,COmole[j]->label,mole.c[i][j],mole.b[i]);
431 
432  }
433  }
434  }
435 # endif
436  for( i=0; i < mole.num_comole_calc; i++ )
437  {
438  /* fprintf(stderr,"%ld [%d] %g\n",i,mole.num_comole_calc,mole.b[i]); */
439  if( mole.b[i] < 0. )
440  {
441 
442  /* >>chng 03 sep 11, the following*/
443  /* comparison between the solution vector from 32 bit machines vs
444  * maple on a pc shows that very small negative numbers are produced by
445  * the linear algebra package used by the code, but is positive with
446  * the math package. allow very small negative numbers to occur,
447  * and silently reset to a postive number.
448  *
449  * >>chng 04 feb 25
450  * Here we want to check if the negative abundance is a significant
451  * fraction of any of the species in the network. The code below
452  * checks to see which elements the molecule in question is made of,
453  * then finds which one has the least abundance. The one with
454  * the least abundance will then be used as the divisor in checking
455  * if the negative abundance is too large to be ignored.*/
456 
457  /* >>chng 04 apr 21, when an element in the chemical network is
458  * turned off, the molecular abundance of a species containing that
459  * element was not zero but rather a small number of order 1e-20. When
460  * these abundances went negative, the code thought the abundances were
461  * important because it checks the ratio of the species to its gas phase
462  * abundance. When the gas phase abundance is zero, the denominator is set
463  * to SMALLFLOAT, which is ~1e-35. This led to a ratio of molecule to gas
464  * phase of ~1e15, clearly unphysical. Here the value of co_denominator
465  * will be set to a high value if the gas phase abundance is zero.*/
466 
467  if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
468  {
469  co_denominator = dense.gas_phase[COmole[i]->nelem_hevmol];
470  }
471  else
472  {
473  /* >>chng 04 apr 20, set to zero if element is off */
474  co_denominator = 1e10;
475  }
476 
477 
478  /* we must have a positive value by now, unlesss element is turned off */
479  ASSERT( co_denominator > SMALLFLOAT || !dense.lgElmtOn[COmole[i]->nelem_hevmol] );
480 
481  /* >>chng 04 feb 28, threshold from -1e-10 to -1e-6, the
482  * level of roundoff in a realnum */
483  /*if( mole.b[i] / MAX2 (co_denominator, SMALLFLOAT) > -1e-10) */
484  /*if( mole.b[i] / SDIV(co_denominator) > -1e-6 ) */
485  /* >>chng 04 oct 31, change limit from -1e-6 to -1e-5, produced only
486  * a few comments in func_map.in */
487  if( mole.b[i] / SDIV(co_denominator) > -1e-5 )
488  {
489  /* this case, it is only slightly negative relative to
490  * the parent species - multiply by -1 to make positive
491  * and press on */
492  mole.b[i] *= -1.;
493  }
494  else
495  {
496  /*>>chng 04 mar 06 press on */
497  *lgNegPop = true;
498  /* 05 jul 16, there had been a bug such that CO was always evaluated as long as
499  * the temperature was below 20,000K. Once fixed, one sim turned on CO mid way
500  * into the H+ zone and had neg pops during first trys - change check on whether it
501  * is to be commented on only if this is not first zone with CO soln */
502  /*if( nzone>0 )*/
503  if( nzone>co.co_nzone )
504  {
505  static long int nzFail=-2;
506  long int nzFail2 = (long)fnzone*100;
507  /* during a map we may be in search phase but not in first zone */
508  if( !conv.lgSearch )
509  fprintf(ioQQQ," PROBLEM ");
510  fprintf(ioQQQ,
511  " CO_solve neg pop for species %li %s, value is %.2e rel value is %.2e zone %.2f Te %.4e Search?%c\n",
512  i ,
513  COmole[i]->label ,
514  /* abs value */
515  mole.b[i],
516  /* relative value */
517  mole.b[i] / SDIV(co_denominator) ,
518  fnzone,
519  phycon.te,TorF( conv.lgSearch ) );
520  /* if well beyond search phase and still having problems, announce them
521  * only announce one failure per sweep through solvers by ConvFail */
522  if( nzFail2 !=nzFail )
523  {
524  nzFail = nzFail2;
525  ConvFail( "pops" , "CO");
526  }
527  }
528  mole.b[i] *= -1;
529  /*>>chng 04 jan 24 set to zero instead of *-1,
530  * h2_cr.in co density went to infinite, with some negative
531  * var getting to - inf
532  *>>chng 04 nov 03, remove this - not needed h2_cr works without
533  mole.b[i] = SMALLFLOAT;*/
534  }
535  }
536  else if( mole.b[i] == 0. )
537  {
538  /* this is not used for anything in calling routine and
539  * could be cleaned up */
540  *lgZerPop = true;
541  /* >>chng 04 feb 28, zero pop not really a problem in very deep neutral
542  * gas - this happens */
543 # if 0
544  if( /*nzone>0 ||*/ CODEBUG>1 )
545  {
546  fprintf(ioQQQ," PROBLEM ");
547  fprintf(ioQQQ,
548  " CO_solve zero pop for species %li %s, value is %.2e zone %li\n",
549  i , COmole[i]->label , mole.b[i], nzone);
550  }
551 # endif
552  mole.b[i] = SMALLFLOAT;
553  }
554  COmole[i]->hevmol = (realnum)mole.b[i];
555  /* check for NaN */
556  ASSERT( !isnan( COmole[i]->hevmol ) );
557  }
558  /* >>chng 04 mar 06 pass negative pop as a problem, but not a fatal one,
559  * also do not call this during 0th zone when search for conditions
560  * is underway */
561  /* 05 jul 16, there had been a bug such that CO was always evaluated as long as
562  * the temperature was below 20,000K. Once fixed, on sim turned on CO mid way
563  * into the H+ zone and had neg pops during first trys - change check on whether it
564  * is to be commented on only if this is not first zone with CO soln */
565  /*if( *lgNegPop && (nzone>0 &&!conv.lgSearch) )*/
566  if( *lgNegPop && (nzone>co.co_nzone &&!conv.lgSearch) )
567  {
568  conv.lgConvPops = false;
569  fprintf(ioQQQ," CO network negative population occurred, Te=%.4e, calling ConvFail. ",
570  phycon.te);
571  fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
572  fprintf( ioQQQ, "\n" );
573  /*ConvFail("pops", "CO");*/
574  *lgNegPop = false;
575  }
576  /* if negative pops were present need to renorm b to get
577  * proper sum rule */
578  if( 0 && *lgNegPop )
579  {
580 
581  for(j=0; j<2; ++j )
582  {
583  double sumcar = 0.;
584  double sumoxy = 0.;
585  double renorm;
586  for( i=0; i < mole.num_comole_calc; i++ )
587  {
588  /* this case, different molecules */
589  sumcar += COmole[i]->hevmol*COmole[i]->nElem[ipCARBON];
590  sumoxy += COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN];
591  }
592  renorm = (tot_ion[0] + tot_ion[1]) / SDIV( sumcar+sumoxy);
593  if(j)
594  fprintf(ioQQQ,"\t%f\n", renorm);
595  else
596  fprintf(ioQQQ,"renormco\t%.2f\t%f", fnzone, renorm);
597  for( i=0; i < mole.num_comole_calc; i++ )
598  {
599  COmole[i]->hevmol *= (realnum)renorm;
600  /* check for NaN */
601  ASSERT( !isnan( COmole[i]->hevmol ) );
602  }
603  }
604  }
605 
606  if( merror != 0 )
607  {
608  fprintf( ioQQQ, " COMOLE matrix inversion error, MERROR=%5ld zone=%5ld\n",
609  (long)merror, nzone );
610  ShowMe();
611  fprintf( ioQQQ, " Product matrix\n " );
612 
613  for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
614  {
615  fprintf( ioQQQ, "%4.4s", COmole[i]->label );
616  }
617  fprintf( ioQQQ, " \n" );
618 
619  for( j=0; j < mole.num_comole_calc; j++ )
620  {
621  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
622  fprintf( ioQQQ, " " );
623 
624  for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
625  {
626  fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
627  COmole[i]->hevmol );
628  }
629  fprintf( ioQQQ, "\n" );
630  }
631 
632  if( mole.num_comole_calc >= 9 )
633  {
634  fprintf( ioQQQ, " COMOLE matrix\n " );
635  for( i=8; i < mole.num_comole_calc; i++ )
636  {
637  fprintf( ioQQQ, "%4.4s", COmole[i]->label );
638  }
639  fprintf( ioQQQ, " \n" );
640 
641  for( j=0; j < mole.num_comole_calc; j++ )
642  {
643  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
644  fprintf( ioQQQ, " " );
645  for( i=8; i < mole.num_comole_calc; i++ )
646  {
647  fprintf( ioQQQ, "%9.1e",
648  mole.amat[i][j]* COmole[i]->hevmol );
649  }
650  fprintf( ioQQQ, "\n" );
651  }
652  }
653 
654  fprintf( ioQQQ, " Mole dens:" );
655  for( j=0; j < mole.num_comole_calc; j++ )
656  {
657  fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label
658  , COmole[j]->hevmol );
659  }
660  fprintf( ioQQQ, " \n" );
661 
662  cdEXIT(EXIT_FAILURE);
663  }
664 
666  {
667  fprintf( ioQQQ, " Product matrix\n " );
668 
669  for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
670  {
671  fprintf( ioQQQ, "%4.4s", COmole[i]->label );
672  }
673  fprintf( ioQQQ, " \n" );
674 
675  for( j=0; j < mole.num_comole_calc; j++ )
676  {
677  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
678  fprintf( ioQQQ, " " );
679  for( i=0; i < MIN2(mole.num_comole_calc,8); i++ )
680  {
681  fprintf( ioQQQ, "%9.1e", mole.amat[i][j]*
682  COmole[i]->hevmol );
683  }
684  fprintf( ioQQQ, "\n" );
685 
686  }
687 
688  if( mole.num_comole_calc >= 9 )
689  {
690  fprintf( ioQQQ, " COMOLE matrix\n " );
691  for( i=8; i < mole.num_comole_calc; i++ )
692  {
693  fprintf( ioQQQ, "%4.4s", COmole[i]->label );
694  }
695  fprintf( ioQQQ, " \n" );
696 
697  for( j=0; j < mole.num_comole_calc; j++ )
698  {
699  fprintf( ioQQQ, " %4.4s", COmole[j]->label );
700  fprintf( ioQQQ, " " );
701 
702  for( i=8; i < mole.num_comole_calc; i++ )
703  {
704  fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* COmole[i]->hevmol );
705  }
706  fprintf( ioQQQ, "\n" );
707  }
708  }
709 
710  fprintf( ioQQQ, " Mole dens:" );
711  for( j=0; j < mole.num_comole_calc; j++ )
712  {
713  fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label
714  , COmole[j]->hevmol );
715  }
716  fprintf( ioQQQ, " \n" );
717  }
718 
719  /* heating due to CO photodissociation */
720  co.CODissHeat = (realnum)(CO_findrate("PHOTON,CO=>C,O")*1e-12);
721 
722  thermal.heating[0][9] = co.CODissHeat;
723 
724  /* the real multi-level model molecule */
725  abundan = findspecies("CO")->hevmol;
726  /* IonStg and nelem were set to 2, 0 in makelevlines */
727  ASSERT( C12O16Rotate[0].Hi->IonStg < LIMELM+2 );
728  ASSERT( C12O16Rotate[0].Hi->nelem-1 < LIMELM+2 );
729  dense.xIonDense[ C12O16Rotate[0].Hi->nelem-1][C12O16Rotate[0].Hi->IonStg-1] = abundan;
730 
731  CO_PopsEmisCool(&C12O16Rotate, nCORotate,abundan , "12CO",
733  /*if( nzone>400 )
734  fprintf(ioQQQ,"DEBUG co cool\t%.2f\t%.4e\t%li\t%.4e\t%.4e\n",
735  fnzone,CoolHeavy.C12O16Rot,nCORotate,abundan,phycon.te );*/
736 
737  abundan = findspecies("CO")->hevmol/co.C12_C13_isotope_ratio;
738  /* IonStg and nelem were set to 3, 0 in makelevlines */
739  ASSERT( C13O16Rotate[0].Hi->IonStg < LIMELM+2 );
740  ASSERT( C13O16Rotate[0].Hi->nelem-1 < LIMELM+2 );
741  dense.xIonDense[ C13O16Rotate[0].Hi->nelem-1][C13O16Rotate[0].Hi->IonStg-1] = abundan;
742 
743  CO_PopsEmisCool(&C13O16Rotate, nCORotate,abundan ,"13CO",
745 
746  /* now set total density of each element locked in gas phase */
747  for( i=0;i<LIMELM; ++i )
748  {
749  dense.xMolecules[i] = 0.;
750  }
751 
752  /* total number of H per unit vol in molecules,
753  * of course not including H0/H+ */
754  for(i=0;i<N_H_MOLEC;i++)
755  {
757  }
759 
760  /* >>chng 02 sep 05, dense.xMolecules[[ipHYDROGEN] is the part of H
761  * that is not computed in hmole
762  * add in gas phase abundances locked in molecules
763  * H is special since treated in ion_solver with all hmole molecules
764  * xMolecules are the densities of these species that are done in co,
765  * this sum is only up to mole.num_heavy_molec and so does not include the atoms/ions */
766  dense.H_sum_in_CO = 0;
767  for(i=0; i<mole.num_comole_calc; ++i)
768  {
769  if(COmole[i]->n_nuclei <= 1)
770  continue;
771  /* special sum of H in CO network */
773 
774  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
775  {
776  if( mole.lgElem_in_chemistry[nelem] )
777  {
778  if( COmole[i]->hevmol > BIGFLOAT )
779  {
780  fprintf(ioQQQ, "PROBLEM DISASTER mole_co_solve has found "
781  "a CO-network molecule with an insane abundance.\n");
782  fprintf(ioQQQ, "Species number %li with abundance %.2e\n",
783  i , COmole[i]->hevmol);
784  TotalInsanity();
785  }
786  else if( conv.lgSearch &&
787  COmole[i]->nElem[nelem]*COmole[i]->hevmol > dense.gas_phase[nelem])
788  {
789  if( !conv.lgSearch )
790  fprintf(ioQQQ,
791  "PROBLEM COmole limit trip call %li Seatch? %c nMole %li "
792  "n(mole) %.2e n(gas) %.2e)\n",
794  TorF(conv.lgSearch),i ,
795  COmole[i]->nElem[nelem]*COmole[i]->hevmol,
796  dense.gas_phase[nelem]);
797  COmole[i]->hevmol = dense.gas_phase[nelem]/COmole[i]->nElem[nelem]/2.f;
798  }
799  dense.xMolecules[nelem] +=COmole[i]->nElem[nelem]*COmole[i]->hevmol;
800  }
801  }
802  }
803 
804  /* This takes the average of an element's atomic and singly ionized density from ion_solver
805  and comole and sets the solution to the updated COmole[i]->hevmol value. The macro mole.num_heavy_molec
806  is the number of heavy element molecules in the network. In addition there are 2*mole.num_elements
807  in the network, half atomic elements and half ionized elements. Finally, the total number of species
808  in the network, including elements, equals mole.num_comole_calc. The following will first take the average of
809  the ionized species (between mole.num_heavy_molec and mole.num_heavy_molec+mole.num_elements) and the other statement will
810  take the average of the atomic species (between mole.num_heavy_molec+mole.num_elements and mole.num_comole_calc).
811  This is generalized, so that if one wants to insert, for example, a sodium chemistry, once the macro's
812  are updated this if statement immediately works */
813 
814  /* >> chng 2006 Sep 02 rjrw: use vectors rather than offsets above to index, for ease of generalization */
815 
816  for(i=0;i<mole.num_elements; ++i)
817  {
818  element = chem_element[i];
819  /*printf("MOL\t%3e\tION\t%e\tSPECIES\t%s\n", COmole[i]->hevmol,dense.xIonDense[COmole[i]->nelem_hevmol][1],COmole[i].label);*/
820  COmole[element->ipMlP]->hevmol = ((dense.xIonDense[element->ipCl][1]+ COmole[element->ipMlP]->hevmol)/2)*dense.lgElmtOn[element->ipCl];
821  /*printf("MOL\t%3e\tION\t%e\tSPECIES\t%s\n", COmole[i]->hevmol,dense.xIonDense[COmole[i]->nelem_hevmol][0],COmole[i].label);*/
822  COmole[element->ipMl]->hevmol = ((dense.xIonDense[element->ipCl][0]+ COmole[element->ipMl]->hevmol)/2)*dense.lgElmtOn[element->ipCl];
823 
824  /* check for NaN */
825  ASSERT( !isnan( COmole[element->ipMlP]->hevmol ) );
826  ASSERT( !isnan( COmole[element->ipMl]->hevmol ) );
827  }
828 
829  /* check whether ion and chem solvers agree yet */
830  if( dense.lgElmtOn[ipCARBON] &&
831  fabs(dense.xIonDense[ipCARBON][0]- findspecies("C")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
832  {
833  conv.lgConvIoniz = false;
834  sprintf( conv.chConvIoniz, "CO C0 con");
836  conv.BadConvIoniz[1] = findspecies("C")->hevmol;
837  }
838  else if( dense.lgElmtOn[ipCARBON] &&
839  fabs(dense.xIonDense[ipCARBON][1]- findspecies("C+")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
840  {
841  conv.lgConvIoniz = false;
842  sprintf( conv.chConvIoniz, "CO C1 con");
844  conv.BadConvIoniz[1] = findspecies("C+")->hevmol;
845  }
846  else if( dense.lgElmtOn[ipOXYGEN] &&
847  fabs(dense.xIonDense[ipOXYGEN][0]- findspecies("O")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
848  {
849  conv.lgConvIoniz = false;
850  sprintf( conv.chConvIoniz, "CO O0 con");
852  conv.BadConvIoniz[1] = findspecies("O")->hevmol;
853  }
854  else if( dense.lgElmtOn[ipOXYGEN] &&
855  fabs(dense.xIonDense[ipOXYGEN][1]- findspecies("O+")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
856  {
857  conv.lgConvIoniz = false;
858  sprintf( conv.chConvIoniz, "CO O1 con");
860  conv.BadConvIoniz[1] = findspecies("O+")->hevmol;
861  }
862 
863  /* now update ionization distribution of the elements we just did */
864  for(i=0;i<mole.num_elements; ++i)
865  {
866  element = chem_element[i];
867  dense.xIonDense[element->ipCl][0] = COmole[element->ipMl]->hevmol*dense.lgElmtOn[element->ipCl];
868  dense.xIonDense[element->ipCl][1] = COmole[element->ipMlP]->hevmol*dense.lgElmtOn[element->ipCl];
869  dense.IonLow[element->ipCl] = 0;
870  dense.IonHigh[element->ipCl] = MAX2( dense.IonHigh[element->ipCl] , 1 );
871  }
872 
873  /* if populations not conserved then not converged */
874 # define EPS_MOLE 0.1
875  if( dense.lgElmtOn[ipHYDROGEN] &&
877  {
878  conv.lgConvIoniz = false;
879  sprintf( conv.chConvIoniz, "COcon%2i",ipHYDROGEN );
882  }
883  else if( dense.lgElmtOn[ipCARBON] &&
885  {
886  conv.lgConvIoniz = false;
887  sprintf( conv.chConvIoniz, "COcon%2i",ipCARBON );
890  }
891  else if( dense.lgElmtOn[ipOXYGEN] &&
893  {
894  conv.lgConvIoniz = false;
895  sprintf( conv.chConvIoniz, "COcon%2i",ipOXYGEN );
898  }
899  else if( dense.lgElmtOn[ipSILICON] &&
901  {
902  conv.lgConvIoniz = false;
903  sprintf( conv.chConvIoniz, "COcon%2i",ipSILICON );
906  }
907  else if( dense.lgElmtOn[ipSULPHUR] &&
909  {
910  conv.lgConvIoniz = false;
911  sprintf( conv.chConvIoniz, "COcon%2i",ipSULPHUR );
914  }
915  else if( dense.lgElmtOn[ipNITROGEN] &&
917  {
918  conv.lgConvIoniz = false;
919  sprintf( conv.chConvIoniz, "COcon%2i",ipNITROGEN );
922  }
923 
924  else if( dense.lgElmtOn[ipCHLORINE] &&
926  {
927  conv.lgConvIoniz = false;
928  sprintf( conv.chConvIoniz, "COcon%2i",ipCHLORINE );
931  }
932 
933  /* the total N photo dissociation rate, cm-3 s-1 */
935  return;
936 
937 /*lint +e550 */
938 /*lint +e778 constant expression evaluates to 0 in operation '-' */
939 }
940 
941 # undef EPS_MOLE

Generated for cloudy by doxygen 1.8.3.1