dense.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 */
#include "cddefines.h"

#include "dense.h"

#include "abund.h"
#include "colden.h"
#include "conv.h"
#include "dynamics.h"
#include "elementnames.h"
#include "deuterium.h"
#include "hmi.h"
#include "h2.h"
#include "mole.h"
#include "mole_priv.h"
#include "phycon.h"
#include "prt.h"
#include "radius.h"
#include "struc.h"
#include "thermal.h"
#include "trace.h"

t_dense dense;

void t_dense::updateXMolecules()
{
        total_molecule_elems(m_xMolecules);
}

void t_dense::zero()
{
        for (long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem)
                m_xMolecules[nelem] = 0.f;
}

void ScaleAllDensities(realnum factor)
{
        double edensave = dense.eden;

        for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
        {
                if (dense.lgElmtOn[nelem])
                {
                        ScaleIonDensities( nelem, factor );
                        dense.SetGasPhaseDensity( nelem, dense.gas_phase[nelem] * factor);
                }
        }

        EdenChange( dense.eden * factor );
        
        if( trace.lgTrace && trace.lgNeBug )
        {
                fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n", 
                                        edensave, dense.eden, edensave/dense.eden );
        }

        hmi.H2_total *= factor;
        h2.ortho_density *= factor;
        h2.para_density *= factor;
        for( long mol=0; mol < mole_global.num_total; mol++ )
        {
                mole.species[mol].den *= factor;
        }
        dense.updateXMolecules();

        ASSERT(lgElemsConserved());
}

void ScaleIonDensities( const long nelem, const realnum factor )
{
        double renorm;
        for( long ion=0; ion < (nelem + 2); ion++ )
        {
                dense.xIonDense[nelem][ion] *= factor;
                if (nelem-ion >= 0 && nelem-ion < NISO)
                        iso_renorm(nelem, nelem-ion, renorm);
        }
        
        if( nelem==ipHYDROGEN && deut.lgElmtOn )
                ScaleDensitiesDeuterium( factor );

        return;
}

void t_dense::SetGasPhaseDensity( const long nelem, const realnum density )
{
        gas_phase[nelem] = density;

        if( nelem==ipHYDROGEN && deut.lgElmtOn )
        {
                // fractionation is taken care of inside called routine
                SetGasPhaseDeuterium( density );
        }

        return;
}

bool lgElemsConserved (void)
{
        bool lgOK=true;

        /* this checks all molecules and ions add up to abundance */
        for( ChemAtomList::iterator atom = atom_list.begin(); atom != atom_list.end(); ++atom )
        {
                long nelem = (*atom)->el->Z-1;
                /* check that species add up if element is turned on */
                if( dense.lgElmtOn[nelem] )
                {
                        /* this sum of over the molecules did not include the atom and first
                         * ion that is calculated in the molecular solver */
                        double sum_monatomic = 0.;
                        for( long ion=0; ion<nelem+2; ++ion )
                        {
                                sum_monatomic += dense.xIonDense[nelem][ion] * (*atom)->frac;
                        }
                        realnum sum_molecules = dense.xMolecules(nelem) * (*atom)->frac;
                        realnum sum_gas_phase = dense.gas_phase[nelem] * (*atom)->frac;
                        if( sum_monatomic + sum_molecules <= SMALLFLOAT && 
                                sum_gas_phase > SMALLFLOAT)
                        {
                                fprintf(ioQQQ,"PROBLEM non-conservation of nuclei %s\tions %g moles %g error %g of %g\n",
                                        (*atom)->label().c_str(),
                                        sum_monatomic,
                                        sum_molecules,
                                        sum_monatomic + sum_molecules-sum_gas_phase,
                                        sum_gas_phase );
                                lgOK=false;
                        }

                        /* check that sum agrees with total gas phase abundance */
                        if( fabs( sum_monatomic + sum_molecules - sum_gas_phase ) > 
                                conv.GasPhaseAbundErrorAllowed * sum_gas_phase )
                        {
                                fprintf(ioQQQ,"PROBLEM non-conservation of nuclei %s\t nzone %li atoms %.12e moles %.12e sum %.12e tot gas %.12e rel err %.3e\n",
                                        (*atom)->label().c_str(),
                                        nzone,
                                        sum_monatomic,
                                        sum_molecules,
                                        sum_monatomic + sum_molecules,
                                        sum_gas_phase,
                                        (sum_monatomic + sum_molecules - sum_gas_phase)/sum_gas_phase );
                                lgOK=false;
                        }
                
                        if( 0 )
                        {
                                // print reactions involving the neutral 
                                mole_print_species_reactions( findspecies( elementnames.chElementSym[nelem] ) );
                        }
                }
        }

        return lgOK;
}

