41 double *
const ervals,
double *
const amat,
const bool lgJac,
bool *lgConserve);
47#define SMALLABUND 1e-24
51 int nBad, nPrevBad, i;
54 valarray<double> oldmols(
mole_global.num_compacted),
60 if (
hmi.H2_frac_abund_set>0.)
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);
72 const double BIGERROR = 1e-8;
75 const int LIM_LOOP = 100;
82 double rlimit=-1., rmax=0.0;
85 for(i=0; i < LIM_LOOP;i++)
90 enum {DEBUG_LOC=
false};
92 if( DEBUG_LOC && (
nzone>140) )
94 fprintf(
ioQQQ,
"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
99 fprintf(
ioQQQ,
"\t%.2e",
mole.species[k].den );
100 fprintf(
ioQQQ,
"\n" );
127 for(
long mol=0; mol <
mole_global.num_compacted; mol++ )
129 if( newmols[mol] < 0. )
131 if(-newmols[mol] > fracneg*
SDIV(oldmols[mol]))
133 fracneg = -newmols[mol]/
SDIV(oldmols[mol]);
148 fprintf(
ioQQQ,
"-ve density in inner chemistry loop, worst is %ld\n",iworst);
167 if ( error < 0.01*eqerror )
173 if (eqerror < BIGERROR && nBad == 0 && nPrevBad == 0)
179 if( (i == LIM_LOOP && eqerror > BIGERROR) || nBad != 0 )
181 conv.setConvIonizFail(
"Chem conv", eqerror, nBad);
188 realnum delta_threshold = 0.2*
conv.IonizErrorAllowed;
189 if (effect_delta > delta_threshold)
191 conv.setConvIonizFail(
"chem eft chg", effect_delta, 0.0);
197 fprintf(
ioQQQ,
"Internal error %15.8g nBad %4d loops %3d\n",
199 fprintf(
ioQQQ,
"Scaling delta on ions %15.8g; threshold %15.8g\n",
200 effect_delta, delta_threshold);
208 enum {DEBUG_LOC=
false};
212 fprintf(
ioQQQ,
"DEBUG mole out\t%i\t%.2f\t%.5e\t%.5e",
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",
231 fprintf(
ioQQQ,
"\n" );
249 conv.setConvIonizFail(
"CO C0 con",
256 conv.setConvIonizFail(
"CO C1 con",
263 conv.setConvIonizFail(
270 conv.setConvIonizFail(
289 if(
hmi.H2_frac_abund_set>0.)
293 mole.species[mol].den = 0.;
305 h2.ortho_density = 0.75*
hmi.H2_total;
306 h2.para_density = 0.25*
hmi.H2_total;
314 hmi.h2plus_heat = 0.;
317 hmi.HeatH2Dish_TH85 = 0.;
318 hmi.HeatH2Dexc_TH85 = 0.;
319 hmi.deriv_HeatH2Dexc_TH85 = 0.;
322 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
324 gv.bin[nd]->rate_h2_form_grains_used = 0.;
336# define MAT(a,I_,J_) ((a)[(I_)*(mole_global.num_compacted)+(J_)])
343 double *
const amat,
const bool lgJac,
bool *lgConserve)
370 if(printsol || (
trace.lgTrace &&
trace.lgTraceMole ))
373 fprintf(
ioQQQ,
" ");
378 fprintf(
ioQQQ,
" \n" );
380 fprintf(
ioQQQ,
" MOLE old abundances\t%.2f",
fnzone);
382 fprintf(
ioQQQ,
"\t%.2e",
mole.species[i].den );
383 fprintf(
ioQQQ,
"\n" );
390 fprintf(
ioQQQ,
"%10.2e", c[j][i] );
392 fprintf(
ioQQQ,
"%10.2e", b[i] );
393 fprintf(
ioQQQ,
"\n" );
401 for (
unsigned long j=0; j<
atom_list.size(); ++j )
409 for (
unsigned long ion=0;ion<
atom_list[j]->ipMl.size();ion++)
420 for (
unsigned long j=0; j<
atom_list.size(); ++j )
428 for (
unsigned long ion=0;ion<
atom_list[j]->ipMl.size();ion++)
441 for (
unsigned long j=0; j<
atom_list.size(); ++j )
449 for (
unsigned long ion=0;ion<
atom_list[j]->ipMl.size();ion++)
471 for (
unsigned long j = 0; j <
atom_list.size(); ++j )
477 double scale=1.0/MoleMap.
molElems[j];
490 valarray<double> molnow(
atom_list.size());
492 for (
unsigned long j = 0; j <
atom_list.size(); ++j )
498 double scale = c[ncons][ncons];
499 b[ncons] = (molnow[j]-MoleMap.
molElems[j])*scale;
502 fprintf(
ioQQQ,
"Cons %s err %g rel %g\n",
509 for(
long i=0; i <
mole_global.num_compacted; i++ )
517 for(
long i=0; i <
mole_global.num_compacted; i++ )
519 for(
long j=0; j <
mole_global.num_compacted; j++ )
546 map<chem_atom*, long> atom_to_index;
547 for (
unsigned long j=0; j<
atom_list.size(); ++j )
550 atom_to_index[
atom_list[j].get_ptr()] = j;
552 for(
long i=0; i <
mole_global.num_compacted; i++ )
557 for (molecule::nAtomsMap::const_iterator el =
groupspecies[i]->nAtom.begin();
560 mole_elems[atom_to_index[el->first.get_ptr()]] += bvec[i]*el->second;
572 calcv[i] =
mole.species[i].den;
575 for (
unsigned long j=0; j<
atom_list.size(); ++j )
580 for (
unsigned long ion=0; ion<
atom_list[j]->ipMl.size(); ion++)
587 double factor = 1./sum;
588 for (
unsigned long ion=0; ion<
atom_list[j]->ipMl.size(); ion++)
599 for (
unsigned long ion=0; ion<
atom_list[j]->ipMl.size(); ion++)
601 if (
atom_list[j]->ipMl[ion] != -1 && !lgSet)
614 for (
unsigned long ion=0; ion<
atom_list[j]->ipMl.size(); ion++)
628 for(
long i=0; i <
mole_global.num_compacted; i++ )
635 for (
unsigned long i = 0; i<
atom_list.size(); ++i)
637 double densAtom = 0.;
642 densAtom =
deut.gas_phase;
645 else if( !
atom_list[i]->lgMeanAbundance() )
650 densAtom =
dense.gas_phase[nelem];
654 ( fabs(
molElems[i]- densAtom) <=
conv.GasPhaseAbundErrorAllowed*densAtom );
656 fprintf(
ioQQQ,
"PROBLEM molElems BAD %li\t%s\t%.12e\t%.12e\t%.12e\n",
671 mole.species[mol].den = 0.;
673 for (
long mol=0;mol<
mole_global.num_compacted;mol++)
684 long parentIndex =
mole_global.list[mol]->parentIndex;
685 mole.species[mol].den =
mole.species[parentIndex].den;
688 if( !it->first->lgMeanAbundance() )
690 mole.species[mol].den *= pow( it->first->frac, it->second );
697 for (
unsigned long j=0; j<
atom_list.size(); ++j )
703 for (
unsigned long ion=0;ion<
atom_list[j]->ipMl.size();ion++)
711 ASSERT(fabs(sum-grouptot) <= 1e-10 * fabs(grouptot));
715 mole.set_isotope_abundances();
754 long nelem = atom->el->Z-1;
758 for (
long ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
valarray< double > molElems
multi_arr< double, 2 > fion
void updateMolecules(const valarray< double > &b2)
void setup(double *b0vec)
UNUSED const realnum BIGFLOAT
bool lgElemsConserved(void)
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_mole_global mole_global
molezone * findspecieslocal(const char buf[])
ChemAtomList unresolved_atom_list
molecule::nAtomsMap::iterator nAtoms_i
void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr< double, 2 > &c)
realnum mole_return_cached_species(const GroupMap &MoleMap)
void check_co_ion_converge(void)
STATIC void funjac(GroupMap &MoleMap, const valarray< double > &b2vec, double *const ervals, double *const amat, const bool lgJac, bool *lgConserve)
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))