cloudy trunk
Loading...
Searching...
No Matches
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 * */
38void check_co_ion_converge(void);
39
40STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec,
41 double * const ervals, double * const amat, const bool lgJac, bool *lgConserve);
42
43STATIC void mole_h_fixup(void);
44
45STATIC void grouped_elems(const double bvec[], double mole_elems[]);
46
47#define SMALLABUND 1e-24
48
49double mole_solve()
50{
51 int nBad, nPrevBad, i;
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 {
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
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
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
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
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
339STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c);
340enum {PRINTSOL = false};
341
342STATIC 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
544STATIC 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
565void 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
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
665void 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
720STATIC 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}
static double b2[63]
static double * amat
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipCARBON
Definition cddefines.h:310
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
float realnum
Definition cddefines.h:103
T * get_ptr(T *v)
Definition cddefines.h:1079
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
valarray< double > molElems
Definition mole_priv.h:24
multi_arr< double, 2 > fion
Definition mole_priv.h:23
void updateMolecules(const valarray< double > &b2)
void setup(double *b0vec)
double den
Definition mole.h:358
t_conv conv
Definition conv.cpp:5
@ MOLE_SOLVE
Definition conv.h:71
@ MOLE_SOLVE_STEPS
Definition conv.h:72
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const realnum SMALLFLOAT
Definition cpu.h:191
bool lgElemsConserved(void)
Definition dense.cpp:99
t_dense dense
Definition dense.cpp:24
t_deuterium deut
Definition deuterium.cpp:8
t_dynamics dynamics
Definition dynamics.cpp:44
GrainVar gv
Definition grainvar.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hmi hmi
Definition hmi.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
#define SMALLABUND
Definition mole.h:13
ChemAtomList atom_list
molezone * findspecieslocal(const char buf[])
ChemAtomList unresolved_atom_list
molecule::nAtomsMap::iterator nAtoms_i
Definition mole.h:214
void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
molecule ** groupspecies
realnum mole_return_cached_species(const GroupMap &MoleMap)
#define MAT(a, I_, J_)
void check_co_ion_converge(void)
@ PRINTSOL
STATIC void funjac(GroupMap &MoleMap, const valarray< double > &b2vec, double *const ervals, double *const amat, const bool lgJac, bool *lgConserve)
double mole_solve()
STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
STATIC void grouped_elems(const double bvec[], double mole_elems[])
STATIC void mole_h_fixup(void)
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))
t_phycon phycon
Definition phycon.cpp:6
static double * source
Definition species2.cpp:28
t_trace trace
Definition trace.cpp:5