Changeset 1686 for branches/newmole/source

Timestamp:

12/16/07 13:43:02 (11 months ago)

Author:

rjrw

Message:

More renaming of variables and moving other stuff to more logical places.

Promote struct mole_s to class Mole.

Location:

branches/newmole/source

Files:

• conv_base.cpp (modified) (4 diffs)
• dynamics.cpp (modified) (2 diffs)
• eden_sum.cpp (modified) (6 diffs)
• h2.h (modified) (1 diff)
• hcmap.cpp (modified) (1 diff)
• heat_sum.cpp (modified) (1 diff)
• init_coreload.cpp (modified) (2 diffs)
• init_defaults_preparse.cpp (modified) (2 diffs)
• init_sim_postparse.cpp (modified) (1 diff)
• mole.h (modified) (4 diffs)
• mole_co_drive.cpp (modified) (7 diffs)
• mole_co_etc.cpp (modified) (16 diffs)
• mole_eval_balance.cpp (modified) (5 diffs)
• mole_h2.cpp (modified) (41 diffs)
• mole_h2_coll.cpp (modified) (20 diffs)
• mole_h2_create.cpp (modified) (21 diffs)
• mole_h2_etc.cpp (modified) (12 diffs)
• mole_h2_io.cpp (modified) (9 diffs)
• mole_newton_step.cpp (modified) (5 diffs)
• mole_reactions.cpp (modified) (5 diffs)
• parse_atom_h2.cpp (modified) (19 diffs)
• parse_set.cpp (modified) (4 diffs)
• prt_comment.cpp (modified) (1 diff)
• radius_next.cpp (modified) (1 diff)

branches/newmole/source/conv_base.cpp

@@ -199 +199 @@
          * will confirm that these are converged */
         EdenTrue_old = dense.EdenTrue;
-        EdenFromMolecOld = co.comole_eden;
+        EdenFromMolecOld = mole.elec;
         EdenFromGrainsOld = gv.TotalEden;
         HeatingOld = thermal.htot;
@@ -867 +867 @@
 
         /* check on molecular electron den */
-        if( fabs(EdenFromMolecOld - co.comole_eden)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
+        if( fabs(EdenFromMolecOld - mole.elec)/dense.EdenTrue > conv.EdenErrorAllowed/2. )
         {
                 conv.lgConvIoniz = false;
@@ -905 +905 @@
 
         /* check on sum of grain and molecular electron den - often two large numbers that cancel */
-        if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (co.comole_eden-gv.TotalEden) )/dense.eden > 
+        if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (mole.elec-gv.TotalEden) )/dense.eden > 
                 conv.EdenErrorAllowed/4. )
         {
@@ -911 +911 @@
                 sprintf( conv.chConvIoniz, "edn grn e" );
                 conv.BadConvIoniz[0] = (EdenFromMolecOld-EdenFromGrainsOld);
-                conv.BadConvIoniz[1] = (co.comole_eden-gv.TotalEden);
+                conv.BadConvIoniz[1] = (mole.elec-gv.TotalEden);
         }
 

branches/newmole/source/dynamics.cpp

@@ -951 +951 @@
                                 frac_next;
                         /* Atoms and positive ions are already counted in UpstreamElem */
-                        if(COmole[mol]->n_nuclei != 1 || COmole[mol]->nElec < 0) 
+                        if(COmole[mol]->n_nuclei != 1 || COmole[mol]->charge < 0) 
                         {
                                 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]->nElec < 0)
+                        if(COmole[mol]->n_nuclei > 1 || COmole[mol]->charge < 0)
                         {
                                 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )

branches/newmole/source/eden_sum.cpp

