Show
Ignore:
Timestamp:
12/16/07 16:54:58 (11 months ago)
Author:
rjrw
Message:

More chemical network rationalization and removal of dead experimental code.

Location:
branches/newmole/source
Files:
1 added
22 modified
3 moved

Legend:

Unmodified
Added
Removed

  • branches/newmole/source/atmdat_char_tran.cpp

    r1610 r1687  
    102102                        He+ + H => He + H+ as defined in the UMIST database.  This is only 
    103103                        used if the "set UMIST rates" command is used */ 
    104                 if(!co.lgUMISTrates) 
     104                if(!mole.lgUMISTrates) 
    105105                { 
    106106                        atmdat.HCharExcRecTo[ipHELIUM][0] = 4.85e-15*pow(phycon.te/300, 0.18); 
    … …  
    180180                /* use UMIST rates for charge transfer if UMIST command is used - disagree 
    181181                 * by 20% at 5000K and diverge at high temperature */ 
    182                 if(!co.lgUMISTrates) 
     182                if(!mole.lgUMISTrates) 
    183183                { 
    184184                        atmdat.HCharExcIonOf[ipOXYGEN][0] = HMRATE(7.31e-10,0.23,225.9); 
    … …  
    525525                TeUsed = phycon.te; 
    526526 
    527                 if(!co.lgUMISTrates) 
     527                if(!mole.lgUMISTrates) 
    528528                { 
    529529                        /* Set all charge transfer rates equal to zero that do not appear 

  • branches/newmole/source/cdinit.cpp

    r1610 r1687  
    4040* to verify that cdInit was called first */ 
    4141bool lgcdInitCalled=false; 
     42#include "co.h" 
    4243#include "colden.h" 
    4344#include "conv.h" 

  • branches/newmole/source/cool_carb.cpp

    r1610 r1687  
    1010#include "hmi.h" 
    1111#include "h2.h" 
     12#include "co.h" 
     13#include "ligbar.h" 
    1214#include "mole.h" 
    13 #include "ligbar.h" 
    1415#include "thermal.h" 
    1516#include "colden.h" 

  • branches/newmole/source/cool_eval.cpp

    r1643 r1687  
    140140 
    141141        /* molecular molecules molecule cooling */ 
    142         if( co.lgNoMole ) 
     142        if( mole.lgNoMole ) 
    143143        { 
    144144                /* this branch - do not include molecules */ 

  • branches/newmole/source/eden_sum.cpp

    r1686 r1687  
    6363        if( (-mole.elec) > dense.EdenTrue/4. && conv.lgSearch ) 
    6464        { 
    65                 /*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,hmole_eden);*/ 
     65                /*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,mole.elec);*/ 
    6666                dense.EdenTrue *= 0.9; 
    6767        } 
    … …  
    7171                 * negative electron density.  occurred in pdr 
    7272                 * with large h2, after hmole failure */ 
    73                 fprintf(ioQQQ," PROBLEM sum eden from hmole too neg, set to  limt.  EdenTrue:%.3e hmole_eden:%.3e \n", 
     73                fprintf(ioQQQ," PROBLEM sum eden from mole too neg, set to  limt.  EdenTrue:%.3e mole.elec:%.3e \n", 
    7474                        dense.EdenTrue,  
    7575                        mole.elec); 

  • branches/newmole/source/init_sim_postparse.cpp

    r1686 r1687  
    2525 
    2626        mole.init(); 
    27         mole_zero(); 
     27        mole.zero(); 
    2828 
    2929        /* >> chng 06 Nov 24 rjrw: Move reaction definitions here so they can be controlled by parsed commands */ 

  • branches/newmole/source/ion_carbo.cpp

    r1610 r1687  
    5757        /* >>chng 05 aug 10, add Leiden hack here to get same C0 photo rate 
    5858         * as UMIST - negates difference in grain opacities */ 
    59         if(!co.lgUMISTrates) 
     59        if(!mole.lgUMISTrates) 
    6060        { 
    6161                int nelem=ipCARBON , ion=0 , ns=2; 

  • branches/newmole/source/ion_trim.cpp

    r1610 r1687  
    106106        /* logic for PDRs, for elements included in chemistry, need stable solutions,  
    107107         * keep 3 ion stages in most cases - added by NA to do HII/PDR sims */ 
    108         if( !co.lgNoMole ) 
     108        if( !mole.lgNoMole ) 
    109109        { 
    110110                trimlo = SMALLFLOAT; 

  • branches/newmole/source/iter_startend.cpp

    r1653 r1687  
    1919#include "geometry.h" 
    2020#include "h2.h" 
     21#include "co.h" 
    2122#include "he.h" 
    2223#include "grains.h" 

  • branches/newmole/source/mole.h

    r1686 r1687  
    99enum mole_state {MOLE_NULL, MOLE_PASSIVE, MOLE_ACTIVE}; 
    1010 
    11 /**mole_drive main driver for heavy molecular equilibrium routines */ 
     11/**mole_drive main driver for chemical equilibrium routines */ 
    1212extern void mole_drive(void); 
    13  
    14 /**mole_zero allocate + initialize workspace */ 
    15 extern void mole_zero(void); 
    1613 
    1714/**mole_create_react build reaction structures */ 
    1815extern void mole_create_react(void); 
    19  
    20 /** called from cdInit to initialized co routines */ 
    21 extern void mole_Init(void); 
    22 /**mole_update_rks update rate coefficients, only temp part */ 
    23 extern void mole_update_rks( void ); 
    24  
    25 extern void mole_update_species_cache(void); 
    26  
    27 extern void mole_return_cached_species(void); 
    2816 
    2917extern double mole_sink_rate(const char chSpecies[]); 
    … …  
    4028 
    4129extern double mole_findrate(const char buf[]); 
    42  
    43 /** Take one Newton step of the chemical network  
    44 \param *nBad 
    45 \param *error 
    46 */ 
    47 void mole_newton_step(int *nBad, realnum *error, realnum numelem[]); 
    4830 
    4931/** \verbatim >>chng 03 feb 09, rm ipH3P_hev, since not used, and decrement NUM_HEAVY_MOLEC to 17  
    … …  
    7961   can be important to the ionization balance */ 
    8062 
    81 EXTERN struct t_co { 
     63class Mole { 
    8264 
    83         /** CODissHeat is CO Photodissociation heating */ 
    84         realnum CODissHeat,  
    85           /**< largest fraction of total heating */ 
    86           codfrc,  
    87           /**< total heating integrated over cloud */ 
    88           codtot; 
    89  
    90         /** flag saying what fraction of total cooling was carried by CO */ 
    91         double COCoolBigFrac; 
    92  
    93         /** flag saying CO is important but lines capped by highest level */ 
    94         bool lgCOCoolCaped; 
     65        public: 
     66        void init(void); 
     67                 
     68  /**mole_zero allocate + initialize workspace */ 
     69        void zero(void); 
    9570 
    9671        /** flag to turn off all molecules, set with no molecules command */ 
    … …  
    9974        /** flag to turn off heavy molecules, set with no heavy molecules command */ 
    10075        bool lgNoHeavyMole; 
    101  
    102         /** C12/C13 isotope ratio, sets the ratio of C12O16 to C13O16 and  
    103         * C13/C12 1909, initialized in zero.c  */ 
    104         realnum C12_C13_isotope_ratio; 
    10576 
    10677        /** flag set true if H2O destruction rate went to zero */ 
    … …  
    12798        bool lgNeutrals; 
    12899 
    129         double h2lim; 
    130  
    131 }       co; 
    132  
    133 class Mole { 
    134  
    135         public: 
    136         void init(void); 
    137                  
    138100        long nzone , iteration; 
    139101 
    140102         /** do we include capture of molecules onto grain surfaces?  default is true, 
    141103          * turned off with NO GRAIN MOLECULES command */ 
    142          bool lgGrain_mole_deplete; 
     104        bool lgGrain_mole_deplete; 
    143105 
    144106        /** flag saying whether an element is in the chemistry network */ 

  • branches/newmole/source/mole_co_atom.cpp

    r1610 r1687  
    1515#include "lines_service.h" 
    1616#include "cddrive.h" 
    17 #include "mole.h" 
     17#include "co.h" 
    1818/*lint -e778 constant expression evaluatess to 0 in operation '-' */ 
    1919/*=================================================================*/ 

  • branches/newmole/source/mole_drive.cpp

    r1686 r1687  
    2727#include "dynamics.h" 
    2828#include "h2.h" 
     29#include "co.h" 
    2930 
    3031/* says whether the co network is currently set to zero */ 
    3132static bool lgMoleZeroed=true; 
    3233 
    33 /* the limit to H2/H where we will solve for CO */ 
    34 static double h2lim; 
    35  
    3634/* this is the error tolerance for reporting non-convergence */ 
    37 static const double COTOLER = 0.10; 
    38  
    39 /* flag to control debug statement, if 0 then none, just loops when 1, more when 2 */ 
    40 static const int CODEBUG = 0; 
     35static const double MOLETOLER = 0.10; 
    4136 
    4237/*lgMolecSetup use previous iteration's converged solution to initialize if available */ 
    43 STATIC void MolecSetup ( 
    44         char *chReason ); 
     38STATIC void MolecSetup (void); 
    4539STATIC void mole_effects(void); 
    4640 
    47 static void mole_null_results(void); 
     41STATIC void mole_null_results(void); 
    4842STATIC void mole_h_rate_diagnostics(void); 
    49 void mole_h_fixup(void); 
    5043 
    5144struct molecule **COmole; 
    … …  
    5548void mole_drive(void) 
    5649{ 
    57         bool lgNegPop,  
    58           lgZerPop, 
    59           lgFirstTry; 
    6050        bool lgConverged; 
    61         long int i; 
    62  
    63         /* this will give reason CO not converged */ 
    64         char chReason[100]; 
     51 
    6552        double solomon_assumed, 
    6653                error; 
    … …  
    6855        DEBUG_ENTRY( "mole_drive()" ); 
    6956 
    70         /* h2lim is smallest ratio of H2/HDEN for which we will 
    71          * even try to invert heavy element network */ 
    72         /* this limit is important since the co molecular network is first turned 
    73          * on when this is hit.  It's conservation law will only ever include the 
    74          * initial O, O+, and C, C+, so must not start before nearly all 
    75          * C and O is in these forms */ 
    76         h2lim = 1e-15; 
    77         /* >>chng 05 jul 15, from 1e-15 to 1e-12, see whether results are stable, 
    78          * this does include CO in the H+ zone in orion_hii_pdr 
    79          * a problem is that, during search phase, the first temp is 4000K and the 
    80          * H2 abundance is larger than it will be at the illuminated face.  try to 
    81          * avoid turning H2 on too soon */ 
    82         h2lim = 1e-12; 
    83  
    84         if( conv.nTotalIoniz==0 && iteration==0 ) 
    85         { 
    86                 /* we should have an H0 soln at this point  
    87                          -- all species for network will be sourced from ionization solvers  
    88                          in mole_update_species_cache below */ 
    89                 ASSERT( dense.xIonDense[ipHYDROGEN][0]>SMALLFLOAT ); 
    90         } 
    91  
    9257        /* this branch we will not do the CO equilibrium, set some variables that 
    9358         * would be calculated by the chemistry, zero others, and return */ 
    94         if( 1 ) 
    95         { 
    96                 /* do we need to zero out the arrays and vars? */ 
    97                 if( !lgMoleZeroed ) 
    98                         for( i=0; i<mole.num_calc; ++i ) 
    99                         { 
    100                                 COmole[i]->hevmol = 0.;  
    101                         } 
    102                          
     59        mole_null_results(); 
     60 
     61        MolecSetup(); 
     62 
     63        mole_update_species_cache();  /* Update densities of species controlled outside the chemical network */ 
     64         
     65        mole_update_rks(); 
     66         
     67        solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/SDIV(hmi.H2_rate_destroy); 
     68         
     69        error = mole_solve();    
     70         
     71        lgConverged = (error < MOLETOLER);  
     72        if( lgConverged && h2.lgH2ON  && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 ) 
     73        { 
     74                /* >>chng 05 mar 11, go from "total" to H2g for solomon rate */ 
     75                if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) >  
     76                                conv.EdenErrorAllowed/5.) 
     77                { 
     78                        lgConverged = false; 
     79                } 
     80        } 
     81 
     82        mole_effects(); 
     83         
     84        return; 
     85} 
     86 
     87STATIC void MolecSetup(void) 
     88{ 
     89        long int i; 
     90        /* this will be used to rescale old saved abundances in constant pressure model */ 
     91        static realnum hden_save_prev_call; 
     92        const bool DEBUG_MOLECAVER = false; 
     93 
     94        DEBUG_ENTRY( "MolecSetup()" ); 
     95 
     96        /* mole_Init set this to -2 when code initialized, so negative 
     97         * number shows very first call in this model */ 
     98        /* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration 
     99         * and we are below limit where it should be evaluated */ 
     100        if( mole.iteration < 0 || lgMoleZeroed ) 
     101        { 
     102 
     103                /* very first attempt at ever obtaining a solution - 
     104                 * called one time per model since co.iteration_co set negative 
     105                 * when cdInit called */ 
     106                 
     107                /* >>chng 05 jun 24, during map phase do not reset molecules to zero 
     108                 * since we probably have a better estimate right now */ 
     109                 
     110                /* we should have a neutral carbon solution at this point */ 
     111                //ASSERT( !dense.lgElmtOn[ipCARBON] || dense.xIonDense[ipCARBON][0]>SMALLFLOAT ); 
     112                                 
     113                /* set iteration flag */ 
     114                mole.iteration = iteration; 
     115                mole.nzone = nzone; 
     116                 
     117                /* for constant pressure models when molecules are reset on second 
     118                 * and higher iterations, total density will be different, so 
     119                 * must take this into account when abundances are reset */ 
     120                hden_save_prev_call = dense.gas_phase[ipHYDROGEN]; 
     121                 
     122                if( DEBUG_MOLECAVER ) 
     123                        fprintf(ioQQQ," DEBUG MolecSetup iteration %li zone %li zeroing since very first call. H2/H=%.2e\n",  
     124                                                        iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 
     125        } 
     126        else if( iteration > mole.iteration ) 
     127        { 
     128                realnum ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_call; 
     129                 
     130                /* this is first call on new iteration, reset 
     131                 * to first initial abundances from previous iteration */ 
     132                for( i=0; i < mole.num_calc; i++ ) 
     133                { 
     134                        COmole[i]->hevmol = COmole[i]->hev_reinit*ratio; 
     135                } 
     136 
     137                mole.iteration = iteration; 
     138                mole.nzone = nzone; 
     139                 
     140                if( DEBUG_MOLECAVER ) 
     141                        fprintf(ioQQQ," DEBUG MolecSetup iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n",  
     142                                                        iteration, 
     143                                                        nzone, 
     144                                                        hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 
     145        } 
     146        else if( iteration == mole.iteration && nzone==mole.nzone+1 ) 
     147        { 
     148                /* this branch, second zone with solution, so save previous 
     149                 * zone's solution to reset things in next iteration */ 
     150                for( i=0; i < mole.num_calc; i++ ) 
     151                { 
     152                        COmole[i]->hev_reinit = (realnum) COmole[i]->hevmol; 
     153                } 
     154                 
     155                mole.nzone = -2; 
     156                hden_save_prev_call = dense.gas_phase[ipHYDROGEN]; 
     157                 
     158                if( DEBUG_MOLECAVER ) 
     159                        fprintf(ioQQQ,"DEBUG MolecSetup iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone); 
     160        } 
     161         
     162        return; 
     163} 
     164 
     165static void mole_null_results(void) 
     166{ 
     167        DEBUG_ENTRY( "mole_null_results()" ); 
     168        /* these are heavy - heavy charge transfer rate that will still be needed 
     169         * for recombination of Si+, S+, and C+ */ 
     170        long int nelem , ion, i; 
     171 
     172        /* heating due to CO photodissociation */ 
     173        /* do we need to zero out the arrays and vars? */ 
     174        if( !lgMoleZeroed ) 
     175        { 
    103176                lgMoleZeroed = true; 
     177                for( i=0; i<mole.num_calc; ++i ) 
     178                { 
     179                        COmole[i]->hevmol = 0.; 
     180                } 
    104181                /* heating due to CO photodissociation */ 
    105182                thermal.heating[0][9] = 0.; 
    … …  
    110187                CoolHeavy.dC13O16Rot = 0.; 
    111188                co.CODissHeat = 0.; 
    112                  
     189 
    113190                /** \todo       2       use TransitionZero here? */ 
    114191                for( i=0; i < nCORotate; i++ ) 
    … …  
    134211                        C13O16Rotate[i].Emis->phots = 0.; 
    135212                } 
    136         } 
    137  
    138         MolecSetup(chReason); 
    139  
    140         mole_update_species_cache();  /* Update densities of species controlled outside the chemical network */ 
    141          
    142         mole_update_rks(); 
    143          
    144         solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/SDIV(hmi.H2_rate_destroy); 
    145          
    146         /* this flag is used to stop calculation due to negative abundances */ 
    147         lgNegPop = false; 
    148          
    149         error = mole_solve(&lgNegPop,&lgZerPop );        
    150          
    151         lgConverged = (error < COTOLER);  
    152         if( lgConverged && h2.lgH2ON  && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 ) 
    153         { 
    154                 /* >>chng 05 mar 11, go from "total" to H2g for solomon rate */ 
    155                 if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) >  
    156                                 conv.EdenErrorAllowed/5.) 
    157                 { 
    158                         lgConverged = false; 
    159                 } 
    160         } 
    161         if( CODEBUG ) 
    162                 fprintf(ioQQQ,"codrivbug %.2f Neg?:%c\tZero?:%c\tOH new\t%.3e\tCO\t%.3e\tTe\t%.3e\tH2/H\t%.2e\n", 
    163                                                 fnzone , 
    164                                                 TorF(lgNegPop) ,  
    165                                                 TorF(lgZerPop) ,  
    166                                                 findspecies("OH")->hevmol , 
    167                                                 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 
    168                                                 phycon.te, 
    169                                                 hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 
    170  
    171         mole_effects(); 
    172         mole_h_rate_diagnostics(); 
    173          
    174         /* say that we have found a solution before */ 
    175         /* lgMoleZeroed = false; */ 
    176          
    177         /* this flag says whether this is the first zone we are trying 
    178          * to converge the molecules - there will be problems during initial 
    179          * search so do not print anything in this case */ 
    180         lgFirstTry = (nzone==mole.nzone && iteration==mole.iteration); 
    181  
    182         /* did the molecule network have negative pops? */ 
    183         if( lgNegPop ) 
    184         { 
    185                 if( conv.lgSearch && (hmi.H2_total/dense.gas_phase[ipHYDROGEN] <h2lim) && 
    186                         (iteration==1) ) 
    187                 { 
    188                         /* we are in search phase during the first iteration,  
    189                          * and the H2/H ratio has fallen 
    190                          * substantially below the threshold for solving the CO network. 
    191                          * Turn it off */ 
    192                         /* >> chng 07 jan 10 rjrw: this was mole_Init(), but the comment suggests 
    193                          * it should really be mole_zero */ 
    194                         mole_zero(); 
    195                 } 
    196                 else 
    197                 { 
    198                         if( called.lgTalk && !lgFirstTry ) 
    199                         { 
    200                                 conv.lgConvPops = false; 
    201                                 fprintf( ioQQQ,  
    202                                                 " PROBLEM mole_drive failed to converge1 due to negative pops, zone=%.2f,  CO/C=%.2e  OH/C=%.2e H2/H=%.2e\n",  
    203                                                 fnzone,  
    204                                                 findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 
    205                                                 findspecies("OH")->hevmol/SDIV(dense.gas_phase[ipCARBON]), 
    206                                                 hmi.H2_total/dense.gas_phase[ipHYDROGEN]); 
    207                                 ConvFail( "pops" , "CO" ); 
    208                         } 
    209                 } 
    210         } 
    211          
    212  
    213         /* make sure we have not used more than all the CO */ 
    214         ASSERT(conv.lgSearch || findspecies("CO")->hevmol/SDIV(dense.gas_phase[ipCARBON]) <= 1.01 ); 
    215         /*fprintf(ioQQQ,"ratioo o\t%c\t%.2f\t%f\n", TorF(conv.lgSearch), 
    216                 fnzone , co.hevmol[ipCO]/dense.gas_phase[ipCARBON] );*/ 
    217  
    218         /* these are the elements whose converge we check */ 
    219         /* this is a self consistency check, for converged solution */ 
    220         /* >>chng 04 dec 02, this test is turn  back on - don't know why it was turned off */ 
    221         if(0) { 
    222                 double sum_ion , sum_mol; 
    223                 sum_ion = dense.xIonDense[ipCARBON][0] + dense.xIonDense[ipCARBON][1]; 
    224                 sum_mol = findspecies("C")->hevmol + findspecies("C+")->hevmol; 
    225                 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 ) 
    226                 { 
    227                         /*fprintf(ioQQQ, 
    228                                 "non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n", 
    229                                 5, fnzone, sum_ion , sum_mol);*/ 
    230                         conv.lgConvIoniz = false; 
    231                         strcpy( conv.chConvIoniz, "CO con2" ); 
    232                         conv.BadConvIoniz[0] = sum_ion; 
    233                         conv.BadConvIoniz[1] = sum_mol; 
    234                         if( CODEBUG ) 
    235                                 fprintf(ioQQQ,"mole_drive not converged\n"); 
    236                 } 
    237                 sum_ion = dense.xIonDense[ipOXYGEN][0] + dense.xIonDense[ipOXYGEN][1]; 
    238                 sum_mol = findspecies("O")->hevmol + findspecies("O+")->hevmol; 
    239                 if( fabs(sum_ion-sum_mol)/SDIV(sum_mol) > 1e-2 ) 
    240                 { 
    241                         /*fprintf(ioQQQ, 
    242                                 "non conservation element %li zone %.2f sum_ion %.3e sum_mol %.3e\n", 
    243                                 7, fnzone, sum_ion , sum_mol);*/ 
    244                         conv.lgConvIoniz = false; 
    245                         strcpy( conv.chConvIoniz, "CO con3" ); 
    246                         conv.BadConvIoniz[0] = sum_ion; 
    247                         conv.BadConvIoniz[1] = sum_mol; 
    248                         if( CODEBUG ) 
    249                                 fprintf(ioQQQ,"mole_drive not converged\n"); 
    250                 } 
    251         } 
    252         return; 
    253 } 
    254  
    255 STATIC void MolecSetup( 
    256         char *chReason ) 
    257 { 
    258         long int i; 
    259         /* this will be used to rescale old saved abundances in constant pressure model */ 
    260         static realnum hden_save_prev_call; 
    261         const bool DEBUG_MOLECAVER = false; 
    262  
    263         DEBUG_ENTRY( "MolecSetup()" ); 
    264  
    265         /* this will become reason not converged */ 
    266         strcpy( chReason , "none given" ); 
    267  
    268          
    269         /* mole_Init set this to -2 when code initialized, so negative 
    270          * number shows very first call in this model */ 
    271