cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ConvFail handle convergence failure */ 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 /*ConvFail handle convergence failure - sets lgAbort if too many failures occur */ 00018 void ConvFail( 00019 /* chMode is one of "pres", "eden", "ioni", "pops", "grai", "temp" */ 00020 const char chMode[], /* chMode[5] */ 00021 /* chDetail - string giving details about the convergence failure */ 00022 const char chDetail[] ) 00023 { 00024 double relerror; 00025 00026 DEBUG_ENTRY( "ConvFail()" ); 00027 00028 /* >>chng 02 jun 15, add this branch */ 00029 /* do not announce a convergence failure - this was an abort before 00030 * convergence was achieved */ 00031 if( lgAbort ) 00032 { 00033 /* an abort is not one of the failures we deal with - simply return and 00034 * let something else handle this */ 00035 /*fprintf( ioQQQ, " ConvFail - abort has been set.\n");*/ 00036 return; 00037 } 00038 00039 /* pressure failure */ 00040 if( strcmp( chMode , "pres" )==0 ) 00041 { 00042 /* record number of pressure failures */ 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 /* this identifies new dynamics that failed near the sonic point */ 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 /* electron density failure */ 00069 else if( strcmp( chMode, "eden" ) == 0 ) 00070 { 00071 /* record number of electron density failures */ 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 /* some extra information that may be printed */ 00088 /* heating cooling failure */ 00089 if( !conv.lgConvTemp ) 00090 { 00091 fprintf( ioQQQ, " Temperature failure also." ); 00092 } 00093 00094 /* heating cooling failure */ 00095 if( !conv.lgConvIoniz ) 00096 { 00097 fprintf( ioQQQ, " Ionization failure also." ); 00098 } 00099 00100 /* electron density is oscillating */ 00101 if( conv.lgEdenOscl ) 00102 { 00103 fprintf( ioQQQ, " Electron density oscillating." ); 00104 } 00105 00106 /* heating cooling oscillating */ 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 /* ionization failure */ 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 /* populations failure */ 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 /* ionization failure */ 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 /* rest of routine is temperature failure */ 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 /* not really a temperature failure, but something else */ 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 /* increment total number of failures */ 00214 ++conv.nTotalFailures; 00215 00216 /* now see how many total failures we have, and if it is time to abort */ 00217 /* remember which zone this is */ 00218 conv.ifailz[MIN2(conv.nTotalFailures,10)-1] = nzone; 00219 00220 /* remember the relative error 00221 * convert to single precision for following max, abs (vax failed here) */ 00222 relerror = fabs((thermal.htot-thermal.ctot)/ thermal.htot); 00223 00224 conv.failmx = MAX2(conv.failmx,(realnum)relerror); 00225 00226 /* this branch is non-abort exit - we have not exceeded the limit to the number of failures */ 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 /* check whether went into cold neutral gas without cosmic rays */ 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 /* if due to pressure failures then recommend looking at pressure map */ 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 /* punt */ 00256 if( conv.lgMap ) 00257 { 00258 /* only do map if requested */ 00259 /* adjust range of punting map */ 00260 hcmap.RangeMap[0] = (realnum)(phycon.te/100.); 00261 hcmap.RangeMap[1] = (realnum)MIN2(phycon.te*100.,9e9); 00262 /* need to make printout out now, before disturbing solution with map */ 00263 PrtZone(); 00264 hcmap.lgMapBeingDone = true; 00265 map_do(ioQQQ,"punt"); 00266 } 00267 00268 /* return out from here and let iter_end_check catch lgAbort set, 00269 * and generate normal output there */ 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 }