Changeset 1691 for branches/newmole/source

Timestamp:

12/17/07 17:55:05 (11 months ago)

Author:

rjrw

Message:

Refactor testing of whether the network holds the species density to
allow for Cgrn species.

Location:

branches/newmole/source

Files:

• dynamics.cpp (modified) (2 diffs)
• heat_sum.cpp (modified) (1 diff)
• mole_newton_step.cpp (modified) (1 diff)
• mole_species.cpp (modified) (6 diffs)
• prt_comment.cpp (modified) (1 diff)
• radius_next.cpp (modified) (1 diff)

branches/newmole/source/dynamics.cpp

@@ -950 +950 @@
                                  Old_molecules[ipUpstream][mol]/Old_density[ipUpstream])*
                                 frac_next;
-                        /* Atoms and positive ions are already counted in UpstreamElem */
-                        if(COmole[mol]->n_nuclei != 1 || COmole[mol]->charge < 0) 
+                        /* Externally represented species are already counted in UpstreamElem */
+                        if(COmole[mol]->location == NULL)
                         {
                                 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
@@ -983 +983 @@
                 {
                         Upstream_molecules[mol] = Old_molecules[ipUpstream][mol]/Old_density[ipUpstream];
-                        if(COmole[mol]->n_nuclei > 1 || COmole[mol]->charge < 0)
+                        if(COmole[mol]->location == NULL)
                         {
                                 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )

branches/newmole/source/heat_sum.cpp

@@ -118 +118 @@
         for( i=0; i < mole.num_calc; i++ )
         {
-                if(COmole[i]->n_nuclei > 1 || COmole[i]->charge < 0)
+                if(COmole[i]->location == NULL)
                         AtomicCollidDensity += COmole[i]->hevmol;
         }

branches/newmole/source/mole_newton_step.cpp

@@ -554 +554 @@
                         b[i] -= COmole[i]->hevmol*dynamics.Rate;
 
-                        if (COmole[i]->n_nuclei != 1 || COmole[i]->charge < 0)
+                        if (COmole[i]->n_nuclei != 1 || COmole[i]->charge < 0 || ! COmole[i]->lgGas_Phase)
                         {
                                 b[i] += dynamics.molecules[i]*dynamics.Rate;

branches/newmole/source/mole_species.cpp

@@ -585 +585 @@
         for (i=0;i<mole.num_total;i++) 
         {
-                if(COmole[i]->location) 
-                {
+                if(COmole[i]->location != NULL) 
                         COmole[i]->hevmol = *(COmole[i]->location);
-                }
         }
 }
@@ -606 +604 @@
         {
                 if (COmole[i]->location == NULL)
-                {
                         mole.elec += COmole[i]->hevmol*COmole[i]->charge;
-                }
         }
 
@@ -687 +683 @@
         for (i=0;i<mole.num_calc;++i) 
         {
-                if (COmole[i]->n_nuclei == 1 && COmole[i]->charge >= 0)
-                        continue;
-                for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
-                {
-                        total[nelem] += COmole[i]->hevmol*COmole[i]->nElem[nelem];
+                if (COmole[i]->location == NULL)
+                {
+                        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
+                        {
+                                total[nelem] += COmole[i]->hevmol*COmole[i]->nElem[nelem];
+                        }
                 }
         }
@@ -705 +702 @@
         for (i=0;i<mole.num_calc;++i) 
         {
-                if (COmole[i]->n_nuclei == 1 && COmole[i]->charge >= 0)
-                        continue;
-                total += (realnum) COmole[i]->hevmol;
+                if (COmole[i]->location == NULL)
+                        total += (realnum) COmole[i]->hevmol;
         }
         return total;
@@ -721 +717 @@
         for (i=0;i<mole.num_calc;++i) 
         {
-                if (COmole[i]->lgGas_Phase && (COmole[i]->n_nuclei != 1 || COmole[i]->charge < 0) )
+                if (COmole[i]->lgGas_Phase && COmole[i]->location == NULL )
                         total += (realnum) COmole[i]->hevmol;
         }
@@ -745 +741 @@
         for (i=0,j=0;i<mole.num_calc;i++) 
         {
-                if (COmole[i]->n_nuclei != 1 || COmole[i]->charge <= 0) 
+                if (COmole[i]->n_nuclei != 1 || COmole[i]->charge <= 0 || ! COmole[i]->lgGas_Phase)
                 {
                         /* Compound molecules and negative ions are represented individually */

branches/newmole/source/prt_comment.cpp

@@ -2316 +2316 @@
         lgLots_of_moles = false;
         lgLotsSolids = false;
-        /* largest fraction in any heavy element molecule */
+        /* largest fraction in any molecule */
         for( i=0; i<mole.num_calc; ++i )
         {
-                if( COmole[i]->n_nuclei == 1 && COmole[i]->charge >= 0)
-                        continue;
-
-                if( COmole[i]->xMoleFracMax > 0.1 )
-                {
-                        sprintf( chLine, "  !The fraction of %s in %s reached %.1f%% at some point in the cloud.", 
-                                elementnames.chElementSym[COmole[i]->xMoleFracMaxElem],
-                                COmole[i]->label,
-                                COmole[i]->xMoleFracMax*100. );
-                        bangin(chLine);
-                        lgLots_of_moles = true;
-                        /* check whether molecules are on grains */
-                        if( !COmole[i]->lgGas_Phase )
-                                lgLotsSolids = true;
-                }
-                else if( COmole[i]->xMoleFracMax>0.01 )
-                {
-                        sprintf( chLine, "   The fraction of %s in %s reached %.2f%% at some point in the cloud.", 
-                                elementnames.chElementSym[COmole[i]->xMoleFracMaxElem],
-                                COmole[i]->label,
-                                COmole[i]->xMoleFracMax*100. );
-                        notein(chLine);
-                        lgLots_of_moles = true;
-                        /* check whether molecules are on grains */
-                        if( !COmole[i]->lgGas_Phase )
-                                lgLotsSolids = true;
-                }
-                else if( COmole[i]->xMoleFracMax > 1e-3 )
-                {
-                        sprintf( chLine, "   The fraction of %s in %s reached %.3f%% at some point in the cloud.", 
-                                elementnames.chElementSym[COmole[i]->xMoleFracMaxElem],
-                                COmole[i]->label,
-                                COmole[i]->xMoleFracMax*100. );
-                        notein(chLine);
-                        /* check whether molecules are on grains */
-                        if( !COmole[i]->lgGas_Phase )
-                                lgLotsSolids = true;
+                if( COmole[i]->location == NULL )
+                {
+                        if( COmole[i]->xMoleFracMax > 0.1 )
+                        {
+                                sprintf( chLine, "  !The fraction of %s in %s reached %.1f%% at some point in the cloud.", 
+                                                                 elementnames.chElementSym[COmole[i]->xMoleFracMaxElem],
+                                                                 COmole[i]->label,
+                                                                 COmole[i]->xMoleFracMax*100. );
+                                bangin(chLine);
+                                lgLots_of_moles = true;
+                                /* check whether molecules are on grains */
+                                if( !COmole[i]->lgGas_Phase )
+                                        lgLotsSolids = true;
+                        }
+                        else if( COmole[i]->xMoleFracMax>0.01 )
+                        {
+                                sprintf( chLine, "   The fraction of %s in %s reached %.2f%% at some point in the cloud.", 
+                                                                 elementnames.chElementSym[COmole[i]->xMoleFracMaxElem],
+                                                                 COmole[i]->label,
+                                                                 COmole[i]->xMoleFracMax*100. );
+                                notein(chLine);
+                                lgLots_of_moles = true;
+                                /* check whether molecules are on grains */
+                                if( !COmole[i]->lgGas_Phase )
+                                        lgLotsSolids = true;
+                        }
+                        else if( COmole[i]->xMoleFracMax > 1e-3 )
+                        {
+                                sprintf( chLine, "   The fraction of %s in %s reached %.3f%% at some point in the cloud.", 
+                                                                 elementnames.chElementSym[COmole[i]->xMoleFracMaxElem],
+                                                                 COmole[i]->label,
+                                                                 COmole[i]->xMoleFracMax*100. );
+                                notein(chLine);
+                                /* check whether molecules are on grains */
+                                if( !COmole[i]->lgGas_Phase )
+                                        lgLotsSolids = true;
+                        }
                 }
         }

branches/newmole/source/radius_next.cpp

@@ -885 +885 @@
                                 abund1,
                                 abund_limit;
-                        if( COmole[i]->n_nuclei <= 1 && COmole[i]->charge >= 0)
-                                continue;
-                        /* >>chng 03 sep 21, add CO logic */
-                        /* >>chng 04 mar 30, generalize to any molecule at all */
-                        /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
-                         * to 0.001 to pick up approach to molecular fronts */
-                        abund = 0.;
-                        for (nelem=ipHYDROGEN;nelem < LIMELM; nelem++) 
-                        {
-                                if (dense.lgElmtOn[nelem])
+                        if( COmole[i]->location == NULL)
+                        {
+                                /* >>chng 03 sep 21, add CO logic */
+                                /* >>chng 04 mar 30, generalize to any molecule at all */
+                                /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
+                                 * to 0.001 to pick up approach to molecular fronts */
+                                abund = 0.;
+                                for (nelem=ipHYDROGEN;nelem < LIMELM; nelem++) 
                                 {
-                                        abund1 = (realnum)COmole[i]->hevmol*COmole[i]->nElem[nelem]/SDIV(dense.gas_phase[nelem]);
-                                        if (abund1 > abund)
+                                        if (dense.lgElmtOn[nelem])
                                         {
-                                                abund = abund1;
+                                                abund1 = (realnum)COmole[i]->hevmol*COmole[i]->nElem[nelem]/SDIV(dense.gas_phase[nelem]);
+                                                if (abund1 > abund)
+                                                {
+                                                        abund = abund1;
+                                                }
                                         }
                                 }
-                        }
-                        /* is this an ice?  need special logic for ice since density increases
-                         * exponentially when grain temp falls below sublimation temperature 
-                         * >>chng 05 dec 06 - detect changes for smaller abundances for ices
-                         * due to large changes in ice abundances */
-                        if( COmole[i]->lgGas_Phase )
-                        {
-                                abund_limit = 1e-3f;
-                        }
-                        else
-                        {
-                                /* this is an ice - track its abundance at lower values so that
-                                 * we resolve the sublimation transition region */
-                                abund_limit = 1e-5f;
-                        }
-
-                        if( abund > abund_limit )
-                        {
-                                double drmole_one, relative_change, relative_change_denom;
-                                /* >>chng 05 dec 08, use smaller abundance for the denominator since just taking
-                                 * current abundance will overlook case where current density is vastly large
-                                 * than old density */
-                                if( struc.molecules[i][nzone-3]>SMALLFLOAT )
+                                /* is this an ice?  need special logic for ice since density increases
+                                 * exponentially when grain temp falls below sublimation temperature 
+                                 * >>chng 05 dec 06 - detect changes for smaller abundances for ices
+                                 * due to large changes in ice abundances */
+                                if( COmole[i]->lgGas_Phase )
                                 {
-                                        relative_change_denom = MIN2( COmole[i]->hevmol , struc.molecules[i][nzone-3] );
+                                        abund_limit = 1e-3f;
                                 }
                                 else
                                 {
-                                        relative_change_denom = COmole[i]->hevmol;
+                                        /* this is an ice - track its abundance at lower values so that
+                                         * we resolve the sublimation transition region */
+                                        abund_limit = 1e-5f;
                                 }
-                                /* the relative change in the abundance */
-                                relative_change = 
-                                        /*fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / COmole[i]->hevmol;*/
-                                        /* >>chng 05 dec 08, use smaller abundance */
-                                        fabs( COmole[i]->hevmol - struc.molecules[i][nzone-3] ) / relative_change_denom;
-                                /* >>chng 03 jun 16, change from 0.05 to 0.1, fine resolution of H2/H exposed
-                                 * small numerical oscillations in Solomon process */
-                                /* >>chng 04 jun 02, from 0.1 back to 0.05, more extensive CO etc network
-                                 * caused oscillations in SiO abundance and Si Si+ density. */
-                                /* >>chng 04 aug 03, from 0.05 to 0.035, leiden pdr model v2 had major
-                                 * jump in eden */
-                                /* >>chng 04 oct 18, from 0.035 back to 0.05, leiden pdr v2 actually due to having
-                                 * PAHs in fully molecular limit (??), this caused cool flow pdr grid to trip on
-                                 * too small dr */
-                                relative_change = SDIV(relative_change);
-                                /*>>chng 05 dec 08, relative_change must be less than one - with
-                                 * revised logic above can be bigger than one */
-                                if( relative_change > 1. )
-                                        relative_change = 1./relative_change;
-                                /*drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );*/
-                                /* >>chng 05 aug 23, thermal front allowed dr to become much too small
-                                 * chng from 0.02 to .6 */
-                                /*drmole_one = radius.drad*MAX2(0.2,0.05/relative_change );*/
-                                drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
-                                /* final dr will be the smallest we encounter */
-                                if( drmole_one < dr_mole_abund )
+                                
+                                if( abund > abund_limit )
                                 {
-                                        /* this is the dr used to set next dr - keep track of which moe was changing */
-                                        dr_mole_abund = drmole_one;
-                                        mole_dr_change = i;
-                                        dCO_abund = relative_change;
+                                        double drmole_one, relative_change, relative_change_denom;
+                                        /* >>chng 05 dec 08, use smaller abundance for the denominator since just taking
+                                         * current abundance will overlook case where current density is vastly large
+                                         * than old density */
+                                        if( struc.molecules[i][nzone-3]>SMALLFLOAT )
+                                        {
+                                                relative_change_denom = MIN2( COmole[i]->hevmol , struc.molecules[i][nzone-3] );
+                                        }
+                                        else
+                                        {
+                                                relative_change_denom = COmole[i]->hevmol;
+                                        }
+                                        /* the relative change in the abundance */
+                                        relative_change = 
+                                                /*fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / COmole[i]->hevmol;*/
+                                                /* >>chng 05 dec 08, use smaller abundance */
+                                                fabs( COmole[i]->hevmol - struc.molecules[i][nzone-3] ) / relative_change_denom;
+                                        /* >>chng 03 jun 16, change from 0.05 to 0.1, fine resolution of H2/H exposed
+                                         * small numerical oscillations in Solomon process */
+                                        /* >>chng 04 jun 02, from 0.1 back to 0.05, more extensive CO etc network
+                                         * caused oscillations in SiO abundance and Si Si+ density. */
+                                        /* >>chng 04 aug 03, from 0.05 to 0.035, leiden pdr model v2 had major
+                                         * jump in eden */
+                                        /* >>chng 04 oct 18, from 0.035 back to 0.05, leiden pdr v2 actually due to having
+                                         * PAHs in fully molecular limit (??), this caused cool flow pdr grid to trip on
+                                         * too small dr */
+                                        relative_change = SDIV(relative_change);
+                                        /*>>chng 05 dec 08, relative_change must be less than one - with
+                                         * revised logic above can be bigger than one */
+                                        if( relative_change > 1. )
+                                                relative_change = 1./relative_change;
+                                        /*drmole_one = radius.drad*MAX2(0.2,0.035/relative_change );*/
+                                        /* >>chng 05 aug 23, thermal front allowed dr to become much too small
+                                         * chng from 0.02 to .6 */
+                                        /*drmole_one = radius.drad*MAX2(0.2,0.05/relative_change );*/
+                                        drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
+                                        /* final dr will be the smallest we encounter */
+                                        if( drmole_one < dr_mole_abund )
+                                        {
+                                                /* this is the dr used to set next dr - keep track of which moe was changing */
+                                                dr_mole_abund = drmole_one;
+                                                mole_dr_change = i;
+                                                dCO_abund = relative_change;
+                                        }
                                 }
                         }