cloudy  trunk
 All Data Structures 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-2008 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 */
18 void ConvFail(
19  /* chMode is one of "pres", "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 Corr Pres:%.4e Pra/gas:%.3e\n",
48  conv.nPreFail,
49  iteration,
50  fnzone,
51  phycon.te,
55  pressure.pbeta);
56 
57  /* this identifies new dynamics that failed near the sonic point */
59  ((strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 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  /* electron density is oscillating */
101  if( conv.lgEdenOscl )
102  {
103  fprintf( ioQQQ, " Electron density oscillating." );
104  }
105 
106  /* heating cooling oscillating */
107  if( conv.lgCmHOsc )
108  {
109  fprintf( ioQQQ, " Heating-cooling oscillating." );
110  }
111  }
112  fprintf( ioQQQ, " \n");
113  }
114 
115  else if( strcmp( chMode, "ioni" ) == 0 )
116  {
117  /* ionization failure */
118  ++conv.nIonFail;
119  if( called.lgTalk )
120  {
121  fprintf( ioQQQ, " PROBLEM ConvFail %li, %s ionization not converged iteration %li zone %li fnzone %.2f \n",
122  conv.nIonFail,
123  chDetail ,
124  iteration ,
125  nzone,
126  fnzone );
127  }
128  }
129 
130  else if( strcmp( chMode, "pops" ) == 0 )
131  {
132  /* populations failure */
133  ++conv.nPopFail;
134  conv.lgConvPops = false;
135  if( called.lgTalk )
136  {
137  fprintf( ioQQQ, " PROBLEM ConvFail %li, %s population not converged iteration %li zone %li fnzone %.2f \n",
138  conv.nPopFail,
139  chDetail ,
140  iteration,
141  nzone ,
142  fnzone );
143  }
144  }
145 
146  else if( strcmp( chMode, "grai" ) == 0 )
147  {
148  /* ionization failure */
149  ++conv.nGrainFail;
150  if( called.lgTalk )
151  {
152  fprintf( ioQQQ, " PROBLEM ConvFail %ld, a grain failure occurred iteration %li zone %li fnzone %.2f \n",
153  conv.nGrainFail,
154  iteration ,
155  nzone ,
156  fnzone );
157  }
158  }
159 
160  /* rest of routine is temperature failure */
161  else if( strcmp( chMode, "temp" ) == 0 )
162  {
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 DTe=%.2e"
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.dTemper ,
176  thermal.htot,
177  thermal.ctot,
180 
181  if( conv.lgTOscl && conv.lgCmHOsc )
182  {
183  fprintf( ioQQQ, " Temperature and d(Cool-Heat)/dT were both oscillating.\n" );
184  }
185 
186  else if( conv.lgTOscl )
187  {
188  fprintf( ioQQQ, " Temperature was oscillating.\n" );
189  }
190 
191  else if( conv.lgCmHOsc )
192  {
193  fprintf( ioQQQ, " d(Cool-Heat)/dT was oscillating.\n" );
194  }
195 
196  /* not really a temperature failure, but something else */
197  if( !conv.lgConvIoniz )
198  {
199  fprintf( ioQQQ, " Solution not converged due to %10.10s\n",
200  conv.chConvIoniz );
201  }
202  }
203  }
204  else
205  {
206  fprintf( ioQQQ, " ConvFail called with insane mode %s detail %s\n",
207  chMode ,
208  chDetail );
209  ShowMe();
210  cdEXIT(EXIT_FAILURE);
211  }
212 
213  /* increment total number of failures */
215 
216  /* now see how many total failures we have, and if it is time to abort */
217  /* remember which zone this is */
219 
220  /* remember the relative error
221  * convert to single precision for following max, abs (vax failed here) */
222  relerror = fabs((thermal.htot-thermal.ctot)/ thermal.htot);
223 
224  conv.failmx = MAX2(conv.failmx,(realnum)relerror);
225 
226  /* this branch is non-abort exit - we have not exceeded the limit to the number of failures */
228  {
229  return;
230  }
231 
232  fprintf( ioQQQ, " Stop due to excessive convergence failures - there have been %ld so far. \n",
233  conv.LimFail );
234  fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n" );
235 
236  /* check whether went into cold neutral gas without cosmic rays */
237  if( phycon.te < 1e3 && dense.eden/dense.gas_phase[ipHYDROGEN] < 0.1 &&
238  (hextra.cryden == 0.) )
239  {
240  fprintf( ioQQQ,"\n This problem may be solved by adding cosmic rays.\n");
241  fprintf( ioQQQ,"\n The gas was cold and neutral.\n");
242  fprintf( ioQQQ,"\n The chemistry is not designed to work without a source of ionization.\n");
243  fprintf( ioQQQ, " >>> Add galactic background cosmic rays with the COSMIC RAYS BACKBOUND command and try again.\n\n" );
244  }
245 
246  /* if due to pressure failures then recommend looking at pressure map */
248  {
249  fprintf( ioQQQ, " These were all pressure failures - we may be near an unstable point in the cooling curve. \n");
250  fprintf( ioQQQ, " The PUNCH PRESSURE HISTORY command will show the n-T-P curve, and may be interesting.\n\n");
251  }
252 
253  ShowMe();
254 
255  /* punt */
256  if( conv.lgMap )
257  {
258  /* only do map if requested */
259  /* adjust range of punting map */
260  hcmap.RangeMap[0] = (realnum)(phycon.te/100.);
261  hcmap.RangeMap[1] = (realnum)MIN2(phycon.te*100.,9e9);
262  /* need to make printout out now, before disturbing solution with map */
263  PrtZone();
264  hcmap.lgMapBeingDone = true;
265  map_do(ioQQQ,"punt");
266  }
267 
268  /* return out from here and let iter_end_check catch lgAbort set,
269  * and generate normal output there */
270  lgAbort = true;
271  if( called.lgTalk )
272  {
273  fprintf( ioQQQ, " ConvFail sets lgAbort since nTotalFailures=%ld is >= LimFail=%ld\n",
275  conv.LimFail );
276  fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n");
277  fflush( ioQQQ );
278  }
279  return;
280 }

Generated for cloudy by doxygen 1.8.3.1