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 /*lgMolecAver average old and new molecular equilibrium balance from CO_solve */ 00004 /*CO_drive - public routine, calls CO_solve to converge molecules */ 00005 #include "cddefines.h" 00006 #include "taulines.h" 00007 #include "dense.h" 00008 #include "hmi.h" 00009 #include "conv.h" 00010 #include "phycon.h" 00011 #include "trace.h" 00012 #include "thermal.h" 00013 #include "called.h" 00014 #include "hcmap.h" 00015 #include "coolheavy.h" 00016 #include "mole.h" 00017 #include "mole_co_priv.h" 00018 00019 /* says whether the co network is currently set to zero */ 00020 static bool lgMoleZeroed=true; 00021 00022 /* the limit to H2/H where we will solve for CO */ 00023 static double h2lim; 00024 00025 /* this is the relative change in OH, which is used as a convergence 00026 * criteria in lgMolecAver */ 00027 static const double COTOLER_MOLAV = 0.10; 00028 00029 /* flag to control debug statement, if 0 then none, just loops when 1, more when 2 */ 00030 static const int CODEBUG = 0; 00031 /* this is the maximum number of loops CO_drive will go over while 00032 * trying to converge the molecular network */ 00033 static const int LUPMAX_CODRIV = 100; 00034 00035 /*lgMolecAver average old and new molecular equilibrium balance from CO_solve, 00036 * returns true if converged */ 00037 STATIC bool lgMolecAver( 00038 /* job == "SETUP", set things up during first call, 00039 * job == "AVER" take average of arrays */ 00040 const char chJob[10] , 00041 char *chReason ); 00042 00043 /*CO_drive main driver for heavy molecular equilibrium routines */ 00044 void CO_drive(void) 00045 { 00046 bool lgNegPop, 00047 lgZerPop, 00048 lgFirstTry; 00049 long int i; 00050 bool lgConverged; 00051 /* count how many times through loop */ 00052 long int loop; 00053 long int nelem , ion; 00054 00055 00056 /* this will give reason CO not converged */ 00057 char chReason[100]; 00058 /*static bool lgMoleZeroed=true;*/ 00059 00060 DEBUG_ENTRY( "CO_drive()" ); 00061 00062 /* h2lim is smallest ratio of H2/HDEN for which we will 00063 * even try to invert heavy element network */ 00064 if( hmi.lgNoH2Mole ) 00065 { 00066 /* H2 network is turned off, ignore this limit */ 00067 h2lim = 0.; 00068 } 00069 else 00070 { 00071 /* this limit is important since the co molecular network is first turned 00072 * on when this is hit. It's conservation law will only ever include the 00073 * initial O, O+, and C, C+, so must not start before nearly all 00074 * C and O is in these forms */ 00075 h2lim = 1e-15; 00076 /* >>chng 05 jul 15, from 1e-15 to 1e-12, see whether results are stable, 00077 * this does include CO in the H+ zone in orion_hii_pdr 00078 * a problem is that, during search phase, the first temp is 4000K and the 00079 * H2 abundance is larger than it will be at the illuminated face. try to 00080 * avoid turning H2 on too soon */ 00081 h2lim = 1e-12; 00082 } 00083 00084 /* do we want to try to calculate the CO? 00085 * >>chng 03 sep 07, add first test, so that we do not turn off 00086 * heavy element network if it has ever been turned on 00087 * >>chng 04 apr 23, change logic so never done when temp > 20000K */ 00088 /* make series of tests setting co.lgCODoCalc for simplicity */ 00089 co.lgCODoCalc = true; 00090 /* too hot - don't bother Te > 2e4K */ 00091 if( phycon.te > 20000. ) 00092 co.lgCODoCalc = false; 00093 00094 /* molecules have been turned off */ 00095 if( co.lgNoCOMole ) 00096 co.lgCODoCalc = false; 00097 00098 /* a critical element has been turned off */ 00099 if( dense.lgElmtOn[ipCARBON]*dense.lgElmtOn[ipOXYGEN] == 0 ) 00100 co.lgCODoCalc = false; 00101 00102 /* >>chng 04 jul 20, do not attempt co if no O or C atoms present, 00103 * since co net needs their atomic data */ 00104 if( dense.IonLow[ipCARBON]>0 || dense.IonLow[ipOXYGEN]>0 ) 00105 co.lgCODoCalc = false; 00106 00107 /* do not do if H2/Htot < h2lim which is set to 1e-12 by default above */ 00108 if( iteration!=co.iteration_co && 00109 hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim ) 00110 co.lgCODoCalc = false; 00111 00112 /* >>chng 06 sep 10, do not call CO network if H+/Htot is greater than 0.5 - 00113 * this change suggested by Nick Abel after cooling plasma developed CO 00114 * convergence problems when H+/H = 0.97 */ 00115 if( dense.xIonDense[ipHYDROGEN][1] / dense.gas_phase[ipHYDROGEN] > 0.5 ) 00116 co.lgCODoCalc = false; 00117 00118 if( dense.lgElmtOn[ipSILICON] ) 00119 { 00120 /* >>chng 04 jun 28, add test on atomic Si, some H II region 00121 * models crashed due to matrix failure when co network turned 00122 * on with atomic Si at very small abundance */ 00123 if( dense.xIonDense[ipSILICON][0]/dense.gas_phase[ipSILICON] < 1e-15 ) 00124 co.lgCODoCalc = false; 00125 } 00126 00127 /* this branch we will not do the CO equilibrium, set some variables that 00128 * would be calculated by the chemistry, zero others, and return */ 00129 if( !co.lgCODoCalc ) 00130 { 00131 /* these are heavy - heavy charge transfer rate that will still be needed 00132 * for recombination of Si+, S+, and C+ */ 00133 struct COmole_rate_s *rate; 00134 00135 mole.xMoleChTrRate[ipSILICON][1][0] = 0.; 00136 rate = CO_findrate_s("Si+,Fe=>Si,Fe+"); 00137 mole.xMoleChTrRate[ipSILICON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol); 00138 rate = CO_findrate_s("Si+,Mg=>Si,Mg+"); 00139 mole.xMoleChTrRate[ipSILICON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol); 00140 00141 mole.xMoleChTrRate[ipSULPHUR][1][0] = 0.; 00142 rate = CO_findrate_s("S+,Fe=>S,Fe+"); 00143 mole.xMoleChTrRate[ipSULPHUR][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol); 00144 rate = CO_findrate_s("S+,Mg=>S,Mg+"); 00145 mole.xMoleChTrRate[ipSULPHUR][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol); 00146 00147 mole.xMoleChTrRate[ipCARBON][1][0] = 0.; 00148 rate = CO_findrate_s("C+,Fe=>C,Fe+"); 00149 mole.xMoleChTrRate[ipCARBON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol); 00150 rate = CO_findrate_s("C+,Mg=>C,Mg+"); 00151 mole.xMoleChTrRate[ipCARBON][1][0] += (realnum)(rate->rk*rate->rate_species[1]->hevmol); 00152 00153 #if 0 00154 mole.xMoleChTrRate[ipSILICON][1][0] = (realnum)(dense.xIonDense[ipMAGNESIUM][0]* 00155 COmole_rate[ipR_SiP_Mg_Si_MgP].rk + 00156 dense.xIonDense[ipIRON][0]*COmole_rate[ipR_SiP_Fe_Si_FeP].rk); 00157 00158 mole.xMoleChTrRate[ipSULPHUR][1][0] = (realnum)(dense.xIonDense[ipMAGNESIUM][0]* 00159 COmole_rate[ipR_SP_Mg_S_MgP].rk + 00160 dense.xIonDense[ipIRON][0]*COmole_rate[ipR_SP_Fe_S_FeP].rk); 00161 00162 mole.xMoleChTrRate[ipCARBON][1][0] = (realnum)(dense.xIonDense[ipMAGNESIUM][0]* 00163 COmole_rate[ipR_CP_Mg_C_MgP].rk + 00164 dense.xIonDense[ipIRON][0]*COmole_rate[ipR_CP_Fe_C_FeP].rk); 00165 #endif 00166 00167 /* do we need to zero out the arrays and vars? */ 00168 if( !lgMoleZeroed ) 00169 { 00170 lgMoleZeroed = true; 00171 for( i=0; i<mole.num_comole_calc; ++i ) 00172 { 00173 COmole[i]->hevmol = 0.; 00174 } 00175 /* heating due to CO photodissociation */ 00176 thermal.heating[0][9] = 0.; 00177 /* CO cooling */ 00178 CoolHeavy.C12O16Rot = 0.; 00179 CoolHeavy.dC12O16Rot = 0.; 00180 CoolHeavy.C13O16Rot = 0.; 00181 CoolHeavy.dC13O16Rot = 0.; 00182 co.CODissHeat = 0.; 00183 00185 for( i=0; i < nCORotate; i++ ) 00186 { 00187 C12O16Rotate[i].Lo->Pop = 0.; 00188 C13O16Rotate[i].Lo->Pop = 0.; 00189 /* population of upper level */ 00190 C12O16Rotate[i].Hi->Pop = 0.; 00191 C13O16Rotate[i].Hi->Pop = 0.; 00192 /* population of lower level with correction for stimulated emission */ 00193 C12O16Rotate[i].Emis->PopOpc = 0.; 00194 C13O16Rotate[i].Emis->PopOpc = 0.; 00195 /* following two heat exchange excitation, deexcitation */ 00196 C12O16Rotate[i].Coll.cool = 0.; 00197 C12O16Rotate[i].Coll.heat = 0.; 00198 C13O16Rotate[i].Coll.cool = 0.; 00199 C13O16Rotate[i].Coll.heat = 0.; 00200 /* intensity of line */ 00201 C12O16Rotate[i].Emis->xIntensity = 0.; 00202 C13O16Rotate[i].Emis->xIntensity = 0.; 00203 /* number of photons emitted in line */ 00204 C12O16Rotate[i].Emis->phots = 0.; 00205 C13O16Rotate[i].Emis->phots = 0.; 00206 } 00207 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem ) 00208 { 00209 for( ion=0; ion<nelem+2; ++ion ) 00210 { 00211 mole.source[nelem][ion] = 0.; 00212 mole.sink[nelem][ion] = 0.; 00213 } 00214 } 00215 } 00216 00217 if( trace.nTrConvg>=4 ) 00218 { 00219 /* trace ionization */ 00220 fprintf( ioQQQ, 00221 " CO_drive4 do not evaluate CO chemistry.\n"); 00222 } 00223 return; 00224 } 00225 00226 CO_update_species_cache(); /* Update densities of species controlled outside the network */ 00227 00228 /* Update the rate constants -- final generic version */ 00229 CO_update_rks(); 00230 00231 /* this flag is used to stop calculation due to negative abundances */ 00232 loop = 0; 00233 lgNegPop = false; 00234 lgConverged = lgMolecAver("SETUP",chReason); 00235 do 00236 { 00237 00238 if( trace.nTrConvg>=5 ) 00239 { 00240 /* trace ionization */ 00241 fprintf( ioQQQ, 00242 " CO_drive5 CO chemistry loop %li chReason %s.\n" , loop, chReason ); 00243 } 00244 00245 /*oh_h_o_h2ld = findspecies("OH")->hevmol;*/ 00246 CO_solve(&lgNegPop,&lgZerPop ); 00247 /* this one takes average first, then updates molecular densities in co.hevmol, 00248 * returns true if change in OH density is less than macro COTOLER_MOLAV set above*/ 00249 lgConverged = lgMolecAver("AVER",chReason); 00250 if( CODEBUG ) 00251 fprintf(ioQQQ,"codrivbug %.2f %li Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n", 00252 fnzone , 00253 loop , 00254 TorF(lgNegPop) , 00255 TorF(lgZerPop) , 00256 findspecies("OH")->hevmol , 00257 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 00258 phycon.te, 00259 hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 00260 00261 ++loop; 00262 } while ( 00263 /* these are a series of conditions to stop this loop - 00264 * has the limit to the number of loops been reached? */ 00265 (loop < LUPMAX_CODRIV) && !lgConverged ); 00266 00267 /* say that we have found a solution before */ 00268 lgMoleZeroed = false; 00269 00270 /* check whether we have converged */ 00271 if( loop == LUPMAX_CODRIV) 00272 { 00273 conv.lgConvIoniz = false; 00274 strcpy( conv.chConvIoniz, "CO con1" ); 00275 if( CODEBUG ) 00276 fprintf(ioQQQ,"CO_drive not converged\n"); 00277 } 00278 00279 /* this flag says whether this is the first zone we are trying 00280 * to converge the molecules - there will be problems during initial 00281 * search so do not print anything in this case */ 00282 lgFirstTry = (nzone==co.co_nzone && iteration==co.iteration_co); 00283 00284 /* did the molecule network have negative pops? */ 00285 if( lgNegPop ) 00286 { 00287 if( conv.lgSearch && (hmi.H2_total/dense.gas_phase[ipHYDROGEN] <h2lim) && 00288 (iteration==1) ) 00289 { 00290 /* we are in search phase during the first iteration, 00291 * and the H2/H ratio has fallen 00292 * substantially below the threshold for solving the CO network. 00293 * Turn it off */ 00294 /* >> chng 07 jan 10 rjrw: this was CO_Init(), but the comment suggests 00295 * it should really be CO_zero */ 00296 CO_zero(); 00297 } 00298 else 00299 { 00300 if( called.lgTalk && !lgFirstTry ) 00301 { 00302 conv.lgConvPops = false; 00303 fprintf( ioQQQ, 00304 " PROBLEM CO_drive failed to converge1 due to negative pops, zone=%.2f, CO/C=%.2e OH/C=%.2e H2/H=%.2e\n", 00305 fnzone, 00306 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 00307 findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 00308 hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 00309 ConvFail( "pops" , "CO" ); 00310 } 00311 } 00312 } 00313 00314 /* this test, hit upper limit to number of iterations */ 00315 else if( loop == LUPMAX_CODRIV ) 00316 { 00317 /* do not print this if we are in the first zone where molecules are turned on */ 00318 if( called.lgTalk && !lgFirstTry ) 00319 { 00320 conv.lgConvPops = false; 00321 fprintf( ioQQQ, 00322 " PROBLEM CO_drive failed to converge2 in %li iter, zone=%.2f, CO/C=%.2e negpop?%1c reason:%s\n", 00323 loop, 00324 fnzone, 00325 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 00326 TorF(lgNegPop) , 00327 chReason); 00328 ConvFail( "pops" , "CO" ); 00329 } 00330 } 00331 00332 /* make sure we have not used more than all the CO */ 00333 ASSERT(conv.lgSearch || findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) <= 1.01 || 00334 /* this is true when loop did not converge co */ 00335 loop == LUPMAX_CODRIV ); 00336 /*fprintf(ioQQQ,"ratioo o\t%c\t%.2f\t%f\n", TorF(conv.lgSearch), 00337 fnzone , co.hevmol[ipCO]/dense.gas_phase[ipCARBON] );*/ 00338 00339 /* these are the elements whose converge we check */ 00340 /* this is a self consistency check, for converged solution */ 00341 /* >>chng 04 dec 02, this test is turn back on - don't know why it was turned off */ 00342 if(0) { 00343 double sum_ion , sum_mol; 00344 sum_ion = dense.xIonDense[ipCARBON][0] + dense.xIonDense[ipCARBON][1]; 00345 sum_mol = findspecies("C")->hevmol + findspecies("C+")->hevmol; 00346 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 ) 00347 { 00348 /*fprintf(ioQQQ, 00349 "non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n", 00350 5, fnzone, sum_ion , sum_mol);*/ 00351 conv.lgConvIoniz = false; 00352 strcpy( conv.chConvIoniz, "CO con2" ); 00353 conv.BadConvIoniz[0] = sum_ion; 00354 conv.BadConvIoniz[1] = sum_mol; 00355 if( CODEBUG ) 00356 fprintf(ioQQQ,"CO_drive not converged\n"); 00357 } 00358 sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1]; 00359 sum_mol = findspecies("O")->hevmol + findspecies("O+")->hevmol; 00360 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 ) 00361 { 00362 /*fprintf(ioQQQ, 00363 "non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n", 00364 7, fnzone, sum_ion , sum_mol);*/ 00365 conv.lgConvIoniz = false; 00366 strcpy( conv.chConvIoniz, "CO con3" ); 00367 conv.BadConvIoniz[0] = sum_ion; 00368 conv.BadConvIoniz[1] = sum_mol; 00369 if( CODEBUG ) 00370 fprintf(ioQQQ,"CO_drive not converged\n"); 00371 } 00372 } 00373 return; 00374 } 00375 00376 STATIC bool lgMolecAver( 00377 /* job == "SETUP", set things up during first call, 00378 * job == "AVER" take average of arrays */ 00379 const char chJob[10] , 00380 char *chReason ) 00381 { 00382 long int i; 00383 bool lgReturn; 00384 /* this will be used to rescale old saved abundances in constant pressure model */ 00385 static realnum hden_save_prev_call; 00386 double conv_fac; 00387 struct chem_element_s *element; 00388 00389 DEBUG_ENTRY( "lgMolecAver()" ); 00390 const bool DEBUG_MOLECAVER = false; 00391 00392 /* this will become reason not converged */ 00393 strcpy( chReason , "none given" ); 00394 00395 /* when called with SETUP, just about to enter solver loop */ 00396 if( strcmp( "SETUP" , chJob ) == 0 ) 00397 { 00398 static realnum hden_save_prev_iter; 00399 00400 /* >>chng 04 apr 29, use PDR solution, scaled to density of neutral carbon, as first guess*/ 00401 /* CO_Init set this to -2 when code initialized, so negative 00402 * number shows very first call in this model */ 00403 /* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration 00404 * and we are below limit where it should be evaluated */ 00405 if( iteration!=co.iteration_co && 00406 hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim ) 00407 { 00408 00409 lgReturn = true; 00410 return lgReturn; 00411 } 00412 00413 else if( co.iteration_co < 0 || lgMoleZeroed ) 00414 { 00415 00416 /* very first attempt at ever obtaining a solution - 00417 * called one time per model since co.iteration_co set negative 00418 * when cdInit called */ 00419 00420 /* >>chng 05 jun 24, during map phase do not reset molecules to zero 00421 * since we probably have a better estimate right now */ 00422 if( !hcmap.lgMapBeingDone || lgMoleZeroed ) 00423 { 00424 for( i=0; i < mole.num_comole_calc; i++ ) 00425 { 00426 COmole[i]->hevmol = dense.xIonDense[ipCARBON][0]*COmole[i]->pdr_mole_co; 00427 COmole[i]->co_save = COmole[i]->hevmol; 00428 } 00429 } 00430 00431 /* we should have a neutral carbon solution at this point */ 00432 ASSERT( dense.xIonDense[ipCARBON][0]>SMALLFLOAT ); 00433 00434 /* first guess of the elements and charges come from ionization solvers */ 00435 for(i=0;i<mole.num_elements;i++) { 00436 element = chem_element[i]; 00437 COmole[element->ipMl]->hevmol = dense.xIonDense[element->ipCl][0]; 00438 COmole[element->ipMlP]->hevmol = dense.xIonDense[element->ipCl][1]; 00439 } 00440 00441 /* set iteration flag */ 00442 co.iteration_co = iteration; 00443 co.co_nzone = nzone; 00444 00445 /* for constant pressure models when molecules are reset on second 00446 * and higher iterations, total density will be different, so 00447 * must take this into account when abundances are reset */ 00448 hden_save_prev_iter = dense.gas_phase[ipHYDROGEN]; 00449 hden_save_prev_call = dense.gas_phase[ipHYDROGEN]; 00450 00451 if( DEBUG_MOLECAVER ) 00452 fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li zeroing since very first call. H2/H=%.2e\n", 00453 iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 00454 } 00455 else if( iteration > co.iteration_co ) 00456 { 00457 realnum ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_iter; 00458 00459 /* this is first call on new iteration, reset 00460 * to first initial abundances from previous iteration */ 00461 for( i=0; i < mole.num_comole_calc; i++ ) 00462 { 00463 COmole[i]->hevmol = COmole[i]->hev_reinit*ratio; 00464 COmole[i]->co_save = COmole[i]->hev_reinit*ratio; 00465 } 00466 00467 co.iteration_co = iteration; 00468 co.co_nzone = nzone; 00469 00470 if( DEBUG_MOLECAVER ) 00471 fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n", 00472 iteration, 00473 nzone, 00474 hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 00475 } 00476 else if( iteration == co.iteration_co && nzone==co.co_nzone+1 ) 00477 { 00478 /* this branch, second zone with solution, so save previous 00479 * zone's solution to reset things in next iteration */ 00480 for( i=0; i < mole.num_comole_calc; i++ ) 00481 { 00482 COmole[i]->hev_reinit = COmole[i]->hevmol; 00483 } 00484 00485 co.co_nzone = -2; 00486 hden_save_prev_iter = dense.gas_phase[ipHYDROGEN]; 00487 hden_save_prev_call = dense.gas_phase[ipHYDROGEN]; 00488 00489 if( DEBUG_MOLECAVER ) 00490 fprintf(ioQQQ,"DEBUG lgMolecAver iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone); 00491 } 00492 00493 /* didn't do anything, but say converged */ 00494 lgReturn = true; 00495 } 00496 00497 else if( strcmp( "AVER" , chJob ) == 0 ) 00498 { 00499 /* >>chng 04 jan 22, co.hevmol is current value, we want to save old 00500 * value, which is in co.co_save */ 00501 /*realnum oldoh = findspecies("OH")->hevmol;*/ 00502 realnum oldoh = findspecies("OH")->co_save; 00503 lgReturn = true; 00504 00505 /* get new numbers - take average of old and new */ 00506 /* >>chng 03 sep 16, only use damper for molecular species, 00507 * not ion/atoms */ 00508 for( i=0; i < mole.num_comole_calc; i++ ) 00509 { 00510 /* parameter to mix old and new, 00511 * damper is fraction of old solution to take */ 00512 realnum damper = 0.2f; 00513 00514 /* >>chng 04 sep 11, correct for den chng in var den models */ 00515 /* the ratio of densities is unity for constant density model, but in const pres or other 00516 * models, account for change in entire density scale */ 00517 COmole[i]->co_save *= dense.gas_phase[ipHYDROGEN] / hden_save_prev_call; 00518 00519 /* if co.co_save is zero then don't take average of it and 00520 * new solution */ 00521 if( COmole[i]->co_save< SMALLFLOAT ) 00522 { 00523 COmole[i]->co_save = COmole[i]->hevmol; 00524 } 00525 else 00526 { 00527 /* >>chng 04 jan 24, add logic to check how different old and 00528 * new solutions are. if quite different then take average 00529 * in the log */ 00530 double ratio = (double)COmole[i]->hevmol/(double)COmole[i]->co_save; 00531 ASSERT( COmole[i]->co_save>0. && COmole[i]->hevmol>=0. ); 00532 if( fabs(ratio-1.) < 1.5 || 00533 fabs(ratio-1.) > 0.66 ) 00534 { 00535 /* this branch moderate changes */ 00536 COmole[i]->hevmol = COmole[i]->hevmol*(1.f - damper) + 00537 COmole[i]->co_save*damper; 00538 } 00539 else 00540 { 00541 /* this branch - large changes take average in the log */ 00542 COmole[i]->hevmol = (realnum)pow(10. , 00543 (log10(COmole[i]->hevmol) + log10(COmole[i]->co_save))/2. ); 00544 } 00545 COmole[i]->co_save = COmole[i]->hevmol; 00546 } 00547 } 00548 /* remember the current H density */ 00549 hden_save_prev_call = dense.gas_phase[ipHYDROGEN]; 00550 00551 /* >>chng 05 jul 14, add logic to not be as fine a tolerance when 00552 * only trace amounts of molecules are present */ 00553 if( oldoh/SDIV(dense.gas_phase[ipOXYGEN]) < 1e-10 ) 00554 conv_fac = 3.; 00555 else 00556 conv_fac = 1.; 00557 00558 if(fabs(oldoh/SDIV(findspecies("OH")->hevmol)-1.) < COTOLER_MOLAV*conv_fac ) 00559 { 00560 lgReturn = true; 00561 } 00562 else 00563 { 00564 /* say we are not converged */ 00565 lgReturn = false; 00566 sprintf( chReason , "OH change from %.3e to %.3e", 00567 oldoh , 00568 findspecies("OH")->hevmol); 00569 } 00570 } 00571 else 00572 TotalInsanity(); 00573 00574 /* also not converged if co exceeds c */ 00575 if( findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) > 1.+FLT_EPSILON ) 00576 { 00577 lgReturn = false; 00578 sprintf( chReason , "CO/C >1, value is %e", 00579 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) ); 00580 } 00581 00582 if( (trace.lgTrace && trace.lgTr_CO_Mole) || DEBUG_MOLECAVER ) 00583 { 00584 fprintf( ioQQQ, " DEBUG lgMolecAver converged? %c" ,TorF(lgReturn) ); 00585 fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) ); 00586 fprintf( ioQQQ, " OH/O=%9.1e", findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]) ); 00587 fprintf( ioQQQ, "\n" ); 00588 } 00589 00590 return lgReturn; 00591 }