00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007
00008 #include "phycon.h"
00009 #include "rt.h"
00010 #include "dense.h"
00011 #include "pressure.h"
00012 #include "trace.h"
00013 #include "conv.h"
00014 #include "pressure_change.h"
00015 #include "thermal.h"
00016 #include "geometry.h"
00017 #include "grainvar.h"
00018 #include "grains.h"
00019
00020
00021
00022
00023 int ConvPresTempEdenIoniz(void)
00024 {
00025 DEBUG_ENTRY( "ConvPresTempEdenIoniz()" );
00026
00027
00028
00029
00030 conv.nPres2Ioniz = 0;
00031 conv.lgFirstSweepThisZone = true;
00032 conv.lgLastSweepThisZone = false;
00033
00034
00035 if( nzone != conv.hist_pres_nzone )
00036 {
00037
00038 conv.hist_pres_nzone = nzone;
00039 conv.hist_pres_density.clear();
00040 conv.hist_pres_current.clear();
00041 conv.hist_pres_error.clear();
00042 }
00043
00044
00045 conv.lgConvPops = true;
00046
00047
00048
00049 if( trace.nTrConvg>=1 )
00050 {
00051 fprintf( ioQQQ,
00052 " ConvPresTempEdenIoniz1 entered, will call ConvIoniz to initialize\n");
00053 }
00054
00055
00056
00057
00058
00059
00060
00061 conv.lgConvPres = false;
00062
00063
00064 if( trace.nTrConvg>=1 )
00065 {
00066 fprintf( ioQQQ,
00067 " ConvPresTempEdenIoniz1 ConvIoniz found following converged: Pres;%c, Temp;%c, Eden;%c, Ion:%c, Pops:%c\n",
00068 TorF(conv.lgConvPres),
00069 TorF(conv.lgConvTemp),
00070 TorF(conv.lgConvEden),
00071 TorF(conv.lgConvIoniz()),
00072 TorF(conv.lgConvPops));
00073 }
00074
00075
00076 if( trace.nTrConvg>=1 )
00077 {
00078 fprintf( ioQQQ,
00079 "\n ConvPresTempEdenIoniz1 entering main pressure loop.\n");
00080 }
00081
00082 AbundChange();
00083
00084 if ( strcmp(dense.chDenseLaw,"CPRE") == 0 ||
00085 strcmp(dense.chDenseLaw,"DYNA") == 0 )
00086 {
00087
00088
00089 double TemperatureInitial = phycon.te;
00090
00091
00092
00093 if( ConvTempEdenIoniz() )
00094 {
00095 return 1;
00096 }
00097
00098 PresMode presmode;
00099 presmode.set();
00100 solverState st;
00101 st.press = pressureZone(presmode);
00102
00103
00104 bool lgPresOscil = false;
00105
00106 long nloop_pres_oscil = 0;
00107
00108 double hden_old = dense.gas_phase[ipHYDROGEN];
00109 double hden_chng = 0.;
00110
00111 double dP_chng_factor = 1.;
00112
00113
00114 const int LOOPMAX = 100;
00115
00116 long LoopMax = LOOPMAX;
00117 long loop = 0;
00118
00119
00120
00121 if( nzone <= 1 && pressure.lgPressureInitialSpecified )
00122 LoopMax = 10*LOOPMAX;
00123
00124 while( (loop < LoopMax) && !conv.lgConvPres && !lgAbort )
00125 {
00126
00127
00128
00129 if( fabs( TemperatureInitial - phycon.te )/phycon.te > 0.3 )
00130 LoopMax = 2*LOOPMAX;
00131
00132
00133
00134
00135 hden_old = dense.gas_phase[ipHYDROGEN];
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 conv.incrementCounter(PRES_CHANGES);
00146 bool lgChange = PressureChange(dP_chng_factor, presmode, st);
00147 if( lgChange )
00148 {
00149
00150
00151
00152 int lgTempStatus = ConvTempEdenIoniz();
00153 if( lgTempStatus != 0 )
00154 {
00155 return 1;
00156 }
00157 }
00158
00159
00160 double hden_chng_old = hden_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
00165 {
00166
00167 enum{DEBUG_LOC=false};
00168
00169 if( DEBUG_LOC && nzone > 150 && iteration > 1 )
00170 {
00171 fprintf(ioQQQ,"%li\t%.2e\t%.2e\n",
00172 nzone,
00173 pressure.PresTotlCurr,
00174 pressure.PresTotlError*100.
00175 );
00176 }
00177 }
00178
00179
00180 if( ( ( hden_chng*hden_chng_old < 0. ) ) && loop > 1 )
00181 {
00182
00183
00184 lgPresOscil = true;
00185 nloop_pres_oscil = loop;
00186
00187 dP_chng_factor = MAX2(0.1 , dP_chng_factor * 0.75 );
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 }
00199
00200
00201 if( trace.nTrConvg>=1 )
00202 {
00203 fprintf( ioQQQ,
00204 " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e err:%6.3f%% dP/dn:%.2e Te:%.4e Osc:%c\n",
00205 fnzone,
00206 loop,
00207 dense.gas_phase[ipHYDROGEN],
00208 dense.eden,
00209 pressure.PresTotlCurr,
00210
00211 100.*pressure.PresTotlError,
00212 dP_chng_factor ,
00213 phycon.te,
00214 TorF(lgPresOscil) );
00215 }
00216
00217
00218 ++loop;
00219
00220
00221
00222 if( loop - nloop_pres_oscil > 4 )
00223 lgPresOscil = false;
00224
00225
00226
00227 if( loop == LoopMax && !lgPresOscil )
00228 {
00229 LoopMax = MIN2( LOOPMAX, LoopMax*2 );
00230 }
00231 }
00232 }
00233 else
00234 {
00235 double targetDensity = zoneDensity();
00236 double startingDensity = scalingDensity();
00237 double pdelta = conv.MaxFractionalDensityStepPerIteration;
00238
00239 double logRatio = log(targetDensity/startingDensity);
00240 long nstep = (long) ceil(fabs(logRatio)/pdelta);
00241
00242 if (nstep == 0)
00243 nstep = 1;
00244 double density_change_factor = exp(logRatio/nstep);
00245
00246 for (long i=0; i<nstep; i++)
00247 {
00248 if (i == nstep - 1)
00249 {
00250
00251 density_change_factor = targetDensity/scalingDensity();
00252 }
00253
00254
00255
00256
00257
00258
00259
00260 conv.hist_pres_density.push_back( dense.gas_phase[ipHYDROGEN] );
00261 conv.hist_pres_current.push_back( pressure.PresTotlCurr );
00262 conv.hist_pres_error.push_back( pressure.PresTotlError );
00263
00264 if( trace.lgTrace )
00265 {
00266 fprintf( ioQQQ,
00267 " DensityUpdate called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
00268 dense.gas_phase[ipHYDROGEN], density_change_factor*dense.gas_phase[ipHYDROGEN],
00269 geometry.FillFac );
00270 }
00271
00272 ScaleAllDensities((realnum) density_change_factor);
00273
00274
00275
00276 TempChange(phycon.te , false);
00277
00278
00279
00280
00281 if( ConvTempEdenIoniz() )
00282 {
00283 return 1;
00284 }
00285
00286
00287 if( trace.nTrConvg>=1 )
00288 {
00289 fprintf( ioQQQ,
00290 " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e err:%6.3f%% Te:%.4e Osc:%c\n",
00291 fnzone,
00292 i,
00293 dense.gas_phase[ipHYDROGEN],
00294 dense.eden,
00295 pressure.PresTotlCurr,
00296
00297 100.*pressure.PresTotlError,
00298 phycon.te,
00299 TorF(false) );
00300 }
00301 }
00302
00303 PresTotCurrent();
00304
00305 conv.lgConvPres = true;
00306 }
00307
00308 thermal.thist = max((realnum)phycon.te,thermal.thist);
00309 thermal.tlowst = min((realnum)phycon.te,thermal.tlowst);
00310
00311
00312 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
00313 GrainDrift();
00314
00315
00316
00317
00318
00319
00320 if( !conv.lgConvIoniz() )
00321 ConvFail("ioni","");
00322 else if( !conv.lgConvEden )
00323 ConvFail("eden","");
00324 else if( !conv.lgConvTemp )
00325 ConvFail("temp","");
00326 else if( !conv.lgConvPres )
00327 ConvFail("pres","");
00328
00329
00330
00331
00332 # if !defined(NDEBUG)
00333 RT_OTS_ChkSum(0);
00334 # endif
00335
00336 return 0;
00337 }