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 Error:%.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.PresTotlError*100.,
00055 pressure.pbeta);
00056
00057
00058 if( fabs(pressure.PresGasCurr - pressure.PresRamCurr)/pressure.PresGasCurr < 0.1 &&
00059 strcmp(dense.chDenseLaw,"DYNA") == 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 fprintf( ioQQQ, " \n");
00101 }
00102
00103 else if( strcmp( chMode, "ioni" ) == 0 )
00104 {
00105
00106 ++conv.nIonFail;
00107 if( called.lgTalk )
00108 {
00109 fprintf( ioQQQ, " PROBLEM ConvFail %li, %s ionization not converged"
00110 " iteration %li zone %li fnzone %.2f reason %s BadConvIoniz0:%g [1]=%g\n",
00111 conv.nIonFail,
00112 chDetail ,
00113 iteration ,
00114 nzone,
00115 fnzone ,
00116 conv.chConvIoniz(),
00117 conv.convIonizOldVal(),
00118 conv.convIonizNewVal());
00119 }
00120 }
00121
00122 else if( strcmp( chMode, "pops" ) == 0 )
00123 {
00124
00125 ++conv.nPopFail;
00126 conv.lgConvPops = false;
00127 if( called.lgTalk )
00128 {
00129 fprintf( ioQQQ, " PROBLEM ConvFail %li, %s population not converged"
00130 " iteration %li zone %li fnzone %.2f %s %g %g\n",
00131 conv.nPopFail,
00132 chDetail ,
00133 iteration,
00134 nzone ,
00135 fnzone ,
00136 conv.chConvIoniz(),
00137 conv.convIonizOldVal(),
00138 conv.convIonizNewVal());
00139 }
00140 }
00141
00142 else if( strcmp( chMode, "grai" ) == 0 )
00143 {
00144
00145 ++conv.nGrainFail;
00146 if( called.lgTalk )
00147 {
00148 fprintf( ioQQQ, " PROBLEM ConvFail %ld, a grain failure occurred"
00149 " iteration %li zone %li fnzone %.2f %s %g %g\n",
00150 conv.nGrainFail,
00151 iteration ,
00152 nzone ,
00153 fnzone ,
00154 conv.chConvIoniz(),
00155 conv.convIonizOldVal(),
00156 conv.convIonizNewVal());
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"
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.htot,
00176 thermal.ctot,
00177 (thermal.htot - thermal.ctot)/ thermal.htot,
00178 conv.HeatCoolRelErrorAllowed );
00179
00180
00181 if( !conv.lgConvIoniz() )
00182 {
00183 fprintf( ioQQQ, " Solution not converged due to %10.10s\n",
00184 conv.chConvIoniz() );
00185 }
00186 }
00187 }
00188 else
00189 {
00190 fprintf( ioQQQ, " ConvFail called with insane mode %s detail %s\n",
00191 chMode ,
00192 chDetail );
00193 ShowMe();
00194 cdEXIT(EXIT_FAILURE);
00195 }
00196
00197
00198 ++conv.nTotalFailures;
00199
00200
00201
00202 conv.ifailz[MIN2(conv.nTotalFailures,10)-1] = nzone;
00203
00204
00205
00206 relerror = fabs((thermal.htot-thermal.ctot)/ thermal.htot);
00207
00208 conv.failmx = MAX2(conv.failmx,(realnum)relerror);
00209
00210
00211 if( conv.nTotalFailures < conv.LimFail )
00212 {
00213 return;
00214 }
00215
00216 fprintf( ioQQQ, " Stop due to excessive convergence failures - there have been %ld so far. \n",
00217 conv.LimFail );
00218 fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n" );
00219
00220
00221 if( phycon.te < 1e3 && dense.eden/dense.gas_phase[ipHYDROGEN] < 0.1 &&
00222 (hextra.cryden == 0.) )
00223 {
00224 fprintf( ioQQQ,"\n This problem may be solved by adding cosmic rays.\n");
00225 fprintf( ioQQQ,"\n The gas was cold and neutral.\n");
00226 fprintf( ioQQQ,"\n The chemistry is not designed to work without a source of ionization.\n");
00227 fprintf( ioQQQ, " >>> Add galactic background cosmic rays with the COSMIC RAYS BACKBOUND command and try again.\n\n" );
00228 }
00229
00230
00231 if( conv.nPreFail==conv.nTotalFailures )
00232 {
00233 fprintf( ioQQQ, " These were all pressure failures - we may be near an unstable point in the cooling curve. \n");
00234 fprintf( ioQQQ, " The PUNCH PRESSURE HISTORY command will show the n-T-P curve, and may be interesting.\n\n");
00235 }
00236
00237
00238 if( conv.lgMap )
00239 {
00240
00241
00242 hcmap.RangeMap[0] = (realnum)(phycon.te/100.);
00243 hcmap.RangeMap[1] = (realnum)MIN2(phycon.te*100.,9e9);
00244
00245 PrtZone();
00246 hcmap.lgMapBeingDone = true;
00247 map_do(ioQQQ,"punt");
00248 }
00249
00250
00251
00252 lgAbort = true;
00253 if( called.lgTalk )
00254 {
00255 fprintf( ioQQQ, " ConvFail sets lgAbort since nTotalFailures=%ld is >= LimFail=%ld\n",
00256 conv.nTotalFailures,
00257 conv.LimFail );
00258 fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n");
00259 fflush( ioQQQ );
00260 }
00261 return;
00262 }