@@ -18 +18 @@
 int eden_sum(void)
 {
-        long int i,
+        long int 
           ion, 
           nelem;
@@ -54 +54 @@
         dense.EdenTrue += sum_all_ions;
         ASSERT( dense.EdenTrue >= 0. );
-
-        /* add on molecules */
-        co.comole_eden = 0.;
-        for(i=0;i<mole.num_calc;i++)
-        {
-                if (COmole[i]->n_nuclei > 1 || COmole[i]->nElec < 0)
-                        co.comole_eden += COmole[i]->hevmol*COmole[i]->nElec;
-        }
         
 #if 1
@@ -69 +61 @@
         /* >>chng 04 feb 20, recoil_ion.in shows this important,
          * update test so prevent neg eden */
-        if( (-co.comole_eden) > dense.EdenTrue/4. && conv.lgSearch )
+        if( (-mole.elec) > dense.EdenTrue/4. && conv.lgSearch )
         {
                 /*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,hmole_eden);*/
                 dense.EdenTrue *= 0.9;
         }
-        else if( (-co.comole_eden) > dense.EdenTrue )
+        else if( (-mole.elec) > dense.EdenTrue )
         {
                 /* >>chng 04 mar 14, add this branch to prevent
@@ -81 +73 @@
                 fprintf(ioQQQ," PROBLEM sum eden from hmole too neg, set to  limt.  EdenTrue:%.3e hmole_eden:%.3e \n",
                         dense.EdenTrue, 
-                        co.comole_eden);
+                        mole.elec);
                 dense.EdenTrue = dense.EdenTrue/2.;
         }
@@ -87 +79 @@
         {
                 /* account for electrons on all molecules */
-                dense.EdenTrue += co.comole_eden;
+                dense.EdenTrue += mole.elec;
         }
 #else
-        dense.EdenTrue += co.comole_eden;
+        dense.EdenTrue += mole.elec;
 #endif
 
@@ -149 +141 @@
                   dense.EdenTrue , 
                   sum_all_ions ,
-                  co.comole_eden ,
+                  mole.elec ,
                   gv.TotalEden*gv.lgGrainElectrons,
                   dense.EdenExtra ,
-                  sum_all_ions + co.comole_eden + gv.TotalEden*gv.lgGrainElectrons + dense.EdenExtra,
+                  sum_all_ions + mole.elec + gv.TotalEden*gv.lgGrainElectrons + dense.EdenExtra,
                   rfield.otslin[2182] );
 

branches/newmole/source/h2.h

@@ -159 +159 @@
         long int Jlowest[N_H2_ELEC];
 
+        /** std and mean for the noise, log normal distribution */
+        double xMeanNoise , xSTDNoise;
+
+        /** limit to the ratio H2/Htot - if ratio is below this, large atom is not called */
+        double H2_to_H_limit;
+
+        /** the number of electronic quantum states to include.
+        * To do both Lyman and Werner bands want nelec = 3 */
+        long int n_h2_elec_states;
+
+        /** this is option to use estimates of the collision rates from g-bar approximations */
+        /** turn mole.lgColl_gbar on/off with atom h2 gbar on off */
+        bool lgColl_gbar;
+
+        /** this is option to turn off the calculated collision rates */
+        bool lgColl_deexec_Calc;
+
+        /** this is option to turn off guesses of collisional dissociation rates */
+        bool lgColl_dissoc_coll;
+
         /* true to use 2007 set of H2 - H collision rate, false use 1999 */
         bool lgH2_H_coll_07;
 
+        /** include collision rates that come from real calculations,
+         * off with atom h2 collisions off command */
+        bool lgH2_grain_deexcitation;
+
+        /** flag to force LTE level populations, atom H2 LTE */
+        bool lgH2_LTE;
+
+        /** option to turn off ortho-para collisions, command ATOM H2 COLLISIONS ORTHO PARA OFF */
+        bool lgH2_ortho_para_coll_on;
+
+        /** which set of He - H2 collisions to use? default is ORNL, other
+         * is Le BOURlet */
+        bool lgH2_He_ORNL;
+
+        /** turn on trace information */
+        int nH2_TRACE;
+
+        /** put noise into collision rates */
+        bool lgH2_NOISE ,
+                /** noise for the CR collisions */
+                lgH2_NOISECOSMIC;
+
+        /** this sets how fine a trace we want for atom  h2 trace */
+        int nH2_trace_final , 
+                nH2_trace_iterations , 
+                nH2_trace_full,
+                nH2_trace_matrix;
+        
 }       h2;
 

