00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "taulines.h"
00007 #include "iso.h"
00008 #include "phycon.h"
00009 #include "mole.h"
00010 #include "elementnames.h"
00011 #include "dynamics.h"
00012 #include "stopcalc.h"
00013 #include "dense.h"
00014 #include "iterations.h"
00015 #include "colden.h"
00016 #include "punch.h"
00017 #include "rt.h"
00018 #include "conv.h"
00019
00020
00021
00022 void ConvIterCheck( void )
00023 {
00024 bool lgConverged;
00025 long int nelem,
00026 i,
00027 ipISO,
00028 ipHi, ipLo;
00029
00030 DEBUG_ENTRY( "ConvIterCheck()" );
00031
00032
00033
00034
00035
00036
00037
00038
00039 lgConverged = true;
00040 strcpy( conv.chNotConverged, "Converged!" );
00041 if( iteration > 1 && conv.lgAutoIt )
00042 {
00043 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00044 {
00045 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00046 {
00047 if( dense.lgElmtOn[nelem] )
00048 {
00049
00050
00051
00052 if( true && EmisLines[ipISO][nelem][3][2].TauIn > 0.5 )
00053 {
00054
00055 if( fabs(EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].TauTot/
00056 (EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].TauIn*rt.DoubleTau)-1.) >
00057 conv.autocv )
00058 {
00059
00060 lgConverged = false;
00061
00062
00063
00064 sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] );
00065
00066 if( punch.lgPunConv )
00067 {
00068 fprintf( punch.ipPunConv, " %s-like Lya, nelem= %li iteration %li old %.3e new %.3e \n",
00069 elementnames.chElementSym[ipISO] ,
00070 nelem,
00071 iteration,
00072 EmisLines[ipISO][nelem][ipH2p][ipH1s].TauTot ,
00073 EmisLines[ipISO][nelem][ipH2p][ipH1s].TauIn );
00074 }
00075 }
00076
00077
00078
00079
00080 if(ipISO==ipH_LIKE )
00081 {
00082 ipHi = 3;
00083 ipLo = 2;
00084 }
00085 else if( ipISO==ipHE_LIKE )
00086 {
00087 ipHi = ipHe2p3P2;
00088 ipLo = ipHe2s3S;
00089 }
00090 else
00091
00092 TotalInsanity();
00093
00094 if( fabs(EmisLines[ipISO][nelem][ipHi][ipLo].TauTot/
00095 (EmisLines[ipISO][nelem][ipHi][ipLo].TauIn*rt.DoubleTau)-1.) > conv.autocv )
00096 {
00097
00098 lgConverged = false;
00099
00100
00101
00102 sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] );
00103
00104 if( punch.lgPunConv )
00105 {
00106 fprintf( punch.ipPunConv, " %s-like subord, nelem= %li iteration %li old %.3e new %.3e \n" ,
00107 elementnames.chElementSym[ipISO],
00108 nelem, iteration,
00109 EmisLines[ipISO][nelem][ipHi][ipLo].TauTot ,
00110 EmisLines[ipISO][nelem][ipHi][ipLo].TauIn );
00111 }
00112 }
00113 }
00114 }
00115 }
00116 }
00117
00118
00119
00120 for( i=0; i<NCOLD; ++i )
00121 {
00122 double differ = fabs(colden.colden_old[i]-colden.colden[i]) /
00123 SDIV(colden.colden[i]);
00124
00125
00126 if( (colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5) &&
00127 (differ > conv.autocv) )
00128 {
00129
00130 lgConverged = false;
00131
00132
00133
00134 strcpy( conv.chNotConverged, "H mole col" );
00135
00136 if( punch.lgPunConv )
00137 {
00138 fprintf( punch.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
00139 i,iteration,
00140 colden.colden_old[i],
00141 colden.colden[i],
00142 colden.colden[ipCOL_HTOT] );
00143 }
00144 }
00145 }
00146
00147
00148
00149 for( i=0; i<mole.num_comole_calc; ++i )
00150 {
00151 double differ;
00152 if(COmole[i]->n_nuclei == 1)
00153 continue;
00154
00155
00156 differ = fabs(COmole[i]->hevcol_old-COmole[i]->hevcol) /
00157 SDIV(COmole[i]->hevcol);
00158 if( (COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5) &&
00159 (differ > conv.autocv) )
00160 {
00161
00162 lgConverged = false;
00163
00164
00165
00166 strcpy( conv.chNotConverged, "CO mol col" );
00167
00168
00169
00170 if( punch.lgPunConv )
00171 {
00172 fprintf( punch.ipPunConv, "CO mol col, old:%.3e new:%.3e\n" ,
00173 COmole[i]->hevcol_old ,
00174 COmole[i]->hevcol );
00175 }
00176 }
00177 }
00178
00179
00180 if( dynamics.lgAdvection )
00181 {
00182 double error1 = dynamics.convergence_error / SDIV(dynamics.error_scale2);
00183 double error2 = dynamics.discretization_error / SDIV(dynamics.error_scale2);
00184
00185 if( MAX2(error1/dynamics.convergence_tolerance, error2 ) > conv.autocv )
00186
00187 {
00188 lgConverged = false;
00189
00190
00191 strcpy( conv.chNotConverged, "Dynamics " );
00192 if( punch.lgPunConv )
00193 {
00194 fprintf( punch.ipPunConv, " Dynamics\n" );
00195 }
00196 }
00197 }
00198
00199 if( punch.lgPunConv && lgConverged )
00200 {
00201 fprintf( punch.ipPunConv, " converged\n" );
00202 }
00203
00204
00205 if( lgConverged )
00206 iterations.itermx = MIN2(iterations.itermx,iteration);
00207
00208
00209
00210 if( phycon.te < StopCalc.tend && nzone == 1 )
00211 {
00212 lgConverged = true;
00213 strcpy( conv.chNotConverged, " " );
00214 iterations.itermx = MIN2(iterations.itermx,iteration);
00215 }
00216 }
00217
00218 DEBUG_EXIT( "ConvIterCheck()" );
00219 return;
00220
00221 }