00001
00002
00003
00004 #include "cddefines.h"
00005 #include "called.h"
00006 #include "dense.h"
00007 #include "deuterium.h"
00008 #include "ionbal.h"
00009 #include "thermal.h"
00010 #include "phycon.h"
00011 #include "hmi.h"
00012 #include "dynamics.h"
00013 #include "conv.h"
00014 #include "trace.h"
00015 #include "timesc.h"
00016 #include "mole.h"
00017 #include "mole_priv.h"
00018 #include "grainvar.h"
00019 #include "h2.h"
00020 #include "newton_step.h"
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 void check_co_ion_converge(void);
00039
00040 STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec,
00041 double * const ervals, double * const amat, const bool lgJac, bool *lgConserve);
00042
00043 STATIC void mole_h_fixup(void);
00044
00045 STATIC void grouped_elems(const double bvec[], double mole_elems[]);
00046
00047 #define SMALLABUND 1e-24
00048
00049 double mole_solve()
00050 {
00051 int nBad, nPrevBad, i;
00052 realnum
00053 eqerror, error;
00054 valarray<double> oldmols(mole_global.num_compacted),
00055 newmols(mole_global.num_compacted);
00056 GroupMap MoleMap( atom_list.size() );
00057
00058 DEBUG_ENTRY( "mole_solve()" );
00059
00060 if (hmi.H2_frac_abund_set>0.)
00061 {
00062 mole_h_fixup();
00063 fixit();
00064 fprintf(stderr,"Need to treat hmi.H2_frac_abund_set in revised network\n");
00065 fprintf(stderr,"%g\n",hmi.H2_frac_abund_set);
00066 exit(-1);
00067 }
00068
00069 ASSERT(lgElemsConserved());
00070
00071
00072 const double BIGERROR = 1e-8;
00073
00074
00075 const int LIM_LOOP = 100;
00076
00077
00078 nBad = nPrevBad = 0;
00079 eqerror = -1.;
00080
00081 valarray<double> escale(mole_global.num_compacted);
00082 double rlimit=-1., rmax=0.0;
00083
00084
00085 for(i=0; i < LIM_LOOP;i++)
00086 {
00087 {
00088
00089
00090 enum {DEBUG_LOC=false};
00091
00092 if( DEBUG_LOC && (nzone>140) )
00093 {
00094 fprintf(ioQQQ,"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
00095 fnzone,
00096 phycon.te,
00097 dense.eden);
00098 for( long k=0; k<mole_global.num_calc; k++ )
00099 fprintf(ioQQQ,"\t%.2e", mole.species[k].den );
00100 fprintf(ioQQQ,"\n" );
00101 }
00102 }
00103
00104
00105
00106 nPrevBad = nBad;
00107
00108 MoleMap.setup(get_ptr(oldmols));
00109
00110 long j;
00111
00112
00113
00114
00115 conv.incrementCounter(MOLE_SOLVE);
00116 for (j=0; j<30; j++)
00117 {
00118 conv.incrementCounter(MOLE_SOLVE_STEPS);
00119 bool lgOK = newton_step(MoleMap, oldmols, newmols, &eqerror, &error, mole_global.num_compacted, &rlimit, &rmax, escale, funjac);
00120
00121
00122 if (lgOK)
00123 {
00124 nBad = 0;
00125 double fracneg = 0.;
00126 long iworst = -1;
00127 for( long mol=0; mol < mole_global.num_compacted; mol++ )
00128 {
00129 if( newmols[mol] < 0. )
00130 {
00131 if(-newmols[mol] > fracneg*SDIV(oldmols[mol]))
00132 {
00133 fracneg = -newmols[mol]/SDIV(oldmols[mol]);
00134 iworst = mol;
00135 }
00136 if( newmols[mol] < -SMALLABUND)
00137 {
00138 ++nBad;
00139 }
00140 newmols[mol] = 0.;
00141 }
00142 }
00143 if ( 0 != nBad)
00144 {
00145 lgOK = false;
00146 if (0)
00147 {
00148 fprintf(ioQQQ,"-ve density in inner chemistry loop, worst is %ld\n",iworst);
00149 }
00150 }
00151 }
00152 if ( lgOK )
00153 {
00154
00155 break;
00156 }
00157
00158 ASSERT( rlimit < BIGFLOAT );
00159 ASSERT( rlimit > 0.0 );
00160 rlimit = 10.*rlimit;
00161 }
00162
00163
00164
00165
00166
00167 if ( error < 0.01*eqerror )
00168 rlimit = 0.1*rlimit;
00169
00170 MoleMap.updateMolecules( newmols );
00171
00172
00173 if (eqerror < BIGERROR && nBad == 0 && nPrevBad == 0)
00174 break;
00175 }
00176
00177
00178
00179 if( (i == LIM_LOOP && eqerror > BIGERROR) || nBad != 0 )
00180 {
00181 conv.setConvIonizFail("Chem conv", eqerror, nBad);
00182 }
00183
00184
00185
00186 realnum effect_delta = mole_return_cached_species();
00187
00188 realnum delta_threshold = 0.2*conv.IonizErrorAllowed;
00189 if (effect_delta > delta_threshold)
00190 {
00191 conv.setConvIonizFail("chem eft chg", effect_delta, 0.0);
00192 }
00193
00194 if (trace.lgTrace)
00195 {
00196 fprintf(ioQQQ,"Mole zn %3ld -- %7.2f\n",nzone,fnzone);
00197 fprintf(ioQQQ,"Internal error %15.8g nBad %4d loops %3d\n",
00198 eqerror,nBad,i);
00199 fprintf(ioQQQ,"Scaling delta on ions %15.8g; threshold %15.8g\n",
00200 effect_delta, delta_threshold);
00201 }
00202
00203 check_co_ion_converge();
00204
00205 {
00206
00207
00208 enum {DEBUG_LOC=false};
00209
00210 if( DEBUG_LOC )
00211 {
00212 fprintf(ioQQQ,"DEBUG mole out\t%i\t%.2f\t%.5e\t%.5e",
00213 i,
00214 fnzone,
00215 phycon.te,
00216 dense.eden);
00217 fprintf(ioQQQ,
00218 "\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",
00219 eqerror,
00220 dense.xIonDense[ipHYDROGEN][0],
00221 dense.xIonDense[ipHYDROGEN][1],
00222 mole.sink[ipHYDROGEN][0],
00223 mole.sink[ipHYDROGEN][1],
00224 mole.source[ipHYDROGEN][0] ,
00225 mole.source[ipHYDROGEN][1] ,
00226 ionbal.RateIonizTot(ipHYDROGEN,0),
00227 ionbal.RateRecomTot[ipHYDROGEN][0]);
00228
00229
00230
00231 fprintf(ioQQQ,"\n" );
00232 }
00233 }
00234
00235 return eqerror;
00236
00237
00238
00239 }
00240
00241
00242 void check_co_ion_converge(void)
00243 {
00244 DEBUG_ENTRY("check_co_ion_converge()");
00245
00246 if( dense.lgElmtOn[ipCARBON] &&
00247 fabs(dense.xIonDense[ipCARBON][0]- findspecieslocal("C")->den)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
00248 {
00249 conv.setConvIonizFail("CO C0 con",
00250 dense.xIonDense[ipCARBON][0],
00251 findspecieslocal("C")->den);
00252 }
00253 else if( dense.lgElmtOn[ipCARBON] &&
00254 fabs(dense.xIonDense[ipCARBON][1]- findspecieslocal("C+")->den)/SDIV(dense.gas_phase[ipCARBON]) >0.001 )
00255 {
00256 conv.setConvIonizFail("CO C1 con",
00257 dense.xIonDense[ipCARBON][1],
00258 findspecieslocal("C+")->den);
00259 }
00260 else if( dense.lgElmtOn[ipOXYGEN] &&
00261 fabs(dense.xIonDense[ipOXYGEN][0]- findspecieslocal("O")->den)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
00262 {
00263 conv.setConvIonizFail(
00264 "CO O0 con",dense.xIonDense[ipOXYGEN][0],
00265 findspecieslocal("O")->den);
00266 }
00267 else if( dense.lgElmtOn[ipOXYGEN] &&
00268 fabs(dense.xIonDense[ipOXYGEN][1]- findspecieslocal("O+")->den)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 )
00269 {
00270 conv.setConvIonizFail(
00271 "CO O1 con", dense.xIonDense[ipOXYGEN][1],
00272 findspecieslocal("O+")->den);
00273 }
00274 }
00275
00276 STATIC void mole_h_fixup(void)
00277 {
00278 long int mol;
00279
00280
00281
00282
00283
00284
00285 DEBUG_ENTRY( "mole_h_fixup()" );
00286
00287
00288
00289 if( hmi.H2_frac_abund_set>0.)
00290 {
00291 for(mol=0;mol<mole_global.num_calc;mol++)
00292 {
00293 mole.species[mol].den = 0.;
00294 }
00295
00296
00297 dense.xIonDense[ipHYDROGEN][0] = dense.xIonDense[ipHYDROGEN][1] =
00298 2.f*SMALLFLOAT*dense.gas_phase[ipHYDROGEN];
00299
00300 findspecieslocal("H2")->den = (realnum)(dense.gas_phase[ipHYDROGEN] * hmi.H2_frac_abund_set);
00301 findspecieslocal("H2*")->den = 0.;
00302
00303 hmi.H2_total = findspecieslocal("H2")->den + findspecieslocal("H2*")->den;
00304
00305 h2.ortho_density = 0.75*hmi.H2_total;
00306 h2.para_density = 0.25*hmi.H2_total;
00307 {
00308 hmi.H2_total_f = (realnum)hmi.H2_total;
00309 h2.ortho_density_f = (realnum)h2.ortho_density;
00310 h2.para_density_f = (realnum)h2.para_density;
00311 }
00312
00313 hmi.hmihet = 0.;
00314 hmi.h2plus_heat = 0.;
00315 hmi.H2Opacity = 0.;
00316 hmi.hmicol = 0.;
00317 hmi.HeatH2Dish_TH85 = 0.;
00318 hmi.HeatH2Dexc_TH85 = 0.;
00319 hmi.deriv_HeatH2Dexc_TH85 = 0.;
00320 hmi.hmidep = 1.;
00321
00322 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00323 {
00324 gv.bin[nd]->rate_h2_form_grains_used = 0.;
00325 }
00326
00327 return;
00328 }
00329 }
00330
00331 #define ABSLIM 1e-12
00332 #define ERRLIM 1e-12
00333 # ifdef MAT
00334 # undef MAT
00335 # endif
00336 # define MAT(a,I_,J_) ((a)[(I_)*(mole_global.num_compacted)+(J_)])
00337
00338
00339 STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c);
00340 enum {PRINTSOL = false};
00341
00342 STATIC void funjac(GroupMap &MoleMap, const valarray<double> &b2vec, double * const ervals,
00343 double * const amat, const bool lgJac, bool *lgConserve)
00344 {
00345 static multi_arr<double,2> c(mole_global.num_total,mole_global.num_total);
00346 bool printsol = PRINTSOL;
00347 valarray<double> b(mole_global.num_total);
00348
00349 DEBUG_ENTRY( "funjac()" );
00350
00351 MoleMap.updateMolecules( b2vec );
00352
00353 if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
00354 dynamics.Rate != 0.0)
00355 {
00356 ASSERT(dynamics.Rate > 0.0);
00357 *lgConserve = false;
00358 }
00359 else
00360 {
00361 *lgConserve = true;
00362 }
00363
00364
00365
00366
00367 mole_eval_dynamic_balance(mole_global.num_total,get_ptr(b),lgJac,c);
00368
00369
00370 if(printsol || (trace.lgTrace && trace.lgTraceMole ))
00371 {
00372
00373 fprintf( ioQQQ, " ");
00374 for( long i=0; i < mole_global.num_calc; i++ )
00375 {
00376 fprintf( ioQQQ, " %-4.4s", mole_global.list[i]->label.c_str() );
00377 }
00378 fprintf( ioQQQ, " \n" );
00379
00380 fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone);
00381 for( long i=0; i<mole_global.num_calc; i++ )
00382 fprintf(ioQQQ,"\t%.2e", mole.species[i].den );
00383 fprintf(ioQQQ,"\n" );
00384
00385 for( long i=0; i < mole_global.num_calc; i++ )
00386 {
00387 fprintf( ioQQQ, " MOLE%2ld %-4.4s", i ,mole_global.list[i]->label.c_str());
00388 for( long j=0; j < mole_global.num_calc; j++ )
00389 {
00390 fprintf( ioQQQ, "%10.2e", c[j][i] );
00391 }
00392 fprintf( ioQQQ, "%10.2e", b[i] );
00393 fprintf( ioQQQ, "\n" );
00394 }
00395 }
00396
00397
00398
00399
00400 if (lgJac) {
00401 for (unsigned long j=0; j<atom_list.size(); ++j )
00402 {
00403 if (atom_list[j]->ipMl[0] != -1)
00404 {
00405 for(long i=0;i<mole_global.num_calc;i++)
00406 {
00407 ASSERT(mole_global.list[i]->index == i);
00408 double sum = 0.;
00409 for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
00410 {
00411 if (atom_list[j]->ipMl[ion] != -1)
00412 {
00413 sum += MoleMap.fion[j][ion]*c[atom_list[j]->ipMl[ion]][i];
00414 }
00415 }
00416 c[atom_list[j]->ipMl[0]][i] = sum;
00417 }
00418 }
00419 }
00420 for (unsigned long j=0; j<atom_list.size(); ++j )
00421 {
00422 if (atom_list[j]->ipMl[0] != -1)
00423 {
00424 for(long i=0;i<mole_global.num_calc;i++)
00425 {
00426 double sum = 0.0;
00427
00428 for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
00429 {
00430 if (atom_list[j]->ipMl[ion] != -1)
00431 {
00432 sum += c[i][atom_list[j]->ipMl[ion]];
00433 }
00434 }
00435 c[i][atom_list[j]->ipMl[0]] = sum;
00436 }
00437 }
00438 }
00439 }
00440
00441 for (unsigned long j=0; j<atom_list.size(); ++j )
00442 {
00443
00444
00445
00446 if (atom_list[j]->ipMl[0] != -1)
00447 {
00448 double sum = 0.0;
00449 for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
00450 {
00451 if (atom_list[j]->ipMl[ion] != -1)
00452 {
00453
00454
00455
00456 sum += b[atom_list[j]->ipMl[ion]];
00457 b[atom_list[j]->ipMl[ion]] = 0.0;
00458 }
00459 }
00460 b[atom_list[j]->ipMl[0]] = sum;
00461 }
00462 }
00463
00464
00465
00466 if (*lgConserve)
00467 {
00468
00469 if (lgJac)
00470 {
00471 for ( unsigned long j = 0; j < atom_list.size(); ++j )
00472 {
00473 long ncons = atom_list[j]->ipMl[0];
00474 if (ncons != -1)
00475 {
00476 ASSERT( mole_global.list[ncons]->isMonatomic() );
00477 double scale=1.0;
00478 for( long i=0;i<mole_global.num_compacted;i++)
00479 {
00480 if( groupspecies[i]->nAtom.find(atom_list[j]) != groupspecies[i]->nAtom.end() )
00481 c[groupspecies[i]->index][ncons] = groupspecies[i]->nAtom[atom_list[j]]*scale;
00482 else
00483 c[groupspecies[i]->index][ncons] = 0.;
00484
00485 }
00486 }
00487 }
00488 }
00489
00490 valarray<double> molnow(atom_list.size());
00491 grouped_elems(get_ptr(b2vec), get_ptr(molnow));
00492 for ( unsigned long j = 0; j < atom_list.size(); ++j )
00493 {
00494 long ncons = atom_list[j]->ipMl[0];
00495 if (ncons != -1)
00496 {
00497 ASSERT( mole_global.list[ncons]->isMonatomic() );
00498 double scale = c[ncons][ncons];
00499 b[ncons] = (molnow[j]-MoleMap.molElems[j])*scale;
00500 if (false)
00501 if (b[ncons] != 0.0)
00502 fprintf(ioQQQ,"Cons %s err %g rel %g\n",
00503 atom_list[j]->label().c_str(),b[ncons],
00504 (molnow[j]-MoleMap.molElems[j])/SDIV(MoleMap.molElems[j]));
00505 }
00506 }
00507 }
00508
00509 for( long i=0; i < mole_global.num_compacted; i++ )
00510 {
00511 ervals[i] = b[groupspecies[i]->index];
00512 }
00513
00514 if (lgJac)
00515 {
00516
00517 for( long i=0; i < mole_global.num_compacted; i++ )
00518 {
00519 for( long j=0; j < mole_global.num_compacted; j++ )
00520 {
00521 MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
00522 }
00523 }
00524
00525
00526
00527 for(long i=0;i<mole_global.num_compacted;i++)
00528 {
00529 double sum1 = 0.;
00530 for (long j=0;j<mole_global.num_compacted;j++)
00531 {
00532 sum1 += fabs(MAT(amat,j,i));
00533 }
00534 if (sum1 == 0.0)
00535 {
00536 ASSERT(ervals[i] == 0.0);
00537 MAT(amat,i,i) = 1.0;
00538
00539 }
00540 }
00541 }
00542 }
00543
00544 STATIC void grouped_elems(const double bvec[], double mole_elems[])
00545 {
00546 map<chem_atom*, long> atom_to_index;
00547 for (unsigned long j=0; j<atom_list.size(); ++j )
00548 {
00549 mole_elems[j] = 0.;
00550 atom_to_index[atom_list[j].get_ptr()] = j;
00551 }
00552 for( long i=0; i < mole_global.num_compacted; i++ )
00553 {
00554 if( groupspecies[i]->parentLabel.empty() == false )
00555 continue;
00556
00557 for (molecule::nAtomsMap::const_iterator el = groupspecies[i]->nAtom.begin();
00558 el != groupspecies[i]->nAtom.end(); ++el)
00559 {
00560 mole_elems[atom_to_index[el->first.get_ptr()]] += bvec[i]*el->second;
00561 }
00562 }
00563 }
00564
00565 void GroupMap::setup(double *b0vec)
00566 {
00567 valarray<double> calcv(mole_global.num_total);
00568 bool lgSet;
00569
00570 for( long i=0;i<mole_global.num_total;i++)
00571 {
00572 calcv[i] = mole.species[i].den;
00573 }
00574
00575 for (unsigned long j=0; j<atom_list.size(); ++j )
00576 {
00577 if (atom_list[j]->ipMl[0] != -1)
00578 {
00579 double sum = 0.;
00580 for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
00581 {
00582 if (atom_list[j]->ipMl[ion] != -1)
00583 sum += calcv[atom_list[j]->ipMl[ion]];
00584 }
00585 if (sum > SMALLFLOAT)
00586 {
00587 double factor = 1./sum;
00588 for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
00589 {
00590 if (atom_list[j]->ipMl[ion] != -1)
00591 fion[j][ion] = calcv[atom_list[j]->ipMl[ion]]*factor;
00592 else
00593 fion[j][ion] = 0.;
00594 }
00595 }
00596 else
00597 {
00598 lgSet = false;
00599 for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
00600 {
00601 if (atom_list[j]->ipMl[ion] != -1 && !lgSet)
00602 {
00603 fion[j][ion] = 1.0;
00604 lgSet = true;
00605 }
00606 else
00607 {
00608 fion[j][ion] = 0.;
00609 }
00610 }
00611 }
00612
00613 lgSet = false;
00614 for (unsigned long ion=0; ion<atom_list[j]->ipMl.size(); ion++)
00615 {
00616 if (atom_list[j]->ipMl[ion] != -1)
00617 {
00618 if (!lgSet)
00619 calcv[atom_list[j]->ipMl[ion]] = sum;
00620 else
00621 calcv[atom_list[j]->ipMl[ion]] = 0.;
00622 lgSet = true;
00623 }
00624 }
00625 }
00626 }
00627
00628 for( long i=0; i < mole_global.num_compacted; i++ )
00629 {
00630 b0vec[i] = calcv[groupspecies[i]->index];
00631 }
00632
00633 grouped_elems(get_ptr(b0vec),get_ptr(molElems));
00634
00635 for (unsigned long i = 0; i<atom_list.size(); ++i)
00636 {
00637 double densAtom = 0.;
00638
00639 if( atom_list[i]->el->Z==1 && atom_list[i]->A==2 )
00640 {
00641 ASSERT( deut.lgElmtOn );
00642 densAtom = deut.gas_phase;
00643 }
00644
00645 else if( !atom_list[i]->lgMeanAbundance() )
00646 continue;
00647 else
00648 {
00649 int nelem = atom_list[i]->el->Z-1;
00650 densAtom = dense.gas_phase[nelem];
00651 }
00652 bool lgTest =
00653 ( densAtom < SMALLABUND && molElems[i] < SMALLABUND ) ||
00654 ( fabs(molElems[i]- densAtom) <= conv.GasPhaseAbundErrorAllowed*densAtom );
00655 if( !lgTest )
00656 fprintf( ioQQQ, "PROBLEM molElems BAD %li\t%s\t%.12e\t%.12e\t%.12e\n",
00657 i, atom_list[i]->label().c_str(), atom_list[i]->frac, densAtom, molElems[i] );
00658
00659 ASSERT( lgTest );
00660 molElems[i] = densAtom;
00661 }
00662
00663 }
00664
00665 void GroupMap::updateMolecules(const valarray<double> & b2 )
00666 {
00667 DEBUG_ENTRY("updateMolecules()");
00668
00669 for (long mol=0;mol<mole_global.num_calc;mol++)
00670 {
00671 mole.species[mol].den = 0.;
00672 }
00673 for (long mol=0;mol<mole_global.num_compacted;mol++)
00674 {
00675 mole.species[ groupspecies[mol]->index ].den = b2[mol];
00676 }
00677
00678
00679 for (long mol=0;mol<mole_global.num_calc;mol++)
00680 {
00681 if( mole_global.list[mol]->parentIndex >= 0 )
00682 {
00683 ASSERT( !mole_global.list[mol]->parentLabel.empty() );
00684 long parentIndex = mole_global.list[mol]->parentIndex;
00685 mole.species[mol].den = mole.species[parentIndex].den;
00686 for( nAtoms_i it = mole_global.list[mol]->nAtom.begin(); it != mole_global.list[mol]->nAtom.end(); ++it )
00687 {
00688 if( !it->first->lgMeanAbundance() )
00689 {
00690 mole.species[mol].den *= pow( it->first->frac, it->second );
00691 }
00692 }
00693 }
00694 }
00695
00696
00697 for (unsigned long j=0; j<atom_list.size(); ++j )
00698 {
00699 if (atom_list[j]->ipMl[0] != -1)
00700 {
00701 double grouptot = mole.species[atom_list[j]->ipMl[0]].den;
00702 double sum = 0.0;
00703 for (unsigned long ion=0;ion<atom_list[j]->ipMl.size();ion++)
00704 {
00705 if (atom_list[j]->ipMl[ion] != -1)
00706 {
00707 mole.species[atom_list[j]->ipMl[ion]].den = grouptot * fion[j][ion];
00708 sum += mole.species[atom_list[j]->ipMl[ion]].den;
00709 }
00710 }
00711 ASSERT(fabs(sum-grouptot) <= 1e-10 * fabs(grouptot));
00712 }
00713 }
00714
00715 mole.set_isotope_abundances();
00716
00717 return;
00718 }
00719
00720 STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c)
00721 {
00722 double source;
00723
00724 DEBUG_ENTRY("mole_eval_dynamic_balance()");
00725
00726 mole_eval_balance(num_total, b, lgJac, c);
00727
00728
00729 if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
00730 dynamics.Rate != 0.0 )
00731 {
00732
00733
00734 for( long i=0;i<mole_global.num_calc;i++)
00735 {
00736 if (lgJac)
00737 c[i][i] -= dynamics.Rate;
00738
00739 if( !mole_global.list[i]->parentLabel.empty() )
00740 continue;
00741
00742 b[i] -= mole.species[i].den*dynamics.Rate;
00743
00744 if (!mole_global.list[i]->isMonatomic() || mole_global.list[i]->charge < 0 || ! mole_global.list[i]->lgGas_Phase ||
00745 ( mole_global.list[i]->isMonatomic() && !mole_global.list[i]->nAtom.begin()->first->lgMeanAbundance() ) )
00746 {
00747 b[i] += dynamics.molecules[i]*dynamics.Rate;
00748 }
00749 else if (mole_global.list[i]->charge == 0)
00750 {
00751 ASSERT( mole_global.list[i]->isMonatomic() );
00752 ASSERT( (int)mole_global.list[i]->nAtom.size() == 1 );
00753 count_ptr<chem_atom> atom = mole_global.list[i]->nAtom.begin()->first;
00754 long nelem = atom->el->Z-1;
00755 if( nelem >= LIMELM )
00756 continue;
00757 source = 0.0;
00758 for ( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00759 {
00760 source += dynamics.Source[nelem][ion] * atom->frac;
00761 if (unresolved_atom_list[nelem]->ipMl[ion] == -1)
00762 source -= dense.xIonDense[nelem][ion]*dynamics.Rate * atom->frac;
00763 }
00764 b[i] += source;
00765 }
00766 }
00767 }
00768 }