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