00001
00002
00003
00004 #include "cddefines.h"
00005 #include "cddrive.h"
00006 #include "prt.h"
00007 #include "phycon.h"
00008 #include "hextra.h"
00009 #include "pressure.h"
00010 #include "dense.h"
00011 #include "thermal.h"
00012 #include "called.h"
00013 #include "map.h"
00014 #include "wind.h"
00015 #include "conv.h"
00016
00017
00018 void ConvFail(
00019
00020 const char chMode[],
00021
00022 const char chDetail[] )
00023 {
00024 double relerror;
00025
00026 DEBUG_ENTRY( "ConvFail()" );
00027
00028
00029
00030
00031 if( lgAbort )
00032 {
00033
00034
00035
00036 DEBUG_EXIT( "ConvFail()" );
00037 return;
00038 }
00039
00040
00041 if( strcmp( chMode , "pres" )==0 )
00042 {
00043
00044 ++conv.nPreFail;
00045 if( called.lgTalk )
00046 {
00047 fprintf( ioQQQ,
00048 " PROBLEM ConvFail %li, pressure not converged; itr %li, zone %.2f Te:%.3e Hden:%.4e curr Pres:%.4e Corr Pres:%.4e Pra/gas:%.3e\n",
00049 conv.nPreFail,
00050 iteration,
00051 fnzone,
00052 phycon.te,
00053 dense.gas_phase[ipHYDROGEN],
00054 pressure.PresTotlCurr,
00055 pressure.PresTotlCorrect,
00056 pressure.pbeta);
00057
00058
00059 if( fabs(pressure.PresGasCurr - pressure.PresRamCurr)/pressure.PresGasCurr < 0.1 &&
00060 ((strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. ) )
00061 {
00062 fprintf( ioQQQ,
00063 "\n PROBLEM continued, pressure not converged; we are stuck at the sonic point.\n\n");
00064 pressure.lgSonicPoint = true;
00065 }
00066 }
00067 }
00068
00069
00070 else if( strcmp( chMode, "eden" ) == 0 )
00071 {
00072
00073 ++conv.nNeFail;
00074
00075 if( called.lgTalk )
00076 {
00077 fprintf( ioQQQ,
00078 " PROBLEM ConvFail %li, eden not converged itr %li zone %li fnzone %.2f correct=%.3e "
00079 "assumed=%.3e.",
00080 conv.nNeFail,
00081 iteration ,
00082 nzone ,
00083 fnzone,
00084 dense.EdenTrue,
00085 dense.eden
00086 );
00087
00088
00089
00090 if( !conv.lgConvTemp )
00091 {
00092 fprintf( ioQQQ, " Temperature failure also." );
00093 }
00094
00095
00096 if( !conv.lgConvIoniz )
00097 {
00098 fprintf( ioQQQ, " Ionization failure also." );
00099 }
00100
00101
00102 if( conv.lgEdenOscl )
00103 {
00104 fprintf( ioQQQ, " Electron density oscillating." );
00105 }
00106
00107
00108 if( conv.lgCmHOsc )
00109 {
00110 fprintf( ioQQQ, " Heating-cooling oscillating." );
00111 }
00112 }
00113 fprintf( ioQQQ, " \n");
00114 }
00115
00116 else if( strcmp( chMode, "ioni" ) == 0 )
00117 {
00118
00119 ++conv.nIonFail;
00120 if( called.lgTalk )
00121 {
00122 fprintf( ioQQQ, " PROBLEM ConvFail %li, %s ionization not converged iteration %li zone %li fnzone %.2f \n",
00123 conv.nIonFail,
00124 chDetail ,
00125 iteration ,
00126 nzone,
00127 fnzone );
00128 }
00129 }
00130
00131 else if( strcmp( chMode, "pops" ) == 0 )
00132 {
00133
00134 ++conv.nPopFail;
00135 conv.lgConvPops = false;
00136 if( called.lgTalk )
00137 {
00138 fprintf( ioQQQ, " PROBLEM ConvFail %li, %s population not converged iteration %li zone %li fnzone %.2f \n",
00139 conv.nPopFail,
00140 chDetail ,
00141 iteration,
00142 nzone ,
00143 fnzone );
00144 }
00145 }
00146
00147 else if( strcmp( chMode, "grai" ) == 0 )
00148 {
00149
00150 ++conv.nGrainFail;
00151 if( called.lgTalk )
00152 {
00153 fprintf( ioQQQ, " PROBLEM ConvFail %ld, a grain failure occurred iteration %li zone %li fnzone %.2f \n",
00154 conv.nGrainFail,
00155 iteration ,
00156 nzone ,
00157 fnzone );
00158 }
00159 }
00160
00161
00162 else if( strcmp( chMode, "temp" ) == 0 )
00163 {
00164 ASSERT( fabs((thermal.htot - thermal.ctot)/thermal.htot ) > conv.HeatCoolRelErrorAllowed );
00165 ++conv.nTeFail;
00166 if( called.lgTalk )
00167 {
00168 fprintf( ioQQQ,
00169 " PROBLEM ConvFail %ld, Temp not converged itr %li zone %li fnzone %.2f Te=%.4e DTe=%.2e"
00170 " Htot=%.3e Ctot=%.3e rel err=%.3e rel tol:%.3e\n",
00171 conv.nTeFail,
00172 iteration ,
00173 nzone ,
00174 fnzone,
00175 phycon.te,
00176 thermal.dTemper ,
00177 thermal.htot,
00178 thermal.ctot,
00179 (thermal.htot - thermal.ctot)/ thermal.htot,
00180 conv.HeatCoolRelErrorAllowed );
00181
00182 if( conv.lgTOscl && conv.lgCmHOsc )
00183 {
00184 fprintf( ioQQQ, " Temperature and d(Cool-Heat)/dT were both oscillating.\n" );
00185 }
00186
00187 else if( conv.lgTOscl )
00188 {
00189 fprintf( ioQQQ, " Temperature was oscillating.\n" );
00190 }
00191
00192 else if( conv.lgCmHOsc )
00193 {
00194 fprintf( ioQQQ, " d(Cool-Heat)/dT was oscillating.\n" );
00195 }
00196
00197
00198 if( !conv.lgConvIoniz )
00199 {
00200 fprintf( ioQQQ, " Solution not converged due to %10.10s\n",
00201 conv.chConvIoniz );
00202 }
00203 }
00204 }
00205 else
00206 {
00207 fprintf( ioQQQ, " ConvFail called with insane mode %s detail %s\n",
00208 chMode ,
00209 chDetail );
00210 ShowMe();
00211 puts( "[Stop in ConvFail]" );
00212 cdEXIT(EXIT_FAILURE);
00213 }
00214
00215
00216 ++conv.nTotalFailures;
00217
00218
00219
00220 conv.ifailz[MIN2(conv.nTotalFailures,10)-1] = nzone;
00221
00222
00223
00224 relerror = fabs((thermal.htot-thermal.ctot)/ thermal.htot);
00225
00226 conv.failmx = MAX2(conv.failmx,(float)relerror);
00227
00228
00229
00230
00231 if( conv.nTotalFailures < conv.LimFail )
00232 {
00233 DEBUG_EXIT( "ConvFail()" );
00234 return;
00235 }
00236
00237 fprintf( ioQQQ, " Stop due to excessive convergence failures - there have been %ld so far. \n",
00238 conv.LimFail );
00239 fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n" );
00240
00241
00242 if( phycon.te < 1e3 && dense.eden/dense.gas_phase[ipHYDROGEN] < 0.1 &&
00243 (hextra.cryden == 0.) )
00244 {
00245 fprintf( ioQQQ,"\n This problem may be solved by adding cosmic rays.\n");
00246 fprintf( ioQQQ,"\n The gas was cold and neutral.\n");
00247 fprintf( ioQQQ,"\n The chemistry is not designed to work without a source of ionization.\n");
00248 fprintf( ioQQQ, " >>> Add galactic background cosmic rays with the COSMIC RAYS BACKBOUND command and try again.\n\n" );
00249 }
00250
00251 fprintf( ioQQQ, "\n The input commands that lead to this failure are the following:\n" );
00252
00253 cdPrintCommands(ioQQQ);
00254 fprintf( ioQQQ, "\n" );
00255
00256
00257 if( conv.nPreFail==conv.nTotalFailures )
00258 {
00259 fprintf( ioQQQ, " These were all pressure failures - we may be near an unstable point in the cooling curve. \n");
00260 fprintf( ioQQQ, " The PUNCH PRESSURE HISTORY command will show the n-T-P curve, and may be interesting.\n\n");
00261 }
00262
00263
00264 if( conv.lgMap )
00265 {
00266
00267
00268 map.RangeMap[0] = (float)(phycon.te/100.);
00269 map.RangeMap[1] = (float)MIN2(phycon.te*100.,9e9);
00270
00271 PrtZone();
00272 map.lgMapBeingDone = true;
00273 map_do(ioQQQ,"punt");
00274 }
00275
00276
00277
00278 lgAbort = true;
00279 if( called.lgTalk )
00280 {
00281 fprintf( ioQQQ, " ConvFail sets lgAbort since nTotalFailures=%ld is >= LimFail=%ld\n",
00282 conv.nTotalFailures,
00283 conv.LimFail );
00284 fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n");
00285 fflush( ioQQQ );
00286 }
00287
00288 DEBUG_EXIT( "ConvFail()" );
00289
00290 return;
00291 }