void lgStatesConserved ( long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion )
{
        if( 0 && conv.lgConvIoniz() && 
                dense.lgElmtOn[nelem] && 
                dense.xIonDense[nelem][ionStage]/dense.eden > conv.EdenErrorAllowed/10. &&
                dense.xIonDense[nelem][ionStage] > SMALLFLOAT )
        {
                double abund = 0.;
                for (long ipLev=0; ipLev<numStates; ipLev++)
                {
                        abund += states[ipLev].Pop();
                }

                double rel_err = fabs(abund-dense.xIonDense[nelem][ionStage])/
                        SDIV(dense.xIonDense[nelem][ionStage]);

                //fprintf(ioQQQ,"Iso consistency error %ld %ld %g\n",nelem,ionStage,rel_err);

                if( rel_err > err_tol ) 
                {
                        /* we'll set the flag for non-convergence, but don't bother complaining in printout for first sweeps */
                        if( 0 && nzone > 0 && loop_ion > 0 )
                        {       
                                fprintf(ioQQQ,"PROBLEM Inconsistent states/stage pops nzone %3ld loop_ion %2ld nelem %2ld ion %2ld states = %e stage = %e error = %e\n",
                                        nzone,
                                        loop_ion,
                                        nelem,
                                        ionStage,
                                        abund,
                                        dense.xIonDense[nelem][ionStage],
                                        rel_err );
                        }
                        char chConvIoniz[INPUT_LINE_LENGTH];
                        sprintf( chConvIoniz , "States!=stage pops nelem %ld ion %ld ", nelem, ionStage );
                        conv.setConvIonizFail( chConvIoniz, abund,
                                                                                 dense.xIonDense[nelem][ionStage] );
                }
        }
}

void SumDensities( void )
{
        realnum den_mole,
                DenseAtomsIons;

        /* find total number of atoms and ions */
        DenseAtomsIons = 0.;
        for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
        {
                if( dense.lgElmtOn[nelem]  )
                {
                        for( long ion=0; ion<=nelem+1; ++ion )
                                DenseAtomsIons += dense.xIonDense[nelem][ion];
                }
        }

        /* DensElements is sum of all gas-phase densities of all elements,
         * DenseAtomsIons is sum of density of atoms and ions, does not 
         * include molecules */
        ASSERT( DenseAtomsIons > 0. );

        /* find number of molecules of the heavy elements in the gas phase. 
         * Atoms and ions were counted above.  Do not include ices, solids
         * on grains */
        den_mole = total_molecules_gasphase();
        
        /* total number of atoms, ions, and molecules in gas phase per unit vol */
        dense.xNucleiTotal = den_mole + DenseAtomsIons;
        if( dense.xNucleiTotal > BIGFLOAT )
        {
                fprintf(ioQQQ, "PROBLEM DISASTER SumDensities has found "
                        "dense.xNucleiTotal with an insane density.\n");
                fprintf(ioQQQ, "The density was %.2e\n",
                        dense.xNucleiTotal);
                TotalInsanity();
        }
        ASSERT( dense.xNucleiTotal > 0. );

        /* particle density that enters into the pressure includes electrons
         * cm-3 */
        dense.pden = (realnum)(dense.eden + dense.xNucleiTotal);
        
        /* dense.wmole is molecular weight, AtomicMassUnit per particle */
        dense.wmole = 0.;
        for( long i=0; i < LIMELM; i++ )
        {
                dense.wmole += dense.gas_phase[i]*(realnum)dense.AtomicWeight[i];
        }
        /* dense.wmole is now molecular weight, AtomicMassUnit per particle */
        dense.wmole /= dense.pden;
        ASSERT( dense.wmole > 0. && dense.pden > 0. );
        
        /* xMassDensity is density in grams per cc */
        dense.xMassDensity = (realnum)(dense.wmole*ATOMIC_MASS_UNIT*dense.pden);
        
        /*>>chng 04 may 25, introduce following test on xMassDensity0 
         * WJH 21 may 04, this is the mass density that corresponds to the hden 
         * specified in the init file. It is used to calculate the mass flux in the 
         * dynamics models. It may not necessarily be the same as struc.DenMass[0],
         * since the pressure solver may have jumped to a different density at the
         * illuminated face from that specified.*/   
        if( dense.xMassDensity0 < 0.0 )
                 dense.xMassDensity0 = dense.xMassDensity;

        return;
}

