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