branches/newmole/source/hcmap.cpp

@@ -235 +235 @@
 
                 /* this will force reinitialization of co network */
-                co.iteration_co = -2;
+                mole.iteration = -2;
 
                 /* now get ionization solution for this temperature */

branches/newmole/source/heat_sum.cpp

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

branches/newmole/source/init_coreload.cpp

@@ -11 +11 @@
 #include "date.h"
 #include "version.h"
-#include "mole.h"
+#include "h2.h"
 #include "extinc.h"
 #include "heavy.h"
@@ -344 +344 @@
         /* these are used to set trace levels of output by options on atom h2 trace command 
          * first is minimum level of trace, keyword is FINAL */
-        mole.nH2_trace_final = 1;
+        h2.nH2_trace_final = 1;
         /* intermediate level of trace, info per iteration, key ITERATION */
-        mole.nH2_trace_iterations = 2;
+        h2.nH2_trace_iterations = 2;
         /* full trace, keyword is FULL */
-        mole.nH2_trace_full = 3;
+        h2.nH2_trace_full = 3;
         /* print matrices used in solving X */
-        mole.nH2_trace_matrix = 4;
+        h2.nH2_trace_matrix = 4;
 
         /* initialization for punch files - must call after input commands have

branches/newmole/source/init_defaults_preparse.cpp

@@ -12 +12 @@
 #include "extinc.h"
 #include "h2.h"
-#include "mole.h"
 #include "wind.h"
 #include "thermal.h"
@@ -69 +68 @@
 
         /* holds options for trace set with atom h2 command */
-        mole.nH2_TRACE = false;
+        h2.nH2_TRACE = false;
 
         /* option to scramble collision data */
-        mole.lgH2_NOISE = false;
-        mole.lgH2_NOISECOSMIC = false; 
+        h2.lgH2_NOISE = false;
+        h2.lgH2_NOISECOSMIC = false; 
 
         /* option to turn off or on gbar collisions of the collision rate,
         * default is to have it on */
-        /* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
-        mole.lgColl_gbar = true;
+        /* turn h2.lgColl_gbar on/off with atom h2 gbar on off */
+        h2.lgColl_gbar = true;
 
         /* include collision rates that come from real calculations,
         * off with atom h2 collisions off command */
-        mole.lgColl_deexec_Calc = true;
-
-        mole.lgColl_dissoc_coll = true;
+        h2.lgColl_deexec_Calc = true;
+
+        h2.lgColl_dissoc_coll = true;
 
         /* option to turn off H2 - grain collision & deexcitation,
         * atom h2 grain collision on/off */
-        mole.lgH2_grain_deexcitation = false;
+        h2.lgH2_grain_deexcitation = false;
 
         /* option to turn off ortho-para collisions, command ATOM H2 COLLISIONS ORTHO PARA OFF */
-        mole.lgH2_ortho_para_coll_on = true;
+        h2.lgH2_ortho_para_coll_on = true;
 
         /* which set of H2 - He collision data to use?  default is Le Bourlot Meudon set,
         * set to ORNL with command
         * atom H2 He collisions ORNL */
-        mole.lgH2_He_ORNL = true;
+        h2.lgH2_He_ORNL = true;
 
         /* use 1999 or 2007 H2 - H collision set atomic data mole H2 h */

branches/newmole/source/init_sim_postparse.cpp

@@ -24 +24 @@
         thermal.lgUnstable = false;
 
-        mole_Init();
+        mole.init();
         mole_zero();
 

branches/newmole/source/mole.h

@@ -129 +129 @@
         double h2lim;
 
-        long co_nzone , iteration_co;
-
-        /** total electron density in molecules */
-        double comole_eden;
-
 }       co;
 
