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 "hcmap.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] += (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
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 co.CODissHeat = 0.;
00178
00179 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00180 {
00181 for( ion=0; ion<nelem+2; ++ion )
00182 {
00183 mole.source[nelem][ion] = 0.;
00184 mole.sink[nelem][ion] = 0.;
00185 }
00186 }
00187 }
00188
00189 if( trace.nTrConvg>=4 )
00190 {
00191
00192 fprintf( ioQQQ,
00193 " CO_drive4 do not evaluate CO chemistry.\n");
00194 }
00195 return;
00196 }
00197
00198 CO_update_species_cache();
00199
00200
00201 CO_update_rks();
00202
00203
00204 loop = 0;
00205 lgNegPop = false;
00206 lgConverged = lgMolecAver("SETUP",chReason);
00207 do
00208 {
00209
00210 if( trace.nTrConvg>=5 )
00211 {
00212
00213 fprintf( ioQQQ,
00214 " CO_drive5 CO chemistry loop %li chReason %s.\n" , loop, chReason );
00215 }
00216
00217
00218 CO_solve(&lgNegPop,&lgZerPop );
00219
00220
00221 lgConverged = lgMolecAver("AVER",chReason);
00222 if( CODEBUG )
00223 fprintf(ioQQQ,"codrivbug %.2f %li Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n",
00224 fnzone ,
00225 loop ,
00226 TorF(lgNegPop) ,
00227 TorF(lgZerPop) ,
00228 findspecies("OH")->hevmol ,
00229 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00230 phycon.te,
00231 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00232
00233 ++loop;
00234 } while (
00235
00236
00237 (loop < LUPMAX_CODRIV) && !lgConverged );
00238
00239
00240 lgMoleZeroed = false;
00241
00242
00243 if( loop == LUPMAX_CODRIV)
00244 {
00245 conv.lgConvIoniz = false;
00246 strcpy( conv.chConvIoniz, "CO con1" );
00247 if( CODEBUG )
00248 fprintf(ioQQQ,"CO_drive not converged\n");
00249 }
00250
00251
00252
00253
00254 lgFirstTry = (nzone==co.co_nzone && iteration==co.iteration_co);
00255
00256
00257 if( lgNegPop )
00258 {
00259 if( conv.lgSearch && (hmi.H2_total/dense.gas_phase[ipHYDROGEN] <h2lim) &&
00260 (iteration==1) )
00261 {
00262
00263
00264
00265
00266
00267
00268 CO_zero();
00269 }
00270 else
00271 {
00272 if( called.lgTalk && !lgFirstTry )
00273 {
00274 conv.lgConvPops = false;
00275 fprintf( ioQQQ,
00276 " PROBLEM CO_drive failed to converge1 due to negative pops, zone=%.2f, CO/C=%.2e OH/C=%.2e H2/H=%.2e\n",
00277 fnzone,
00278 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00279 findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00280 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00281 ConvFail( "pops" , "CO" );
00282 }
00283 }
00284 }
00285
00286
00287 else if( loop == LUPMAX_CODRIV )
00288 {
00289
00290 if( called.lgTalk && !lgFirstTry )
00291 {
00292 conv.lgConvPops = false;
00293 fprintf( ioQQQ,
00294 " PROBLEM CO_drive failed to converge2 in %li iter, zone=%.2f, CO/C=%.2e negpop?%1c reason:%s\n",
00295 loop,
00296 fnzone,
00297 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]),
00298 TorF(lgNegPop) ,
00299 chReason);
00300 ConvFail( "pops" , "CO" );
00301 }
00302 }
00303
00304
00305 ASSERT(conv.lgSearch || findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) <= 1.01 ||
00306
00307 loop == LUPMAX_CODRIV );
00308
00309
00310
00311
00312
00313
00314 if(0) {
00315 double sum_ion , sum_mol;
00316 sum_ion = dense.xIonDense[ipCARBON][0] + dense.xIonDense[ipCARBON][1];
00317 sum_mol = findspecies("C")->hevmol + findspecies("C+")->hevmol;
00318 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
00319 {
00320
00321
00322
00323 conv.lgConvIoniz = false;
00324 strcpy( conv.chConvIoniz, "CO con2" );
00325 conv.BadConvIoniz[0] = sum_ion;
00326 conv.BadConvIoniz[1] = sum_mol;
00327 if( CODEBUG )
00328 fprintf(ioQQQ,"CO_drive not converged\n");
00329 }
00330 sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1];
00331 sum_mol = findspecies("O")->hevmol + findspecies("O+")->hevmol;
00332 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 )
00333 {
00334
00335
00336
00337 conv.lgConvIoniz = false;
00338 strcpy( conv.chConvIoniz, "CO con3" );
00339 conv.BadConvIoniz[0] = sum_ion;
00340 conv.BadConvIoniz[1] = sum_mol;
00341 if( CODEBUG )
00342 fprintf(ioQQQ,"CO_drive not converged\n");
00343 }
00344 }
00345 return;
00346 }
00347
00348 STATIC bool lgMolecAver(
00349
00350
00351 const char chJob[10] ,
00352 char *chReason )
00353 {
00354 long int i;
00355 bool lgReturn;
00356
00357 static realnum hden_save_prev_call;
00358 double conv_fac;
00359 struct chem_element_s *element;
00360
00361 DEBUG_ENTRY( "lgMolecAver()" );
00362 const bool DEBUG_MOLECAVER = false;
00363
00364
00365 strcpy( chReason , "none given" );
00366
00367
00368 if( strcmp( "SETUP" , chJob ) == 0 )
00369 {
00370 static realnum hden_save_prev_iter;
00371
00372
00373
00374
00375
00376
00377 if( iteration!=co.iteration_co &&
00378 hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
00379 {
00380
00381 lgReturn = true;
00382 return lgReturn;
00383 }
00384
00385 else if( co.iteration_co < 0 || lgMoleZeroed )
00386 {
00387
00388
00389
00390
00391
00392
00393
00394 if( !hcmap.lgMapBeingDone || lgMoleZeroed )
00395 {
00396 for( i=0; i < mole.num_comole_calc; i++ )
00397 {
00398 COmole[i]->hevmol = dense.xIonDense[ipCARBON][0]*COmole[i]->pdr_mole_co;
00399 COmole[i]->co_save = COmole[i]->hevmol;
00400 }
00401 }
00402
00403
00404 ASSERT( dense.xIonDense[ipCARBON][0]>SMALLFLOAT );
00405
00406
00407 for(i=0;i<mole.num_elements;i++) {
00408 element = chem_element[i];
00409 COmole[element->ipMl]->hevmol = dense.xIonDense[element->ipCl][0];
00410 COmole[element->ipMlP]->hevmol = dense.xIonDense[element->ipCl][1];
00411 }
00412
00413
00414 co.iteration_co = iteration;
00415 co.co_nzone = nzone;
00416
00417
00418
00419
00420 hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
00421 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
00422
00423 if( DEBUG_MOLECAVER )
00424 fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li zeroing since very first call. H2/H=%.2e\n",
00425 iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00426 }
00427 else if( iteration > co.iteration_co )
00428 {
00429 realnum ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_iter;
00430
00431
00432
00433 for( i=0; i < mole.num_comole_calc; i++ )
00434 {
00435 COmole[i]->hevmol = COmole[i]->hev_reinit*ratio;
00436 COmole[i]->co_save = COmole[i]->hev_reinit*ratio;
00437 }
00438
00439 co.iteration_co = iteration;
00440 co.co_nzone = nzone;
00441
00442 if( DEBUG_MOLECAVER )
00443 fprintf(ioQQQ," DEBUG lgMolecAver iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n",
00444 iteration,
00445 nzone,
00446 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
00447 }
00448 else if( iteration == co.iteration_co && nzone==co.co_nzone+1 )
00449 {
00450
00451
00452 for( i=0; i < mole.num_comole_calc; i++ )
00453 {
00454 COmole[i]->hev_reinit = COmole[i]->hevmol;
00455 }
00456
00457 co.co_nzone = -2;
00458 hden_save_prev_iter = dense.gas_phase[ipHYDROGEN];
00459 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
00460
00461 if( DEBUG_MOLECAVER )
00462 fprintf(ioQQQ,"DEBUG lgMolecAver iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone);
00463 }
00464
00465
00466 lgReturn = true;
00467 }
00468
00469 else if( strcmp( "AVER" , chJob ) == 0 )
00470 {
00471
00472
00473
00474 realnum oldoh = findspecies("OH")->co_save;
00475 lgReturn = true;
00476
00477
00478
00479
00480 for( i=0; i < mole.num_comole_calc; i++ )
00481 {
00482
00483
00484 realnum damper = 0.2f;
00485
00486
00487
00488
00489 COmole[i]->co_save *= dense.gas_phase[ipHYDROGEN] / hden_save_prev_call;
00490
00491
00492
00493 if( COmole[i]->co_save< SMALLFLOAT )
00494 {
00495 COmole[i]->co_save = COmole[i]->hevmol;
00496 }
00497 else
00498 {
00499
00500
00501
00502 double ratio = (double)COmole[i]->hevmol/(double)COmole[i]->co_save;
00503 ASSERT( COmole[i]->co_save>0. && COmole[i]->hevmol>=0. );
00504 if( fabs(ratio-1.) < 1.5 ||
00505 fabs(ratio-1.) > 0.66 )
00506 {
00507
00508 COmole[i]->hevmol = COmole[i]->hevmol*(1.f - damper) +
00509 COmole[i]->co_save*damper;
00510 }
00511 else
00512 {
00513
00514 COmole[i]->hevmol = (realnum)pow(10. ,
00515 (log10(COmole[i]->hevmol) + log10(COmole[i]->co_save))/2. );
00516 }
00517 COmole[i]->co_save = COmole[i]->hevmol;
00518 }
00519 }
00520
00521 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
00522
00523
00524
00525 if( oldoh/SDIV(dense.gas_phase[ipOXYGEN]) < 1e-10 )
00526 conv_fac = 3.;
00527 else
00528 conv_fac = 1.;
00529
00530 if(fabs(oldoh/SDIV(findspecies("OH")->hevmol)-1.) < COTOLER_MOLAV*conv_fac )
00531 {
00532 lgReturn = true;
00533 }
00534 else
00535 {
00536
00537 lgReturn = false;
00538 sprintf( chReason , "OH change from %.3e to %.3e",
00539 oldoh ,
00540 findspecies("OH")->hevmol);
00541 }
00542 }
00543 else
00544 TotalInsanity();
00545
00546
00547 if( findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) > 1.+FLT_EPSILON )
00548 {
00549 lgReturn = false;
00550 sprintf( chReason , "CO/C >1, value is %e",
00551 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
00552 }
00553
00554 if( (trace.lgTrace && trace.lgTr_CO_Mole) || DEBUG_MOLECAVER )
00555 {
00556 fprintf( ioQQQ, " DEBUG lgMolecAver converged? %c" ,TorF(lgReturn) );
00557 fprintf( ioQQQ, " CO/C=%9.1e", findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) );
00558 fprintf( ioQQQ, " OH/O=%9.1e", findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipOXYGEN]) );
00559 fprintf( ioQQQ, "\n" );
00560 }
00561
00562 return lgReturn;
00563 }