00001
00002
00003 #include "cddefines.h"
00004
00005 #include "dense.h"
00006
00007 #include "abund.h"
00008 #include "colden.h"
00009 #include "conv.h"
00010 #include "dynamics.h"
00011 #include "elementnames.h"
00012 #include "deuterium.h"
00013 #include "hmi.h"
00014 #include "h2.h"
00015 #include "mole.h"
00016 #include "mole_priv.h"
00017 #include "phycon.h"
00018 #include "prt.h"
00019 #include "radius.h"
00020 #include "struc.h"
00021 #include "thermal.h"
00022 #include "trace.h"
00023
00024 t_dense dense;
00025
00026 void t_dense::updateXMolecules()
00027 {
00028 total_molecule_elems(m_xMolecules);
00029 }
00030
00031 void t_dense::zero()
00032 {
00033 for (long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem)
00034 m_xMolecules[nelem] = 0.f;
00035 }
00036
00037 void ScaleAllDensities(realnum factor)
00038 {
00039 double edensave = dense.eden;
00040
00041 for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
00042 {
00043 if (dense.lgElmtOn[nelem])
00044 {
00045 ScaleIonDensities( nelem, factor );
00046 dense.SetGasPhaseDensity( nelem, dense.gas_phase[nelem] * factor);
00047 }
00048 }
00049
00050 EdenChange( dense.eden * factor );
00051
00052 if( trace.lgTrace && trace.lgNeBug )
00053 {
00054 fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n",
00055 edensave, dense.eden, edensave/dense.eden );
00056 }
00057
00058 hmi.H2_total *= factor;
00059 h2.ortho_density *= factor;
00060 h2.para_density *= factor;
00061 for( long mol=0; mol < mole_global.num_total; mol++ )
00062 {
00063 mole.species[mol].den *= factor;
00064 }
00065 dense.updateXMolecules();
00066
00067 ASSERT(lgElemsConserved());
00068 }
00069
00070 void ScaleIonDensities( const long nelem, const realnum factor )
00071 {
00072 double renorm;
00073 for( long ion=0; ion < (nelem + 2); ion++ )
00074 {
00075 dense.xIonDense[nelem][ion] *= factor;
00076 if (nelem-ion >= 0 && nelem-ion < NISO)
00077 iso_renorm(nelem, nelem-ion, renorm);
00078 }
00079
00080 if( nelem==ipHYDROGEN && deut.lgElmtOn )
00081 ScaleDensitiesDeuterium( factor );
00082
00083 return;
00084 }
00085
00086 void t_dense::SetGasPhaseDensity( const long nelem, const realnum density )
00087 {
00088 gas_phase[nelem] = density;
00089
00090 if( nelem==ipHYDROGEN && deut.lgElmtOn )
00091 {
00092
00093 SetGasPhaseDeuterium( density );
00094 }
00095
00096 return;
00097 }
00098
00099 bool lgElemsConserved (void)
00100 {
00101 bool lgOK=true;
00102
00103
00104 for( ChemAtomList::iterator atom = atom_list.begin(); atom != atom_list.end(); ++atom )
00105 {
00106 long nelem = (*atom)->el->Z-1;
00107
00108 if( dense.lgElmtOn[nelem] )
00109 {
00110
00111
00112 double sum_monatomic = 0.;
00113 for( long ion=0; ion<nelem+2; ++ion )
00114 {
00115 sum_monatomic += dense.xIonDense[nelem][ion] * (*atom)->frac;
00116 }
00117 realnum sum_molecules = dense.xMolecules(nelem) * (*atom)->frac;
00118 realnum sum_gas_phase = dense.gas_phase[nelem] * (*atom)->frac;
00119 if( sum_monatomic + sum_molecules <= SMALLFLOAT &&
00120 sum_gas_phase > SMALLFLOAT)
00121 {
00122 fprintf(ioQQQ,"PROBLEM non-conservation of nuclei %s\tions %g moles %g error %g of %g\n",
00123 (*atom)->label().c_str(),
00124 sum_monatomic,
00125 sum_molecules,
00126 sum_monatomic + sum_molecules-sum_gas_phase,
00127 sum_gas_phase );
00128 lgOK=false;
00129 }
00130
00131
00135 if( fabs( sum_monatomic + sum_molecules - sum_gas_phase ) >
00136 conv.GasPhaseAbundErrorAllowed * sum_gas_phase )
00137 {
00138 fprintf(ioQQQ,"PROBLEM non-conservation of nuclei %s\t nzone %li atoms %.12e moles %.12e sum %.12e tot gas %.12e rel err %.3e\n",
00139 (*atom)->label().c_str(),
00140 nzone,
00141 sum_monatomic,
00142 sum_molecules,
00143 sum_monatomic + sum_molecules,
00144 sum_gas_phase,
00145 (sum_monatomic + sum_molecules - sum_gas_phase)/sum_gas_phase );
00146 lgOK=false;
00147 }
00148
00149 if( 0 )
00150 {
00151
00152 mole_print_species_reactions( findspecies( elementnames.chElementSym[nelem] ) );
00153 }
00154 }
00155 }
00156
00157 return lgOK;
00158 }
00159
00160 void lgStatesConserved ( long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion )
00161 {
00162 if( 0 && conv.lgConvIoniz() &&
00163 dense.lgElmtOn[nelem] &&
00164 dense.xIonDense[nelem][ionStage]/dense.eden > conv.EdenErrorAllowed/10. &&
00165 dense.xIonDense[nelem][ionStage] > SMALLFLOAT )
00166 {
00167 double abund = 0.;
00168 for (long ipLev=0; ipLev<numStates; ipLev++)
00169 {
00170 abund += states[ipLev].Pop();
00171 }
00172
00173 double rel_err = fabs(abund-dense.xIonDense[nelem][ionStage])/
00174 SDIV(dense.xIonDense[nelem][ionStage]);
00175
00176
00177
00178 if( rel_err > err_tol )
00179 {
00180
00181 if( 0 && nzone > 0 && loop_ion > 0 )
00182 {
00183 fprintf(ioQQQ,"PROBLEM Inconsistent states/stage pops nzone %3ld loop_ion %2ld nelem %2ld ion %2ld states = %e stage = %e error = %e\n",
00184 nzone,
00185 loop_ion,
00186 nelem,
00187 ionStage,
00188 abund,
00189 dense.xIonDense[nelem][ionStage],
00190 rel_err );
00191 }
00192 char chConvIoniz[INPUT_LINE_LENGTH];
00193 sprintf( chConvIoniz , "States!=stage pops nelem %ld ion %ld ", nelem, ionStage );
00194 conv.setConvIonizFail( chConvIoniz, abund,
00195 dense.xIonDense[nelem][ionStage] );
00196 }
00197 }
00198 }
00199
00200 void SumDensities( void )
00201 {
00202 realnum den_mole,
00203 DenseAtomsIons;
00204
00205
00206 DenseAtomsIons = 0.;
00207 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00208 {
00209 if( dense.lgElmtOn[nelem] )
00210 {
00211 for( long ion=0; ion<=nelem+1; ++ion )
00212 DenseAtomsIons += dense.xIonDense[nelem][ion];
00213 }
00214 }
00215
00216
00217
00218
00219 ASSERT( DenseAtomsIons > 0. );
00220
00221
00222
00223
00224 den_mole = total_molecules_gasphase();
00225
00226
00227 dense.xNucleiTotal = den_mole + DenseAtomsIons;
00228 if( dense.xNucleiTotal > BIGFLOAT )
00229 {
00230 fprintf(ioQQQ, "PROBLEM DISASTER SumDensities has found "
00231 "dense.xNucleiTotal with an insane density.\n");
00232 fprintf(ioQQQ, "The density was %.2e\n",
00233 dense.xNucleiTotal);
00234 TotalInsanity();
00235 }
00236 ASSERT( dense.xNucleiTotal > 0. );
00237
00238
00239
00240 dense.pden = (realnum)(dense.eden + dense.xNucleiTotal);
00241
00242
00243 dense.wmole = 0.;
00244 for( long i=0; i < LIMELM; i++ )
00245 {
00246 dense.wmole += dense.gas_phase[i]*(realnum)dense.AtomicWeight[i];
00247 }
00248
00249 dense.wmole /= dense.pden;
00250 ASSERT( dense.wmole > 0. && dense.pden > 0. );
00251
00252
00255 dense.xMassDensity = (realnum)(dense.wmole*ATOMIC_MASS_UNIT*dense.pden);
00256
00257
00258
00259
00260
00261
00262
00263 if( dense.xMassDensity0 < 0.0 )
00264 dense.xMassDensity0 = dense.xMassDensity;
00265
00266 return;
00267 }
00268
00269 bool AbundChange( )
00270 {
00271 DEBUG_ENTRY( "AbundChange()" );
00272
00273 fixit();
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284 bool lgChange = false;
00285
00286
00287
00288 if( abund.lgAbTaON )
00289 {
00290 lgChange = true;
00291 for( long nelem=1; nelem < LIMELM; nelem++ )
00292 {
00293 if( abund.lgAbunTabl[nelem] )
00294 {
00295 double abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
00296 dense.gas_phase[ipHYDROGEN];
00297
00298 double hold = abun/dense.gas_phase[nelem];
00299 dense.SetGasPhaseDensity( nelem, (realnum)abun );
00300
00301 for( long ion=0; ion < (nelem + 2); ion++ )
00302 {
00303 dense.xIonDense[nelem][ion] *= (realnum)hold;
00304 }
00305 }
00306 }
00307 }
00308
00309
00310
00311
00312 double FacAbun = 1.;
00313 if( !dense.lgDenFlucOn )
00314 {
00315 static double FacAbunSav;
00316 double OldAbun = 0.0;
00317
00318 lgChange = true;
00319 if( nzone > 1 )
00320 {
00321 OldAbun = FacAbunSav;
00322 }
00323
00324
00325 if( dense.lgDenFlucRadius )
00326 {
00327
00328 FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
00329 dense.flcPhase) + dense.csecnd;
00330 }
00331 else
00332 {
00333
00334 FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
00335 dense.flcPhase) + dense.csecnd;
00336 }
00337
00338 if( nzone > 1 )
00339 {
00340 FacAbun = FacAbunSav/OldAbun;
00341 }
00342 }
00343
00344 if( FacAbun != 1. )
00345 {
00346 ASSERT(!dynamics.lgAdvection);
00347
00348
00349 for( long nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
00350 {
00351 if (dense.lgElmtOn[nelem])
00352 {
00353 dense.SetGasPhaseDensity(nelem, dense.gas_phase[nelem] * (realnum) FacAbun);
00354 ScaleIonDensities( nelem, (realnum)(FacAbun) );
00355 }
00356 }
00357
00358 for( long mol=0; mol < mole_global.num_total; mol++ )
00359 {
00360 mole.species[mol].den *= (realnum)(FacAbun);
00361 }
00362 }
00363
00364 if (lgChange)
00365 {
00366
00367
00368 TempChange(phycon.te , false);
00369 }
00370
00371 return lgChange;
00372 }
00373
00374 namespace {
00375 const int SCALE_NEW = 1;
00376 }
00377
00378 realnum scalingDensity(void)
00379 {
00380 if (SCALE_NEW)
00381 return dense.xMassDensity/(realnum)ATOMIC_MASS_UNIT;
00382 else
00383 return dense.gas_phase[ipHYDROGEN];
00384 }
00385 realnum scalingZoneDensity(long i)
00386 {
00387 if (SCALE_NEW)
00388 return struc.DenMass[i]/(realnum)ATOMIC_MASS_UNIT;
00389 else
00390 return struc.hden[i];
00391 }