cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
conv_itercheck.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 /*ConvIterCheck check whether model has converged or whether more iterations
4  * are needed - implements the iter to converg comnd */
5 #include "cddefines.h"
6 #include "taulines.h"
7 #include "iso.h"
8 #include "phycon.h"
9 #include "mole.h"
10 #include "elementnames.h"
11 #include "dynamics.h"
12 #include "stopcalc.h"
13 #include "dense.h"
14 #include "iterations.h"
15 #include "colden.h"
16 #include "punch.h"
17 #include "rt.h"
18 #include "conv.h"
19 
20 /*ConvIterCheck check whether model has converged or whether more iterations
21  * are needed - implements the iterate to convergence command */
22 void ConvIterCheck( void )
23 {
24  bool lgConverged;
25  long int nelem,
26  i,
27  ipISO,
28  ipHi, ipLo;
29 
30  DEBUG_ENTRY( "ConvIterCheck()" );
31 
32  /* =======================================================================*/
33  /* this is an option to keep iterating until it converges
34  * iterate to convergence option
35  * autocv is percentage difference in optical depths allowed,
36  * =0.20 in block data
37  * checking on Ly and Balmer lines */
38  /*>>chng 04 oct 19, promote loop to do all iso-electronic series */
39  lgConverged = true;
40  strcpy( conv.chNotConverged, "Converged!" );
41  if( iteration > 1 && conv.lgAutoIt )
42  {
43  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
44  {
45  for( nelem=ipISO; nelem < LIMELM; nelem++ )
46  {
47  if( dense.lgElmtOn[nelem] )
48  {
49  /* now check if major subordinate line is converged - for H-like this will
50  * be Ha, and for He-like, the 23P - 23S transition - this will not work for
51  * NISO > 2 so must check against this */
52  if(ipISO==ipH_LIKE )
53  {
54  ipHi = ipH3p;
55  ipLo = ipH2s;
56  }
57  else if( ipISO==ipHE_LIKE )
58  {
59  ipHi = ipHe2p3P2;
60  ipLo = ipHe2s3S;
61  }
62  else
63  /* fails when NISO increased, add more sequences */
64  TotalInsanity();
65 
66  /* check both H-alpha and Ly-alpha for all species -
67  * only if Balmer lines thick
68  * so check if Ha optical depth significant */
69  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn > 0.5 )
70  {
71  /* test if Lya converged, nLyaLevel is upper level of Lya for iso seq */
72  if( fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot -
73  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) >
74  conv.autocv*fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) )
75  {
76  /* not converged to within AUTOCV, normally 15 percent */
77  lgConverged = false;
78 
79  /* for iterate to convergence, print reason why it was not converged
80  * on 3rd and higher iterations */
81  sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] );
82 
83  if( punch.lgPunConv )
84  {
85  fprintf( punch.ipPunConv, " %s-like Lya thick, "
86  "nelem= %s iteration %li old %.3e new %.3e \n",
87  elementnames.chElementSym[ipISO] ,
88  elementnames.chElementSym[nelem],
89  iteration,
90  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot ,
91  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
92  }
93  }
94 
95  if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot -
96  Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) >
97  conv.autocv*fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) )
98  {
99  /* not converged to within AUTOCV, normally 15 percent */
100  lgConverged = false;
101 
102  /* for iterate to convergence, print reason why it was not converged
103  * on 3rd and higher iterations */
104  sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] );
105 
106  if( punch.lgPunConv )
107  {
108  fprintf( punch.ipPunConv, " %s-like subord, nelem= %s iteration %li old %.3e new %.3e \n" ,
109  elementnames.chElementSym[ipISO],
110  elementnames.chElementSym[nelem],
111  iteration,
112  Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot ,
113  Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn );
114  }
115  }
116  }
117  }
118  }
119  }
120 
121  /* >>chng 03 sep 07, add this test */
122  /* check on changes in major column densities */
123  for( i=0; i<NCOLD; ++i )
124  {
125  /* was the species column density significant relative to
126  * the total H column density, and was its abundance changing? */
127  if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 &&
129  {
130  /* not converged to within conv.autocv, normally 20 percent */
131  lgConverged = false;
132 
133  /* for iterate to convergence, print reason why it was not converged
134  * on 3rd and higher iterations */
135  strcpy( conv.chNotConverged, "H mole col" );
136 
137  if( punch.lgPunConv )
138  {
139  fprintf( punch.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
140  i,iteration,
141  colden.colden_old[i],
142  colden.colden[i],
144  }
145  }
146  }
147 
148  /* >>chng 03 sep 07, add this test */
149  /* check on changes in major column densities */
150  for( i=0; i<mole.num_comole_calc; ++i )
151  {
152  if(COmole[i]->n_nuclei == 1)
153  continue;
154 
155  /* was the species abundance and changing? */
156  if( COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5 &&
157  fabs(COmole[i]->hevcol_old-COmole[i]->hevcol) > conv.autocv*COmole[i]->hevcol )
158  {
159  /* not converged to within conv.autocv, normally 20 percent */
160  lgConverged = false;
161 
162  /* for iterate to convergence, print reason why it was not converged
163  * on 3rd and higher iterations */
164  strcpy( conv.chNotConverged, "CO mol col" );
165  /*fprintf(ioQQQ,"debugggreset\t CO mole %li %li %.2e %.2e\n",
166  i,iteration,COmole[i]->hevcol_old,COmole[i]->hevcol);*/
167 
168  if( punch.lgPunConv )
169  {
170  fprintf( punch.ipPunConv, "CO mol col, old:%.3e new:%.3e\n" ,
171  COmole[i]->hevcol_old ,
172  COmole[i]->hevcol );
173  }
174  }
175  }
176 
177  /* check on dynamical convergence in wind model with negative velocity */
178  if( dynamics.lgAdvection )
179  {
180  /* >>chng 02 nov 29, as per Will Henney email */
183  {
184  lgConverged = false;
185  /* for iterate to convergence, print reason why it was not converged
186  * on 3rd and higher iterations */
187  strcpy( conv.chNotConverged, "Dynamics " );
188  if( punch.lgPunConv )
189  {
190  fprintf( punch.ipPunConv, " Dynamics\n" );
191  }
192  }
193  }
194 
195  if( punch.lgPunConv && lgConverged )
196  {
197  fprintf( punch.ipPunConv, " converged\n" );
198  }
199 
200  /* lower limit to number of iterations if converged */
201  if( lgConverged )
203 
204  /* >>chng 96 dec 20, moved following to within if on lgAutoIt
205  * this is test for stopping on first zone */
206  if( phycon.te < StopCalc.tend && nzone == 1 )
207  {
208  lgConverged = true;
209  strcpy( conv.chNotConverged, " " );
211  }
212  }
213  return;
214 }

Generated for cloudy by doxygen 1.8.1.2