conv_pres_temp_eden_ioniz.cpp

Go to the documentation of this file.

/* 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 */
/*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIonize,
 * called by cloudy */
/*ConvFail handle convergence failure */
#include "cddefines.h"
 
#include "phycon.h"
#include "rt.h"
#include "dense.h"
#include "pressure.h"
#include "trace.h"
#include "conv.h"
#include "pressure_change.h"
#include "thermal.h"
#include "geometry.h" 
#include "grainvar.h"
#include "grains.h"
 
/*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIoniz,
 * called by cloudy 
 * returns 0 if ok, 1 if disaster */
int ConvPresTempEdenIoniz(void)
{
        DEBUG_ENTRY( "ConvPresTempEdenIoniz()" );
 
        /* this will count number of times we call ConvBase in this zone,
         * counter is incremented there 
         * zero indicates first pass through solvers on this zone */
        conv.nPres2Ioniz = 0;
        conv.lgFirstSweepThisZone = true;
        conv.lgLastSweepThisZone = false;
        /* this will save the  history of density - pressure relationship
         * for the current zone */
        if( nzone != conv.hist_pres_nzone )
        {
                /* first time in this zone - reset history */
                conv.hist_pres_nzone = nzone;
                conv.hist_pres_density.clear();
                conv.hist_pres_current.clear();
                conv.hist_pres_error.clear();
        }
 
        /* should still be true at end */
        conv.lgConvPops = true;
 
        /*rel_slope = 0.;*/
 
        if( trace.nTrConvg>=1  )
        {
                fprintf( ioQQQ, 
                        " ConvPresTempEdenIoniz1 entered, will call ConvIoniz to initialize\n");
        }
 
        /* converge the ionization first, so that we know where we are, and have
         * a valid foundation to begin the search */
        /* the true electron density dense.EdenTrue is set in eden_sum called by ConvBase */
 
        /* this evaluates current pressure, and returns whether or not 
         * it is within tolerance of correct pressure */
        conv.lgConvPres = false; 
 
        /* convergence trace at this level */
        if( trace.nTrConvg>=1  )
        {
                fprintf( ioQQQ, 
                        " ConvPresTempEdenIoniz1 ConvIoniz found following converged: Pres;%c, Temp;%c, Eden;%c, Ion:%c, Pops:%c\n", 
                        TorF(conv.lgConvPres), 
                        TorF(conv.lgConvTemp),
                        TorF(conv.lgConvEden),
                        TorF(conv.lgConvIoniz()),
                        TorF(conv.lgConvPops));
        }
 
        /* trace convergence print at this level */
        if( trace.nTrConvg>=1  )
        {
                fprintf( ioQQQ, 
                        "\n ConvPresTempEdenIoniz1 entering main pressure loop.\n");
        }
 
        AbundChange();
 
        if ( strcmp(dense.chDenseLaw,"CPRE") == 0 || 
                  strcmp(dense.chDenseLaw,"DYNA") == 0 ) 
        {
                /* set the initial temperature to the current value, so we will know
                 * if we are trying to jump over a thermal front */
                double TemperatureInitial = phycon.te;
 
                /* chng 02 dec 11 rjrw -- ConvIoniz => ConvTempEdenIoniz() here for consistency with inside loop */
                /* ConvIoniz; */
                if( ConvTempEdenIoniz() )
                {
                        return 1;
                }
 
                PresMode presmode;
                presmode.set();
                solverState st;
                st.press = pressureZone(presmode);
 
                /* this will be flag to check for pressure oscillations */
                bool lgPresOscil = false;
                /* this is loop where it happened */
                long nloop_pres_oscil = 0;
                /* we will use these to check whether hden oscillating - would need to decrease step size */
                double hden_old = dense.gas_phase[ipHYDROGEN];
                double hden_chng = 0.;
                /* this is dP_chng_factor, cut in half when pressure error changes sign */
                double dP_chng_factor = 1.;
 
                /* the limit to the number of loops */
                const int LOOPMAX = 100;
                /* this will be the limit, which we will increase if no oscillations occur */
                long LoopMax = LOOPMAX;
                long loop = 0;
 
                /* if start of calculation and we are solving for set pressure,
                 * allow a lot more iterations */
                if( nzone <= 1 && pressure.lgPressureInitialSpecified )
                        LoopMax = 10*LOOPMAX;
                
                while( (loop < LoopMax) && !conv.lgConvPres &&  !lgAbort )
                {
                        /* there can be a pressure or density oscillation early in the search - if not persistent
                         * ok to clear flag */
                        /* >>chng 01 aug 24, if large change in temperature allow lots more loops */
                        if( fabs( TemperatureInitial - phycon.te )/phycon.te > 0.3 )
                                LoopMax = 2*LOOPMAX;
                        
                        /* change current densities of all constituents if necessary, 
                         * PressureChange evaluates lgPresOK, true if pressure is now ok
                         * sets CurrentPressure and CorrectPressure */
                        hden_old = dense.gas_phase[ipHYDROGEN];
                        /*pres_old = pressure.PresTotlCurr;*/
                        
                        /* this will evaluate current pressure, update the densities, 
                         * determine the wind velocity, and set conv.lgConvPres,
                         * return value is true if density was changed, false if no changes were necessary 
                         * if density changes then we must redo the temperature and ionization 
                         * PressureChange contains the logic that determines how to change the density to get
                         * the right pressure */
 
                        conv.incrementCounter(PRES_CHANGES);                    
                        bool lgChange = PressureChange(dP_chng_factor, presmode, st);
                        if( lgChange )
                        {
                                /* heating cooling balance while doing ionization,
                                 * this is where the heavy lifting is done, this calls PresTotCurrent,
                                 * which sets pressure.PresTotlCurr */
                                int lgTempStatus = ConvTempEdenIoniz();
                                if( lgTempStatus != 0 )
                                {
                                        return 1;
                                }
                        }
                        
                        /* if product of these two is negative then hden is oscillating */
                        double hden_chng_old = hden_chng;
                        hden_chng = dense.gas_phase[ipHYDROGEN] - hden_old;
                        if( fabs(hden_chng) < SMALLFLOAT )
                                hden_chng = sign( (double)SMALLFLOAT, hden_chng );
                        
                        {
                                /*@-redef@*/
                                enum{DEBUG_LOC=false};
                                /*@+redef@*/
                                if( DEBUG_LOC && nzone > 150 && iteration > 1 )
                                {
                                        fprintf(ioQQQ,"%li\t%.2e\t%.2e\n", 
                                                          nzone,
                                                          pressure.PresTotlCurr, 
                                                          pressure.PresTotlError*100.
                                                );
                                }
                        }
                        
                        /* check whether pressure is oscillating */
                        if( ( ( hden_chng*hden_chng_old < 0. ) ) && loop > 1 )
                        {
                                /*fprintf(ioQQQ,"DEBUG\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                  fnzone,pres_chng,pres_chng_old ,hden_chng,hden_chng_old);*/
                                lgPresOscil = true;
                                nloop_pres_oscil = loop;
                                /* >>chng 04 dec 09, go to this factor becoming smaller every time oscillation occurs */
                                dP_chng_factor = MAX2(0.1 , dP_chng_factor * 0.75 );
                                /* dP_chng_factor is how pressure changes with density - pass this to
                                 * changing routine if it is stable */
                                
                                /*fprintf(ioQQQ,"oscilll %li %.2e %.2e %.2e %.2e dP_chng_factor %.2e\n", 
                                  loop ,
                                  pres_chng, 
                                  pres_chng_old,
                                  hden_chng , 
                                  hden_chng_old ,
                                  rel_slope);*/
                        }
                        
                        /* convergence trace at this level */
                        if( trace.nTrConvg>=1  )
                        {
                                fprintf( ioQQQ, 
                                                        " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e err:%6.3f%% dP/dn:%.2e Te:%.4e Osc:%c\n", 
                                                        fnzone,
                                                        loop, 
                                                        dense.gas_phase[ipHYDROGEN], 
                                                        dense.eden,
                                                        pressure.PresTotlCurr, 
                                                        /* this is percentage error */
                                                        100.*pressure.PresTotlError,
                                                        dP_chng_factor ,
                                                        phycon.te,
                                                        TorF(lgPresOscil)  );
                        }
                        
                        /* increment loop counter */
                        ++loop;
                        
                        /* there can be a pressure or density oscillation early in the search - if not persistent
                         * ok to clear flag */
                        if( loop - nloop_pres_oscil > 4 )
                                lgPresOscil = false;
                        
                        /* if we hit limit of loop, but no oscillations have happened, then we are
                         * making progress, and can keep going */
                        if( loop == LoopMax && !lgPresOscil )
                        {
                                LoopMax = MIN2( LOOPMAX, LoopMax*2 );
                        }
                }
        }
        else
        {
                double targetDensity = zoneDensity();
                double startingDensity = scalingDensity();
                double pdelta = conv.MaxFractionalDensityStepPerIteration;
 
                double logRatio = log(targetDensity/startingDensity);
                long nstep = (long) ceil(fabs(logRatio)/pdelta);
                // Ensure at least one update per zone
                if (nstep == 0)
                        nstep = 1;
                double density_change_factor = exp(logRatio/nstep);
 
                for (long i=0; i<nstep; i++)
                {
                        if (i == nstep - 1)
                        {
                                // Try to ensure exact answer at last iteration
                                density_change_factor = targetDensity/scalingDensity();
                        }
 
                        /* this will evaluate current pressure, update the densities, 
                         * determine the wind velocity, and set conv.lgConvPres,
                         * return value is true if density was changed, false if no changes were necessary 
                         * if density changes then we must redo the temperature and ionization  */
 
                        /* >>chng 04 feb 11, add option to remember current density and pressure */
                        conv.hist_pres_density.push_back( dense.gas_phase[ipHYDROGEN] );
                        conv.hist_pres_current.push_back( pressure.PresTotlCurr );
                        conv.hist_pres_error.push_back( pressure.PresTotlError );
                        
                        if( trace.lgTrace )
                        {
                                fprintf( ioQQQ, 
                                                        " DensityUpdate called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n", 
                                                        dense.gas_phase[ipHYDROGEN], density_change_factor*dense.gas_phase[ipHYDROGEN], 
                                                        geometry.FillFac );
                        }
                        
                        ScaleAllDensities((realnum) density_change_factor);
                                
                        /* must call TempChange since ionization has changed, there are some
                         * terms that affect collision rates (H0 term in electron collision) */
                        TempChange(phycon.te , false);
                        
                        /* heating cooling balance while doing ionization, this is
                         * where the heavy lifting is done, this calls
                         * PresTotCurrent, which sets pressure.PresTotlCurr */
                        if( ConvTempEdenIoniz() )
                        {
                                return 1;
                        }
                        
                        /* convergence trace at this level */
                        if( trace.nTrConvg>=1  )
                        {
                                fprintf( ioQQQ, 
                                                        " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e err:%6.3f%% Te:%.4e Osc:%c\n", 
                                                        fnzone,
                                                        i, 
                                                        dense.gas_phase[ipHYDROGEN], 
                                                        dense.eden,
                                                        pressure.PresTotlCurr, 
                                                        /* this is percentage error */
                                                        100.*pressure.PresTotlError,
                                                        phycon.te,
                                                        TorF(false)  );
                        }
                }
 
                PresTotCurrent();
                        
                conv.lgConvPres = true;
        }
        /* keep track of minimum and maximum temperature */
        thermal.thist = max((realnum)phycon.te,thermal.thist);
        thermal.tlowst = min((realnum)phycon.te,thermal.tlowst);
 
        /* >>chng 04 jan 31, now that all of the physics is converged, determine grain drift velocity */
        if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
                GrainDrift();
 
        /* >>chng 01 mar 14, all ConvFail one more time, no matter how
         * many failures occurred below.  Had been series of if, so multiple
         * calls per failure possible. */
        /* >>chng 04 au 07, only announce pres fail here,
         * we did not converge the pressure */
        if( !conv.lgConvIoniz() )
                ConvFail("ioni","");
        else if( !conv.lgConvEden )
                ConvFail("eden","");
        else if( !conv.lgConvTemp )
                ConvFail("temp","");
        else if( !conv.lgConvPres )
                ConvFail("pres","");
 
        /* this is only a sanity check that the summed continua accurately
         * reflect all of the individual components.  Only include this
         * when NDEBUG is not set, we are in not debug compile */
#       if !defined(NDEBUG)
        RT_OTS_ChkSum(0);
#       endif
 
        return 0;
}