cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "cddrive.h"
10 #include "mole.h"
11 #include "elementnames.h"
12 #include "dynamics.h"
13 #include "stopcalc.h"
14 #include "dense.h"
15 #include "iterations.h"
16 #include "colden.h"
17 #include "save.h"
18 #include "rt.h"
19 #include "conv.h"
20 
21 /*ConvIterCheck check whether model has converged or whether more iterations
22  * are needed - implements the iterate to convergence command */
23 void ConvIterCheck( void )
24 {
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 */
40  strcpy( conv.chNotConverged, "Converged!" );
41 
42  // set up intensities used to converge outward intensity of Hb
43  static double HbFracOutOld=-1. , HbFracOutNew=-1.;
44  HbFracOutOld = HbFracOutNew;
45 
46  double a, total, BeamedIn;
47  long int ipTotal = cdLine( "H 1" , 4861.33f , &a , &total );
48  long int ipInwd = cdLine( "INWD" , 4861.33f , &a , &BeamedIn );
49 
50  /* 2014 aug 23, mchatzikos
51  * when cdLine returned -37 for log of zero intensity, the ratio below evaluated to 0;
52  * preserve this behaviour now that cdLine returns linear intensities */
53  HbFracOutNew = 0.;
54  if( total > 0. )
55  HbFracOutNew = 1. - BeamedIn / total;
56 
57  ASSERT( iteration == 1 || (HbFracOutNew>=0 && HbFracOutNew<=1.) );
58  // this disables the test on the outward Hb
59  ipInwd = -1;
60 
61  if( save.lgPunConv )
62  {
63  fprintf( save.ipPunConv, " iteration %li of %li\n" ,
65  }
66 
67  bool lgReasonGiven = false;
68  if( iteration > 1 && conv.lgAutoIt )
69  {
70  if( nzone>3 && ipInwd>=0 && ipTotal>=0 )
71  {
72  // check whether outward intensity of Hb has converged
73  if( fabs(HbFracOutNew-HbFracOutOld)/HbFracOutNew> conv.autocv )
74  {
76  sprintf( conv.chNotConverged, "change in outward Hb");
77  if( save.lgPunConv )
78  {
79  lgReasonGiven = true;
80  fprintf( save.ipPunConv, " Change in outward Hb, "
81  "old=%.3e new=%.3e \n",
82  HbFracOutOld , HbFracOutNew);
83  }
84  }
85  }
86  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
87  {
88  for( nelem=ipISO; nelem < LIMELM; nelem++ )
89  {
90  if( dense.lgElmtOn[nelem] )
91  {
92  /* now check if major subordinate line is converged - for H-like this will
93  * be Ha, and for He-like, the 23P - 23S transition - this will not work for
94  * NISO > 2 so must check against this */
95  if(ipISO==ipH_LIKE )
96  {
97  ipHi = ipH3p;
98  ipLo = ipH2s;
99  }
100  else if( ipISO==ipHE_LIKE )
101  {
102  ipHi = ipHe2p3P2;
103  ipLo = ipHe2s3S;
104  }
105  else
106  /* fails when NISO increased, add more sequences */
107  TotalInsanity();
108 
109  /* check both H-alpha and Ly-alpha for all species -
110  * only if Balmer lines thick
111  * so check if Ha optical depth significant */
112  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() > 0.5 )
113  {
114  /* test if Lya converged, nLyaLevel is upper level of Lya for iso seq */
115  if( fabs(iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() -
116  iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn()*rt.DoubleTau) >
117  conv.autocv*fabs(iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn()*rt.DoubleTau) )
118  {
119  /* not converged to within AUTOCV, normally 15 percent */
121 
122  /* for iterate to convergence, print reason why it was not converged
123  * on 3rd and higher iterations */
124  sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] );
125 
126  if( save.lgPunConv )
127  {
128  lgReasonGiven = true;
129  fprintf( save.ipPunConv, " %s-like Lya thick, "
130  "nelem= %s iteration %li old %.3e new %.3e \n",
131  elementnames.chElementSym[ipISO] ,
132  elementnames.chElementSym[nelem],
133  iteration,
134  iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
135  iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
136  }
137  }
138 
139  if( fabs(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
140  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()*rt.DoubleTau) >
141  conv.autocv*fabs(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()*rt.DoubleTau) )
142  {
143  /* not converged to within AUTOCV, normally 15 percent */
145 
146  /* for iterate to convergence, print reason why it was not converged
147  * on 3rd and higher iterations */
148  sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] );
149 
150  if( save.lgPunConv )
151  {
152  lgReasonGiven = true;
153  fprintf( save.ipPunConv, " %s-like subord, nelem= %s iteration %li old %.3e new %.3e \n" ,
154  elementnames.chElementSym[ipISO],
155  elementnames.chElementSym[nelem],
156  iteration,
157  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() ,
158  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() );
159  }
160  }
161  }
162  }
163  }
164  }
165 
166  if(0)
167  {
168  // all database lines
169  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
170  {
171  if( dBaseSpecies[ipSpecies].lgActive )
172  {
173  for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
174  tr != dBaseTrans[ipSpecies].end(); ++tr)
175  {
176  if( (*tr).Emis().TauIn() > 1. &&
177  fabs((*tr).Emis().TauTot() - (*tr).Emis().TauIn()*rt.DoubleTau) >
178  conv.autocv*fabs((*tr).Emis().TauIn()*rt.DoubleTau) )
179  {
180  /* not converged to within AUTOCV, normally 15 percent */
182 
183  /* for iterate to convergence, print reason why it was not converged
184  * on 3rd and higher iterations */
185  sprintf( conv.chNotConverged, "database %s %.1f",
186  dBaseSpecies[ipSpecies].chLabel ,
187  (*tr).WLAng() );
188 
189  if( save.lgPunConv )
190  {
191  lgReasonGiven = true;
192  fprintf( save.ipPunConv, " database %s %.1f iteration %li old %.3e new %.3e \n",
193  dBaseSpecies[ipSpecies].chLabel ,
194  (*tr).WLAng(),
195  iteration,
196  (*tr).Emis().TauTot() ,
197  (*tr).Emis().TauIn());
198  }
199  }
200  }
201  }
202  }
203  }
204 
205  if( conv.lgAllTransitions )
206  {
207  realnum baddiff=0., badscale=0.;
208  bool lgTauConv=true;
210 
211  // all lines
212  for( vector<TransitionList>::iterator it = AllTransitions.begin(); it != AllTransitions.end(); ++it )
213  {
214  for (TransitionList::iterator tr=it->begin();
215  tr != it->end(); ++tr)
216  {
217  if ( tr->ipCont() <= 0 )
218  continue;
219  realnum diff = fabs((*tr).Emis().TauTot() - (*tr).Emis().TauIn()*rt.DoubleTau);
220  realnum scale = fabs((*tr).Emis().TauIn()*rt.DoubleTau);
221  if( (*tr).Emis().TauIn() > 1. && diff > conv.autocv*scale )
222  {
223  /* not converged to within AUTOCV, normally 15 percent */
225  if ( lgTauConv || diff*badscale > scale*baddiff )
226  {
227  lgTauConv = false;
228  badscale = scale;
229  baddiff = diff;
230  badtr = tr;
231  }
232  if( save.lgPunConv )
233  {
234  lgReasonGiven = true;
235  fprintf( save.ipPunConv, " database %s line %s iteration %li old %.3e new %.3e \n",
236  (*it).chLabel().c_str(),
237  chLineLbl(*tr).c_str() ,
238  iteration,
239  (*tr).Emis().TauTot() ,
240  (*tr).Emis().TauIn()*rt.DoubleTau);
241  }
242  }
243  }
244  }
245  if (!lgTauConv)
246  {
247  /* for iterate to convergence, print reason why it was not converged
248  * on 3rd and higher iterations */
249  sprintf( conv.chNotConverged, "%s line '%s' %.3e=>%.3e",
250  (*badtr).system().chLabel.c_str(),
251  chLineLbl(*badtr).c_str(),
252  (*badtr).Emis().TauTot(), (*badtr).Emis().TauIn()*rt.DoubleTau);
253  }
254  }
255 
256  /* >>chng 03 sep 07, add this test */
257  /* check on changes in major column densities */
258  for( i=0; i<NCOLD; ++i )
259  {
260  /* was the species column density significant relative to
261  * the total H column density, and was its abundance changing? */
262  if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 &&
264  {
265  /* not converged to within conv.autocv, normally 20 percent */
267 
268  /* for iterate to convergence, print reason why it was not converged
269  * on 3rd and higher iterations */
270  strcpy( conv.chNotConverged, "H mole col" );
271 
272  if( save.lgPunConv )
273  {
274  lgReasonGiven = true;
275  fprintf( save.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
276  i,iteration,
277  colden.colden_old[i],
278  colden.colden[i],
280  }
281  }
282  }
283 
284  double biggestDiffer = 0.;
285  /* >>chng 03 sep 07, add this test */
286  /* check on changes in major column densities */
287  for( i=0; i<mole_global.num_calc; ++i )
288  {
289  if(mole_global.list[i]->isMonatomic())
290  continue;
291 
292  /* was the species abundance and changing? */
293  double differ = (double)fabs(mole.species[i].column_old-mole.species[i].column) ;
294  if( (mole.species[i].column/colden.colden[ipCOL_HTOT] > 1e-5) &&
295  (differ > conv.autocv*mole.species[i].column) )
296  {
297  /* not converged to within conv.autocv, normally 20 percent */
299 
300  /* for iterate to convergence, print reason why it was not converged
301  * on 3rd and higher iterations */
302  if( differ > biggestDiffer )
303  {
304  strcpy( conv.chNotConverged, mole_global.list[i]->label.c_str() );
305  strcat( conv.chNotConverged, " column" );
306  /*fprintf(ioQQQ,"debugggreset\t CO mole %li %li %.2e %.2e\n",
307  i,iteration,mole.species[i].column_old,mole.species[i].column);*/
308  biggestDiffer = differ;
309  }
310 
311  if( save.lgPunConv )
312  {
313  lgReasonGiven = true;
314  fprintf( save.ipPunConv, "%s, old:%.3e new:%.3e\n" ,
315  mole_global.list[i]->label.c_str(),
316  mole.species[i].column_old ,
317  mole.species[i].column );
318  }
319  }
320  }
321 
322  /* check on dynamical convergence in wind model with negative velocity */
323  if( dynamics.lgAdvection )
324  {
325  /* >>chng 02 nov 29, as per Will Henney email */
329  {
331  /* for iterate to convergence, print reason why it was not converged
332  * on 3rd and higher iterations */
333  strcpy( conv.chNotConverged, "Dynamics " );
334  if( save.lgPunConv )
335  {
336  lgReasonGiven = true;
337  fprintf( save.ipPunConv, " Dynamics\n" );
338  }
339  }
340  }
341 
342  if( save.lgPunConv )
343  {
345  fprintf( save.ipPunConv, " conv_itercheck exits converged\n" );
346  else
347  fprintf( save.ipPunConv, " conv_itercheck exits NOT converged\n" );
348  }
349 
350  /* lower limit to number of iterations if converged */
353 
354  /* test for stopping on first zone due to too-low temperature */
355  if( phycon.te < StopCalc.TempLoStopZone && nzone == 1 )
356  {
358  strcpy( conv.chNotConverged, " " );
360  }
361 
362  /* Fails if we have not fully implemented save convergence reason -
363  * should generate string stating reason, and set this flat,
364  * if lgOpticalDepthonverged is set true.
365  * These tests are only done when save output is requested
366  */
367  if( save.lgPunConv )
368  if( !iterations.lgOpticalDepthonverged && !lgReasonGiven )
369  TotalInsanity();
370  }
371  return;
372 }
t_mole_global mole_global
Definition: mole.cpp:7
t_colden colden
Definition: colden.cpp:5
const int ipHE_LIKE
Definition: iso.h:65
double convergence_tolerance
Definition: dynamics.h:162
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
int num_calc
Definition: mole.h:362
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:311
char chNotConverged[INPUT_LINE_LENGTH]
Definition: conv.h:128
const int ipHe2s3S
Definition: iso.h:46
bool lgAdvection
Definition: dynamics.h:66
realnum & TauTot() const
Definition: emission.h:490
t_StopCalc StopCalc
Definition: stopcalc.cpp:7
t_conv conv
Definition: conv.cpp:5
t_phycon phycon
Definition: phycon.cpp:6
realnum colden_old[NCOLD]
Definition: colden.h:32
bool lgAllTransitions
Definition: conv.h:252
long int nzone
Definition: cddefines.cpp:14
t_dynamics dynamics
Definition: dynamics.cpp:42
#define MIN2(a, b)
Definition: cddefines.h:803
long int nSpecies
Definition: taulines.cpp:22
t_dense dense
Definition: global.cpp:15
double error_scale2
Definition: dynamics.h:168
realnum autocv
Definition: conv.h:258
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
FILE * ipPunConv
Definition: save.h:454
long int iteration
Definition: cddefines.cpp:16
bool lgOpticalDepthonverged
Definition: iterations.h:57
bool lgPunConv
Definition: save.h:453
EmissionList::reference Emis() const
Definition: transition.h:447
t_mole_local mole
Definition: mole.cpp:8
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
double discretization_error
Definition: dynamics.h:165
t_iterations iterations
Definition: iterations.cpp:6
Definition: colden.h:17
realnum TempLoStopZone
Definition: stopcalc.h:42
bool lgElmtOn[LIMELM]
Definition: dense.h:160
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
int nLyaLevel[NISO]
Definition: iso.h:397
long int itermx
Definition: iterations.h:37
#define ASSERT(exp)
Definition: cddefines.h:613
const int ipH2s
Definition: iso.h:30
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1067
void ConvIterCheck(void)
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
const int ipHe2p3P2
Definition: iso.h:50
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
bool lgAutoIt
Definition: conv.h:249
vector< species > dBaseSpecies
Definition: taulines.cpp:15
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
MoleculeList list
Definition: mole.h:365
long int n_initial_relax
Definition: dynamics.h:132
realnum DoubleTau
Definition: rt.h:178
vector< TransitionList > AllTransitions
Definition: taulines.cpp:9
double convergence_error
Definition: dynamics.h:159
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
realnum colden[NCOLD]
Definition: colden.h:32
const int ipH3p
Definition: iso.h:33
realnum & TauIn() const
Definition: emission.h:470
t_rt rt
Definition: rt.cpp:5