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 /*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIonize, 00004 * called by cloudy */ 00005 /*ConvFail handle conergece failure */ 00006 #include "cddefines.h" 00007 #include "phycon.h" 00008 #include "rt.h" 00009 #include "dense.h" 00010 #include "pressure.h" 00011 #include "trace.h" 00012 #include "conv.h" 00013 #include "grainvar.h" 00014 #include "grains.h" 00015 00016 /* the limit to the number of loops */ 00017 /* >>chng 02 jun 13, from 40 to 50 */ 00018 static const int LOOPMAX = 100; 00019 00020 /*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIoniz, 00021 * called by cloudy 00022 * returns 0 if ok, 1 if disaster */ 00023 int ConvPresTempEdenIoniz(void) 00024 { 00025 long int loop, 00026 LoopMax=LOOPMAX; 00027 double 00028 hden_old , 00029 hden_chng_old , 00030 hden_chng, 00031 /*pres_old ,*/ 00032 /*pres_chng_old ,*/ 00033 /*pres_chng,*/ 00034 /*old_slope,*/ 00035 /*rel_slope,*/ 00036 dP_chng_factor; 00037 bool lgPresOscil; 00038 long int nloop_pres_oscil; 00039 double TemperatureInitial; 00040 00041 DEBUG_ENTRY( "ConvPresTempEdenIoniz()" ); 00042 00043 /* this will count number of times we call ConvBase in this zone, 00044 * counter is incremented there */ 00045 conv.nPres2Ioniz = 0; 00046 loop = 0; 00047 00048 /* this will be the limit, which we will increase if no oscillations occur */ 00049 LoopMax = LOOPMAX; 00050 /* set the initial temperature to the current value, so we will know 00051 * if we are trying to jump over a thermal front */ 00052 TemperatureInitial = phycon.te; 00053 00054 /* this will be flag to check for pressure oscillations */ 00055 lgPresOscil = false; 00056 /* this is loop where it happened */ 00057 nloop_pres_oscil = 0; 00058 /* should still be true at end */ 00059 conv.lgConvPops = true; 00060 00061 /* we will use these to check whether hden oscillating - would need to decrease step size */ 00062 hden_old = dense.gas_phase[ipHYDROGEN]; 00063 hden_chng = 0.; 00064 /*pres_old = pressure.PresTotlCurr;*/ 00065 /*pres_chng = 0.;*/ 00066 00067 /* this is dP_chng_factor 00068 * cut in half when pressure error changes sign */ 00069 dP_chng_factor = 1.; 00070 /*rel_slope = 0.;*/ 00071 00072 if( trace.nTrConvg>=1 ) 00073 { 00074 fprintf( ioQQQ, 00075 " ConvPresTempEdenIoniz1 entered, will call ConvIoniz to initialize\n"); 00076 } 00077 00078 /* converge the ionization first, so that we know where we are, and have 00079 * a valid foundation to begin the search */ 00080 /* the true electron density dense.EdenTrue is set in eden_sum called by ConvBase */ 00081 00082 /* chng 02 dec 11 rjrw -- ConvIoniz => ConvTempEdenIoniz() here for consistency with inside loop */ 00083 /* ConvIoniz; */ 00084 if( ConvTempEdenIoniz() ) 00085 { 00086 00087 return 1; 00088 } 00089 00090 /* this evaluates current pressure, and returns whether or not 00091 * it is within tolerance of correct pressure */ 00092 conv.lgConvPres = false; /* lgConvPres(); */ 00093 00094 /* convergence trace at this level */ 00095 if( trace.nTrConvg>=1 ) 00096 { 00097 fprintf( ioQQQ, 00098 " ConvPresTempEdenIoniz1 ConvIoniz found following converged: Pres;%c, Eden;%c, Temp;%c, Ion:%c Pops:%c\n", 00099 TorF(conv.lgConvPres) , 00100 TorF(lgConvEden() ), 00101 TorF(lgConvTemp() ) , 00102 TorF(conv.lgConvIoniz) , 00103 TorF(conv.lgConvPops)); 00104 } 00105 00106 /* >>chng 01 apr 01, add test for at least 2 loops to get better pressure convergence */ 00107 /* >>chng 01 oct 31, add test for number of times converged, for constant 00108 * pressure will get two valid solutions */ 00109 /*while( (loop < LoopMax) && !(conv.lgConvPres &&nPresConverged > 1) && !lgAbort )*/ 00110 /* >>chng 02 dec 12, do not demand two constant pressure being valid - should not be 00111 * necessary if first one really is valid, and this caused unneeded second evaluation 00112 * in the constant density cases */ 00113 00114 /* trace convergence print at this level */ 00115 if( trace.nTrConvg>=1 ) 00116 { 00117 fprintf( ioQQQ, 00118 "\n ConvPresTempEdenIoniz1 entering main pressure loop.\n"); 00119 } 00120 while( (loop < LoopMax) && !conv.lgConvPres && !lgAbort ) 00121 { 00122 /* there can be a pressure or density oscillation early in the search - if not persistent 00123 * ok to clear flag */ 00124 /* >>chng 01 aug 24, if large change in temperature allow lots more loops */ 00125 if( fabs( TemperatureInitial - phycon.te )/phycon.te > 0.3 ) 00126 LoopMax = 2*LOOPMAX; 00127 /* if start of calculation and we are solving for set pressure, 00128 * allow a lot more iterations */ 00129 if( nzone <= 1 && pressure.lgPressureInitialSpecified ) 00130 LoopMax = 10*LOOPMAX; 00131 00132 /* change current densities of all constituents if necessary, 00133 * PressureChange evaluates lgPresOK, true if pressure is now ok 00134 * sets CurrentPressure and CorrectPressure */ 00135 hden_old = dense.gas_phase[ipHYDROGEN]; 00136 /*pres_old = pressure.PresTotlCurr;*/ 00137 00138 /* this will evaluate current pressure, update the densities, 00139 * determine the wind velocity, and set conv.lgConvPres, 00140 * return value is true if density was changed, false if no changes were necessary 00141 * if density changes then we must redo the temperature and ionization 00142 * PressureChange contains the logic that determines how to change the density to get 00143 * the right pressure */ 00144 if( PressureChange( dP_chng_factor ) ) 00145 { 00146 /* heating cooling balance while doing ionization, 00147 * this is where the heavy lifting is done, this calls PresTotCurrent, 00148 * which sets pressure.PresTotlCurr */ 00149 if( ConvTempEdenIoniz() ) 00150 { 00151 return 1; 00152 } 00153 } 00154 00155 /* if product of these two is negative then hden is oscillating */ 00156 hden_chng_old = hden_chng; 00157 /*pres_chng_old = pres_chng;*/ 00158 hden_chng = dense.gas_phase[ipHYDROGEN] - hden_old; 00159 if( fabs(hden_chng) < SMALLFLOAT ) 00160 hden_chng = sign( (double)SMALLFLOAT, hden_chng ); 00161 /*pres_chng = pressure.PresTotlCurr - pres_old;*/ 00162 /*old_slope = rel_slope;*/ 00163 /*rel_slope = (pres_chng/pressure.PresTotlCurr) / (hden_chng/dense.gas_phase[ipHYDROGEN]);*/ 00164 00165 { 00166 /*@-redef@*/ 00167 enum{DEBUG_LOC=false}; 00168 /*@+redef@*/ 00169 if( DEBUG_LOC && nzone > 150 && iteration > 1 ) 00170 { 00171 fprintf(ioQQQ,"%li\t%.2e\t%.2e\t%.2e\n", 00172 nzone, 00173 pressure.PresTotlCurr, 00174 pressure.PresTotlCorrect, 00175 (pressure.PresTotlCorrect - pressure.PresTotlCurr)*100./pressure.PresTotlCorrect 00176 ); 00177 } 00178 } 00179 00180 /* check whether pressure is oscillating */ 00181 /* >>chng 02 may 31, add check on sign on hden changes */ 00182 /* >>chng 02 jun 05, add loop > 1 so that don't trigger off old vals that have 00183 * not stabilized yet */ 00184 /* >>chng 04 dec 11, rm test of pressure changing, only check whether hden is changing 00185 * very small changes in hden resulted in changes in pressure with some noise due to 00186 * finite precision in te solver - that makes very small oscillations that are not 00187 * actually oscillations, only noise */ 00188 if( ( /*( pres_chng*pres_chng_old < 0. )||*/ ( hden_chng*hden_chng_old < 0. ) ) && 00189 loop > 1 ) 00190 { 00191 /*fprintf(ioQQQ,"DEBUG\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\n", 00192 fnzone,pres_chng,pres_chng_old ,hden_chng,hden_chng_old);*/ 00193 /* the sign of the change in pressure has changed, so things 00194 * are oscillating. This would be a problem */ 00195 lgPresOscil = true; 00196 nloop_pres_oscil = loop; 00197 /* >>chng 04 dec 09, go to this factor becoming smaller every time oscillation occurs */ 00198 dP_chng_factor = MAX2(0.1 , dP_chng_factor * 0.75 ); 00199 /* dP_chng_factor is how pressure changes with density - pass this to 00200 * changing routine if it is stable */ 00201 /* >>chng 04 dec 09, take average of old and new dP_chng_factor to avoid pressure 00202 * failures in orion_hii_pdr_pp.in 00203 if( loop > 4 && old_slope*rel_slope > 0. ) 00204 dP_chng_factor = (1.-FRACNEW)*dP_chng_factor + FRACNEW*rel_slope;*/ 00205 00206 /*fprintf(ioQQQ,"oscilll %li %.2e %.2e %.2e %.2e dP_chng_factor %.2e\n", 00207 loop , 00208 pres_chng, 00209 pres_chng_old, 00210 hden_chng , 00211 hden_chng_old , 00212 rel_slope);*/ 00213 } 00214 00215 /* convergence trace at this level */ 00216 if( trace.nTrConvg>=1 ) 00217 { 00218 fprintf( ioQQQ, 00219 " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e PCorct:%.4e err:%6.3f%% dP/dn:%.2e Te:%.4e Osc:%c\n", 00220 fnzone, 00221 loop, 00222 dense.gas_phase[ipHYDROGEN], 00223 dense.eden, 00224 pressure.PresTotlCurr, 00225 pressure.PresTotlCorrect, 00226 /* this is percentage error */ 00227 100.*(pressure.PresTotlCurr - pressure.PresTotlCorrect )/pressure.PresTotlCorrect, 00228 dP_chng_factor , 00229 phycon.te, 00230 TorF(lgPresOscil) ); 00231 } 00232 00233 /* increment loop counter */ 00234 ++loop; 00235 00236 /* there can be a pressure or density oscillation early in the search - if not persistent 00237 * ok to clear flag 00238 * >>chng 04 sep 22, add this logic */ 00239 if( loop - nloop_pres_oscil > 4 ) 00240 lgPresOscil = false; 00241 00242 /* if we hit limit of loop, but no oscillations have happened, then we are 00243 * making progress, and can keep going */ 00244 if( loop == LoopMax && !lgPresOscil ) 00245 { 00246 LoopMax = MIN2( 100 , LoopMax*2 ); 00247 } 00248 } 00249 00250 /* >>chng 04 jan 31, now that all of the physics is converged, determine grain drift velocity */ 00251 if( gv.lgDustOn && gv.lgGrainPhysicsOn ) 00252 GrainDrift(); 00253 00254 /* >>chng 01 mar 14, all ConvFail one more time, no matter how 00255 * many failures occurred below. Had been series of if, so multiple 00256 * calls per failure possible. */ 00257 /* >>chng 04 au 07, only announce pres fail here, 00258 * we did not converge the pressure */ 00259 if( !conv.lgConvIoniz ) 00260 ConvFail("ioni",""); 00261 else if( !conv.lgConvEden ) 00262 ConvFail("eden",""); 00263 else if( !conv.lgConvTemp ) 00264 ConvFail("temp",""); 00265 else if( !conv.lgConvPres ) 00266 ConvFail("pres",""); 00267 00268 /* this is only a sanity check that the summed continua accurately reflect 00269 * all of the individual components. Only include this when NDEBUG is not set, 00270 * we are in not debug compile */ 00271 # if !defined(NDEBUG) 00272 RT_OTS_ChkSum(0); 00273 # endif 00274 00275 return 0; 00276 }
1.4.7