/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and * others. For conditions of distribution and use see copyright notice in license.txt */ /*ConvInitSolution drive search for initial temperature, for illuminated face */ /*lgTempTooHigh returns true if temperature is too high */ /*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */ #include "cddefines.h" #include "physconst.h" #include "trace.h" #include "struc.h" #include "rfield.h" #include "mole.h" #include "dense.h" #include "stopcalc.h" #include "heavy.h" #include "wind.h" #include "geometry.h" #include "thermal.h" #include "radius.h" #include "phycon.h" #include "pressure.h" #include "conv.h" #include "hmi.h" #include "dynamics.h" /* derivative of net cooling wrt temperature to check on sign oscillations */ static double dCoolNetDTOld = 0; static double OxyInGrains , FracMoleMax; /* determine whether chemistry is important - what is the largest * fraction of an atom in molecules? Also is ice formation * important * sets above two file static variables */ /*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */ STATIC bool lgCoolHeatCheckConverge( double *CoolNet, bool lgReset ) { static double HeatOld=-1. , CoolOld=-1.; DEBUG_ENTRY( "lgCoolHeatCheckConverge()" ); if( lgReset ) { HeatOld = -1; CoolOld = -1; } double HeatChange = thermal.htot - HeatOld; double CoolChange = thermal.ctot - CoolOld; /* will use larger of heating or cooling in tests for relative convergence * because in constant temperature models one may have trivially small value */ double HeatCoolMax = MAX2( thermal.htot , thermal.ctot ); /* is the heating / cooling converged? * get max relative change, determines whether converged heating/cooling */ double HeatRelChange = fabs(HeatChange)/SDIV( HeatCoolMax ); double CoolRelChange = fabs(CoolChange)/SDIV( HeatCoolMax ); bool lgConverged = true; if( MAX2( HeatRelChange , CoolRelChange ) > conv.HeatCoolRelErrorAllowed ) lgConverged = false; *CoolNet = thermal.ctot - thermal.htot; HeatOld = thermal.htot; CoolOld = thermal.ctot; return lgConverged; } /* call lgCoolHeatCheckConverge until cooling and heating are converged */ STATIC bool lgCoolNetConverge( double *CoolNet , double *dCoolNetDT, bool lgReset ) { const long int LOOP_MAX=20; long int loop = 0; bool lgConvergeCoolHeat = false; DEBUG_ENTRY( "lgCoolNetConverge()" ); while( loop < LOOP_MAX && !lgConvergeCoolHeat && !lgAbort ) { if( ConvEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; } lgConvergeCoolHeat = lgCoolHeatCheckConverge( CoolNet, lgReset && loop==0 ); if( trace.lgTrace || trace.nTrConvg>=3 ) fprintf(ioQQQ," lgCoolNetConverge %li calls to ConvEdenIoniz, converged? %c\n", loop , TorF( lgConvergeCoolHeat ) ); ++loop; } if( thermal.ConstTemp > 0 ) { /* constant temperature model - this is trick so that temperature * updater uses derivative to go to the set constant temperature */ *CoolNet = phycon.te - thermal.ConstTemp; *dCoolNetDT = 1.f; } else if( phycon.te < 4000.f ) { /* at low temperatures the analytical cooling-heating derivative * is often of no value - the species densities change when the * temperature changes. Use simple dCnet/dT=(C-H)/T - this is * usually too large and results is smaller than optical dT, but * we do eventually converge */ *dCoolNetDT = thermal.ctot / phycon.te; } else if( thermal.htot / thermal.ctot > 2.f ) { /* we are far from the solution and the temperature is much lower than * equilibrium - analytical derivative from cooling evaluation is probably * wrong since coolants are not correct */ *dCoolNetDT = thermal.ctot / phycon.te; } else { *dCoolNetDT = thermal.dCooldT - thermal.dHeatdT; if( dCoolNetDTOld * *dCoolNetDT < 0. ) { /* derivative is oscillating */ fprintf(ioQQQ,"PROBLEM CoolNet derivative oscillating in initial solution "\ "Te=%.3e dCoolNetDT=%.3e CoolNet=%.3e dCooldT=%.3e dHeatdT=%.3e\n", phycon.te , *dCoolNetDT , *CoolNet, thermal.dCooldT , thermal.dHeatdT); *dCoolNetDT = *CoolNet / phycon.te; } } return lgConvergeCoolHeat; } /*ChemImportance find fraction of heavy elements in molecules, O in ices */ STATIC void ChemImportance( void ) { DEBUG_ENTRY( "ChemImportance()" ); FracMoleMax = 0.; long int ipMax = -1; for( long i=0; iisMonatomic()) { for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin(); atom != mole_global.list[i]->nAtom.end(); ++atom) { long nelem = atom->first->el->Z-1; if( dense.lgElmtOn[nelem] ) { realnum dens_elemsp = (realnum) mole.species[i].den*atom->second; if ( FracMoleMax*dense.gas_phase[nelem] < dens_elemsp ) { FracMoleMax = dens_elemsp/dense.gas_phase[nelem]; ipMax = i; } } } } } /* fraction of all oxygen in ices on grains */ OxyInGrains = 0; count_ptr elOxygen = unresolved_atom_list[ipOXYGEN]; for(long i=0;ilgGas_Phase && mole_global.list[i]->parentLabel.empty() && mole_global.list[i]->nAtom.find(elOxygen) != mole_global.list[i]->nAtom.end() ) OxyInGrains += mole.species[i].den*mole_global.list[i]->nAtom[elOxygen]; } /* this is now fraction of O in ices */ OxyInGrains /= SDIV(dense.gas_phase[ipOXYGEN]); { /* option to print out entire matrix */ enum {DEBUG_LOC=false}; if( DEBUG_LOC ) { fprintf(ioQQQ,"DEBUG fraction molecular %.2e species %li O ices %.2e\n" , FracMoleMax , ipMax ,OxyInGrains ); } } return; } double FindTempChangeFactor(void) { double TeChangeFactor; DEBUG_ENTRY( "FindTempChangeFactor()" ); /* find fraction of atoms in molecules - need small changes * in temperature if fully molecular since chemistry * network is complex and large changes in solution would * cause linearization to fail */ /* this evaluates the global variables FracMoleMax and * OxyInGrains */ ChemImportance(); /* Following are series of chemical * tests that determine the factor by which * which can change the temperature. must be VERY small when * ice formation on grains is important due to exponential sublimation * rate. Also small when molecules are important, to keep chemistry * within limits of linearized solver * the complete linearization that is implicit in all these solvers * will not work when large Te jumps occur when molecules/ices * are important - solvers can't return to solution if too far away * keep temp steps small when mole/ice is important */ if( OxyInGrains > 0.99 ) { TeChangeFactor = 0.999; } else if( OxyInGrains > 0.9 ) { TeChangeFactor = 0.99; } /* >>chng 97 mar 3, make TeChangeFactor small when close to lower bound */ else if( phycon.te < 5. || /*>>chng 06 jul 30, or if molecules/ices are important */ OxyInGrains > 0.1 ) { TeChangeFactor = 0.97; } /*>>chng 07 feb 23, add test of chemistry being very important */ else if( phycon.te < 20. || FracMoleMax > 0.1 ) { /* >>chng 07 feb 24, changed from 0,9 to 0.95 to converge * pdr_coolbb.in test */ TeChangeFactor = 0.95; } else { /* this is the default largest step */ TeChangeFactor = 0.8; } return TeChangeFactor; } /*ConvInitSolution drive search for initial temperature, for illuminated face */ bool ConvInitSolution() { long int i, ionlup, nelem , nflux_old, nelem_reflux , ion_reflux; /* current value of Cooling - Heating */ bool lgConvergeCoolHeat; double TeChangeFactor, TeBoundLow, TeBoundHigh; DEBUG_ENTRY( "ConvInitSolution()" ); /* set NaN for safety */ set_NaN( TeBoundLow ); set_NaN( TeBoundHigh ); /* this counts number of times ConvBase is called by PressureChange, in current zone */ conv.nPres2Ioniz = 0; /* this counts how many times ConvBase is called in this iteration * and is flag used by various routines to understand they are * being called the first time*/ conv.nTotalIoniz = 0; conv.hist_pres_nzone = -1; conv.hist_temp_nzone = -1; /* ots rates not oscillating */ conv.lgOscilOTS = false; lgAbort = false; dense.lgEdenBad = false; dense.nzEdenBad = 0; /* these are variables to remember the biggest error in the * electron density, and the zone in which it occurred */ conv.BigEdenError = 0.; conv.AverEdenError = 0.; conv.BigHeatCoolError = 0.; conv.AverHeatCoolError = 0.; conv.BigPressError = 0.; conv.AverPressError = 0.; /* these are equal if set dr was used, and we know the first dr */ if( fp_equal( radius.sdrmin, radius.sdrmax ) ) { // lgSdrmaxRel true if sdrmax is relative to current radius, false if limit in cm double rfac = radius.lgSdrmaxRel ? radius.Radius : 1.; radius.drad = MIN2( rfac*radius.sdrmax, radius.drad ); radius.drad_mid_zone = radius.drad/2.; radius.drad_x_fillfac = radius.drad * geometry.FillFac; } /* the H+ density is zero in sims with no H-ionizing radiation and no cosmic rays, * the code does work in this limit. The H0 density can be zero in limit * of very high ionization where H0 underflows to zero. */ ASSERT( dense.xIonDense[ipHYDROGEN][0] >=0 && dense.xIonDense[ipHYDROGEN][1]>= 0.); if( trace.nTrConvg ) fprintf( ioQQQ, "\nConvInitSolution entered \n" ); /******************************************************************** * * this is second or higher iteration, reestablish original temperature * *********************************************************************/ if( iteration != 1 ) { /* this is second or higher iteration on multi-iteration model */ if( trace.lgTrace || trace.nTrConvg ) { fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n", iteration, struc.testr[0] ); } if( trace.lgTrace || trace.nTrConvg ) { fprintf( ioQQQ, " search set true\n" ); } /* search phase must be turned on so that variables such as the ots rates, * secondary ionization, and auger yields, can converge more quickly to * proper values */ conv.lgSearch = true; conv.lgFirstSweepThisZone = true; conv.lgLastSweepThisZone = false; /* reset to the zone one temperature from the previous iteration */ TempChange( struc.testr[0] , false); /* find current pressure - sets pressure.PresTotlCurr */ PresTotCurrent(); /* get new initial temperature and pressure */ if( ConvPresTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } if( trace.lgTrace || trace.nTrConvg ) { fprintf( ioQQQ, " ========================================================================================================================\n"); fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te); fprintf( ioQQQ, " ========================================================================================================================\n"); } conv.lgSearch = false; /* get temperature a second time */ if( ConvTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } /* the initial pressure should now be valid * this sets values of pressure.PresTotlCurr */ PresTotCurrent(); } else { /******************************************************************** * * do first te from scratch * *********************************************************************/ // Set to false to switch to using only standard temperature convergence method const bool lgDoInitConv = true; /* say that we are in search phase */ conv.lgSearch = true; conv.lgFirstSweepThisZone = true; conv.lgLastSweepThisZone = false; if( trace.lgTrace ) { fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" ); } /* coming up to here Te is either 4,000 K (usually) or 10^6 */ TeBoundLow = phycon.TEMP_LIMIT_LOW/1.001; /* set temperature floor option - StopCalc.TeFloor is usually * zero but reset this this command - will go over to constant * temperature calculation if temperature falls below floor */ TeBoundLow = MAX2( TeBoundLow , StopCalc.TeFloor ); /* highest possible temperature - must not get this high since * parts of code will abort if value is reached. * divide by 1.2 to prevent checking on temp > 1e10 */ TeBoundHigh = phycon.TEMP_LIMIT_HIGH/1.2; /* set initial temperature, options are constant temperature, * or approach equilibrium from either high or low TE */ double TeFirst; if( thermal.ConstTemp > 0 ) { /* constant temperature , set 4000 K then relax to desired value * for cold temperatures, to slowly approach fully molecular solution. * if constant temperature is higher than 4000., go right to it */ /* true allow ionization stage trimming, false blocks it */ TeFirst = thermal.ConstTemp ; if (lgDoInitConv) TeFirst = MAX2( 4000. , TeFirst ); /* override this if molecule deliberately disabled */ if( mole_global.lgNoMole ) TeFirst = thermal.ConstTemp; } else if( thermal.lgTeHigh ) { /* approach from high TE */ TeFirst = MIN2( 1e6 , TeBoundHigh ); } else { /* approach from low TE */ TeFirst = MAX2( 4000., TeBoundLow ); } // initial kinetic temperature should be at or above the local // energy density temperature if the radiation field is well // coupled to the matter, but never let this overrule // constant temperature command if( !thermal.lgTemperatureConstant ) TeFirst = max( TeFirst , phycon.TEnerDen ); TempChange(TeFirst , false); if( trace.lgTrace || trace.nTrConvg>=2 ) fprintf(ioQQQ,"ConvInitSolution doing initial solution with temp=%.2e\n", phycon.te); if (lgDoInitConv) { /* this sets values of pressure.PresTotlCurr */ PresTotCurrent(); thermal.htot = 1.; thermal.ctot = 1.; /* call cooling, heating, opacity, loop to convergence * this is very first call to it, by default is at 4000K */ double CoolNet=0, dCoolNetDT=0; /* do ionization twice to get stable solution, evaluating initial dR each time */ lgConvergeCoolHeat = false; for( ionlup=0; ionlup<2; ++ionlup ) { if( trace.lgTrace || trace.nTrConvg>=2 ) fprintf( ioQQQ, "ConvInitSolution calling lgCoolNetConverge " "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n", ionlup,phycon.te , thermal.ctot , thermal.htot,dCoolNetDT ); /* do not flag oscillating d(C-H)/dT until stable */ dCoolNetDTOld = 0.; lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , &dCoolNetDT, true ); if( lgAbort ) break; /* set thickness of first zone, this affects the pumping rates due * to correction for attenuation across zone, so must be known * for ConvEdenIoniz to get valid solution - this is why it * is in this loop */ radius_first(); } if( !lgConvergeCoolHeat ) fprintf(ioQQQ," PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n"); if( lgAbort ) { /* we hit an abort */ fprintf( ioQQQ, " DISASTER PROBLEM ConvInitSolution: Search for " "initial conditions aborted - lgAbort set true.\n" ); ShowMe(); /* this is an error return, calculation will immediately stop */ return lgAbort; } /* we now have error in C-H and its derivative - following is better * derivative for global case where we may be quite far from C==H */ if( thermal.ConstTemp<=0 ) dCoolNetDT = thermal.ctot / phycon.te; bool lgConvergedLoop = false; long int LoopThermal = 0; /* initial temperature is usually 4000K, if the dTe scale factor is 0.95 * it will take 140 integrations to lower temperature to 3K, * and many more if ices are important. */ const long int LIMIT_THERMAL_LOOP=300; double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,CoolSave=-1; while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP ) { /* change temperature until sign of C-H changes */ CoolMHeatSave = CoolNet; dCoolNetDTOld = dCoolNetDT; CoolSave = thermal.ctot; TempSave = phycon.te; /* find temperature change factor, a number less than 1*/ TeChangeFactor = FindTempChangeFactor(); ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. ); TeNew = phycon.te - CoolNet / dCoolNetDT; TeNew = MAX2( phycon.te*TeChangeFactor , TeNew ); TeNew = MIN2( phycon.te/TeChangeFactor , TeNew ); /* change temperature */ TempChange(TeNew , true); lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , & dCoolNetDT, false ); if( trace.lgTrace || trace.nTrConvg>=2 ) fprintf(ioQQQ,"ConvInitSolution %4li TeNnew=%.3e " "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e" " Ograins %.2e FracMoleMax %.2e\n", LoopThermal , TeNew , TeChangeFactor , CoolNet/SDIV(thermal.htot) , dCoolNetDT, dense.eden , OxyInGrains , FracMoleMax ); if( lgAbort ) return lgAbort; /* keep changing temperature until sign of CoolNet changes * in constant temperature case CoolNet is delta Temperature * so is zero for last pass in this loop * this is for constant temperature case */ if( fabs(CoolNet)< SMALLDOUBLE ) /* CoolNet is zero */ lgConvergedLoop = true; else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0) /* change in sign of net cooling */ lgConvergedLoop = true; else if( phycon.te <= TeBoundLow || phycon.te>=TeBoundHigh) lgConvergedLoop = true; ++LoopThermal; } if( !lgConvergedLoop ) fprintf(ioQQQ,"PROBLEM ConvInitSolution: temperature solution not " "found in initial search, final Te=%.2e\n", phycon.te ); /* interpolate temperature where C=H if not constant temperature sim * will have set constant temperature mode above if we hit temperature * floor */ if( thermal.ConstTemp<=0 && ! fp_equal( TempSave , phycon.te ) ) { if( trace.lgTrace || trace.nTrConvg>=2 ) fprintf(ioQQQ," Change in sign of C-H found, Te brackets %.3e " "%.3e Cool brackets %.3e %.3e ", TempSave , phycon.te , CoolMHeatSave , CoolNet); /* interpolate new temperature assuming Cool = a T^alpha, * first find alpha */ double alpha = log(thermal.ctot/CoolSave) / log( phycon.te/TempSave); if( fabs(alpha) < SMALLFLOAT ) /* alpha close to 0 means constant temperature */ TeNew = phycon.te; else { /* next find log a - work in logs to avoid infinities */ if( thermal.ctot<0. || thermal.htot<0. ) TotalInsanity(); double Alog = log( thermal.ctot ) - alpha * log( phycon.te ); /* the interpolated temperature where heating equals cooling */ double TeLn = (log( thermal.htot ) - Alog) / alpha ; /* a sanity check */ if( TeLn < log( MIN2(phycon.te , TempSave )) ) TeNew = MIN2(phycon.te , TempSave ); else if( TeLn > log( MAX2(phycon.te , TempSave )) ) TeNew = MAX2(phycon.te , TempSave ); else TeNew = pow( EE , TeLn ); } ASSERT( TeNew>=MIN2 ( TempSave , phycon.te ) || TeNew<=MAX2 ( TempSave , phycon.te )); if( trace.lgTrace || trace.nTrConvg>=2 ) fprintf(ioQQQ," interpolating to Te %.3e \n\n", TeNew); /* change temperature */ TempChange(TeNew , true); } } if( ConvTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } conv.lgSearch = false; // Solve again without search phase lo-fi physics before starting on // anything which might have a long-term impact if( ConvTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } /* this sets values of pressure.PresTotlCurr */ PresTotCurrent(); if( trace.lgTrace || trace.nTrConvg ) { fprintf( ioQQQ, "\nConvInitSolution: end 1st iteration search phase " "finding Te=%.3e\n\n" , phycon.te); } if( trace.lgTrace ) { fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n", phycon.te ); } } /* we now have a fairly good temperature and ionization * iteration is 1 on first iteration - remember current pressure * on first iteration so we can hold this constant in constant * pressure simulation. * * flag dense.lgDenseInitConstant false if * *constant pressure reset* is used - * default is true, after first iteration initial density is used for * first zone no matter what pressure results. Small changes occur due * to radiative transfer convergence, * when set false with *reset* option the density on later iterations * can change to keep the pressure constant */ if( iteration==1 || dense.lgDenseInitConstant ) { double PresNew = pressure.PresTotlCurr; /* option to specify the pressure as option on constant pressure command */ if( pressure.lgPressureInitialSpecified ) /* this is log of nT product - if not present then set zero */ PresNew = pressure.PressureInitialSpecified; if( trace.lgTrace ) { fprintf( ioQQQ, " PresTotCurrent 1st zn old values of PresTotlInit:%.3e" " to %.3e \n", pressure.PresTotlInit, PresNew); } pressure.PresTotlInit = PresNew; } if( dynamics.lgTimeDependentStatic ) { // some sort of time dependent sim static double PresTotalInitTime; if( iteration <= dynamics.n_initial_relax ) { PresTotalInitTime = pressure.PresTotlInit; } else if (dense.lgPressureVaryTime) { pressure.PresTotlInit = PresTotalInitTime * pow( 1.+(dynamics.time_elapsed/dense.PressureVaryTimeTimescale) , dense.PressureVaryTimeIndex); } fprintf(ioQQQ,"DEBUG conv_init_solution time dependent time=%.2e sets " "pressure to %.2e\n", dynamics.time_elapsed , pressure.PresTotlInit); } /* now bring current pressure to the correct pressure - must do this * before initial values are saved in iter start/end */ ConvPresTempEdenIoniz(); /* this counts number of times ConvBase is called by PressureChange, in * current zone these are reset here, so that we count from first zone * not search */ conv.nPres2Ioniz = 0; dense.lgEdenBad = false; dense.nzEdenBad = 0; /* save counter upon exit so we can see how many iterations were * needed to do true zones */ conv.nTotalIoniz_start = conv.nTotalIoniz; /* these are variables to remember the biggest error in the * electron density, and the zone in which it occurred */ conv.BigEdenError = 0.; conv.AverEdenError = 0.; conv.BigHeatCoolError = 0.; conv.AverHeatCoolError = 0.; conv.BigPressError = 0.; conv.AverPressError = 0.; /*>>chng 06 jun 09, only reset on first try - otherwise have stable solution? */ if( iteration == 1 ) { /* now remember some things we may need even in first zone, these * are normally set towards end of zone calculation in RT_tau_inc */ struc.testr[0] = (realnum)phycon.te; /* >>chng 02 May 2001 rjrw: add hden for dilution */ struc.hden[0] = dense.gas_phase[ipHYDROGEN]; /* pden is the total number of particles per unit vol */ struc.DenParticles[0] = dense.pden; struc.heatstr[0] = thermal.htot; struc.coolstr[0] = thermal.ctot; struc.volstr[0] = (realnum)radius.dVeffAper; struc.drad_x_fillfac[0] = (realnum)radius.drad_x_fillfac; struc.histr[0] = dense.xIonDense[ipHYDROGEN][0]; struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1]; struc.ednstr[0] = (realnum)dense.eden; struc.o3str[0] = dense.xIonDense[ipOXYGEN][2]; struc.DenMass[0] = dense.xMassDensity; struc.drad[0] = (realnum)radius.drad; } /* check that nflux extends above IP of highest ionization species present. * for collisional case it is possible for species to exist that are higher * IP than the limit to the continuum. Need continuum to encompass all * possible emission - to account for diffuse emission * NB * on second iteration of multi-iteration model this may result in rfield.nflux increasing * which can change the solution */ nflux_old = rfield.nflux; nelem_reflux = -1; ion_reflux = -1; for( nelem=2; nelem < LIMELM; nelem++ ) { /* do not include hydrogenic species in following */ for( i=0; i < nelem; i++ ) { if( dense.xIonDense[nelem][i+1] > 0. ) { if( Heavy.ipHeavy[nelem][i] > rfield.nflux ) { rfield.nflux = Heavy.ipHeavy[nelem][i]; nelem_reflux = nelem; ion_reflux = i; } } } } /* was the upper limit to the continuum updated? if so need to define * continuum variables to this new range */ if( nflux_old != rfield.nflux ) { /* zero out parts of rfield arrays that were previously undefined */ rfield_opac_zero( nflux_old-1 , rfield.nflux ); /* if continuum reset up, we need to define gaunt factors through high end */ /*tfidle(false);*/ /* this calls tfidle, among other things */ /* this sets values of pressure.PresTotlCurr */ PresTotCurrent(); /* redo ionization and update opacities */ if( ConvBase(1) ) { /* this is catastrophic failure */ lgAbort = true; return lgAbort; } /* need to update continuum opacities */ if( trace.lgTrace ) { fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n", nflux_old , rfield.nflux , rfield.anu[nflux_old] , rfield.anu[rfield.nflux] ); fprintf(ioQQQ," caused by element %li ion %li \n", nelem_reflux ,ion_reflux ); } } /* dynamics solution - change density slightly * and call pressure solver to see if it returns to original density */ if( strcmp(dense.chDenseLaw,"DYNA") == 0 ) { /* >>chng 04 apr 23, change pressure and let solver come back to correct * pressure. This trick makes sure we are correctly either super or sub sonic * since solver will jump from one branch to the other if necessary */ static const double PCHNG = 0.98; /* this ignores return condition, assume must be ok if got this far */ if( ConvPresTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } pressure.PresTotlInit *= PCHNG; if( ConvPresTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } pressure.PresTotlInit /= PCHNG; if( ConvPresTempEdenIoniz() ) { /* this is an error return, calculation will immediately stop */ lgAbort = true; return lgAbort; } # undef PCHNG } /* this is the only valid return and lgAbort should be false*/ return lgAbort; }