cloudy trunk
Loading...
Searching...
No Matches
conv_fail.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*ConvFail handle convergence failure */
4#include "cddefines.h"
5#include "cddrive.h"
6#include "prt.h"
7#include "phycon.h"
8#include "hextra.h"
9#include "pressure.h"
10#include "dense.h"
11#include "thermal.h"
12#include "called.h"
13#include "hcmap.h"
14#include "wind.h"
15#include "conv.h"
16
17/*ConvFail handle convergence failure - sets lgAbort if too many failures occur */
19 /* chMode is one of "pres", "chem", "eden", "ioni", "pops", "grai", "temp" */
20 const char chMode[], /* chMode[5] */
21 /* chDetail - string giving details about the convergence failure */
22 const char chDetail[] )
23{
24 double relerror;
25
26 DEBUG_ENTRY( "ConvFail()" );
27
28 /* >>chng 02 jun 15, add this branch */
29 /* do not announce a convergence failure - this was an abort before
30 * convergence was achieved */
31 if( lgAbort )
32 {
33 /* an abort is not one of the failures we deal with - simply return and
34 * let something else handle this */
35 /*fprintf( ioQQQ, " ConvFail - abort has been set.\n");*/
36 return;
37 }
38
39 /* pressure failure */
40 if( strcmp( chMode , "pres" )==0 )
41 {
42 /* record number of pressure failures */
43 ++conv.nPreFail;
44 if( called.lgTalk )
45 {
46 fprintf( ioQQQ,
47 " PROBLEM ConvFail %li, pressure not converged; itr %li, zone %.2f Te:%.3e Hden:%.4e curr Pres:%.4e Error:%.4e%% Pra/gas:%.3e\n",
48 conv.nPreFail,
50 fnzone,
51 phycon.te,
52 dense.gas_phase[ipHYDROGEN],
53 pressure.PresTotlCurr,
54 pressure.PresTotlError*100.,
55 pressure.pbeta);
56
57 /* this identifies new dynamics that failed near the sonic point */
58 if( fabs(pressure.PresGasCurr - pressure.PresRamCurr)/pressure.PresGasCurr < 0.1 &&
59 strcmp(dense.chDenseLaw,"DYNA") == 0 )
60 {
61 fprintf( ioQQQ,
62 "\n PROBLEM continued, pressure not converged; we are stuck at the sonic point.\n\n");
63 pressure.lgSonicPoint = true;
64 }
65 }
66 }
67
68 /* electron density failure */
69 else if( strcmp( chMode, "eden" ) == 0 )
70 {
71 /* record number of electron density failures */
72 ++conv.nNeFail;
73
74 if( called.lgTalk )
75 {
76 fprintf( ioQQQ,
77 " PROBLEM ConvFail %li, eden not converged itr %li zone %li fnzone %.2f correct=%.3e "
78 "assumed=%.3e.",
79 conv.nNeFail,
80 iteration ,
81 nzone ,
82 fnzone,
83 dense.EdenTrue,
84 dense.eden
85 );
86
87 /* some extra information that may be printed */
88 /* heating cooling failure */
89 if( !conv.lgConvTemp )
90 {
91 fprintf( ioQQQ, " Temperature failure also." );
92 }
93
94 /* heating cooling failure */
95 if( !conv.lgConvIoniz() )
96 {
97 fprintf( ioQQQ, " Ionization failure also." );
98 }
99 }
100 fprintf( ioQQQ, " \n");
101 }
102
103 else if( strcmp( chMode, "ioni" ) == 0 )
104 {
105 /* ionization failure */
106 ++conv.nIonFail;
107 if( called.lgTalk )
108 {
109 fprintf( ioQQQ, " PROBLEM ConvFail %li, %s ionization not converged"
110 " iteration %li zone %li fnzone %.2f reason %s BadConvIoniz0:%g [1]=%g\n",
111 conv.nIonFail,
112 chDetail ,
113 iteration ,
114 nzone,
115 fnzone ,
116 conv.chConvIoniz(),
117 conv.convIonizOldVal(),
118 conv.convIonizNewVal());
119 }
120 }
121
122 else if( strcmp( chMode, "pops" ) == 0 )
123 {
124 /* populations failure */
125 ++conv.nPopFail;
126 conv.lgConvPops = false;
127 if( called.lgTalk )
128 {
129 fprintf( ioQQQ, " PROBLEM ConvFail %li, %s population not converged"
130 " iteration %li zone %li fnzone %.2f %s %g %g\n",
131 conv.nPopFail,
132 chDetail ,
133 iteration,
134 nzone ,
135 fnzone ,
136 conv.chConvIoniz(),
137 conv.convIonizOldVal(),
138 conv.convIonizNewVal());
139 }
140 }
141
142 else if( strcmp( chMode, "grai" ) == 0 )
143 {
144 /* ionization failure */
145 ++conv.nGrainFail;
146 if( called.lgTalk )
147 {
148 fprintf( ioQQQ, " PROBLEM ConvFail %ld, a grain failure occurred"
149 " iteration %li zone %li fnzone %.2f %s %g %g\n",
150 conv.nGrainFail,
151 iteration ,
152 nzone ,
153 fnzone ,
154 conv.chConvIoniz(),
155 conv.convIonizOldVal(),
156 conv.convIonizNewVal());
157 }
158 }
159
160 /* rest of routine is temperature failure */
161 else if( strcmp( chMode, "temp" ) == 0 )
162 {
163 ASSERT( fabs((thermal.htot - thermal.ctot)/thermal.htot ) > conv.HeatCoolRelErrorAllowed );
164 ++conv.nTeFail;
165 if( called.lgTalk )
166 {
167 fprintf( ioQQQ,
168 " PROBLEM ConvFail %ld, Temp not converged itr %li zone %li fnzone %.2f Te=%.4e"
169 " Htot=%.3e Ctot=%.3e rel err=%.3e rel tol:%.3e\n",
170 conv.nTeFail,
171 iteration ,
172 nzone ,
173 fnzone,
174 phycon.te,
175 thermal.htot,
176 thermal.ctot,
177 (thermal.htot - thermal.ctot)/ thermal.htot,
178 conv.HeatCoolRelErrorAllowed );
179
180 /* not really a temperature failure, but something else */
181 if( !conv.lgConvIoniz() )
182 {
183 fprintf( ioQQQ, " Solution not converged due to %10.10s\n",
184 conv.chConvIoniz() );
185 }
186 }
187 }
188 else
189 {
190 fprintf( ioQQQ, " ConvFail called with insane mode %s detail %s\n",
191 chMode ,
192 chDetail );
193 ShowMe();
195 }
196
197 /* increment total number of failures */
198 ++conv.nTotalFailures;
199
200 /* now see how many total failures we have, and if it is time to abort */
201 /* remember which zone this is */
202 conv.ifailz[MIN2(conv.nTotalFailures,10)-1] = nzone;
203
204 /* remember the relative error
205 * convert to single precision for following max, abs (vax failed here) */
206 relerror = fabs((thermal.htot-thermal.ctot)/ thermal.htot);
207
208 conv.failmx = MAX2(conv.failmx,(realnum)relerror);
209
210 /* this branch is non-abort exit - we have not exceeded the limit to the number of failures */
211 if( conv.nTotalFailures < conv.LimFail )
212 {
213 return;
214 }
215
216 fprintf( ioQQQ, " Stop due to excessive convergence failures - there have been %ld so far. \n",
217 conv.LimFail );
218 fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n" );
219
220 /* check whether went into cold neutral gas without cosmic rays */
221 if( phycon.te < 1e3 && dense.eden/dense.gas_phase[ipHYDROGEN] < 0.1 &&
222 (hextra.cryden == 0.) )
223 {
224 fprintf( ioQQQ,"\n This problem may be solved by adding cosmic rays.\n");
225 fprintf( ioQQQ,"\n The gas was cold and neutral.\n");
226 fprintf( ioQQQ,"\n The chemistry is not designed to work without a source of ionization.\n");
227 fprintf( ioQQQ, " >>> Add galactic background cosmic rays with the COSMIC RAYS BACKBOUND command and try again.\n\n" );
228 }
229
230 /* if due to pressure failures then recommend looking at pressure map */
231 if( conv.nPreFail==conv.nTotalFailures )
232 {
233 fprintf( ioQQQ, " These were all pressure failures - we may be near an unstable point in the cooling curve. \n");
234 fprintf( ioQQQ, " The PUNCH PRESSURE HISTORY command will show the n-T-P curve, and may be interesting.\n\n");
235 }
236
237 /* punt */
238 if( conv.lgMap )
239 {
240 /* only do map if requested */
241 /* adjust range of punting map */
242 hcmap.RangeMap[0] = (realnum)(phycon.te/100.);
243 hcmap.RangeMap[1] = (realnum)MIN2(phycon.te*100.,9e9);
244 /* need to make printout out now, before disturbing solution with map */
245 PrtZone();
246 hcmap.lgMapBeingDone = true;
247 map_do(ioQQQ,"punt");
248 }
249
250 /* return out from here and let iter_end_check catch lgAbort set,
251 * and generate normal output there */
252 lgAbort = true;
253 if( called.lgTalk )
254 {
255 fprintf( ioQQQ, " ConvFail sets lgAbort since nTotalFailures=%ld is >= LimFail=%ld\n",
256 conv.nTotalFailures,
257 conv.LimFail );
258 fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n");
259 fflush( ioQQQ );
260 }
261 return;
262}
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define MIN2
Definition cddefines.h:761
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
t_conv conv
Definition conv.cpp:5
void ConvFail(const char chMode[], const char chDetail[])
Definition conv_fail.cpp:18
t_dense dense
Definition dense.cpp:24
t_hcmap hcmap
Definition hcmap.cpp:21
void map_do(FILE *io, const char *chType)
Definition hcmap.cpp:23
t_hextra hextra
Definition hextra.cpp:5
t_phycon phycon
Definition phycon.cpp:6
t_pressure pressure
Definition pressure.cpp:5
void PrtZone(void)
Definition prt_zone.cpp:36
t_thermal thermal
Definition thermal.cpp:5