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( Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn > 0.5 )
00100 {
00101
00102 if( fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot -
00103 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) >
00104 conv.autocv*fabs(Transitions[ipISO][nelem][iso.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 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot ,
00122 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
00123 }
00124 }
00125
00126 if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot -
00127 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) >
00128 conv.autocv*fabs(Transitions[ipISO][nelem][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 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot ,
00145 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn );
00146 }
00147 }
00148 }
00149 }
00150 }
00151 }
00152
00153
00154
00155 for( i=0; i<NCOLD; ++i )
00156 {
00157
00158
00159 if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 &&
00160 fabs(colden.colden_old[i]-colden.colden[i]) > conv.autocv*colden.colden[i] )
00161 {
00162
00163 lgConverged = false;
00164
00165
00166
00167 strcpy( conv.chNotConverged, "H mole col" );
00168
00169 if( save.lgPunConv )
00170 {
00171 lgReasonGiven = true;
00172 fprintf( save.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
00173 i,iteration,
00174 colden.colden_old[i],
00175 colden.colden[i],
00176 colden.colden[ipCOL_HTOT] );
00177 }
00178 }
00179 }
00180
00181
00182
00183 for( i=0; i<mole.num_comole_calc; ++i )
00184 {
00185 if(COmole[i]->n_nuclei == 1)
00186 continue;
00187
00188
00189 if( COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5 &&
00190 fabs(COmole[i]->hevcol_old-COmole[i]->hevcol) > conv.autocv*COmole[i]->hevcol )
00191 {
00192
00193 lgConverged = false;
00194
00195
00196
00197 strcpy( conv.chNotConverged, "CO mol col" );
00198
00199
00200
00201 if( save.lgPunConv )
00202 {
00203 lgReasonGiven = true;
00204 fprintf( save.ipPunConv, "CO mol col, old:%.3e new:%.3e\n" ,
00205 COmole[i]->hevcol_old ,
00206 COmole[i]->hevcol );
00207 }
00208 }
00209 }
00210
00211
00212 if( dynamics.lgAdvection )
00213 {
00214
00215 if( dynamics.convergence_error > conv.autocv*dynamics.error_scale2*dynamics.convergence_tolerance ||
00216 dynamics.discretization_error > conv.autocv*dynamics.error_scale2 )
00217 {
00218 lgConverged = false;
00219
00220
00221 strcpy( conv.chNotConverged, "Dynamics " );
00222 if( save.lgPunConv )
00223 {
00224 lgReasonGiven = true;
00225 fprintf( save.ipPunConv, " Dynamics\n" );
00226 }
00227 }
00228 }
00229
00230 if( save.lgPunConv && lgConverged )
00231 {
00232 lgReasonGiven = true;
00233 fprintf( save.ipPunConv, " converged\n" );
00234 }
00235
00236
00237 if( lgConverged )
00238 iterations.itermx = MIN2(iterations.itermx,iteration);
00239
00240
00241
00242 if( phycon.te < StopCalc.TempLoStopZone && nzone == 1 )
00243 {
00244 lgConverged = true;
00245 strcpy( conv.chNotConverged, " " );
00246 iterations.itermx = MIN2(iterations.itermx,iteration);
00247 }
00248
00249 ASSERT( !save.lgPunConv || lgReasonGiven );
00250 }
00251 return;
00252 }