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