cloudy trunk
Loading...
Searching...
No Matches
conv_itercheck.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/*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 */
23void ConvIterCheck( void )
24{
25 bool lgConverged;
26 long int nelem,
27 i,
28 ipISO,
29 ipHi, ipLo;
30
31 DEBUG_ENTRY( "ConvIterCheck()" );
32
33 /* =======================================================================*/
34 /* this is an option to keep iterating until it converges
35 * iterate to convergence option
36 * autocv is percentage difference in optical depths allowed,
37 * =0.20 in block data
38 * checking on Ly and Balmer lines */
39 /*>>chng 04 oct 19, promote loop to do all iso-electronic series */
40 lgConverged = true;
41 strcpy( conv.chNotConverged, "Converged!" );
42
43 // set up intensities used to converge outward intensity of Hb
44 static double HbFracOutOld=-1. , HbFracOutNew=-1.;
45 HbFracOutOld = HbFracOutNew;
46 double a, total, BeamedIn;
47 long int ipTotal = cdLine( "TOTL" , 4861.36f , &a , &total );
48 long int ipInwd = cdLine( "INWD" , 4861.36f , &a , &BeamedIn );
49 HbFracOutNew = 1. - pow(10. , (BeamedIn-total));
50 ASSERT( iteration == 1 || (HbFracOutNew>=0 && HbFracOutNew<=1.) );
51 // this disables the test on the outward Hb
52 ipInwd = -1;
53
54 bool lgReasonGiven = false;
55 if( iteration > 1 && conv.lgAutoIt )
56 {
57 if( nzone>3 && ipInwd>=0 && ipTotal>=0 )
58 {
59 // check whether outward intensity of Hb has converged
60 if( fabs(HbFracOutNew-HbFracOutOld)/HbFracOutNew> conv.autocv )
61 {
62 lgConverged = false;
63 sprintf( conv.chNotConverged, "change in outward Hb");
64 if( save.lgPunConv )
65 {
66 lgReasonGiven = true;
67 fprintf( save.ipPunConv, " Change in outward Hb, "
68 "old=%.3e new=%.3e \n",
69 HbFracOutOld , HbFracOutNew);
70 }
71 }
72 }
73 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
74 {
75 for( nelem=ipISO; nelem < LIMELM; nelem++ )
76 {
77 if( dense.lgElmtOn[nelem] )
78 {
79 /* now check if major subordinate line is converged - for H-like this will
80 * be Ha, and for He-like, the 23P - 23S transition - this will not work for
81 * NISO > 2 so must check against this */
82 if(ipISO==ipH_LIKE )
83 {
84 ipHi = ipH3p;
85 ipLo = ipH2s;
86 }
87 else if( ipISO==ipHE_LIKE )
88 {
89 ipHi = ipHe2p3P2;
90 ipLo = ipHe2s3S;
91 }
92 else
93 /* fails when NISO increased, add more sequences */
95
96 /* check both H-alpha and Ly-alpha for all species -
97 * only if Balmer lines thick
98 * so check if Ha optical depth significant */
99 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() > 0.5 )
100 {
101 /* test if Lya converged, nLyaLevel is upper level of Lya for iso seq */
102 if( fabs(iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() -
103 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn()*rt.DoubleTau) >
104 conv.autocv*fabs(iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn()*rt.DoubleTau) )
105 {
106 /* not converged to within AUTOCV, normally 15 percent */
107 lgConverged = false;
108
109 /* for iterate to convergence, print reason why it was not converged
110 * on 3rd and higher iterations */
111 sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] );
112
113 if( save.lgPunConv )
114 {
115 lgReasonGiven = true;
116 fprintf( save.ipPunConv, " %s-like Lya thick, "
117 "nelem= %s iteration %li old %.3e new %.3e \n",
118 elementnames.chElementSym[ipISO] ,
119 elementnames.chElementSym[nelem],
120 iteration,
121 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
122 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
123 }
124 }
125
126 if( fabs(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
127 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()*rt.DoubleTau) >
128 conv.autocv*fabs(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()*rt.DoubleTau) )
129 {
130 /* not converged to within AUTOCV, normally 15 percent */
131 lgConverged = false;
132
133 /* for iterate to convergence, print reason why it was not converged
134 * on 3rd and higher iterations */
135 sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] );
136
137 if( save.lgPunConv )
138 {
139 lgReasonGiven = true;
140 fprintf( save.ipPunConv, " %s-like subord, nelem= %s iteration %li old %.3e new %.3e \n" ,
141 elementnames.chElementSym[ipISO],
142 elementnames.chElementSym[nelem],
143 iteration,
144 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() ,
145 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() );
146 }
147 }
148 }
149 }
150 }
151 }
152
153 if(0)
154 {
155 // check convergence of all level 1 lines
156 for( i=1; i <=nLevel1; i++ )
157 {
158 if( TauLines[i].Emis().TauIn() > 1. &&
159 fabs(TauLines[i].Emis().TauTot() - TauLines[i].Emis().TauIn()*rt.DoubleTau) >
160 conv.autocv*fabs(TauLines[i].Emis().TauIn()*rt.DoubleTau) )
161 {
162 /* not converged to within AUTOCV, normally 15 percent */
163 lgConverged = false;
164
165 /* for iterate to convergence, print reason why it was not converged
166 * on 3rd and higher iterations */
167 sprintf( conv.chNotConverged, "level 1 line %.1f",TauLines[i].WLAng() );
168
169 if( save.lgPunConv )
170 {
171 lgReasonGiven = true;
172 fprintf( save.ipPunConv, " level 1 line %.1f iteration %li old %.3e new %.3e \n",
173 TauLines[i].WLAng(),
174 iteration,
175 TauLines[i].Emis().TauTot() ,
176 TauLines[i].Emis().TauIn());
177 }
178 }
179 }
180 }
181
182 /* >>chng 03 sep 07, add this test */
183 /* check on changes in major column densities */
184 for( i=0; i<NCOLD; ++i )
185 {
186 /* was the species column density significant relative to
187 * the total H column density, and was its abundance changing? */
188 if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 &&
189 fabs(colden.colden_old[i]-colden.colden[i]) > conv.autocv*colden.colden[i] )
190 {
191 /* not converged to within conv.autocv, normally 20 percent */
192 lgConverged = false;
193
194 /* for iterate to convergence, print reason why it was not converged
195 * on 3rd and higher iterations */
196 strcpy( conv.chNotConverged, "H mole col" );
197
198 if( save.lgPunConv )
199 {
200 lgReasonGiven = true;
201 fprintf( save.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
202 i,iteration,
203 colden.colden_old[i],
204 colden.colden[i],
205 colden.colden[ipCOL_HTOT] );
206 }
207 }
208 }
209
210 double biggestDiffer = 0.;
211 /* >>chng 03 sep 07, add this test */
212 /* check on changes in major column densities */
213 for( i=0; i<mole_global.num_calc; ++i )
214 {
215 if(mole_global.list[i]->isMonatomic())
216 continue;
217
218 /* was the species abundance and changing? */
219 double differ = (double)fabs(mole.species[i].column_old-mole.species[i].column) ;
220 if( (mole.species[i].column/colden.colden[ipCOL_HTOT] > 1e-5) &&
221 (differ > conv.autocv*mole.species[i].column) )
222 {
223 /* not converged to within conv.autocv, normally 20 percent */
224 lgConverged = false;
225
226 /* for iterate to convergence, print reason why it was not converged
227 * on 3rd and higher iterations */
228 if( differ > biggestDiffer )
229 {
230 strcpy( conv.chNotConverged, mole_global.list[i]->label.c_str() );
231 strcat( conv.chNotConverged, " column" );
232 /*fprintf(ioQQQ,"debugggreset\t CO mole %li %li %.2e %.2e\n",
233 i,iteration,mole.species[i].column_old,mole.species[i].column);*/
234 biggestDiffer = differ;
235 }
236
237 if( save.lgPunConv )
238 {
239 lgReasonGiven = true;
240 fprintf( save.ipPunConv, "%s, old:%.3e new:%.3e\n" ,
241 mole_global.list[i]->label.c_str(),
242 mole.species[i].column_old ,
243 mole.species[i].column );
244 }
245 }
246 }
247
248 /* check on dynamical convergence in wind model with negative velocity */
249 if( dynamics.lgAdvection )
250 {
251 /* >>chng 02 nov 29, as per Will Henney email */
252 if( dynamics.convergence_error > conv.autocv*dynamics.error_scale2*dynamics.convergence_tolerance ||
253 dynamics.discretization_error > conv.autocv*dynamics.error_scale2 )
254 {
255 lgConverged = false;
256 /* for iterate to convergence, print reason why it was not converged
257 * on 3rd and higher iterations */
258 strcpy( conv.chNotConverged, "Dynamics " );
259 if( save.lgPunConv )
260 {
261 lgReasonGiven = true;
262 fprintf( save.ipPunConv, " Dynamics\n" );
263 }
264 }
265 }
266
267 if( save.lgPunConv && lgConverged )
268 {
269 lgReasonGiven = true;
270 fprintf( save.ipPunConv, " converged\n" );
271 }
272
273 /* lower limit to number of iterations if converged */
274 if( lgConverged )
275 iterations.itermx = MIN2(iterations.itermx,iteration);
276
277 /* >>chng 96 dec 20, moved following to within if on lgAutoIt
278 * this is test for stopping on first zone */
279 if( phycon.te < StopCalc.TempLoStopZone && nzone == 1 )
280 {
281 lgConverged = true;
282 strcpy( conv.chNotConverged, " " );
283 iterations.itermx = MIN2(iterations.itermx,iteration);
284 }
285 // fails if we have not fully implemented save convergence reason
286 ASSERT( !save.lgPunConv || lgReasonGiven );
287 }
288 return;
289}
long int nzone
Definition cddefines.cpp:14
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
t_colden colden
Definition colden.cpp:5
#define ipCOL_HTOT
Definition colden.h:12
#define NCOLD
Definition colden.h:9
t_conv conv
Definition conv.cpp:5
void ConvIterCheck(void)
t_dense dense
Definition dense.cpp:24
t_dynamics dynamics
Definition dynamics.cpp:44
t_elementnames elementnames
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH3p
Definition iso.h:31
const int ipHe2s3S
Definition iso.h:44
const int ipHE_LIKE
Definition iso.h:63
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
const int ipHe2p3P2
Definition iso.h:48
t_iterations iterations
Definition iterations.cpp:5
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
t_phycon phycon
Definition phycon.cpp:6
t_rt rt
Definition rt.cpp:5
t_save save
Definition save.cpp:5
t_StopCalc StopCalc
Definition stopcalc.cpp:5
long int nLevel1
Definition taulines.cpp:28
TransitionList TauLines("TauLines", &AnonStates)