00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "taulines.h"
00007 #include "iso.h"
00008 #include "phycon.h"
00009 #include "cddrive.h"
00010 #include "mole.h"
00011 #include "elementnames.h"
00012 #include "dynamics.h"
00013 #include "stopcalc.h"
00014 #include "dense.h"
00015 #include "iterations.h"
00016 #include "colden.h"
00017 #include "save.h"
00018 #include "rt.h"
00019 #include "conv.h"
00020
00021
00022
00023 void ConvIterCheck( void )
00024 {
00025 bool lgConverged;
00026 long int nelem,
00027 i,
00028 ipISO,
00029 ipHi, ipLo;
00030
00031 DEBUG_ENTRY( "ConvIterCheck()" );
00032
00033
00034
00035
00036
00037
00038
00039
00040 lgConverged = true;
00041 strcpy( conv.chNotConverged, "Converged!" );
00042
00043
00044 static double HbFracOutOld=-1. , HbFracOutNew=-1.;
00045 HbFracOutOld = HbFracOutNew;
00046 double a, total, BeamedIn;
00047 long int ipTotal = cdLine( "TOTL" , 4861.36f , &a , &total );
00048 long int ipInwd = cdLine( "INWD" , 4861.36f , &a , &BeamedIn );
00049 HbFracOutNew = 1. - pow(10. , (BeamedIn-total));
00050 ASSERT( iteration == 1 || (HbFracOutNew>=0 && HbFracOutNew<=1.) );
00051
00052 ipInwd = -1;
00053
00054 bool lgReasonGiven = false;
00055 if( iteration > 1 && conv.lgAutoIt )
00056 {
00057 if( nzone>3 && ipInwd>=0 && ipTotal>=0 )
00058 {
00059
00060 if( fabs(HbFracOutNew-HbFracOutOld)/HbFracOutNew> conv.autocv )
00061 {
00062 lgConverged = false;
00063 sprintf( conv.chNotConverged, "change in outward Hb");
00064 if( save.lgPunConv )
00065 {
00066 lgReasonGiven = true;
00067 fprintf( save.ipPunConv, " Change in outward Hb, "
00068 "old=%.3e new=%.3e \n",
00069 HbFracOutOld , HbFracOutNew);
00070 }
00071 }
00072 }
00073 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00074 {
00075 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00076 {
00077 if( dense.lgElmtOn[nelem] )
00078 {
00079
00080
00081
00082 if(ipISO==ipH_LIKE )
00083 {
00084 ipHi = ipH3p;
00085 ipLo = ipH2s;
00086 }
00087 else if( ipISO==ipHE_LIKE )
00088 {
00089 ipHi = ipHe2p3P2;
00090 ipLo = ipHe2s3S;
00091 }
00092 else
00093
00094 TotalInsanity();
00095
00096
00097
00098
00099 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() > 0.5 )
00100 {
00101
00102 if( fabs(iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() -
00103 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn()*rt.DoubleTau) >
00104 conv.autocv*fabs(iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn()*rt.DoubleTau) )
00105 {
00106
00107 lgConverged = false;
00108
00109
00110
00111 sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] );
00112
00113 if( save.lgPunConv )
00114 {
00115 lgReasonGiven = true;
00116 fprintf( save.ipPunConv, " %s-like Lya thick, "
00117 "nelem= %s iteration %li old %.3e new %.3e \n",
00118 elementnames.chElementSym[ipISO] ,
00119 elementnames.chElementSym[nelem],
00120 iteration,
00121 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
00122 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
00123 }
00124 }
00125
00126 if( fabs(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
00127 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()*rt.DoubleTau) >
00128 conv.autocv*fabs(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()*rt.DoubleTau) )
00129 {
00130
00131 lgConverged = false;
00132
00133
00134
00135 sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] );
00136
00137 if( save.lgPunConv )
00138 {
00139 lgReasonGiven = true;
00140 fprintf( save.ipPunConv, " %s-like subord, nelem= %s iteration %li old %.3e new %.3e \n" ,
00141 elementnames.chElementSym[ipISO],
00142 elementnames.chElementSym[nelem],
00143 iteration,
00144 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() ,
00145 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn() );
00146 }
00147 }
00148 }
00149 }
00150 }
00151 }
00152
00153 if(0)
00154 {
00155
00156 for( i=1; i <=nLevel1; i++ )
00157 {
00158 if( TauLines[i].Emis().TauIn() > 1. &&
00159 fabs(TauLines[i].Emis().TauTot() - TauLines[i].Emis().TauIn()*rt.DoubleTau) >
00160 conv.autocv*fabs(TauLines[i].Emis().TauIn()*rt.DoubleTau) )
00161 {
00162
00163 lgConverged = false;
00164
00165
00166
00167 sprintf( conv.chNotConverged, "level 1 line %.1f",TauLines[i].WLAng() );
00168
00169 if( save.lgPunConv )
00170 {
00171 lgReasonGiven = true;
00172 fprintf( save.ipPunConv, " level 1 line %.1f iteration %li old %.3e new %.3e \n",
00173 TauLines[i].WLAng(),
00174 iteration,
00175 TauLines[i].Emis().TauTot() ,
00176 TauLines[i].Emis().TauIn());
00177 }
00178 }
00179 }
00180 }
00181
00182
00183
00184 for( i=0; i<NCOLD; ++i )
00185 {
00186
00187
00188 if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 &&
00189 fabs(colden.colden_old[i]-colden.colden[i]) > conv.autocv*colden.colden[i] )
00190 {
00191
00192 lgConverged = false;
00193
00194
00195
00196 strcpy( conv.chNotConverged, "H mole col" );
00197
00198 if( save.lgPunConv )
00199 {
00200 lgReasonGiven = true;
00201 fprintf( save.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
00202 i,iteration,
00203 colden.colden_old[i],
00204 colden.colden[i],
00205 colden.colden[ipCOL_HTOT] );
00206 }
00207 }
00208 }
00209
00210 double biggestDiffer = 0.;
00211
00212
00213 for( i=0; i<mole_global.num_calc; ++i )
00214 {
00215 if(mole_global.list[i]->isMonatomic())
00216 continue;
00217
00218
00219 double differ = (double)fabs(mole.species[i].column_old-mole.species[i].column) ;
00220 if( (mole.species[i].column/colden.colden[ipCOL_HTOT] > 1e-5) &&
00221 (differ > conv.autocv*mole.species[i].column) )
00222 {
00223
00224 lgConverged = false;
00225
00226
00227
00228 if( differ > biggestDiffer )
00229 {
00230 strcpy( conv.chNotConverged, mole_global.list[i]->label.c_str() );
00231 strcat( conv.chNotConverged, " column" );
00232
00233
00234 biggestDiffer = differ;
00235 }
00236
00237 if( save.lgPunConv )
00238 {
00239 lgReasonGiven = true;
00240 fprintf( save.ipPunConv, "%s, old:%.3e new:%.3e\n" ,
00241 mole_global.list[i]->label.c_str(),
00242 mole.species[i].column_old ,
00243 mole.species[i].column );
00244 }
00245 }
00246 }
00247
00248
00249 if( dynamics.lgAdvection )
00250 {
00251
00252 if( dynamics.convergence_error > conv.autocv*dynamics.error_scale2*dynamics.convergence_tolerance ||
00253 dynamics.discretization_error > conv.autocv*dynamics.error_scale2 )
00254 {
00255 lgConverged = false;
00256
00257
00258 strcpy( conv.chNotConverged, "Dynamics " );
00259 if( save.lgPunConv )
00260 {
00261 lgReasonGiven = true;
00262 fprintf( save.ipPunConv, " Dynamics\n" );
00263 }
00264 }
00265 }
00266
00267 if( save.lgPunConv && lgConverged )
00268 {
00269 lgReasonGiven = true;
00270 fprintf( save.ipPunConv, " converged\n" );
00271 }
00272
00273
00274 if( lgConverged )
00275 iterations.itermx = MIN2(iterations.itermx,iteration);
00276
00277
00278
00279 if( phycon.te < StopCalc.TempLoStopZone && nzone == 1 )
00280 {
00281 lgConverged = true;
00282 strcpy( conv.chNotConverged, " " );
00283 iterations.itermx = MIN2(iterations.itermx,iteration);
00284 }
00285
00286 ASSERT( !save.lgPunConv || lgReasonGiven );
00287 }
00288 return;
00289 }