bool AbundChange( )
{
        DEBUG_ENTRY( "AbundChange()" );

        fixit(); // Breaks conservation if molecular abundance is finite

        // Suggested improved procedure:
        // 1. Scale all molecules by the scaling for their most restrictive constituent
        // 2. dense.updateXMolecules();
        // 3. Scale ionic species so as to maintain conservation
        // This will probably lead to material being dumped into the ionic species,
        // but this should slosh back at the next molecular network solve.  This
        // procedure is at least a) conservative; b) consistent with the pure-ion case.

        /* this will be set true if density or abundances change in this zone */
        bool lgChange = false;

        /* >>chng 97 jun 03, added variable abundances for element table command
         * option to get abundances off a table with element read command */
        if( abund.lgAbTaON )
        {
                lgChange = true;
                for( long nelem=1; nelem < LIMELM; nelem++ )
                {
                        if( abund.lgAbunTabl[nelem] )
                        {
                                double abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
                                  dense.gas_phase[ipHYDROGEN];

                                double hold = abun/dense.gas_phase[nelem];
                                dense.SetGasPhaseDensity( nelem, (realnum)abun );

                                for( long ion=0; ion < (nelem + 2); ion++ )
                                {
                                        dense.xIonDense[nelem][ion] *= (realnum)hold;
                                }
                        }
                }
        }

        /* this is set false if fluctuations abundances command entered,
         * when density variations are in place this is true,
         * and dense.chDenseLaw is "SINE" */
        double FacAbun = 1.;
        if( !dense.lgDenFlucOn )
        {
                static double FacAbunSav;
                double OldAbun = 0.0;
                /* abundance variations are in place */
                lgChange = true;
                if( nzone > 1 )
                {
                        OldAbun = FacAbunSav;
                }

                /* rapid abundances fluctuation */
                if( dense.lgDenFlucRadius )
                {
                        /* cycle over radius */
                        FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
                                dense.flcPhase) + dense.csecnd;
                }
                else
                {
                        /* cycle over column density */
                        FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
                                dense.flcPhase) + dense.csecnd;
                }

                if( nzone > 1 )
                {
                        FacAbun = FacAbunSav/OldAbun;
                }
        }

        if( FacAbun != 1. )
        {
                ASSERT(!dynamics.lgAdvection); // Local abundance change doesn't make sense with advective transport

                /* Li on up is affect by abundance variations */
                for( long nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
                {
                        if (dense.lgElmtOn[nelem])
                        {
                                dense.SetGasPhaseDensity(nelem, dense.gas_phase[nelem] * (realnum) FacAbun);
                                ScaleIonDensities( nelem, (realnum)(FacAbun) );
                        }
                }

                for( long mol=0; mol < mole_global.num_total; mol++ )
                {
                        mole.species[mol].den *= (realnum)(FacAbun);
                }
        }

        if (lgChange)
        {
                /* must call TempChange since ionization has changed, there are some
                 * terms that affect collision rates (H0 term in electron collision) */
                TempChange(phycon.te , false);
        }

        return lgChange;
}

namespace {
        const int SCALE_NEW = 1;
}

realnum scalingDensity(void)
{
        if (SCALE_NEW)
                return dense.xMassDensity/(realnum)ATOMIC_MASS_UNIT;
        else
                return dense.gas_phase[ipHYDROGEN];
}
realnum scalingZoneDensity(long i)
{
        if (SCALE_NEW)
                return struc.DenMass[i]/(realnum)ATOMIC_MASS_UNIT;
        else
                return struc.hden[i];
}