-EXTERN struct t_mole {
-
+class Mole {
+
+        public:
+        void init(void);
                 
-        /** limit to the ratio H2/Htot - if ratio is below this, large atom is not called */
-        double H2_to_H_limit;
-
-        /** the number of electronic quantum states to include.
-        * To do both Lyman and Werner bands want nelec = 3 */
-        long int n_h2_elec_states;
-
-        /** this is option to use estimates of the collision rates from g-bar approximations */
-        /** turn mole.lgColl_gbar on/off with atom h2 gbar on off */
-        bool lgColl_gbar;
-
-        /** this is option to turn off the calculated collision rates */
-        bool lgColl_deexec_Calc;
-
-        /** this is option to turn off guesses of collisional dissociation rates */
-        bool lgColl_dissoc_coll;
-
-        /** include collision rates that come from real calculations,
-         * off with atom h2 collisions off command */
-        bool lgH2_grain_deexcitation;
-
-        /** flag to force LTE level populations, atom H2 LTE */
-        bool lgH2_LTE;
-
-        /** option to turn off ortho-para collisions, command ATOM H2 COLLISIONS ORTHO PARA OFF */
-        bool lgH2_ortho_para_coll_on;
-
-        /** which set of He - H2 collisions to use? default is ORNL, other
-         * is Le BOURlet */
-        bool lgH2_He_ORNL;
-
-        /** turn on trace information */
-        int nH2_TRACE;
-
-        /** put noise into collision rates */
-        bool lgH2_NOISE ,
-                /** noise for the CR collisions */
-                lgH2_NOISECOSMIC;
-
-        /** this sets how fine a trace we want for atom  h2 trace */
-        int nH2_trace_final , 
-                nH2_trace_iterations , 
-                nH2_trace_full,
-                nH2_trace_matrix;
-        
+        long nzone , iteration;
+
          /** do we include capture of molecules onto grain surfaces?  default is true,
           * turned off with NO GRAIN MOLECULES command */
          bool lgGrain_mole_deplete;
-
-        /** std and mean for the noise, log normal distribution */
-        double xMeanNoise , xSTDNoise;
 
         /** flag saying whether an element is in the chemistry network */
@@ -195 +146 @@
         
         realnum eden_f, grain_area;
+
+        /** total charge in molecules */
+        double elec;
 
         /** these are source and sink terms for the ionization ladder, for chemical
@@ -205 +159 @@
                 *snk;
         realnum ***xMoleChTrRate;/***[LIMELM][LIMELM+1][LIMELM+1];*/
-}       mole;
+};
+
+extern Mole mole;
 
 enum {CHARS_SPECIES=7};
@@ -214 +170 @@
         char  label[CHARS_SPECIES];      /** molecule name */
         int   nElem[LIMELM];    /** number of each element in molecule */
-        int   nElec;          /* Charge on species/number of e- liberated by formation */
+        int   charge;          /* Charge on species/number of e- liberated by formation */
         int   Excit;          /* Is species excited (e.g. H2*) */
         bool  lgGas_Phase;    /** Solid or gas phase? */

branches/newmole/source/mole_co_drive.cpp

@@ -178 +178 @@
          * to converge the molecules - there will be problems during initial
          * search so do not print anything in this case */
-        lgFirstTry = (nzone==co.co_nzone && iteration==co.iteration_co);
+        lgFirstTry = (nzone==mole.nzone && iteration==mole.iteration);
 
         /* did the molecule network have negative pops? */
@@ -271 +271 @@
         /* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration
          * and we are below limit where it should be evaluated */
