00001
00002
00003
00004 #include "cddefines.h"
00005 #include "mole.h"
00006 #include "mole_priv.h"
00007 #include "save.h"
00008 #include "dense.h"
00009 #include "atmdat.h"
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 STATIC bool lgNucleiConserved(const multi_arr<double,2> &c);
00035
00036 void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c)
00037 {
00038 long int i, j;
00039 mole_reaction *rate;
00040 double rate_tot, rate_deriv[MAXREACTANTS], rk;
00041 molecule *sp;
00042
00043 DEBUG_ENTRY("mole_eval_balance()");
00044
00045 for( i=0; i < num_total; i++ )
00046 {
00047 b[i] = 0.;
00048 }
00049 if (lgJac)
00050 c.zero();
00051
00052 for(mole_reaction_i p=mole_priv::reactab.begin();
00053 p != mole_priv::reactab.end(); ++p)
00054 {
00055 rate = &(*p->second);
00056 rk = mole.reaction_rks[ rate->index ];
00057
00058 rate_tot = rk;
00059 for(i=0;i<rate->nreactants;i++)
00060 {
00061 rate_tot *= mole.species[ rate->reactants[i]->index ].den;
00062 }
00063
00064 for(i=0;i<rate->nreactants;i++)
00065 {
00066 sp = rate->reactants[i];
00067 if (rate->rvector[i] == NULL)
00068 {
00069 b[sp->index] -= rate_tot;
00070 }
00071 }
00072
00073 for(i=0;i<rate->nproducts;i++)
00074 {
00075 sp = rate->products[i];
00076 if (rate->pvector[i] == NULL)
00077 {
00078 b[sp->index] += rate_tot;
00079 }
00080 }
00081
00082 if (lgJac)
00083 {
00084 for(i=0;i<rate->nreactants;i++)
00085 {
00086 rate_deriv[i] = rk;
00087 for(j=0;j<rate->nreactants;j++)
00088 {
00089 if(i!=j)
00090 {
00091 rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
00092 }
00093 }
00094 }
00095 for(j=0;j<rate->nreactants;j++)
00096 {
00097 sp = rate->reactants[j];
00098 const double rated = rate_deriv[j];
00099 for(i=0;i<rate->nreactants;i++)
00100 {
00101 if (rate->rvector[i] == NULL)
00102 c[sp->index][rate->reactants[i]->index] -= rated;
00103 }
00104 for(i=0;i<rate->nproducts;i++)
00105 {
00106 if (rate->pvector[i] == NULL)
00107 c[sp->index][rate->products[i]->index] += rated;
00108 }
00109 }
00110 }
00111 }
00112
00113 if (lgJac)
00114 {
00115 ASSERT( lgNucleiConserved(c) );
00116 }
00117
00118
00119 return;
00120 }
00121
00122 void mole_eval_sources(long int num_total)
00123 {
00124 long int i, j, nelem, ion, ion2;
00125 mole_reaction *rate;
00126 double rate_tot, rate_deriv[MAXREACTANTS], rk;
00127 molecule *sp;
00128
00129 DEBUG_ENTRY("mole_eval_sources()");
00130
00131 for( i=0; i < num_total; i++ )
00132 {
00133 mole.species[i].src = mole.species[i].snk = 0.;
00134 }
00135
00136 for( nelem=0; nelem< LIMELM; ++nelem )
00137 {
00138
00139 for( ion=0; ion<nelem+2; ++ion )
00140 {
00141
00142 for( ion2=0; ion2<nelem+2; ++ion2 )
00143 {
00144 mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
00145 }
00146 }
00147 }
00148
00149 for(mole_reaction_i p=mole_priv::reactab.begin();
00150 p != mole_priv::reactab.end(); ++p)
00151 {
00152 rate = &(*p->second);
00153 rk = mole.reaction_rks[ rate->index ];
00154
00155 for(i=0;i<rate->nreactants;i++)
00156 {
00157 rate_deriv[i] = rk;
00158 for(j=0;j<rate->nreactants;j++)
00159 {
00160 if(i!=j)
00161 {
00162 rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
00163 }
00164 }
00165 }
00166
00167 rate_tot = rate_deriv[0] * mole.species[ rate->reactants[0]->index ].den;
00168
00169 for(i=0;i<rate->nreactants;i++)
00170 {
00171 sp = rate->reactants[i];
00172 if (rate->rvector[i] == NULL)
00173 {
00174 mole.species[sp->index].snk += rate_deriv[i];
00175 }
00176 else
00177 {
00178
00179 if( rate->nreactants==2 )
00180 {
00181 long otherIndex = 1-i;
00182 molecule *sp2 = rate->reactants[otherIndex];
00183 if( sp->isMonatomic() && sp2->isMonatomic() && sp->nAtom.begin()->first->el == sp2->nAtom.begin()->first->el )
00184 continue;
00185 }
00186
00187 if ( atmdat.lgCTOn )
00188 {
00189 for( ChemAtomList::iterator atom = unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
00190 {
00191 nelem = (*atom)->el->Z-1;
00192 if (sp->nAtom.find(*atom) != sp->nAtom.end() && sp->nAtom[*atom] != 0 && rate->rvector[i]->charge != sp->charge)
00193 {
00194 mole.xMoleChTrRate[nelem][sp->charge][rate->rvector[i]->charge] +=
00195 (realnum) rate_deriv[i];
00196 break;
00197 }
00198 }
00199 }
00200 }
00201 }
00202
00203 for(i=0;i<rate->nproducts;i++)
00204 {
00205 sp = rate->products[i];
00206 if (rate->pvector[i] == NULL)
00207 {
00208 mole.species[sp->index].src += rate_tot;
00209 }
00210 }
00211 }
00212
00213 for (ChemAtomList::iterator atom = unresolved_atom_list.begin();
00214 atom != unresolved_atom_list.end(); ++atom)
00215 {
00216 const long int nelem=(*atom)->el->Z-1;
00217 if( !dense.lgElmtOn[nelem] )
00218 continue;
00219
00220 for (long int ion=0;ion<nelem+2;ion++)
00221 {
00222 if ((*atom)->ipMl[ion] != -1)
00223 {
00224 mole.source[nelem][ion] = mole.species[(*atom)->ipMl[ion]].src;
00225 mole.sink[nelem][ion] = mole.species[(*atom)->ipMl[ion]].snk;
00226 }
00227 else
00228 {
00229 mole.source[nelem][ion] = 0.0;
00230 mole.sink[nelem][ion] = 0.0;
00231 }
00232 }
00233 }
00234
00235
00236 return;
00237 }
00238
00239 STATIC bool lgNucleiConserved(const multi_arr<double,2> &c)
00240 {
00241 DEBUG_ENTRY ("lgNucleiConserved()");
00242 bool checkAllOK = true;
00243 multi_arr<double,2> test(atom_list.size(),mole_global.num_calc),
00244 tot(atom_list.size(),mole_global.num_calc);
00245 vector<double> ccache(mole_global.num_calc);
00246 vector<long> ncache(mole_global.num_calc);
00247 test.zero();
00248 tot.zero();
00249
00250 map<chem_atom*, long> atom_to_index;
00251 for (unsigned long j=0; j<atom_list.size(); ++j )
00252 {
00253 atom_to_index[atom_list[j].get_ptr()] = j;
00254 }
00255 for (long j=0;j<mole_global.num_calc;j++)
00256 {
00257 long nc = 0;
00258 for (long i=0;i<mole_global.num_calc;i++)
00259 {
00260 if (c[i][j] != 0.)
00261 {
00262 ccache[nc] = c[i][j];
00263 ncache[nc] = i;
00264 ++nc;
00265 }
00266 }
00267 if (nc > 0)
00268 {
00269 for (molecule::nAtomsMap::const_iterator el = mole_global.list[j]->nAtom.begin();
00270 el != mole_global.list[j]->nAtom.end(); ++el)
00271 {
00272 long natom = atom_to_index[el->first.get_ptr()];
00273 const int nAtomj = el->second;
00274 for (long i=0;i<nc;i++)
00275 {
00276 const double term = ccache[i] * nAtomj;
00277 test[natom][ncache[i]] += term;
00278 tot[natom][ncache[i]] += fabs(term);
00279 }
00280 }
00281 }
00282 }
00283
00284 for( unsigned long natom=0; natom < atom_list.size(); ++natom)
00285 {
00286 for (long i=0;i<mole_global.num_calc;i++)
00287 {
00288 const bool checkOK =
00289 ( fabs(test[natom][i]) <= MAX2(1e-10*tot[natom][i], 1e10*DBL_MIN) );
00290 if ( UNLIKELY(!checkOK) )
00291 {
00292 chem_atom *atom = atom_list[natom].get_ptr();
00293 fprintf(stdout,"Network conservation error %s %s %g %g %g %g\n",
00294 atom->label().c_str(),
00295 mole_global.list[i]->label.c_str(),
00296 test[natom][i],
00297 test[natom][i]/tot[natom][i],
00298 mole.species[atom->ipMl[0]].den,
00299 mole.species[atom->ipMl[1]].den);
00300
00301 checkAllOK = false;
00302 }
00303 }
00304 }
00305 return checkAllOK;
00306 }
00307
00308 void mole_dominant_rates( const molecule *debug_species, FILE *ioOut )
00309 {
00310 long int i, j;
00311 mole_reaction *rate, *ratesnk=NULL, *ratesrc=NULL;
00312 double rate_tot, rate_deriv[MAXREACTANTS], rk;
00313 molecule *sp;
00314 double snkx=0.,srcx=0.;
00315
00316 for(mole_reaction_i p=mole_priv::reactab.begin();
00317 p != mole_priv::reactab.end(); ++p)
00318 {
00319 rate = &(*p->second);
00320 rk = mole.reaction_rks[ rate->index ];
00321
00322 for(i=0;i<rate->nreactants;i++)
00323 {
00324 rate_deriv[i] = rk;
00325 for(j=0;j<rate->nreactants;j++)
00326 {
00327 if(i!=j)
00328 {
00329 rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
00330 }
00331 }
00332 }
00333
00334 rate_tot = rate_deriv[0] * mole.species[ rate->reactants[0]->index ].den;
00335
00336 if (debug_species != null_mole)
00337 {
00338 for(i=0;i<rate->nproducts;++i)
00339 {
00340 sp = rate->products[i];
00341 if (sp == debug_species && rate->pvector[i] == NULL)
00342 {
00343 if (fabs(rate_tot) > srcx)
00344 {
00345 srcx = rate_tot;
00346 ratesrc = rate;
00347 }
00348 }
00349 }
00350 for(i=0;i<rate->nreactants;++i)
00351 {
00352 sp = rate->reactants[i];
00353 if (sp == debug_species && rate->rvector[i] == NULL)
00354 {
00355 if (fabs(rate_deriv[i]) > snkx)
00356 {
00357 snkx = rate_deriv[i];
00358 ratesnk = rate;
00359 }
00360 }
00361 }
00362 }
00363 }
00364
00365
00366 if (debug_species != null_mole)
00367 {
00368 if (ratesrc)
00369 {
00370 fprintf( ioOut, "%20.20s src %13.7g of %13.7g [",
00371 ratesrc->label.c_str(),srcx,mole.species[debug_species->index].src);
00372 for (j=0;j<ratesrc->nreactants;j++)
00373 {
00374 if (j)
00375 {
00376 fprintf( ioOut, "," );
00377 }
00378 fprintf( ioOut, "%-6.6s %13.7g",
00379 ratesrc->reactants[j]->label.c_str(),
00380 mole.species[ ratesrc->reactants[j]->index ].den);
00381 }
00382 fprintf( ioOut, "]" );
00383 }
00384 if (ratesnk)
00385 {
00386 fprintf( ioOut, "%20.20s snk %13.7g of %13.7g [",
00387 ratesnk->label.c_str(), snkx * mole.species[debug_species->index].den,
00388 mole.species[debug_species->index].snk * mole.species[debug_species->index].den);
00389 for (j=0;j<ratesnk->nreactants;j++)
00390 {
00391 if (j)
00392 {
00393 fprintf( ioOut, "," );
00394 }
00395 fprintf( ioOut, "%-6.6s %13.7g",
00396 ratesnk->reactants[j]->label.c_str(),
00397 mole.species[ ratesnk->reactants[j]->index ].den);
00398 }
00399 fprintf( ioOut, "]" );
00400 }
00401 }
00402 fprintf( ioOut, "\n" );
00403
00404 return;
00405 }