46 for( i=0; i < num_total; i++ )
57 fprintf(
ioQQQ,
" TrChemSp %20.13g %s\n",
mole.species[mol].den,
65 fprintf(
ioQQQ,
" TrChem %20.13g %s\n",rk,rate->
label.c_str());
84 fprintf(
ioQQQ,
" TrChemTot %20.13g %20.13g %s\n",rate_tot,rk,rate->
label.c_str());
96 b[sp->
index] -= rate_tot;
105 b[sp->
index] += rate_tot;
125 const double rated = rate_deriv[j];
151 long int i, j, nelem, ion, ion2;
158 for( i=0; i < num_total; i++ )
160 mole.species[i].src =
mole.species[i].snk = 0.;
163 for( nelem=0; nelem<
LIMELM; ++nelem )
166 for( ion=0; ion<nelem+2; ++ion )
169 for( ion2=0; ion2<nelem+2; ++ion2 )
171 mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
179 rate = &(*p->second);
201 mole.species[sp->
index].snk += rate_deriv[i];
208 long otherIndex = 1-i;
218 nelem = (*atom)->el->Z-1;
243 const long int nelem=(*atom)->el->Z-1;
244 if( !
dense.lgElmtOn[nelem] )
247 for (
long int ion=0;ion<nelem+2;ion++)
249 if ((*atom)->ipMl[ion] != -1)
251 mole.source[nelem][ion] =
mole.species[(*atom)->ipMl[ion]].src;
252 mole.sink[nelem][ion] =
mole.species[(*atom)->ipMl[ion]].snk;
256 mole.source[nelem][ion] = 0.0;
257 mole.sink[nelem][ion] = 0.0;
269 bool checkAllOK =
true;
277 map<chem_atom*, long> atom_to_index;
278 for (
unsigned long j=0; j<
atom_list.size(); ++j )
280 atom_to_index[
atom_list[j].get_ptr()] = j;
289 ccache[nc] = c[i][j];
296 for (molecule::nAtomsMap::const_iterator el =
mole_global.list[j]->nAtom.begin();
299 long natom = atom_to_index[el->first.get_ptr()];
300 const int nAtomj = el->second;
301 for (
long i=0;i<nc;i++)
303 const double term = ccache[i] * nAtomj;
304 test[natom][ncache[i]] += term;
305 tot[natom][ncache[i]] += fabs(term);
311 for(
unsigned long natom=0; natom <
atom_list.size(); ++natom)
316 ( fabs(test[natom][i]) <=
MAX2(1e-10*tot[natom][i], 1e10*DBL_MIN) );
320 fprintf(stdout,
"Network conservation error %s %s %g %g %g %g\n",
321 atom->
label().c_str(),
324 test[natom][i]/tot[natom][i],
341 double snkx=0.,srcx=0.;
346 rate = &(*p->second);
368 if (sp == debug_species && rate->
pvector[i] == NULL)
370 if (fabs(rate_tot) > srcx)
380 if (sp == debug_species && rate->
rvector[i] == NULL)
382 if (fabs(rate_deriv[i]) > snkx)
384 snkx = rate_deriv[i];
397 fprintf( ioOut,
"%20.20s src %13.7g of %13.7g [",
398 ratesrc->label.c_str(),srcx,
mole.species[debug_species->
index].src);
399 for (j=0;j<ratesrc->nreactants;j++)
403 fprintf( ioOut,
"," );
405 fprintf( ioOut,
"%-6.6s %13.7g",
406 ratesrc->reactants[j]->label.c_str(),
407 mole.species[ ratesrc->reactants[j]->index ].den);
409 fprintf( ioOut,
"]" );
413 fprintf( ioOut,
"%20.20s snk %13.7g of %13.7g [",
414 ratesnk->
label.c_str(), snkx *
mole.species[debug_species->
index].den,
420 fprintf( ioOut,
"," );
422 fprintf( ioOut,
"%-6.6s %13.7g",
426 fprintf( ioOut,
"]" );
429 fprintf( ioOut,
"\n" );
#define DEBUG_ENTRY(funcname)
molecule * pvector[MAXPRODUCTS]
molecule * reactants[MAXREACTANTS]
molecule * rvector[MAXREACTANTS]
molecule * products[MAXPRODUCTS]
bool isMonatomic(void) const
t_mole_global mole_global
ChemAtomList unresolved_atom_list
void mole_eval_sources(long int num_total)
void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
void mole_dominant_rates(const molecule *debug_species, FILE *ioOut)
STATIC bool lgNucleiConserved(const multi_arr< double, 2 > &c)
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
map< string, count_ptr< mole_reaction > > reactab