00001
00002
00003
00004
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 "map.h"
00015 #include "coolheavy.h"
00016 #include "mole.h"
00017 #include "mole_co_priv.h"
00018
00019
00020 static bool lgMoleZeroed=true;
00021
00022
00023 static double h2lim;
00024
00025
00026
00027 static const double COTOLER_MOLAV = 0.10;
00028
00029
00030 static const int CODEBUG = 0;
00031
00032
00033 static const int LUPMAX_CODRIV = 100;
00034
00035
00036
00037 static bool lgMolecAver(
00038
00039
00040 const char chJob[10] ,
00041 char *chReason );
00042
00043
00044 void CO_drive(void)
00045 {
00046 bool lgNegPop,
00047 lgZerPop,
00048 lgFirstTry;
00049 long int i;
00050 bool lgConverged;
00051
00052 long int loop;
00053 long int nelem , ion;
00054
00055
00056
00057 char chReason[100];
00058
00059
00060 DEBUG_ENTRY( "CO_drive()" );
00061
00062
00063
00064 if( hmi.lgNoH2Mole )
00065 {
00066
00067 h2lim = 0.;
00068 }
00069 else
00070 {
00071
00072
00073
00074
00075 h2lim = 1e-15;
00076
00077
00078
00079
00080
00081 h2lim = 1e-12;
00082 }
00083
00084
00085
00086
00087
00088
00089 co.lgCODoCalc = true;
00090
00091 if( phycon.te > 20000. )
00092 co.lgCODoCalc = false;
00093
00094
00095 if( co.lgNoCOMole )
00096 co.lgCODoCalc = false;
00097
00098
00099 if( dense.lgElmtOn[ipCARBON]*dense.lgElmtOn[ipOXYGEN] == 0 )
00100 co.lgCODoCalc = false;
00101
00102
00103
00104 if( dense.IonLow[ipCARBON]>0 || dense.IonLow[ipOXYGEN]>0 )
00105 co.lgCODoCalc = false;
00106
00107
00108 if( iteration!=co.iteration_co &&
00109 hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
00110 co.lgCODoCalc = false;
00111
00112
00113
00114
00115 if( dense.xIonDense[ipHYDROGEN][1] / dense.gas_phase[ipHYDROGEN] > 0.5 )
00116 co.lgCODoCalc = false;
00117
00118 if( dense.lgElmtOn[ipSILICON] )
00119 {
00120
00121
00122
00123 if( dense.xIonDense[ipSILICON][0]/dense.gas_phase[ipSILICON] < 1e-15 )
00124 co.lgCODoCalc = false;
00125 }
00126
00127
00128
00129 if( !co.lgCODoCalc )
00130 {
00131
00132
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] += (float)(rate->rk*rate->rate_species[1]->hevmol);
00138 rate = CO_findrate_s("Si+,Mg=>Si,Mg+");
00139 mole.xMoleChTrRate[ipSILICON][1][0] += (float)(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] += (float)(rate->rk*rate->rate_species[1]->hevmol);
00144 rate = CO_findrate_s("S+,Mg=>S,Mg+");
00145 mole.xMoleChTrRate[ipSULPHUR][1][0] += (float)(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] += (float)(rate->rk*rate->rate_species[1]->hevmol);
00150 rate = CO_findrate_s("C+,Mg=>C,Mg+");
00151 mole.xMoleChTrRate[ipCARBON][1][0] += (float)(rate->rk*rate->rate_species[1]->hevmol);
00152
00153 #if 0
00154 mole.xMoleChTrRate[ipSILICON][1][0] = (float)(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] = (float)(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] = (float)(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
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
00176 thermal.heating[0][9] = 0.;
00177
00178 CoolHeavy.C12O16Rot = 0.;
00179 CoolHeavy.dC12O16Rot = 0.;
00180 CoolHeavy.C13O16Rot = 0.;
00181 CoolHeavy.dC13O16Rot = 0.;
00182 co.CODissHeat = 0.;
00183
00184 for( i=0; i < nCORotate; i++ )
00185 {
00186 C12O16Rotate[i].PopLo = 0.;
00187 C13O16Rotate[i].PopLo = 0.;
00188
00189 C12O16Rotate[i].PopHi = 0.;
00190 C13O16Rotate[i].PopHi = 0.;
00191
00192 C12O16Rotate[i].PopOpc = 0.;
00193 C13O16Rotate[i].PopOpc = 0.;
00194
00195 C12O16Rotate[i].cool = 0.;
00196 C12O16Rotate[i].heat = 0.;
00197 C13O16Rotate[i].cool = 0.;
00198 C13O16Rotate[i].heat = 0.;
00199
00200 C12O16Rotate[i].xIntensity = 0.;
00201 C13O16Rotate[i].xIntensity = 0.;
00202
00203 C12O16Rotate[i].phots = 0.;
00204 C13O16Rotate[i].phots = 0.;
00205 }
00206 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00207 {
00208 for( ion=0; ion<nelem+2; ++ion )
00209 {
00210 mole.source[nelem][ion] = 0.;
00211 mole.sink[nelem][ion] = 0.;
00212 }
00213 }
00214 }
00215
00216 if( trace.nTrConvg>=4 )
00217 {
00218
00219 fprintf( ioQQQ,
00220 " CO_drive4 do not evaluate CO chemistry.\n");
00221 }
00222
00223 DEBUG_EXIT( "CO_drive()" );
00224 return;
00225 }
00226
00227 CO_update_species_cache();
00228
00229
00230 CO_update_rks();
00231
00232
00233 loop = 0;
00234 lgNegPop = false;
00235 lgConverged = lgMolecAver("SETUP",chReason);
00236 do
00237 {
00238
00239 if( trace.nTrConvg>=5 )
00240 {
00241
00242 fprintf( ioQQQ,
00243 " CO_drive5 CO chemistry loop %li chReason %s.\n" , loop, chReason );
00244 }
00245
00246
00247 CO_solve(&lgNegPop,&lgZerPop );
00248
00249
00250 lgConverged = lgMolecAver("AVER",chReason);
00251 if( CODEBUG )
00252 fprintf(ioQQQ,"codrivbug %.2f %li Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n",
00253 fnzone ,
00254 loop ,
00255 TorF(lgNegPop) ,
00256 TorF(lgZerPop) ,
00257 findspecies("OH")->hevmol ,
00258 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00259 phycon.te,
00260 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00261
00262 ++loop;
00263 } while (
00264
00265
00266 (loop < LUPMAX_CODRIV) && !lgConverged );
00267
00268
00269 lgMoleZeroed = false;
00270
00271
00272 if( loop == LUPMAX_CODRIV)
00273 {
00274 conv.lgConvIoniz = false;
00275 strcpy( conv.chConvIoniz, "CO con1" );
00276 if( CODEBUG )
00277 fprintf(ioQQQ,"CO_drive not converged\n");
00278 }
00279
00280
00281
00282
00283 lgFirstTry = (nzone==co.co_nzone && iteration==co.iteration_co);
00284
00285
00286 if( lgNegPop )
00287 {
00288 if( conv.lgSearch && (hmi.H2_total/dense.gas_phase[ipHYDROGEN] <h2lim) &&
00289 (iteration==1) )
00290 {
00291
00292
00293
00294
00295
00296
00297 CO_zero();
00298 }
00299 else
00300 {
00301 if( called.lgTalk && !lgFirstTry )
00302 {
00303 conv.lgConvPops = false;
00304 fprintf( ioQQQ,
00305 " PROBLEM CO_drive failed to converge1 due to negative pops, zone=%.2f, CO/C=%.2e OH/C=%.2e H2/H=%.2e\n",
00306 fnzone,
00307 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00308 findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00309 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00310 ConvFail( "pops" , "CO" );
00311 }
00312 }
00313 }
00314
00315
00316 else if( loop == LUPMAX_CODRIV )
00317 {
00318
00319 if( called.lgTalk && !lgFirstTry )
00320 {
00321 conv.lgConvPops = false;
00322 fprintf( ioQQQ,
00323 " PROBLEM CO_drive failed to converge2 in %li iter, zone=%.2f, CO/C=%.2e negpop?%1c reason:%s\n",
00324 loop,
00325 fnzone,
00326 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00327 TorF(lgNegPop) ,
00328 chReason);
00329 ConvFail( "pops" , "CO" );
00330 }
00331 }
00332
00333
00334 ASSERT(conv.lgSearch || findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) <= 1.01 ||
00335
00336 loop == LUPMAX_CODRIV );
00337
00338
00339
00340
00341
00342
00343 if(0) {
00344 double sum_ion , sum_mol;
00345 sum_ion = dense.xIonDense[ipCARBON][0] + dense.xIonDense[ipCARBON][1];
00346 sum_mol = findspecies("C")->hevmol + findspecies("C+")->hevmol;
00347 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
00348 {
00349
00350
00351
00352 conv.lgConvIoniz = false;
00353 strcpy( conv.chConvIoniz, "CO con2" );
00354 conv.BadConvIoniz[0] = sum_ion;
00355 conv.BadConvIoniz[1] = sum_mol;
00356 if( CODEBUG )
00357 fprintf(ioQQQ,"CO_drive not converged\n");
00358 }
00359 sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1];
00360 sum_mol = findspecies("O")->hevmol + findspecies("O+")->hevmol;
00361 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
00362 {
00363
00364
00365
00366 conv.lgConvIoniz = false;
00367 strcpy( conv.chConvIoniz, "CO con3" );
00368 conv.BadConvIoniz[0] = sum_ion;
00369 conv.BadConvIoniz[1] = sum_mol;
00370 if( CODEBUG )
00371 fprintf(ioQQQ,"CO_drive not converged\n");
00372 }
00373 }
00374
00375 DEBUG_EXIT( "CO_drive()" );
00376 return;
00377 }
00378
00379 static bool lgMolecAver(
00380
00381
00382 const char chJob[10] ,
00383 char *chReason )
00384 {
00385 long int i;
00386 bool lgReturn;
00387
00388 static float hden_save_prev_call;
00389 double conv_fac;
00390 struct chem_element_s *element;
00391
00392 DEBUG_ENTRY( "lgMolecAver()" );
00393 const bool DEBUG_MOLECAVER = false;
00394
00395
00396 strcpy( chReason , "none given" );
00397
00398
00399 if( strcmp( "SETUP" , chJob ) == 0 )
00400 {
00401 static float hden_save_prev_iter;
00402
00403
00404
00405
00406
00407
00408 if( iteration!=co.iteration_co &&
00409 hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
00410 {
00411
00412 DEBUG_EXIT( "lgMolecAver()" );
00413
00414 lgReturn = true;
00415 return lgReturn;
00416 }
00417
00418 else if( co.iteration_co < 0 || lgMoleZeroed )
00419 {
00420
00421
00422
00423
00424
00425
00426
00427 if( !map.lgMapBeingDone || lgMoleZeroed )
00428 {
00429 for( i=0; i < mole.num_comole_calc; i++ )
00430 {
00431 COmole[i]->hevmol = dense.xIonDense[ipCARBON][0]*COmole[i]->pdr_mole_co;
00432 COmole[i]->co_save = COmole[i]->hevmol;
00433 }
00434 }
00435
00436
00437 ASSERT( dense.xIonDense[ipCARBON][0]>SMALLFLOAT );
00438
00439
00440 for(i=0;i<mole.num_elements;i++) {
00441 element = chem_element[i];
00442 COmole[element->ipMl]->hevmol = dense.xIonDense[element->ipCl][0];
00443 COmole[element->ipMlP]->hevmol = dense.xIonDense[element->ipCl][1];
00444 }
00445
00446
00447 co.iteration_co = iteration;
00448 co.co_nzone = nzone;
00449
00450
00451
00452
00453 hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
00454 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
00455
00456 if( DEBUG_MOLECAVER )
00457 fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li zeroing since very first call. H2/H=%.2e\n",
00458 iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00459 }
00460 else if( iteration > co.iteration_co )
00461 {
00462 float ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_iter;
00463
00464
00465
00466 for( i=0; i < mole.num_comole_calc; i++ )
00467 {
00468 COmole[i]->hevmol = COmole[i]->hev_reinit*ratio;
00469 COmole[i]->co_save = COmole[i]->hev_reinit*ratio;
00470 }
00471
00472 co.iteration_co = iteration;
00473 co.co_nzone = nzone;
00474
00475 if( DEBUG_MOLECAVER )
00476 fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n",
00477 iteration,
00478 nzone,
00479 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00480 }
00481 else if( iteration == co.iteration_co && nzone==co.co_nzone+1 )
00482 {
00483
00484
00485 for( i=0; i < mole.num_comole_calc; i++ )
00486 {
00487 COmole[i]->hev_reinit = COmole[i]->hevmol;
00488 }
00489
00490 co.co_nzone = -2;
00491 hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
00492 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
00493
00494 if( DEBUG_MOLECAVER )
00495 fprintf(ioQQQ,"DEBUG lgMolecAver iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone);
00496 }
00497
00498
00499 lgReturn = true;
00500 }
00501
00502 else if( strcmp( "AVER" , chJob ) == 0 )
00503 {
00504
00505
00506
00507 float oldoh = findspecies("OH")->co_save;
00508 lgReturn = true;
00509
00510
00511
00512
00513 for( i=0; i < mole.num_comole_calc; i++ )
00514 {
00515
00516
00517 float damper = 0.2f;
00518
00519
00520
00521
00522 COmole[i]->co_save *= dense.gas_phase[ipHYDROGEN] / hden_save_prev_call;
00523
00524
00525
00526 if( COmole[i]->co_save< SMALLFLOAT )
00527 {
00528 COmole[i]->co_save = COmole[i]->hevmol;
00529 }
00530 else
00531 {
00532
00533
00534
00535 double ratio = (double)COmole[i]->hevmol/(double)COmole[i]->co_save;
00536 ASSERT( COmole[i]->co_save>0. && COmole[i]->hevmol>=0. );
00537 if( fabs(ratio-1.) < 1.5 ||
00538 fabs(ratio-1.) > 0.66 )
00539 {
00540
00541 COmole[i]->hevmol = COmole[i]->hevmol*(1.f - damper) +
00542 COmole[i]->co_save*damper;
00543 }
00544 else
00545 {
00546
00547 COmole[i]->hevmol = (float)pow(10. ,
00548 (log10(COmole[i]->hevmol) + log10(COmole[i]->co_save))/2. );
00549 }
00550 COmole[i]->co_save = COmole[i]->hevmol;
00551 }
00552 }
00553
00554 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
00555
00556
00557
00558 if( oldoh/SDIV(dense.gas_phase[ipOXYGEN]) < 1e-10 )
00559 conv_fac = 3.;
00560 else
00561 conv_fac = 1.;
00562
00563 if(fabs(oldoh/SDIV(findspecies("OH")->hevmol)-1.) < COTOLER_MOLAV*conv_fac )
00564 {
00565 lgReturn = true;
00566 }
00567 else
00568 {
00569
00570 lgReturn = false;
00571 sprintf( chReason , "OH change from %.3e to %.3e",
00572 oldoh ,
00573 findspecies("OH")->hevmol);
00574 }
00575 }
00576 else
00577 TotalInsanity();
00578
00579
00580 if( findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) > 1.+FLT_EPSILON )
00581 {
00582 lgReturn = false;
00583 sprintf( chReason , "CO/C >1, value is %e",
00584 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
00585 }
00586
00587 if( (trace.lgTrace && trace.lgTr_CO_Mole) || DEBUG_MOLECAVER )
00588 {
00589 fprintf( ioQQQ, " DEBUG lgMolecAver converged? %c" ,TorF(lgReturn) );
00590 fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
00591 fprintf( ioQQQ, " OH/O=%9.1e", findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]) );
00592 fprintf( ioQQQ, "\n" );
00593 }
00594
00595 DEBUG_EXIT( "lgMolecAver()" );
00596
00597 return lgReturn;
00598 }