cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ConvIterCheck check whether model has converged or whether more iterations 00004 * are needed - implements the iter to converg comnd */ 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 /*ConvIterCheck check whether model has converged or whether more iterations 00021 * are needed - implements the iterate to convergence command */ 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 /* this is an option to keep iterating until it converges 00034 * iterate to convergence option 00035 * autocv is percentage difference in optical depths allowed, 00036 * =0.20 in block data 00037 * checking on Ly and Balmer lines */ 00038 /*>>chng 04 oct 19, promote loop to do all iso-electronic series */ 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 /* now check if major subordinate line is converged - for H-like this will 00050 * be Ha, and for He-like, the 23P - 23S transition - this will not work for 00051 * NISO > 2 so must check against this */ 00052 if(ipISO==ipH_LIKE ) 00053 { 00054 ipHi = ipH3p; 00055 ipLo = ipH2s; 00056 } 00057 else if( ipISO==ipHE_LIKE ) 00058 { 00059 ipHi = ipHe2p3P2; 00060 ipLo = ipHe2s3S; 00061 } 00062 else 00063 /* fails when NISO increased, add more sequences */ 00064 TotalInsanity(); 00065 00066 /* check both H-alpha and Ly-alpha for all species - 00067 * only if Balmer lines thick 00068 * so check if Ha optical depth significant */ 00069 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn > 0.5 ) 00070 { 00071 /* test if Lya converged, nLyaLevel is upper level of Lya for iso seq */ 00072 if( fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot - 00073 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) > 00074 conv.autocv*fabs(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn*rt.DoubleTau) ) 00075 { 00076 /* not converged to within AUTOCV, normally 15 percent */ 00077 lgConverged = false; 00078 00079 /* for iterate to convergence, print reason why it was not converged 00080 * on 3rd and higher iterations */ 00081 sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] ); 00082 00083 if( punch.lgPunConv ) 00084 { 00085 fprintf( punch.ipPunConv, " %s-like Lya thick, " 00086 "nelem= %s iteration %li old %.3e new %.3e \n", 00087 elementnames.chElementSym[ipISO] , 00088 elementnames.chElementSym[nelem], 00089 iteration, 00090 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot , 00091 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn); 00092 } 00093 } 00094 00095 if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot - 00096 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) > 00097 conv.autocv*fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn*rt.DoubleTau) ) 00098 { 00099 /* not converged to within AUTOCV, normally 15 percent */ 00100 lgConverged = false; 00101 00102 /* for iterate to convergence, print reason why it was not converged 00103 * on 3rd and higher iterations */ 00104 sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] ); 00105 00106 if( punch.lgPunConv ) 00107 { 00108 fprintf( punch.ipPunConv, " %s-like subord, nelem= %s iteration %li old %.3e new %.3e \n" , 00109 elementnames.chElementSym[ipISO], 00110 elementnames.chElementSym[nelem], 00111 iteration, 00112 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot , 00113 Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn ); 00114 } 00115 } 00116 } 00117 } 00118 } 00119 } 00120 00121 /* >>chng 03 sep 07, add this test */ 00122 /* check on changes in major column densities */ 00123 for( i=0; i<NCOLD; ++i ) 00124 { 00125 /* was the species column density significant relative to 00126 * the total H column density, and was its abundance changing? */ 00127 if( colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5 && 00128 fabs(colden.colden_old[i]-colden.colden[i]) > conv.autocv*colden.colden[i] ) 00129 { 00130 /* not converged to within conv.autocv, normally 20 percent */ 00131 lgConverged = false; 00132 00133 /* for iterate to convergence, print reason why it was not converged 00134 * on 3rd and higher iterations */ 00135 strcpy( conv.chNotConverged, "H mole col" ); 00136 00137 if( punch.lgPunConv ) 00138 { 00139 fprintf( punch.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n", 00140 i,iteration, 00141 colden.colden_old[i], 00142 colden.colden[i], 00143 colden.colden[ipCOL_HTOT] ); 00144 } 00145 } 00146 } 00147 00148 /* >>chng 03 sep 07, add this test */ 00149 /* check on changes in major column densities */ 00150 for( i=0; i<mole.num_comole_calc; ++i ) 00151 { 00152 if(COmole[i]->n_nuclei == 1) 00153 continue; 00154 00155 /* was the species abundance and changing? */ 00156 if( COmole[i]->hevcol/colden.colden[ipCOL_HTOT] > 1e-5 && 00157 fabs(COmole[i]->hevcol_old-COmole[i]->hevcol) > conv.autocv*COmole[i]->hevcol ) 00158 { 00159 /* not converged to within conv.autocv, normally 20 percent */ 00160 lgConverged = false; 00161 00162 /* for iterate to convergence, print reason why it was not converged 00163 * on 3rd and higher iterations */ 00164 strcpy( conv.chNotConverged, "CO mol col" ); 00165 /*fprintf(ioQQQ,"debugggreset\t CO mole %li %li %.2e %.2e\n", 00166 i,iteration,COmole[i]->hevcol_old,COmole[i]->hevcol);*/ 00167 00168 if( punch.lgPunConv ) 00169 { 00170 fprintf( punch.ipPunConv, "CO mol col, old:%.3e new:%.3e\n" , 00171 COmole[i]->hevcol_old , 00172 COmole[i]->hevcol ); 00173 } 00174 } 00175 } 00176 00177 /* check on dynamical convergence in wind model with negative velocity */ 00178 if( dynamics.lgAdvection ) 00179 { 00180 /* >>chng 02 nov 29, as per Will Henney email */ 00181 if( dynamics.convergence_error > conv.autocv*dynamics.error_scale2*dynamics.convergence_tolerance || 00182 dynamics.discretization_error > conv.autocv*dynamics.error_scale2 ) 00183 { 00184 lgConverged = false; 00185 /* for iterate to convergence, print reason why it was not converged 00186 * on 3rd and higher iterations */ 00187 strcpy( conv.chNotConverged, "Dynamics " ); 00188 if( punch.lgPunConv ) 00189 { 00190 fprintf( punch.ipPunConv, " Dynamics\n" ); 00191 } 00192 } 00193 } 00194 00195 if( punch.lgPunConv && lgConverged ) 00196 { 00197 fprintf( punch.ipPunConv, " converged\n" ); 00198 } 00199 00200 /* lower limit to number of iterations if converged */ 00201 if( lgConverged ) 00202 iterations.itermx = MIN2(iterations.itermx,iteration); 00203 00204 /* >>chng 96 dec 20, moved following to within if on lgAutoIt 00205 * this is test for stopping on first zone */ 00206 if( phycon.te < StopCalc.tend && nzone == 1 ) 00207 { 00208 lgConverged = true; 00209 strcpy( conv.chNotConverged, " " ); 00210 iterations.itermx = MIN2(iterations.itermx,iteration); 00211 } 00212 } 00213 return; 00214 }