cloudy  trunk
mole_solve.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 "called.h"
6 #include "dense.h"
7 #include "deuterium.h"
8 #include "ionbal.h"
9 #include "thermal.h"
10 #include "phycon.h"
11 #include "hmi.h"
12 #include "dynamics.h"
13 #include "conv.h"
14 #include "trace.h"
15 #include "timesc.h"
16 #include "mole.h"
17 #include "mole_priv.h"
18 #include "grainvar.h"
19 #include "h2.h"
20 #include "newton_step.h"
21 /* Nick Abel between July and October of 2003 assisted Dr. Ferland in
22  * improving the heavy element molecular network in Cloudy. Before
23  * this routine would predict negative abundances if the fraction of
24  * carbon in the form of molecules came close to 100%. A reorganizing
25  * of the reaction network detected several bugs. Treatment of
26  * "coupled reactions", in which both densities in the reaction rate
27  * were being predicted by Cloudy, were also added. Due to these
28  * improvements, Cloudy can now perform calculations where 100% of the
29  * carbon is in the form of CO without predicting negative abundances
30  *
31  * Additional changes were made in November of 2003 so that our reaction
32  * network would include all reactions from the TH85 paper. This involved
33  * adding silicon to the chemical network. Also the reaction rates were
34  * labeled to make identification with the reaction easier and the matrix
35  * elements of atomic C, O, and Si are now done in a loop, which makes
36  * the addition of future chemical species (like N or S) easy.
37  * */
38 void check_co_ion_converge(void);
39 
40 STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec,
41  double * const ervals, double * const amat, const bool lgJac, bool *lgConserve);
42 
43 STATIC void mole_h_fixup(void);
44 
45 STATIC void grouped_elems(const double bvec[], double mole_elems[]);
46 
47 #define SMALLABUND 1e-24
48 
49 double mole_solve()
50 {
51  int nBad, nPrevBad, i;
52  realnum
53  eqerror, error;
54  valarray<double> oldmols(mole_global.num_compacted),
55  newmols(mole_global.num_compacted);
56  GroupMap MoleMap( atom_list.size() );
57 
58  DEBUG_ENTRY( "mole_solve()" );
59 
60  if (hmi.H2_frac_abund_set>0.)
61  {
62  mole_h_fixup();
63  fixit();
64  fprintf(stderr,"Need to treat hmi.H2_frac_abund_set in revised network\n");
65  fprintf(stderr,"%g\n",hmi.H2_frac_abund_set);
66  exit(-1);
67  }
68 
69  ASSERT(lgElemsConserved());
70 
71  /* will test against this error for quitting */
72  const double BIGERROR = 1e-8;
73  /* >>chng 04 jul 20, upper limit had been 6, why? change to 20
74  * to get closer to soln */
75  const int LIM_LOOP = 100;
76  /* loop until run to limit, or no fix ups needed and error is small */
77  /* will be used to keep track of neg solns */
78  nBad = nPrevBad = 0;
79  eqerror = -1.;
80 
81  valarray<double> escale(mole_global.num_compacted);
82  double rlimit=-1., rmax=0.0;
83 
84  // Run enough times through to converge nonlinear system.
85  for(i=0; i < LIM_LOOP;i++)
86  {
87  {
88  /* option to print deduced abundances */
89  /*@-redef@*/
90  enum {DEBUG_LOC=false};
91  /*@+redef@*/
92  if( DEBUG_LOC && (nzone>140) )
93  {
94  fprintf(ioQQQ,"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
95  fnzone,
96  phycon.te,
97  dense.eden);
98  for( long k=0; k<mole_global.num_calc; k++ )
99  fprintf(ioQQQ,"\t%.2e", mole.species[k].den );
100  fprintf(ioQQQ,"\n" );
101  }
102  }
103  /* nBad returns number of times negative abundances were fixed
104  * for latest step, should be zero at return for valid soln */
105 
106  nPrevBad = nBad;
107 
108  MoleMap.setup(get_ptr(oldmols));
109 
110  long j;
111 
112  // Need to allow enough loops for j that rlimit increases until
113  // no species get to -ve densities. Often, though, this will
114  // just be the first time through.
115  conv.incrementCounter(MOLE_SOLVE);
116  for (j=0; j<30; j++)
117  {
118  conv.incrementCounter(MOLE_SOLVE_STEPS);
119  bool lgOK = newton_step(MoleMap, oldmols, newmols, &eqerror, &error, mole_global.num_compacted, &rlimit, &rmax, escale, funjac);
120 
121  /* check for negative populations */
122  if (lgOK)
123  {
124  nBad = 0;
125  double fracneg = 0.;
126  long iworst = -1;
127  for( long mol=0; mol < mole_global.num_compacted; mol++ )
128  {
129  if( newmols[mol] < 0. )
130  {
131  if(-newmols[mol] > fracneg*SDIV(oldmols[mol]))
132  {
133  fracneg = -newmols[mol]/SDIV(oldmols[mol]);
134  iworst = mol;
135  }
136  if( newmols[mol] < -SMALLABUND)
137  {
138  ++nBad;
139  }
140  newmols[mol] = 0.;
141  }
142  }
143  if ( 0 != nBad)
144  {
145  lgOK = false;
146  if (0)
147  {
148  fprintf(ioQQQ,"-ve density in inner chemistry loop, worst is %ld\n",iworst);
149  }
150  }
151  }
152  if ( lgOK )
153  {
154  //fprintf(ioQQQ,"OK at z %ld l %d j %ld t %g e %g\n",nzone,i,j,1./rlimit,eqerror);
155  break;
156  }
157 
158  ASSERT( rlimit < BIGFLOAT );
159  ASSERT( rlimit > 0.0 );
160  rlimit = 10.*rlimit;
161  }
162 
163  //fprintf(stdout,"Mole zone %ld -- %7.2f error %15.8g eqerror %15.8g rlimit %15.8g nBad %d ninner %ld loop %d\n",nzone,fnzone,error,eqerror,rlimit,nBad,j+1,i);
164 
165  // If the accuracy of the solution obtained was strongly limited
166  // by the pseudo-timestep limit, extend it
167  if ( error < 0.01*eqerror )
168  rlimit = 0.1*rlimit;
169 
170  MoleMap.updateMolecules( newmols );
171 
172  /* Stop early if good enough */
173  if (eqerror < BIGERROR && nBad == 0 && nPrevBad == 0)
174  break;
175  }
176 
177  //fprintf(stdout,"Mole: zn %ld -- %7.2f err %13.6g rlimit %13.6g nBad %d loop %d\n",nzone,fnzone,eqerror,rlimit,nBad,i);
178 
179  if( (i == LIM_LOOP && eqerror > BIGERROR) || nBad != 0 )
180  {
181  conv.setConvIonizFail("Chem conv", eqerror, nBad);
182  }
183 
184  // Effect_delta determines whether changes due to the molecular
185  // network rescale values enough to invalidate other solutions
186  realnum effect_delta = mole_return_cached_species(MoleMap);
187 
188  realnum delta_threshold = 0.2*conv.IonizErrorAllowed;
189  if (effect_delta > delta_threshold)
190  {
191  conv.setConvIonizFail("chem eft chg", effect_delta, 0.0);
192  }
193 
194  if (trace.lgTrace)
195  {
196  fprintf(ioQQQ,"Mole zn %3ld -- %7.2f\n",nzone,fnzone);
197  fprintf(ioQQQ,"Internal error %15.8g nBad %4d loops %3d\n",
198  eqerror,nBad,i);
199  fprintf(ioQQQ,"Scaling delta on ions %15.8g; threshold %15.8g\n",
200  effect_delta, delta_threshold);
201  }
202 
203  check_co_ion_converge();
204 
205  {
206  /* option to print deduced abundances */
207  /*@-redef@*/
208  enum {DEBUG_LOC=false};
209  /*@+redef@*/
210  if( DEBUG_LOC /*&& (nzone>68)*/ )
211  {
212  fprintf(ioQQQ,"DEBUG mole out\t%i\t%.2f\t%.5e\t%.5e",
213  i,
214  fnzone,
215  phycon.te,
216  dense.eden);
217  fprintf(ioQQQ,
218  "\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e",
219  eqerror,
220  dense.xIonDense[ipHYDROGEN][0],
221  dense.xIonDense[ipHYDROGEN][1],
222  mole.sink[ipHYDROGEN][0],
223  mole.sink[ipHYDROGEN][1],
224  mole.source[ipHYDROGEN][0] ,
225  mole.source[ipHYDROGEN][1] ,
226  ionbal.RateIonizTot(ipHYDROGEN,0),
227  ionbal.RateRecomTot[ipHYDROGEN][0]);
228 
229  /*for( j=0; j<mole_global.num_calc; j++ )
230  fprintf(ioQQQ,"\t%.4e", HMOLEC(j) );*/
231  fprintf(ioQQQ,"\n" );
232  }
233  }
234 
235  return eqerror;
236 
237 /*lint +e550 */
238 /*lint +e778 constant expression evaluates to 0 in operation '-' */
239 }
240 
241 
242 void check_co_ion_converge(void)
243 {
244  DEBUG_ENTRY("check_co_ion_converge()");
245  /* check whether ion and chem solvers agree yet */
246  if( dense.lgElmtOn[ipCARBON] &&
247  fabs(dense.xIonDense[ipCARBON][0]- findspecieslocal("C")->den)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
248  {
249  conv.setConvIonizFail("CO C0 con",
250  dense.xIonDense[ipCARBON][0],
251  findspecieslocal("C")->den);
252  }
253  else if( dense.lgElmtOn[ipCARBON] &&
254  fabs(dense.xIonDense[ipCARBON][1]- findspecieslocal("C+")->den)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
255  {
256  conv.setConvIonizFail("CO C1 con",
257  dense.xIonDense[ipCARBON][1],
258  findspecieslocal("C+")->den);
259  }
260  else if( dense.lgElmtOn[ipOXYGEN] &&
261  fabs(dense.xIonDense[ipOXYGEN][0]- findspecieslocal("O")->den)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
262  {
263  conv.setConvIonizFail(
264  "CO O0 con",dense.xIonDense[ipOXYGEN][0],
265  findspecieslocal("O")->den);
266  }
267  else if( dense.lgElmtOn[ipOXYGEN] &&
268  fabs(dense.xIonDense[ipOXYGEN][1]- findspecieslocal("O+")->den)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
269  {
270  conv.setConvIonizFail(
271  "CO O1 con", dense.xIonDense[ipOXYGEN][1],
272  findspecieslocal("O+")->den);
273  }
274 }
275 
276 STATIC void mole_h_fixup(void)
277 {
278  long int mol;
279 
280  /* there are two "no molecules" options, the no co, which turns off everything,
281  * and the no n2, which only turns off the h2. in order to not kill the co
282  * part we still need to compute the hydrogen network here, and then set h2 to
283  * small values */
284  /* >> chng 03 jan 15 rjrw -- suddenly switching off molecules confuses the solvers... */
285  DEBUG_ENTRY( "mole_h_fixup()" );
286 
287  /* >>chng 02 jun 19, add option to force H2 abundance, for testing h2 molecules,
288  * hmi.H2_frac_abund_set is fraction in molecules that is set by set h2 fraction command */
289  if( hmi.H2_frac_abund_set>0.)
290  {
291  for(mol=0;mol<mole_global.num_calc;mol++)
292  {
293  mole.species[mol].den = 0.;
294  }
295  /* >>chng 03 jul 19, from 0 to SMALLFLOAT, to pass asserts in ConvBase,
296  * problem is that ion range has not been reset for hydrogen */
297  dense.xIonDense[ipHYDROGEN][0] = dense.xIonDense[ipHYDROGEN][1] =
298  2.f*SMALLFLOAT*dense.gas_phase[ipHYDROGEN];
299  /* put it all in the ground state */
300  findspecieslocal("H2")->den = (realnum)(dense.gas_phase[ipHYDROGEN] * hmi.H2_frac_abund_set);
301  findspecieslocal("H2*")->den = 0.;
302 
303  hmi.H2_total = findspecieslocal("H2")->den + findspecieslocal("H2*")->den;
304  /* first guess at ortho and para densities */
305  h2.ortho_density = 0.75*hmi.H2_total;
306  h2.para_density = 0.25*hmi.H2_total;
307  {
308  hmi.H2_total_f = (realnum)hmi.H2_total;
309  h2.ortho_density_f = (realnum)h2.ortho_density;
310  h2.para_density_f = (realnum)h2.para_density;
311  }
312 
313  hmi.hmihet = 0.;
314  hmi.h2plus_heat = 0.;
315  hmi.H2Opacity = 0.;
316  hmi.hmicol = 0.;
317  hmi.HeatH2Dish_TH85 = 0.;
318  hmi.HeatH2Dexc_TH85 = 0.;
319  hmi.deriv_HeatH2Dexc_TH85 = 0.;
320  hmi.hmidep = 1.;
321 
322  for( size_t nd=0; nd < gv.bin.size(); nd++ )
323  {
324  gv.bin[nd]->rate_h2_form_grains_used = 0.;
325  }
326 
327  return;
328  }
329 }
330 
331 #define ABSLIM 1e-12
332 #define ERRLIM 1e-12
333 # ifdef MAT
334 # undef MAT
335 # endif
336 # define MAT(a,I_,J_) ((a)[(I_)*(mole_global.num_compacted)+(J_)])
337 
338 
339 STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c);
340 enum {PRINTSOL = false};
341 
342 STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec, double * const ervals,
343  double * const amat, const bool lgJac, bool *lgConserve)
344 {
345  static multi_arr<double,2> c(mole_global.num_total,mole_global.num_total);
346  bool printsol = PRINTSOL;
347  valarray<double> b(mole_global.num_total);
348 
349  DEBUG_ENTRY( "funjac()" );
350 
351  MoleMap.updateMolecules( b2vec );
352 
353  if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
354  dynamics.Rate != 0.0)
355  {
356  ASSERT(dynamics.Rate > 0.0);
357  *lgConserve = false;
358  }
359  else
360  {
361  *lgConserve = true;
362  }
363 
364  /* Generate chemical balance vector (mole.b[]) and Jacobian array
365  (c[][], first iteration) from reaction list */
366 
367  mole_eval_dynamic_balance(mole_global.num_total,get_ptr(b),lgJac,c);
368 
369  /*------------------------------------------------------------------ */
370  if(printsol || (trace.lgTrace && trace.lgTraceMole ))
371  {
372  /* print the full matrix */
373  fprintf( ioQQQ, " ");
374  for( long i=0; i < mole_global.num_calc; i++ )
375  {
376  fprintf( ioQQQ, " %-4.4s", mole_global.list[i]->label.c_str() );
377  }
378  fprintf( ioQQQ, " \n" );
379 
380  fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone);
381  for( long i=0; i<mole_global.num_calc; i++ )
382  fprintf(ioQQQ,"\t%.2e", mole.species[i].den );
383  fprintf(ioQQQ,"\n" );
384 
385  for( long i=0; i < mole_global.num_calc; i++ )
386  {
387  fprintf( ioQQQ, " MOLE%2ld %-4.4s", i ,mole_global.list[i]->label.c_str());
388  for( long j=0; j < mole_global.num_calc; j++ )
389  {
390  fprintf( ioQQQ, "%10.2e", c[j][i] );
391  }
392  fprintf( ioQQQ, "%10.2e", b[i] );
393  fprintf( ioQQQ, "\n" );
394  }
395  }
396 
397  /* add positive ions and neutral atoms: ratios are set by ion_solver,
398  * we determine abundance of the group as a whole here */
399 
400  if (lgJac) {
401  for (unsigned long j=0; j<atom_list.size(); ++j )
402  {
403  if (atom_list[j]->ipMl[0] != -1)
404  {
405  for(long i=0;i<mole_global.num_calc;i++)
406  {
407  ASSERT(mole_global.list[i]->index == i);
408  double sum = 0.;
409  for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
410  {
411  if (atom_list[j]->ipMl[ion] != -1)
412  {
413  sum += MoleMap.fion[j][ion]*c[atom_list[j]->ipMl[ion]][i];
414  }
415  }
416  c[atom_list[j]->ipMl[0]][i] = sum;
417  }
418  }
419  }
420  for (unsigned long j=0; j<atom_list.size(); ++j )
421  {
422  if (atom_list[j]->ipMl[0] != -1)
423  {
424  for(long i=0;i<mole_global.num_calc;i++)
425  {
426  double sum = 0.0;
427 
428  for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
429  {
430  if (atom_list[j]->ipMl[ion] != -1)
431  {
432  sum += c[i][atom_list[j]->ipMl[ion]];
433  }
434  }
435  c[i][atom_list[j]->ipMl[0]] = sum;
436  }
437  }
438  }
439  }
440 
441  for (unsigned long j=0; j<atom_list.size(); ++j )
442  {
443  //if (j == 8)
444  // fprintf(ioQQQ,"Nsum1 %s %ld %d %d\n",
445  // atom_list[j]->label().c_str(), j, groupspecies[29]->index,atom_list[j]->ipMl[0]);
446  if (atom_list[j]->ipMl[0] != -1)
447  {
448  double sum = 0.0;
449  for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
450  {
451  if (atom_list[j]->ipMl[ion] != -1)
452  {
453  // if (j == 8)
454  // fprintf(ioQQQ,"Nsum %d %15.8g\n",atom_list[j]->ipMl[ion],
455  // b[atom_list[j]->ipMl[ion]]);
456  sum += b[atom_list[j]->ipMl[ion]];
457  b[atom_list[j]->ipMl[ion]] = 0.0;
458  }
459  }
460  b[atom_list[j]->ipMl[0]] = sum;
461  }
462  }
463 
464  /* Species are now grouped -- only mole_global.num_compacted elements now active */
465 
466  if (*lgConserve)
467  {
468  // Replace atom rows with conservation constraints
469  if (lgJac)
470  {
471  for ( unsigned long j = 0; j < atom_list.size(); ++j )
472  {
473  long ncons = atom_list[j]->ipMl[0];
474  if (ncons != -1)
475  {
476  ASSERT( mole_global.list[ncons]->isMonatomic() );
477  double scale=1.0/MoleMap.molElems[j]; //fabs(c[ncons][ncons]);//
478  for( long i=0;i<mole_global.num_compacted;i++)
479  {
480  if( groupspecies[i]->nAtom.find(atom_list[j]) != groupspecies[i]->nAtom.end() )
481  c[groupspecies[i]->index][ncons] = groupspecies[i]->nAtom[atom_list[j]]*scale;
482  else
483  c[groupspecies[i]->index][ncons] = 0.;
484 
485  }
486  }
487  }
488  }
489 
490  valarray<double> molnow(atom_list.size());
491  grouped_elems(get_ptr(b2vec), get_ptr(molnow));
492  for ( unsigned long j = 0; j < atom_list.size(); ++j )
493  {
494  long ncons = atom_list[j]->ipMl[0];
495  if (ncons != -1)
496  {
497  ASSERT( mole_global.list[ncons]->isMonatomic() );
498  double scale = c[ncons][ncons];
499  b[ncons] = (molnow[j]-MoleMap.molElems[j])*scale;
500  if (false)
501  if (b[ncons] != 0.0)
502  fprintf(ioQQQ,"Cons %s err %g rel %g\n",
503  atom_list[j]->label().c_str(),b[ncons],
504  (molnow[j]-MoleMap.molElems[j])/SDIV(MoleMap.molElems[j]));
505  }
506  }
507  }
508 
509  for( long i=0; i < mole_global.num_compacted; i++ )
510  {
511  ervals[i] = b[groupspecies[i]->index];
512  }
513 
514  if (lgJac)
515  {
516 
517  for( long i=0; i < mole_global.num_compacted; i++ )
518  {
519  for( long j=0; j < mole_global.num_compacted; j++ )
520  {
521  MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
522  }
523  }
524 
525  // Replace rows for species with no sources and sinks with
526  // an identity
527  for(long i=0;i<mole_global.num_compacted;i++)
528  {
529  double sum1 = 0.;
530  for (long j=0;j<mole_global.num_compacted;j++)
531  {
532  sum1 += fabs(MAT(amat,j,i));
533  }
534  if (sum1 == 0.0)
535  {
536  ASSERT(ervals[i] == 0.0);
537  MAT(amat,i,i) = 1.0;
538  // fprintf(ioQQQ,"Fixing %ld %s\n",i,groupspecies[i]->label.c_str());
539  }
540  }
541  }
542 }
543 
544 STATIC void grouped_elems(const double bvec[], double mole_elems[])
545 {
546  map<chem_atom*, long> atom_to_index;
547  for (unsigned long j=0; j<atom_list.size(); ++j )
548  {
549  mole_elems[j] = 0.;
550  atom_to_index[atom_list[j].get_ptr()] = j;
551  }
552  for( long i=0; i < mole_global.num_compacted; i++ )
553  {
554  if( groupspecies[i]->parentLabel.empty() == false )
555  continue;
556 
557  for (molecule::nAtomsMap::const_iterator el = groupspecies[i]->nAtom.begin();
558  el != groupspecies[i]->nAtom.end(); ++el)
559  {
560  mole_elems[atom_to_index[el->first.get_ptr()]] += bvec[i]*el->second;
561  }
562  }
563 }
564 
565 void GroupMap::setup(double *b0vec)
566 {
567  valarray<double> calcv(mole_global.num_total);
568  bool lgSet;
569 
570  for( long i=0;i<mole_global.num_total;i++)
571  {
572  calcv[i] = mole.species[i].den;
573  }
574 
575  for (unsigned long j=0; j<atom_list.size(); ++j )
576  {
577  if (atom_list[j]->ipMl[0] != -1)
578  {
579  double sum = 0.;
580  for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
581  {
582  if (atom_list[j]->ipMl[ion] != -1)
583  sum += calcv[atom_list[j]->ipMl[ion]];
584  }
585  if (sum > SMALLFLOAT)
586  {
587  double factor = 1./sum;
588  for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
589  {
590  if (atom_list[j]->ipMl[ion] != -1)
591  fion[j][ion] = calcv[atom_list[j]->ipMl[ion]]*factor;
592  else
593  fion[j][ion] = 0.;
594  }
595  }
596  else
597  {
598  lgSet = false;
599  for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
600  {
601  if (atom_list[j]->ipMl[ion] != -1 && !lgSet)
602  {
603  fion[j][ion] = 1.0;
604  lgSet = true;
605  }
606  else
607  {
608  fion[j][ion] = 0.;
609  }
610  }
611  }
612 
613  lgSet = false;
614  for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
615  {
616  if (atom_list[j]->ipMl[ion] != -1)
617  {
618  if (!lgSet)
619  calcv[atom_list[j]->ipMl[ion]] = sum;
620  else
621  calcv[atom_list[j]->ipMl[ion]] = 0.;
622  lgSet = true;
623  }
624  }
625  }
626  }
627 
628  for( long i=0; i < mole_global.num_compacted; i++ )
629  {
630  b0vec[i] = calcv[groupspecies[i]->index];
631  }
632 
633  grouped_elems(get_ptr(b0vec),get_ptr(molElems));
634 
635  for (unsigned long i = 0; i<atom_list.size(); ++i)
636  {
637  double densAtom = 0.;
638  // deuterium is special case
639  if( atom_list[i]->el->Z==1 && atom_list[i]->A==2 )
640  {
641  ASSERT( deut.lgElmtOn );
642  densAtom = deut.gas_phase;
643  }
644  // skip other isotopes
645  else if( !atom_list[i]->lgMeanAbundance() )
646  continue;
647  else
648  {
649  int nelem = atom_list[i]->el->Z-1;
650  densAtom = dense.gas_phase[nelem];
651  }
652  bool lgTest =
653  ( densAtom < SMALLABUND && molElems[i] < SMALLABUND ) ||
654  ( fabs(molElems[i]- densAtom) <= conv.GasPhaseAbundErrorAllowed*densAtom );
655  if( !lgTest )
656  fprintf( ioQQQ, "PROBLEM molElems BAD %li\t%s\t%.12e\t%.12e\t%.12e\n",
657  i, atom_list[i]->label().c_str(), atom_list[i]->frac, densAtom, molElems[i] );
658 
659  ASSERT( lgTest );
660  molElems[i] = densAtom;
661  }
662 
663 }
664 
665 void GroupMap::updateMolecules(const valarray<double> & b2 )
666 {
667  DEBUG_ENTRY("updateMolecules()");
668 
669  for (long mol=0;mol<mole_global.num_calc;mol++)
670  {
671  mole.species[mol].den = 0.;
672  }
673  for (long mol=0;mol<mole_global.num_compacted;mol++)
674  {
675  mole.species[ groupspecies[mol]->index ].den = b2[mol]; /* put derived abundances back into appropriate molecular species */
676  }
677 
678  // calculate isotopologue densities
679  for (long mol=0;mol<mole_global.num_calc;mol++)
680  {
681  if( mole_global.list[mol]->parentIndex >= 0 )
682  {
683  ASSERT( !mole_global.list[mol]->parentLabel.empty() );
684  long parentIndex = mole_global.list[mol]->parentIndex;
685  mole.species[mol].den = mole.species[parentIndex].den;
686  for( nAtoms_i it = mole_global.list[mol]->nAtom.begin(); it != mole_global.list[mol]->nAtom.end(); ++it )
687  {
688  if( !it->first->lgMeanAbundance() )
689  {
690  mole.species[mol].den *= pow( it->first->frac, it->second );
691  }
692  }
693  }
694  }
695 
696  // calculate densities for monatomic species that had been collapsed into a group
697  for (unsigned long j=0; j<atom_list.size(); ++j )
698  {
699  if (atom_list[j]->ipMl[0] != -1)
700  {
701  double grouptot = mole.species[atom_list[j]->ipMl[0]].den;
702  double sum = 0.0;
703  for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
704  {
705  if (atom_list[j]->ipMl[ion] != -1)
706  {
707  mole.species[atom_list[j]->ipMl[ion]].den = grouptot * fion[j][ion];
708  sum += mole.species[atom_list[j]->ipMl[ion]].den;
709  }
710  }
711  ASSERT(fabs(sum-grouptot) <= 1e-10 * fabs(grouptot));
712  }
713  }
714 
715  mole.set_isotope_abundances();
716 
717  return;
718 }
719 
720 STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c)
721 {
722  double source;
723 
724  DEBUG_ENTRY("mole_eval_dynamic_balance()");
725 
726  mole_eval_balance(num_total, b, lgJac, c);
727 
728  /* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */
729  if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
730  dynamics.Rate != 0.0 )
731  {
732  /* Don't use conservation form in matrix solution -- dynamics rate terms make c[][] non-singular */
733 
734  for( long i=0;i<mole_global.num_calc;i++)
735  {
736  if (lgJac)
737  c[i][i] -= dynamics.Rate;
738 
739  if( !mole_global.list[i]->parentLabel.empty() )
740  continue;
741 
742  b[i] -= mole.species[i].den*dynamics.Rate;
743 
744  if (!mole_global.list[i]->isMonatomic() || mole_global.list[i]->charge < 0 || ! mole_global.list[i]->lgGas_Phase ||
745  ( mole_global.list[i]->isMonatomic() && !mole_global.list[i]->nAtom.begin()->first->lgMeanAbundance() ) )
746  {
747  b[i] += dynamics.molecules[i]*dynamics.Rate;
748  }
749  else if (mole_global.list[i]->charge == 0)
750  {
751  ASSERT( mole_global.list[i]->isMonatomic() );
752  ASSERT( (int)mole_global.list[i]->nAtom.size() == 1 );
753  count_ptr<chem_atom> atom = mole_global.list[i]->nAtom.begin()->first;
754  long nelem = atom->el->Z-1;
755  if( nelem >= LIMELM )
756  continue;
757  source = 0.0;
758  for ( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
759  {
760  source += dynamics.Source[nelem][ion] * atom->frac;
761  if (unresolved_atom_list[nelem]->ipMl[ion] == -1)
762  source -= dense.xIonDense[nelem][ion]*dynamics.Rate * atom->frac;
763  }
764  b[i] += source;
765  }
766  }
767  }
768 }
thermal.h
t_hmi::H2Opacity
realnum H2Opacity
Definition: hmi.h:29
check_co_ion_converge
void check_co_ion_converge(void)
Definition: mole_solve.cpp:242
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
GroupMap::molElems
valarray< double > molElems
Definition: mole_priv.h:24
t_dense::eden
double eden
Definition: dense.h:190
count_ptr
Definition: count_ptr.h:11
dense
t_dense dense
Definition: dense.cpp:24
t_dynamics::lgAdvection
bool lgAdvection
Definition: dynamics.h:60
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
t_deuterium::lgElmtOn
bool lgElmtOn
Definition: deuterium.h:19
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
newton_step.h
realnum
float realnum
Definition: cddefines.h:103
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:198
nAtoms_i
molecule::nAtomsMap::iterator nAtoms_i
Definition: mole.h:214
unresolved_atom_list
ChemAtomList unresolved_atom_list
Definition: mole_species.cpp:70
diatomics::para_density_f
realnum para_density_f
Definition: h2_priv.h:325
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
mole.h
MOLE_SOLVE_STEPS
@ MOLE_SOLVE_STEPS
Definition: conv.h:72
t_conv::IonizErrorAllowed
realnum IonizErrorAllowed
Definition: conv.h:280
molecule::index
int index
Definition: mole.h:169
multi_arr< double, 2 >
t_conv::setConvIonizFail
void setConvIonizFail(const char *reason, double oldval, double newval)
Definition: conv.h:107
t_dynamics::Rate
double Rate
Definition: dynamics.h:71
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
funjac
STATIC void funjac(GroupMap &MoleMap, const valarray< double > &b2vec, double *const ervals, double *const amat, const bool lgJac, bool *lgConserve)
Definition: mole_solve.cpp:342
t_mole_global::list
MoleculeList list
Definition: mole.h:317
t_trace::lgTraceMole
bool lgTraceMole
Definition: trace.h:55
t_conv::GasPhaseAbundErrorAllowed
realnum GasPhaseAbundErrorAllowed
Definition: conv.h:285
mole_eval_dynamic_balance
STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
Definition: mole_solve.cpp:720
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
t_dynamics::Source
double ** Source
Definition: dynamics.h:74
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
t_dynamics::n_initial_relax
long int n_initial_relax
Definition: dynamics.h:126
PRINTSOL
@ PRINTSOL
Definition: mole_solve.cpp:340
t_hmi::hmicol
double hmicol
Definition: hmi.h:26
dynamics.h
diatomics::ortho_density
double ortho_density
Definition: h2_priv.h:319
diatomics::ortho_density_f
realnum ortho_density_f
Definition: h2_priv.h:324
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
t_hmi::H2_total
double H2_total
Definition: hmi.h:16
chem_atom::el
chem_element * el
Definition: mole.h:43
b2
static double b2[63]
Definition: atmdat_3body.cpp:19
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
GroupMap::updateMolecules
void updateMolecules(const valarray< double > &b2)
Definition: mole_solve.cpp:665
t_mole_local::sink
double ** sink
Definition: mole.h:394
nzone
long int nzone
Definition: cddefines.cpp:14
mole_return_cached_species
realnum mole_return_cached_species(const GroupMap &MoleMap)
Definition: mole_species.cpp:935
t_deuterium::gas_phase
realnum gas_phase
Definition: deuterium.h:20
source
static double * source
Definition: species2.cpp:28
GroupMap
Definition: mole_priv.h:21
MOLE_SOLVE
@ MOLE_SOLVE
Definition: conv.h:71
dense.h
mole
t_mole_local mole
Definition: mole.cpp:7
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
GroupMap::fion
multi_arr< double, 2 > fion
Definition: mole_priv.h:23
t_hmi::HeatH2Dish_TH85
double HeatH2Dish_TH85
Definition: hmi.h:130
t_mole_global::num_total
int num_total
Definition: mole.h:314
diatomics::para_density
double para_density
Definition: h2_priv.h:321
t_hmi::hmihet
double hmihet
Definition: hmi.h:24
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
t_hmi::deriv_HeatH2Dexc_TH85
realnum deriv_HeatH2Dexc_TH85
Definition: hmi.h:146
groupspecies
molecule ** groupspecies
Definition: mole_species.cpp:72
deut
t_deuterium deut
Definition: deuterium.cpp:8
hmi.h
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_hmi::hmidep
double hmidep
Definition: hmi.h:33
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
t_mole_local::set_isotope_abundances
void set_isotope_abundances(void)
Definition: mole_species.cpp:989
t_mole_local::source
double ** source
Definition: mole.h:394
BIGFLOAT
const UNUSED realnum BIGFLOAT
Definition: cpu.h:189
grouped_elems
STATIC void grouped_elems(const double bvec[], double mole_elems[])
Definition: mole_solve.cpp:544
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
mole_solve
double mole_solve()
Definition: mole_solve.cpp:49
iteration
long int iteration
Definition: cddefines.cpp:16
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
grainvar.h
timesc.h
SMALLABUND
#define SMALLABUND
Definition: mole_solve.cpp:47
ionbal.h
chem_element::Z
const int Z
Definition: mole.h:27
newton_step
bool newton_step(GroupMap &MoleMap, const valarray< double > &b0vec, valarray< double > &b2vec, realnum *eqerror, realnum *error, const long n, double *rlimit, double *rmax, valarray< double > &escale, void(*jacobn)(GroupMap &MoleMap, const valarray< double > &b2vec, double *const ervals, double *const amat, const bool lgJac, bool *lgConserve))
Definition: newton_step.cpp:30
hmi
t_hmi hmi
Definition: hmi.cpp:5
t_hmi::HeatH2Dexc_TH85
double HeatH2Dexc_TH85
Definition: hmi.h:138
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
t_mole_global::num_compacted
int num_compacted
Definition: mole.h:314
conv
t_conv conv
Definition: conv.cpp:5
t_dynamics::molecules
double * molecules
Definition: dynamics.h:80
gv
GrainVar gv
Definition: grainvar.cpp:5
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion)
Definition: ionbal.h:254
t_hmi::H2_frac_abund_set
double H2_frac_abund_set
Definition: hmi.h:185
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1079
t_conv::incrementCounter
void incrementCounter(const counter_type type)
Definition: conv.h:308
fixit
void fixit(void)
Definition: service.cpp:991
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
molecule::nAtom
nAtomsMap nAtom
Definition: mole.h:143
lgElemsConserved
bool lgElemsConserved(void)
Definition: dense.cpp:99
phycon.h
atom_list
ChemAtomList atom_list
Definition: mole_species.cpp:71
mole_h_fixup
STATIC void mole_h_fixup(void)
Definition: mole_solve.cpp:276
called.h
mole_priv.h
amat
static double * amat
Definition: atom_feii.cpp:173
t_hmi::h2plus_heat
double h2plus_heat
Definition: hmi.h:39
h2.h
GroupMap::setup
void setup(double *b0vec)
Definition: mole_solve.cpp:565
fnzone
double fnzone
Definition: cddefines.cpp:15
molezone::den
double den
Definition: mole.h:358
t_phycon::te
double te
Definition: phycon.h:11
MAT
#define MAT(a, I_, J_)
Definition: mole_solve.cpp:336
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
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
deuterium.h
chem_atom::frac
double frac
Definition: mole.h:47
t_hmi::H2_total_f
realnum H2_total_f
Definition: hmi.h:17
Generated by   doxygen 1.8.17