cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*CO_step fills in matrix for heavy elements molecular routines */ 00004 #include "cddefines.h" 00005 #include "taulines.h" 00006 #include "dense.h" 00007 #include "ionbal.h" 00008 #include "thermal.h" 00009 #include "phycon.h" 00010 #include "hmi.h" 00011 #include "dynamics.h" 00012 #include "conv.h" 00013 #include "trace.h" 00014 #include "coolheavy.h" 00015 #include "timesc.h" 00016 #include "thirdparty.h" 00017 #include "mole.h" 00018 #include "mole_co_priv.h" 00019 #include "mole_co_atom.h" 00020 /* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element 00021 * molecular network in Cloudy. Before this routine would predict negative abundances if 00022 * the fraction of carbon in the form of molecules came close to 100%. A reorganizing of 00023 * the reaction network detected several bugs. Treatment of "coupled reactions", 00024 * in which both densities in the reaction rate were being predicted by Cloudy, were also 00025 * added. Due to these improvements, Cloudy can now perform calculations 00026 * where 100% of the carbon is in the form of CO without predicting negative abundances 00027 * 00028 * Additional changes were made in November of 2003 so that our reaction 00029 * network would include all reactions from the TH85 paper. This involved 00030 * adding silicon to the chemical network. Also the reaction rates were 00031 * labeled to make identification with the reaction easier and the matrix 00032 * elements of atomic C, O, and Si are now done in a loop, which makes 00033 * the addition of future chemical species (like N or S) easy. 00034 * */ 00035 00036 void CO_solve( 00037 /* set true if we found neg pops */ 00038 bool *lgNegPop, 00039 /* set true if we tried to compute the pops, but some were zero */ 00040 bool *lgZerPop ) 00041 { 00042 00043 int32 merror; 00044 long int i, j, k, n, 00045 nelem , ion , ion2; 00046 double 00047 co_denominator; 00048 realnum cartot_mol, oxytot_mol; 00049 realnum abundan; 00050 struct chem_element_s *element; 00051 double sum; 00052 00053 DEBUG_ENTRY( "CO_solve()" ); 00054 00055 CO_step(); /* Calculate the matrix elements */ 00056 00057 /* Ugly hack to deal with species which have become uncoupled */ 00058 for(i=0;i<mole.num_comole_calc;i++) 00059 { 00060 sum = 0.; 00061 for(j=0;j<mole.num_comole_calc;j++) 00062 { 00063 sum = sum+fabs(mole.c[i][j]); 00064 } 00065 if(sum < SMALLFLOAT && fabs(mole.b[i]) < SMALLFLOAT) 00066 { 00067 mole.c[i][i] = 1.; 00068 } 00069 } 00070 00071 cartot_mol = dense.xMolecules[ipCARBON] + findspecies("C")->hevmol + findspecies("C+")->hevmol; 00072 oxytot_mol = dense.xMolecules[ipOXYGEN] + findspecies("O")->hevmol + findspecies("O+")->hevmol; 00073 00074 ASSERT( cartot_mol >= 0. && oxytot_mol >= 0.); 00075 00076 *lgNegPop = false; 00077 *lgZerPop = false; 00078 00079 /* zero out molecular charge transfer rates */ 00080 for(nelem=ipLITHIUM; nelem < LIMELM; ++nelem) 00081 { 00082 for(ion=0; ion < nelem+2; ++ion) 00083 { 00084 /*zero out the arrays */ 00085 mole.sink[nelem][ion] = 0.; 00086 mole.source[nelem][ion] = 0.; 00087 00088 for(ion2=0; ion2< nelem+2; ++ion2) 00089 { 00090 mole.xMoleChTrRate[nelem][ion][ion2] = 0.; 00091 } 00092 } 00093 } 00094 00095 00096 /* >>06 jun 29, add advective terms here */ 00097 /* Add rate terms for dynamics to equilibrium, makes c[][] non-singular 00098 * do before cross talk with heavy elements is done to not double count charge 00099 * transfer */ 00100 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection /*&& radius.depth < dynamics.oldFullDepth*/ ) 00101 { 00102 for(i=0;i<mole.num_comole_calc;i++) 00103 { 00104 if(COmole[i]->n_nuclei > 1) { 00105 mole.c[i][i] -= dynamics.Rate; 00106 /* this is the net gain of the species due to advection - 00107 * in this form we have the net gain as the difference between 00108 * current species flowing out of this region, downstream, and 00109 * new advected material coming into this region from upstream */ 00110 mole.b[i] -= (COmole[i]->hevmol*dynamics.Rate-dynamics.CO_molec[i]); 00111 } 00112 } 00113 } 00114 00115 /* >>chng Oct. 21, 2004 -- Terms that contribute to the ionization balance of C, O, S, Si, and N 00116 will now be inserted directly into the ion solver. Before, we just corrected the recombination rate 00117 by scaling saying that IONIZATION_RATE = RECOMBINATION_RATE * [n(X+)/n(X0)]. This code follows the logic 00118 of hmole_step, written by Robin Williams. */ 00119 00120 00121 /* sink terms should include only terms that either form or destroy an atomic or singly ionized 00122 element in the network, but not both. This means that terms that destroy X at the expense of 00123 forming X+ can't be included. The rate of destruction of X or X+ is equal to the formation of 00124 the inverse. We therefore add this rate to sink, which gets rid of this dependence */ 00125 00126 /* The following code is all generalized. After all the molecules, the total number being mole.num_heavy_molec, 00127 there are 2*mole.num_elements remaining. The first set, equal to mole.num_elements, are the ionized elemental 00128 species. The second set, also equal to mole.num_elements, are the atomic species. They are in order such that 00129 00130 array number for X = array number for (X+) + mole.num_elements 00131 00132 In the future, if one wants to add another element to the network, such as Na the macros must be changed. 00133 Also, the array number for Na and Na+ must obey the equation above. If this is done, then the sources, 00134 sinks, and molecular recombination terms are automatically added!!! */ 00135 00136 /* COmole[i]->nelem_hevmol[j] is just the dominant element for a given molecule. For an element this 00137 string is just the element itself!! */ 00138 00139 /* >> chng 2006 Sep 02 rjrw: Change to using element properties rather than properties of elements as molecules, for ease of 00140 generalization -- ordering of atoms and ions no longer hardwired */ 00141 00142 for(i=0;i<mole.num_elements; ++i) 00143 { 00144 element = chem_element[i]; 00145 mole.sink[element->ipCl][0] -= (mole.c[element->ipMl][element->ipMl] + mole.c[element->ipMl][element->ipMlP])*dense.lgElmtOn[element->ipCl]; 00146 mole.sink[element->ipCl][1] -= (mole.c[element->ipMlP][element->ipMlP] + mole.c[element->ipMlP][element->ipMl] )*dense.lgElmtOn[element->ipCl]; 00147 } 00148 00149 /* source terms must be multiplied by the density of the corresponding matrix element */ 00150 00151 for(j=0;j < mole.num_elements; ++j) 00152 { 00153 element = chem_element[j]; 00154 for(i=0; i < mole.num_comole_calc; ++i) 00155 { 00156 mole.source[element->ipCl][1] += COmole[i]->hevmol*mole.c[i][element->ipMlP]*dense.lgElmtOn[element->ipCl]; 00157 mole.source[element->ipCl][0] += COmole[i]->hevmol*mole.c[i][element->ipMl]*dense.lgElmtOn[element->ipCl]; 00158 } 00159 } 00160 00161 /* subtract out diagonal terms, they are already in sink. Also subtract recombination terms, 00162 as these are done by mole.xMoleChTrRate */ 00163 00164 for(i=0;i < mole.num_elements; ++i) 00165 { 00166 element = chem_element[i]; 00167 mole.source[element->ipCl][1] -= ( mole.c[element->ipMlP][element->ipMlP]*COmole[element->ipMlP]->hevmol + 00168 mole.c[element->ipMl][element->ipMlP]*COmole[element->ipMl]->hevmol)*dense.lgElmtOn[element->ipCl]; 00169 00170 mole.source[element->ipCl][0] -= ( mole.c[element->ipMl][element->ipMl]*COmole[element->ipMl]->hevmol + 00171 mole.c[element->ipMlP][element->ipMl]*COmole[element->ipMlP]->hevmol)*dense.lgElmtOn[element->ipCl]; 00172 00173 /* The source terms, as they are right now, include negative contributions. This is because the 00174 linearization "trick" creates source terms. Take for example the reaction: 00175 C+ H2O => HCO+ H 00176 This reaction destroys C+, but in the act of linearizing, there will be the following term: 00177 mole.c[ipH2O][ipCP] 00178 This is only a numerical "trick" and not a real matrix term. All terms like this 00179 have to be removed to get the "true" source terms. Fortunately, the sum of all these terms 00180 is just mole.b[i]. To remove these terms, we have to add mole.b[i] if the overall contribution from these terms 00181 was negative, subtract if their contribution was positive. This is done by subtracting mole.b[i] 00182 from the current value of source. */ 00183 00184 mole.source[element->ipCl][1] = mole.source[element->ipCl][1] - mole.b[element->ipMlP]; 00185 mole.source[element->ipCl][0] = mole.source[element->ipCl][0] - mole.b[element->ipMl]; 00186 00187 } 00188 00189 /* negative source terms are actually destruction mechanisms, so should go into the sink vector. 00190 source and sinks have different units, so if source is negative divide by the density of 00191 the corresponding element to get same units. For example, if: 00192 00193 mole.source[ipCARBON][0] 00194 00195 is negative, we divide by dense.xIonDense[ipCARBON][0] to get rate in same units as mole.sink */ 00196 00197 for(i=2; i < LIMELM; ++i) 00198 { 00199 for(j=0; j < 2; ++j) 00200 { 00201 /*if source term is negative, make it a sink and then set to zero*/ 00202 if(mole.source[i][j] < 0) 00203 { 00204 mole.sink[i][j] -= (mole.source[i][j]/SDIV(dense.xIonDense[i][j])); 00205 mole.source[i][j] = 0; 00206 } 00207 } 00208 } 00209 /* These are rates of recombination for atomic and singly ionized species that are due to molecular processes 00210 This will be added to ion_solver to get correct ionization balance */ 00211 00212 for(i=0;i < mole.num_elements; ++i) 00213 { 00214 element = chem_element[i]; 00215 mole.xMoleChTrRate[element->ipCl][1][0] = (realnum)(mole.c[element->ipMlP][element->ipMl] - 00216 ionbal.RateRecomTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl]; 00217 00218 mole.xMoleChTrRate[element->ipCl][0][1] = (realnum)(mole.c[element->ipMl][element->ipMlP] - 00219 ionbal.RateIonizTot[element->ipCl][0])*dense.lgElmtOn[element->ipCl]; 00220 } 00221 00222 /* If rate for going from 1-0 or 0-1 is negative then it should be added to the inverse process */ 00223 for(i=2; i < LIMELM; ++i) 00224 { 00225 for(j=0; j < 2; ++j) 00226 { 00227 for(k=0; k< 2; ++k) 00228 { 00229 /*only possible charge transfers are 0-1 or 1-0 */ 00230 if(j != k) 00231 { 00232 if( mole.xMoleChTrRate[i][j][k] < 0. ) 00233 { 00234 mole.xMoleChTrRate[i][k][j] -= mole.xMoleChTrRate[i][j][k]; 00235 mole.xMoleChTrRate[i][j][k] = 0.; 00236 } 00237 } 00238 } 00239 } 00240 } 00241 00242 for(n=0;n<mole.num_elements;n++) 00243 { 00244 tot_ion[n] = dense.gas_phase[chem_element[n]->ipCl]; 00245 for( i=2; i < chem_element[n]->ipCl+2; i++ ) 00246 { 00247 tot_ion[n] -= dense.xIonDense[chem_element[n]->ipCl][i]; 00248 } 00249 } 00250 00251 00252 /* at this point the cartot_mol, sum of molecules and atom/first ion, 00253 * should be equal to the gas_phase minus double and higher ions */ 00254 /*fprintf(ioQQQ," dbuggggas\t %.2f\t%f\t%f\n",fnzone, 00255 cartot_mol/cartot_ion, 00256 oxytot_mol/oxytot_ion);*/ 00257 00258 /* these will be used in population equation in case of homogeneous equation */ 00259 00260 00261 /* rjrw 2006 Aug 08: messing with the matrix like this will likely break advection -- 00262 was done correctly in hmole */ 00263 00264 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection) 00265 fixit(); 00266 00267 for(n=0;n<mole.num_elements;n++) 00268 { 00269 mole.b[chem_element[n]->ipMl] = tot_ion[n]; 00270 } 00271 00272 /* <<chng 03 Nov 11--Nick Abel, Set up the non-zero matrix elements that go into the conservation equations 00273 for atomic C, O, and Si, this is now set up by looping over the atomic species in 00274 co.h and setting the number of atoms of C, O, and Si equal to the variable findspecies("CARBON")->nElem, 00275 findspecies("OXYGEN")->nElem, and findspecies("SILICON")->nElem respectively. For every element (excluding Hydrogen) not in 00276 the network an if statement will be needed */ 00277 00278 /* loop over last mole.num_elements elements in co vector - these are atoms of the mole.num_elements */ 00279 for( j=0; j < mole.num_comole_calc; j++ ) 00280 { 00281 for(i=0;i<mole.num_elements;i++) 00282 { 00283 element = chem_element[i]; 00284 mole.c[j][element->ipMl] = COmole[j]->nElem[element->ipCl]; 00285 } 00286 } 00287 00288 /*-------------------------------------------------------------------- 00289 * */ 00290 /*printf( " Here are all matrix elements\n "); 00291 00292 for(i=0; i < mole.num_comole_calc; i++) 00293 00294 { 00295 00296 00297 for(j=0; j < mole.num_comole_calc; j++) 00298 { 00299 printf( "%4.4s", COmole[i]->label ); 00300 printf( "%4.4s\t", COmole[j]->label ); 00301 printf( "%.8e,\n", mole.c[j][i] ); 00302 } 00303 } 00304 printf( "b's are:\n" ); 00305 for(i=0; i < mole.num_comole_calc; i++) 00306 { 00307 printf( "%.8e,\n", mole.b[i] ); 00308 }*/ 00309 00310 00311 if( trace.lgTrace && trace.lgTr_CO_Mole ) 00312 { 00313 fprintf( ioQQQ, " COMOLE matrix\n " ); 00314 00315 # if 0 00316 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ ) 00317 { 00318 fprintf( ioQQQ, "%4.4s", COmole[i]->label ); 00319 } 00320 fprintf( ioQQQ, " \n" ); 00321 00322 for( j=0; j < mole.num_comole_calc; j++ ) 00323 { 00324 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00325 fprintf( ioQQQ, " " ); 00326 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ ) 00327 { 00328 fprintf( ioQQQ, "%9.1e", mole.c[i][j] ); 00329 } 00330 fprintf( ioQQQ, "\n" ); 00331 } 00332 # endif 00333 00334 fprintf( ioQQQ, " COMOLE matrix\n " ); 00335 00336 for( i=0; i < mole.num_comole_calc; i++ ) 00337 { 00338 fprintf( ioQQQ, "%4.4s", COmole[i]->label ); 00339 } 00340 fprintf( ioQQQ, " \n" ); 00341 00342 for( j=0; j < mole.num_comole_calc; j++ ) 00343 { 00344 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00345 fprintf( ioQQQ, " " ); 00346 for( i=0; i < mole.num_comole_calc; i++ ) 00347 { 00348 /*fprintf( ioQQQ, "%9.1e", mole.c[i][j] );*/ 00349 fprintf( ioQQQ, ", %.15e", mole.c[i][j] ); 00350 } 00351 fprintf( ioQQQ, "\n" ); 00352 } 00353 00354 fprintf( ioQQQ, " COMOLE b\n " ); 00355 for( j=0; j < mole.num_comole_calc; j++ ) 00356 { 00357 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00358 fprintf( ioQQQ, ", %.15e\n",mole.b[j] ); 00359 } 00360 00361 #if 0 00362 fprintf( ioQQQ, 00363 " COMOLE H2 den:%10.3e, H2,3+=%10.2e%10.2e Carb sum=%10.3e Oxy sum=%10.3e\n", 00364 hmi.H2_total, 00365 hmi.Hmolec[ipMH2p], 00366 hmi.Hmolec[ipMH3p], 00367 mole.c[mole.num_heavy_molec][findspecies("C")->index], 00368 mole.c[mole.num_heavy_molec][findspecies("O")->index] ); 00369 00370 #endif 00371 } 00372 00373 /* remember longest molecule destruction timescales */ 00374 for( j=0; j < mole.num_comole_calc; j++ ) 00375 { 00376 if( -mole.c[j][j] > SMALLFLOAT ) 00377 { 00378 /* this is rate CO is destroyed, equal to formation rate in equilibrium */ 00379 timesc.AgeCOMoleDest[j] = -1./mole.c[j][j]; 00380 /* moved to radinc */ 00381 /*timesc.BigCOMoleForm = (realnum)MAX2(timesc.BigCOMoleForm,timesc.AgeCOMoleForm);*/ 00382 } 00383 else 00384 { 00385 timesc.AgeCOMoleDest[j] = 0.; 00386 } 00387 } 00388 00389 /* copy to amat, saving c for later use */ 00390 for( j=0; j < mole.num_comole_calc; j++ ) 00391 { 00392 for( i=0; i < mole.num_comole_calc; i++ ) 00393 { 00394 mole.amat[i][j] = mole.c[i][j]; 00395 } 00396 } 00397 00398 merror = 0; 00399 00400 getrf_wrapper(mole.num_comole_calc,mole.num_comole_calc,mole.amat[0],mole.num_comole_calc, ipiv,&merror); 00401 00402 if( merror != 0 ) 00403 { 00404 fprintf( ioQQQ, " CO_solve getrf_wrapper error\n" ); 00405 cdEXIT(EXIT_FAILURE); 00406 } 00407 00408 getrs_wrapper('N',mole.num_comole_calc,1,mole.amat[0],mole.num_comole_calc,ipiv,mole.b,mole.num_comole_calc,&merror); 00409 00410 if( merror != 0 ) 00411 { 00412 fprintf( ioQQQ, " CO_solve: dgetrs finds singular or ill-conditioned matrix\n" ); 00413 cdEXIT(EXIT_FAILURE); 00414 } 00415 00416 /* check for negative populations, which happens when 100% co */ 00417 *lgNegPop = false; 00418 *lgZerPop = false; 00419 co_denominator = 0; 00420 # if 0 00421 if(fnzone > 1.02) 00422 { 00423 for( i=0; i < mole.num_comole_calc; i++) 00424 { 00425 00426 for( j=61; j < 62; j++ ) 00427 { 00428 00429 printf("ZONE\t%.5e\tSPECIES_1\t%s\tSPECIES_2\t%s\tDEST_RATE\t%.3e\tbvec\t%.3e\n", 00430 fnzone, COmole[i]->label,COmole[j]->label,mole.c[i][j],mole.b[i]); 00431 00432 } 00433 } 00434 } 00435 # endif 00436 for( i=0; i < mole.num_comole_calc; i++ ) 00437 { 00438 /* fprintf(stderr,"%ld [%d] %g\n",i,mole.num_comole_calc,mole.b[i]); */ 00439 if( mole.b[i] < 0. ) 00440 { 00441 00442 /* >>chng 03 sep 11, the following*/ 00443 /* comparison between the solution vector from 32 bit machines vs 00444 * maple on a pc shows that very small negative numbers are produced by 00445 * the linear algebra package used by the code, but is positive with 00446 * the math package. allow very small negative numbers to occur, 00447 * and silently reset to a postive number. 00448 * 00449 * >>chng 04 feb 25 00450 * Here we want to check if the negative abundance is a significant 00451 * fraction of any of the species in the network. The code below 00452 * checks to see which elements the molecule in question is made of, 00453 * then finds which one has the least abundance. The one with 00454 * the least abundance will then be used as the divisor in checking 00455 * if the negative abundance is too large to be ignored.*/ 00456 00457 /* >>chng 04 apr 21, when an element in the chemical network is 00458 * turned off, the molecular abundance of a species containing that 00459 * element was not zero but rather a small number of order 1e-20. When 00460 * these abundances went negative, the code thought the abundances were 00461 * important because it checks the ratio of the species to its gas phase 00462 * abundance. When the gas phase abundance is zero, the denominator is set 00463 * to SMALLFLOAT, which is ~1e-35. This led to a ratio of molecule to gas 00464 * phase of ~1e15, clearly unphysical. Here the value of co_denominator 00465 * will be set to a high value if the gas phase abundance is zero.*/ 00466 00467 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] ) 00468 { 00469 co_denominator = dense.gas_phase[COmole[i]->nelem_hevmol]; 00470 } 00471 else 00472 { 00473 /* >>chng 04 apr 20, set to zero if element is off */ 00474 co_denominator = 1e10; 00475 } 00476 00477 00478 /* we must have a positive value by now, unlesss element is turned off */ 00479 ASSERT( co_denominator > SMALLFLOAT || !dense.lgElmtOn[COmole[i]->nelem_hevmol] ); 00480 00481 /* >>chng 04 feb 28, threshold from -1e-10 to -1e-6, the 00482 * level of roundoff in a realnum */ 00483 /*if( mole.b[i] / MAX2 (co_denominator, SMALLFLOAT) > -1e-10) */ 00484 /*if( mole.b[i] / SDIV(co_denominator) > -1e-6 ) */ 00485 /* >>chng 04 oct 31, change limit from -1e-6 to -1e-5, produced only 00486 * a few comments in func_map.in */ 00487 if( mole.b[i] / SDIV(co_denominator) > -1e-5 ) 00488 { 00489 /* this case, it is only slightly negative relative to 00490 * the parent species - multiply by -1 to make positive 00491 * and press on */ 00492 mole.b[i] *= -1.; 00493 } 00494 else 00495 { 00496 /*>>chng 04 mar 06 press on */ 00497 *lgNegPop = true; 00498 /* 05 jul 16, there had been a bug such that CO was always evaluated as long as 00499 * the temperature was below 20,000K. Once fixed, one sim turned on CO mid way 00500 * into the H+ zone and had neg pops during first trys - change check on whether it 00501 * is to be commented on only if this is not first zone with CO soln */ 00502 /*if( nzone>0 )*/ 00503 if( nzone>co.co_nzone ) 00504 { 00505 static long int nzFail=-2; 00506 long int nzFail2 = (long)fnzone*100; 00507 /* during a map we may be in search phase but not in first zone */ 00508 if( !conv.lgSearch ) 00509 fprintf(ioQQQ," PROBLEM "); 00510 fprintf(ioQQQ, 00511 " CO_solve neg pop for species %li %s, value is %.2e rel value is %.2e zone %.2f Te %.4e Search?%c\n", 00512 i , 00513 COmole[i]->label , 00514 /* abs value */ 00515 mole.b[i], 00516 /* relative value */ 00517 mole.b[i] / SDIV(co_denominator) , 00518 fnzone, 00519 phycon.te,TorF( conv.lgSearch ) ); 00520 /* if well beyond search phase and still having problems, announce them 00521 * only announce one failure per sweep through solvers by ConvFail */ 00522 if( nzFail2 !=nzFail ) 00523 { 00524 nzFail = nzFail2; 00525 ConvFail( "pops" , "CO"); 00526 } 00527 } 00528 mole.b[i] *= -1; 00529 /*>>chng 04 jan 24 set to zero instead of *-1, 00530 * h2_cr.in co density went to infinite, with some negative 00531 * var getting to - inf 00532 *>>chng 04 nov 03, remove this - not needed h2_cr works without 00533 mole.b[i] = SMALLFLOAT;*/ 00534 } 00535 } 00536 else if( mole.b[i] == 0. ) 00537 { 00538 /* this is not used for anything in calling routine and 00539 * could be cleaned up */ 00540 *lgZerPop = true; 00541 /* >>chng 04 feb 28, zero pop not really a problem in very deep neutral 00542 * gas - this happens */ 00543 # if 0 00544 if( /*nzone>0 ||*/ CODEBUG>1 ) 00545 { 00546 fprintf(ioQQQ," PROBLEM "); 00547 fprintf(ioQQQ, 00548 " CO_solve zero pop for species %li %s, value is %.2e zone %li\n", 00549 i , COmole[i]->label , mole.b[i], nzone); 00550 } 00551 # endif 00552 mole.b[i] = SMALLFLOAT; 00553 } 00554 COmole[i]->hevmol = (realnum)mole.b[i]; 00555 /* check for NaN */ 00556 ASSERT( !isnan( COmole[i]->hevmol ) ); 00557 } 00558 /* >>chng 04 mar 06 pass negative pop as a problem, but not a fatal one, 00559 * also do not call this during 0th zone when search for conditions 00560 * is underway */ 00561 /* 05 jul 16, there had been a bug such that CO was always evaluated as long as 00562 * the temperature was below 20,000K. Once fixed, on sim turned on CO mid way 00563 * into the H+ zone and had neg pops during first trys - change check on whether it 00564 * is to be commented on only if this is not first zone with CO soln */ 00565 /*if( *lgNegPop && (nzone>0 &&!conv.lgSearch) )*/ 00566 if( *lgNegPop && (nzone>co.co_nzone &&!conv.lgSearch) ) 00567 { 00568 conv.lgConvPops = false; 00569 fprintf(ioQQQ," CO network negative population occurred, Te=%.4e, calling ConvFail. ", 00570 phycon.te); 00571 fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) ); 00572 fprintf( ioQQQ, "\n" ); 00573 /*ConvFail("pops", "CO");*/ 00574 *lgNegPop = false; 00575 } 00576 /* if negative pops were present need to renorm b to get 00577 * proper sum rule */ 00578 if( 0 && *lgNegPop ) 00579 { 00580 00581 for(j=0; j<2; ++j ) 00582 { 00583 double sumcar = 0.; 00584 double sumoxy = 0.; 00585 double renorm; 00586 for( i=0; i < mole.num_comole_calc; i++ ) 00587 { 00588 /* this case, different molecules */ 00589 sumcar += COmole[i]->hevmol*COmole[i]->nElem[ipCARBON]; 00590 sumoxy += COmole[i]->hevmol*COmole[i]->nElem[ipOXYGEN]; 00591 } 00592 renorm = (tot_ion[0] + tot_ion[1]) / SDIV( sumcar+sumoxy); 00593 if(j) 00594 fprintf(ioQQQ,"\t%f\n", renorm); 00595 else 00596 fprintf(ioQQQ,"renormco\t%.2f\t%f", fnzone, renorm); 00597 for( i=0; i < mole.num_comole_calc; i++ ) 00598 { 00599 COmole[i]->hevmol *= (realnum)renorm; 00600 /* check for NaN */ 00601 ASSERT( !isnan( COmole[i]->hevmol ) ); 00602 } 00603 } 00604 } 00605 00606 if( merror != 0 ) 00607 { 00608 fprintf( ioQQQ, " COMOLE matrix inversion error, MERROR=%5ld zone=%5ld\n", 00609 (long)merror, nzone ); 00610 ShowMe(); 00611 fprintf( ioQQQ, " Product matrix\n " ); 00612 00613 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ ) 00614 { 00615 fprintf( ioQQQ, "%4.4s", COmole[i]->label ); 00616 } 00617 fprintf( ioQQQ, " \n" ); 00618 00619 for( j=0; j < mole.num_comole_calc; j++ ) 00620 { 00621 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00622 fprintf( ioQQQ, " " ); 00623 00624 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ ) 00625 { 00626 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* 00627 COmole[i]->hevmol ); 00628 } 00629 fprintf( ioQQQ, "\n" ); 00630 } 00631 00632 if( mole.num_comole_calc >= 9 ) 00633 { 00634 fprintf( ioQQQ, " COMOLE matrix\n " ); 00635 for( i=8; i < mole.num_comole_calc; i++ ) 00636 { 00637 fprintf( ioQQQ, "%4.4s", COmole[i]->label ); 00638 } 00639 fprintf( ioQQQ, " \n" ); 00640 00641 for( j=0; j < mole.num_comole_calc; j++ ) 00642 { 00643 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00644 fprintf( ioQQQ, " " ); 00645 for( i=8; i < mole.num_comole_calc; i++ ) 00646 { 00647 fprintf( ioQQQ, "%9.1e", 00648 mole.amat[i][j]* COmole[i]->hevmol ); 00649 } 00650 fprintf( ioQQQ, "\n" ); 00651 } 00652 } 00653 00654 fprintf( ioQQQ, " Mole dens:" ); 00655 for( j=0; j < mole.num_comole_calc; j++ ) 00656 { 00657 fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label 00658 , COmole[j]->hevmol ); 00659 } 00660 fprintf( ioQQQ, " \n" ); 00661 00662 cdEXIT(EXIT_FAILURE); 00663 } 00664 00665 if( trace.lgTrace && trace.lgTr_CO_Mole ) 00666 { 00667 fprintf( ioQQQ, " Product matrix\n " ); 00668 00669 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ ) 00670 { 00671 fprintf( ioQQQ, "%4.4s", COmole[i]->label ); 00672 } 00673 fprintf( ioQQQ, " \n" ); 00674 00675 for( j=0; j < mole.num_comole_calc; j++ ) 00676 { 00677 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00678 fprintf( ioQQQ, " " ); 00679 for( i=0; i < MIN2(mole.num_comole_calc,8); i++ ) 00680 { 00681 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* 00682 COmole[i]->hevmol ); 00683 } 00684 fprintf( ioQQQ, "\n" ); 00685 00686 } 00687 00688 if( mole.num_comole_calc >= 9 ) 00689 { 00690 fprintf( ioQQQ, " COMOLE matrix\n " ); 00691 for( i=8; i < mole.num_comole_calc; i++ ) 00692 { 00693 fprintf( ioQQQ, "%4.4s", COmole[i]->label ); 00694 } 00695 fprintf( ioQQQ, " \n" ); 00696 00697 for( j=0; j < mole.num_comole_calc; j++ ) 00698 { 00699 fprintf( ioQQQ, " %4.4s", COmole[j]->label ); 00700 fprintf( ioQQQ, " " ); 00701 00702 for( i=8; i < mole.num_comole_calc; i++ ) 00703 { 00704 fprintf( ioQQQ, "%9.1e", mole.amat[i][j]* COmole[i]->hevmol ); 00705 } 00706 fprintf( ioQQQ, "\n" ); 00707 } 00708 } 00709 00710 fprintf( ioQQQ, " Mole dens:" ); 00711 for( j=0; j < mole.num_comole_calc; j++ ) 00712 { 00713 fprintf( ioQQQ, " %4.4s:%.2e", COmole[j]->label 00714 , COmole[j]->hevmol ); 00715 } 00716 fprintf( ioQQQ, " \n" ); 00717 } 00718 00719 /* heating due to CO photodissociation */ 00720 co.CODissHeat = (realnum)(CO_findrate("PHOTON,CO=>C,O")*1e-12); 00721 00722 thermal.heating[0][9] = co.CODissHeat; 00723 00724 /* the real multi-level model molecule */ 00725 abundan = findspecies("CO")->hevmol; 00726 /* IonStg and nelem were set to 2, 0 in makelevlines */ 00727 ASSERT( C12O16Rotate[0].Hi->IonStg < LIMELM+2 ); 00728 ASSERT( C12O16Rotate[0].Hi->nelem-1 < LIMELM+2 ); 00729 dense.xIonDense[ C12O16Rotate[0].Hi->nelem-1][C12O16Rotate[0].Hi->IonStg-1] = abundan; 00730 00731 CO_PopsEmisCool(&C12O16Rotate, nCORotate,abundan , "12CO", 00732 &CoolHeavy.C12O16Rot,&CoolHeavy.dC12O16Rot ); 00733 /*if( nzone>400 ) 00734 fprintf(ioQQQ,"DEBUG co cool\t%.2f\t%.4e\t%li\t%.4e\t%.4e\n", 00735 fnzone,CoolHeavy.C12O16Rot,nCORotate,abundan,phycon.te );*/ 00736 00737 abundan = findspecies("CO")->hevmol/co.C12_C13_isotope_ratio; 00738 /* IonStg and nelem were set to 3, 0 in makelevlines */ 00739 ASSERT( C13O16Rotate[0].Hi->IonStg < LIMELM+2 ); 00740 ASSERT( C13O16Rotate[0].Hi->nelem-1 < LIMELM+2 ); 00741 dense.xIonDense[ C13O16Rotate[0].Hi->nelem-1][C13O16Rotate[0].Hi->IonStg-1] = abundan; 00742 00743 CO_PopsEmisCool(&C13O16Rotate, nCORotate,abundan ,"13CO", 00744 &CoolHeavy.C13O16Rot,&CoolHeavy.dC13O16Rot ); 00745 00746 /* now set total density of each element locked in gas phase */ 00747 for( i=0;i<LIMELM; ++i ) 00748 { 00749 dense.xMolecules[i] = 0.; 00750 } 00751 00752 /* total number of H per unit vol in molecules, 00753 * of course not including H0/H+ */ 00754 for(i=0;i<N_H_MOLEC;i++) 00755 { 00756 dense.xMolecules[ipHYDROGEN] += hmi.Hmolec[i]*hmi.nProton[i]; 00757 } 00758 dense.xMolecules[ipHYDROGEN] -= (hmi.Hmolec[ipMH] + hmi.Hmolec[ipMHp]); 00759 00760 /* >>chng 02 sep 05, dense.xMolecules[[ipHYDROGEN] is the part of H 00761 * that is not computed in hmole 00762 * add in gas phase abundances locked in molecules 00763 * H is special since treated in ion_solver with all hmole molecules 00764 * xMolecules are the densities of these species that are done in co, 00765 * this sum is only up to mole.num_heavy_molec and so does not include the atoms/ions */ 00766 dense.H_sum_in_CO = 0; 00767 for(i=0; i<mole.num_comole_calc; ++i) 00768 { 00769 if(COmole[i]->n_nuclei <= 1) 00770 continue; 00771 /* special sum of H in CO network */ 00772 dense.H_sum_in_CO += COmole[i]->nElem[ipHYDROGEN]*COmole[i]->hevmol; 00773 00774 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00775 { 00776 if( mole.lgElem_in_chemistry[nelem] ) 00777 { 00778 if( COmole[i]->hevmol > BIGFLOAT ) 00779 { 00780 fprintf(ioQQQ, "PROBLEM DISASTER mole_co_solve has found " 00781 "a CO-network molecule with an insane abundance.\n"); 00782 fprintf(ioQQQ, "Species number %li with abundance %.2e\n", 00783 i , COmole[i]->hevmol); 00784 TotalInsanity(); 00785 } 00786 else if( conv.lgSearch && 00787 COmole[i]->nElem[nelem]*COmole[i]->hevmol > dense.gas_phase[nelem]) 00788 { 00789 if( !conv.lgSearch ) 00790 fprintf(ioQQQ, 00791 "PROBLEM COmole limit trip call %li Seatch? %c nMole %li " 00792 "n(mole) %.2e n(gas) %.2e)\n", 00793 conv.nTotalIoniz, 00794 TorF(conv.lgSearch),i , 00795 COmole[i]->nElem[nelem]*COmole[i]->hevmol, 00796 dense.gas_phase[nelem]); 00797 COmole[i]->hevmol = dense.gas_phase[nelem]/COmole[i]->nElem[nelem]/2.f; 00798 } 00799 dense.xMolecules[nelem] +=COmole[i]->nElem[nelem]*COmole[i]->hevmol; 00800 } 00801 } 00802 } 00803 00804 /* This takes the average of an element's atomic and singly ionized density from ion_solver 00805 and comole and sets the solution to the updated COmole[i]->hevmol value. The macro mole.num_heavy_molec 00806 is the number of heavy element molecules in the network. In addition there are 2*mole.num_elements 00807 in the network, half atomic elements and half ionized elements. Finally, the total number of species 00808 in the network, including elements, equals mole.num_comole_calc. The following will first take the average of 00809 the ionized species (between mole.num_heavy_molec and mole.num_heavy_molec+mole.num_elements) and the other statement will 00810 take the average of the atomic species (between mole.num_heavy_molec+mole.num_elements and mole.num_comole_calc). 00811 This is generalized, so that if one wants to insert, for example, a sodium chemistry, once the macro's 00812 are updated this if statement immediately works */ 00813 00814 /* >> chng 2006 Sep 02 rjrw: use vectors rather than offsets above to index, for ease of generalization */ 00815 00816 for(i=0;i<mole.num_elements; ++i) 00817 { 00818 element = chem_element[i]; 00819 /*printf("MOL\t%3e\tION\t%e\tSPECIES\t%s\n", COmole[i]->hevmol,dense.xIonDense[COmole[i]->nelem_hevmol][1],COmole[i].label);*/ 00820 COmole[element->ipMlP]->hevmol = ((dense.xIonDense[element->ipCl][1]+ COmole[element->ipMlP]->hevmol)/2)*dense.lgElmtOn[element->ipCl]; 00821 /*printf("MOL\t%3e\tION\t%e\tSPECIES\t%s\n", COmole[i]->hevmol,dense.xIonDense[COmole[i]->nelem_hevmol][0],COmole[i].label);*/ 00822 COmole[element->ipMl]->hevmol = ((dense.xIonDense[element->ipCl][0]+ COmole[element->ipMl]->hevmol)/2)*dense.lgElmtOn[element->ipCl]; 00823 00824 /* check for NaN */ 00825 ASSERT( !isnan( COmole[element->ipMlP]->hevmol ) ); 00826 ASSERT( !isnan( COmole[element->ipMl]->hevmol ) ); 00827 } 00828 00829 /* check whether ion and chem solvers agree yet */ 00830 if( dense.lgElmtOn[ipCARBON] && 00831 fabs(dense.xIonDense[ipCARBON][0]- findspecies("C")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 ) 00832 { 00833 conv.lgConvIoniz = false; 00834 sprintf( conv.chConvIoniz, "CO C0 con"); 00835 conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][0]; 00836 conv.BadConvIoniz[1] = findspecies("C")->hevmol; 00837 } 00838 else if( dense.lgElmtOn[ipCARBON] && 00839 fabs(dense.xIonDense[ipCARBON][1]- findspecies("C+")->hevmol)/SDIV(dense.gas_phase[ipCARBON]) >0.001 ) 00840 { 00841 conv.lgConvIoniz = false; 00842 sprintf( conv.chConvIoniz, "CO C1 con"); 00843 conv.BadConvIoniz[0] = dense.xIonDense[ipCARBON][1]; 00844 conv.BadConvIoniz[1] = findspecies("C+")->hevmol; 00845 } 00846 else if( dense.lgElmtOn[ipOXYGEN] && 00847 fabs(dense.xIonDense[ipOXYGEN][0]- findspecies("O")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 ) 00848 { 00849 conv.lgConvIoniz = false; 00850 sprintf( conv.chConvIoniz, "CO O0 con"); 00851 conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][0]; 00852 conv.BadConvIoniz[1] = findspecies("O")->hevmol; 00853 } 00854 else if( dense.lgElmtOn[ipOXYGEN] && 00855 fabs(dense.xIonDense[ipOXYGEN][1]- findspecies("O+")->hevmol)/SDIV(dense.gas_phase[ipOXYGEN]) >0.001 ) 00856 { 00857 conv.lgConvIoniz = false; 00858 sprintf( conv.chConvIoniz, "CO O1 con"); 00859 conv.BadConvIoniz[0] = dense.xIonDense[ipOXYGEN][1]; 00860 conv.BadConvIoniz[1] = findspecies("O+")->hevmol; 00861 } 00862 00863 /* now update ionization distribution of the elements we just did */ 00864 for(i=0;i<mole.num_elements; ++i) 00865 { 00866 element = chem_element[i]; 00867 dense.xIonDense[element->ipCl][0] = COmole[element->ipMl]->hevmol*dense.lgElmtOn[element->ipCl]; 00868 dense.xIonDense[element->ipCl][1] = COmole[element->ipMlP]->hevmol*dense.lgElmtOn[element->ipCl]; 00869 dense.IonLow[element->ipCl] = 0; 00870 dense.IonHigh[element->ipCl] = MAX2( dense.IonHigh[element->ipCl] , 1 ); 00871 } 00872 00873 /* if populations not conserved then not converged */ 00874 # define EPS_MOLE 0.1 00875 if( dense.lgElmtOn[ipHYDROGEN] && 00876 dense.xMolecules[ipHYDROGEN] > dense.gas_phase[ipHYDROGEN]*(1.+EPS_MOLE) ) 00877 { 00878 conv.lgConvIoniz = false; 00879 sprintf( conv.chConvIoniz, "COcon%2i",ipHYDROGEN ); 00880 conv.BadConvIoniz[0] = dense.xMolecules[ipHYDROGEN]; 00881 conv.BadConvIoniz[1] = dense.gas_phase[ipHYDROGEN]; 00882 } 00883 else if( dense.lgElmtOn[ipCARBON] && 00884 dense.xMolecules[ipCARBON] > dense.gas_phase[ipCARBON]*(1.+EPS_MOLE) ) 00885 { 00886 conv.lgConvIoniz = false; 00887 sprintf( conv.chConvIoniz, "COcon%2i",ipCARBON ); 00888 conv.BadConvIoniz[0] = dense.xMolecules[ipCARBON]; 00889 conv.BadConvIoniz[1] = dense.gas_phase[ipCARBON]; 00890 } 00891 else if( dense.lgElmtOn[ipOXYGEN] && 00892 dense.xMolecules[ipOXYGEN] > dense.gas_phase[ipOXYGEN]*(1.+EPS_MOLE) ) 00893 { 00894 conv.lgConvIoniz = false; 00895 sprintf( conv.chConvIoniz, "COcon%2i",ipOXYGEN ); 00896 conv.BadConvIoniz[0] = dense.xMolecules[ipOXYGEN]; 00897 conv.BadConvIoniz[1] = dense.gas_phase[ipOXYGEN]; 00898 } 00899 else if( dense.lgElmtOn[ipSILICON] && 00900 dense.xMolecules[ipSILICON] > dense.gas_phase[ipSILICON]*(1.+EPS_MOLE) ) 00901 { 00902 conv.lgConvIoniz = false; 00903 sprintf( conv.chConvIoniz, "COcon%2i",ipSILICON ); 00904 conv.BadConvIoniz[0] = dense.xMolecules[ipSILICON]; 00905 conv.BadConvIoniz[1] = dense.gas_phase[ipSILICON]; 00906 } 00907 else if( dense.lgElmtOn[ipSULPHUR] && 00908 dense.xMolecules[ipSULPHUR] > dense.gas_phase[ipSULPHUR]*(1.+EPS_MOLE) ) 00909 { 00910 conv.lgConvIoniz = false; 00911 sprintf( conv.chConvIoniz, "COcon%2i",ipSULPHUR ); 00912 conv.BadConvIoniz[0] = dense.xMolecules[ipSULPHUR]; 00913 conv.BadConvIoniz[1] = dense.gas_phase[ipSULPHUR]; 00914 } 00915 else if( dense.lgElmtOn[ipNITROGEN] && 00916 dense.xMolecules[ipNITROGEN] > dense.gas_phase[ipNITROGEN]*(1.+EPS_MOLE) ) 00917 { 00918 conv.lgConvIoniz = false; 00919 sprintf( conv.chConvIoniz, "COcon%2i",ipNITROGEN ); 00920 conv.BadConvIoniz[0] = dense.xMolecules[ipNITROGEN]; 00921 conv.BadConvIoniz[1] = dense.gas_phase[ipNITROGEN]; 00922 } 00923 00924 else if( dense.lgElmtOn[ipCHLORINE] && 00925 dense.xMolecules[ipCHLORINE] > dense.gas_phase[ipCHLORINE]*(1.+EPS_MOLE) ) 00926 { 00927 conv.lgConvIoniz = false; 00928 sprintf( conv.chConvIoniz, "COcon%2i",ipCHLORINE ); 00929 conv.BadConvIoniz[0] = dense.xMolecules[ipCHLORINE]; 00930 conv.BadConvIoniz[1] = dense.gas_phase[ipCHLORINE]; 00931 } 00932 00933 /* the total N photo dissociation rate, cm-3 s-1 */ 00934 co.nitro_dissoc_rate = (realnum)(CO_dissoc_rate("N")); 00935 return; 00936 00937 /*lint +e550 */ 00938 /*lint +e778 constant expression evaluates to 0 in operation '-' */ 00939 } 00940 00941 # undef EPS_MOLE