-        if( iteration!=co.iteration_co && 
+        if( iteration!=mole.iteration && 
                         hmi.H2_total/dense.gas_phase[ipHYDROGEN] < h2lim )
         {
                 return;
         }
-        else if( co.iteration_co < 0 || lgMoleZeroed )
+        else if( mole.iteration < 0 || lgMoleZeroed )
         {
 
@@ -290 +290 @@
                                 
                 /* set iteration flag */
-                co.iteration_co = iteration;
-                co.co_nzone = nzone;
+                mole.iteration = iteration;
+                mole.nzone = nzone;
                 
                 /* for constant pressure models when molecules are reset on second
@@ -302 +302 @@
                                                         iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
         }
-        else if( iteration > co.iteration_co )
+        else if( iteration > mole.iteration )
         {
                 realnum ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_call;
@@ -313 +313 @@
                 }
 
-                co.iteration_co = iteration;
-                co.co_nzone = nzone;
+                mole.iteration = iteration;
+                mole.nzone = nzone;
                 
                 if( DEBUG_MOLECAVER )
@@ -322 +322 @@
                                                         hmi.H2_total/dense.gas_phase[ipHYDROGEN]);
         }
-        else if( iteration == co.iteration_co && nzone==co.co_nzone+1 )
+        else if( iteration == mole.iteration && nzone==mole.nzone+1 )
         {
                 /* this branch, second zone with solution, so save previous
@@ -331 +331 @@
                 }
                 
-                co.co_nzone = -2;
+                mole.nzone = -2;
                 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];
                 

branches/newmole/source/mole_co_etc.cpp

@@ -51 +51 @@
 struct molecule **groupspecies;
 
+Mole mole;
+
 /*=================================================================*/
-/*CO_Init called from cdInit to initialize CO routines */
-void mole_Init(void)
+/*mole_Init called from cdInit to initialize CO routines */
+void Mole::init(void)
 {
         long int i,
@@ -59 +61 @@
         struct molecule *sp;
         static realnum one=1.0;
-        static bool lgCO_Init_called=false;
+        static bool lgmole_Init_called=false;
 
         DEBUG_ENTRY( "mole_Init()" );
@@ -65 +67 @@
         /* these tell the molecular solver what zone and iteration it has
          * been evaluated on */
-        co.co_nzone = -2;
-        co.iteration_co = -2;
+        mole.nzone = -2;
+        mole.iteration = -2;
 
         /* prevent memory leaks */
         /* \todo        this is a temporary fix for PR14. We should improve the overall design
          * of this code to prevent valid pointers being overwritten in a second call to mole_Init */
-        if( lgCO_Init_called )
+        if( lgmole_Init_called )
         {
                 return;
@@ -77 +79 @@
 
         /* say that we have been called */
-        lgCO_Init_called = true;
+        lgmole_Init_called = true;
 
         for(nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
@@ -227 +229 @@
         /* Add passive species to complete network */
         sp = newspecies("e-    ",OTHER,MOLE_PASSIVE,&(mole.eden_f),0.00e+00, 0.0f);
-        sp->nElec = -1; sp->mole_mass = (realnum)ELECTRON_MASS; /* Augment properties for this non-molecular species */
+        sp->charge = -1;        sp->mole_mass = (realnum)ELECTRON_MASS; /* Augment properties for this non-molecular species */
         sp = newspecies("grn   ",OTHER,MOLE_PASSIVE,&(mole.grain_area),0.00e+00, 0.0f);
         sp = newspecies("PHOTON",OTHER,MOLE_PASSIVE,&one,0.00e+00, 0.0f);
@@ -317 +319 @@
                                 }
                         }
-                        if(sp->nElec < N_MOLE_ION && sp->nElec >= 0) 
-                        {
-                                element_list[nelem]->ipMl[sp->nElec] = i;
+                        if(sp->charge < N_MOLE_ION && sp->charge >= 0) 
+                        {
+                                element_list[nelem]->ipMl[sp->charge] = i;
                         }
                 }
@@ -382 +384 @@
                 mol->nElem[nelem] = 0;
         }
-        mol->nElec = mol->Excit = 0;
+        mol->charge = mol->Excit = 0;
         mol->mole_mass = 0.0;
         mol->state = state;
@@ -418 +420 @@
                                 n = 1;
                         if(*s == '+')
-                                mol->nElec = n;
+                                mol->charge = n;
                         else
-                                mol->nElec = -n;
+                                mol->charge = -n;
                         *s = '\0';
                 }
@@ -653 +655 @@
         total_molecule_elems(dense.xMolecules); 
 
+        /* add on molecules */
+        mole.elec = 0.;
+        for(i=0;i<mole.num_calc;i++)
+        {
+                if (COmole[i]->n_nuclei > 1 || COmole[i]->charge < 0)
+                        mole.elec += COmole[i]->hevmol*COmole[i]->charge;
+        }
+