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