Show
Ignore:
Timestamp:
01/12/08 16:50:34 (10 months ago)
Author:
rjrw
Message:

Convert molecule list and network source/sink vectors to valarrays.

Similar job on 1D vectors in nemala2.

Location:
branches/newmole/source
Files:
8 modified

Legend:

Unmodified
Added
Removed

  • branches/newmole/source/mole.h

    r1741 r1742  
    2828 
    2929extern double mole_findrate(const char buf[]); 
     30 
     31extern struct molecule *findspecies(const char buf[]); 
    3032 
    3133/** \verbatim >>chng 03 feb 09, rm ipH3P_hev, since not used, and decrement NUM_HEAVY_MOLEC to 17  
    … …  
    6062   to improve the calculation, as deep in molecular regions reactions with molecules 
    6163   can be important to the ionization balance */ 
    62  
    63 struct molecule; 
    6464 
    6565class Mole { 
    … …  
    118118        double **source , **sink; 
    119119         
    120          /** rate s-1 for molecular charge transfer, nelem from to */ 
    121         double  
    122                 *src, 
    123                 *snk; 
    124120        realnum ***xMoleChTrRate;/***[LIMELM][LIMELM+1][LIMELM+1];*/ 
    125121 
     122        valarray<struct molecule *> list; 
    126123 
    127         struct molecule **list; 
     124        valarray<struct molezone> zone; 
    128125}; 
    129126 
    130127extern Mole mole; 
     128 
     129struct molezone { 
     130        /** rate s-1 for molecular charge transfer, nelem from to */ 
     131        double src, snk; 
     132}; 
    131133 
    132134enum {CHARS_SPECIES=7}; 
    … …  
    149151        int index, groupnum; 
    150152 
    151         /* Current solution data */ 
     153        /* Current zone data */ 
    152154        double  den;       /** density (cm-3) */ 
    153155        realnum column;    /** total column density in this iteration */ 
    … …  
    159161        realnum column_old;   /** total column density in previous iteration */ 
    160162};  
    161  
    162 extern struct molecule *findspecies(const char buf[]); 
    163163 
    164164extern void mole_punch(FILE *punit, const char speciesname[],  

  • branches/newmole/source/mole_drive.cpp

    r1741 r1742  
    923923        if( trace.lgTrace && trace.lgTr_H2_Mole ) 
    924924        { 
    925                 rate = mole.snk[findspecies("H2+")->index]; 
     925                rate = mole.zone[findspecies("H2+")->index].snk; 
    926926                if( rate != 0. ) 
    927927                { 

  • branches/newmole/source/mole_eval_balance.cpp

    r1741 r1742  
    3535        double rate_tot, rate_deriv[MAXREACTANTS], rated, rk; 
    3636        char debug_name[] = ""; /* Insert name of species to print largest rates for, e.g. "CH4" */ 
    37         struct molecule *sp, *debug_species; 
     37        struct molecule *sp; 
     38        const struct molecule *debug_species = findspecies(debug_name); 
    3839        double snkx=0.,srcx=0.; 
    39  
    40         debug_species = findspecies(debug_name); 
    4140 
    4241        DEBUG_ENTRY("mole_eval_balance()"); 
    … …  
    4443        for( i=0; i < num_total; i++ ) 
    4544        { 
    46                 mole.src[i] = mole.snk[i] = 0.; 
     45                mole.zone[i].src = mole.zone[i].snk = 0.; 
    4746                if (c) 
    4847                        for( j=0; j < num_total; j++ ) 
    … …  
    101100                        { 
    102101                                sp = rate->reactants[i]; 
    103                                 mole.snk[sp->index] += rate_deriv[i]; 
     102                                mole.zone[sp->index].snk += rate_deriv[i]; 
    104103                                if (sp == debug_species)  
    105104                                { 
    … …  
    129128                        { 
    130129                                sp = rate->products[i]; 
    131                                 mole.src[sp->index] += rate_tot; 
     130                                mole.zone[sp->index].src += rate_tot; 
    132131                                if (sp == debug_species)  
    133132                                { 
    … …  
    188187        for (i=0;i<num_total;i++) 
    189188        { 
    190                 b[i] = mole.src[i]-mole.snk[i]*mole.list[i]->den; 
     189                b[i] = mole.zone[i].src-mole.zone[i].snk*mole.list[i]->den; 
    191190        } 
    192191         
    … …  
    196195                { 
    197196                        fprintf(stdout,"%20.20s src %13.7g of %13.7g [", 
    198                                                         ratesrc->label,srcx,mole.src[debug_species->index]); 
     197                                                        ratesrc->label,srcx,mole.zone[debug_species->index].src); 
    199198                        for (j=0;j<ratesrc->nreactants;j++)  
    200199                        { 
    … …  
    206205                { 
    207206                        fprintf(stdout,"%20.20s snk %13.7g of %13.7g [", 
    208                                                         ratesnk->label,snkx,mole.snk[debug_species->index]); 
     207                                                        ratesnk->label,snkx,mole.zone[debug_species->index].snk); 
    209208                        for (j=0;j<ratesnk->nreactants;j++)  
    210209                        { 

  • branches/newmole/source/mole_solve.cpp

    r1741 r1742  
    124124                        if (element_list[nelem]->ipMl[ion] != -1) 
    125125                        { 
    126                                 mole.source[nelem][ion] = mole.src[element_list[nelem]->ipMl[ion]]; 
    127                                 mole.sink[nelem][ion] = mole.snk[element_list[nelem]->ipMl[ion]]; 
     126                                mole.source[nelem][ion] = mole.zone[element_list[nelem]->ipMl[ion]].src; 
     127                                mole.sink[nelem][ion] = mole.zone[element_list[nelem]->ipMl[ion]].snk; 
    128128                                if (dense.IonLow[nelem] > ion) 
    129129                                        dense.IonLow[nelem] = ion; 

  • branches/newmole/source/mole_species.cpp

    r1741 r1742  
    4646struct chem_element_s *element_list[LIMELM]; 
    4747struct molecule **groupspecies; 
     48 
     49#include <functional> 
     50 
     51class MoleCmp : public binary_function<const struct molecule *,const struct molecule *,bool> 
     52{ 
     53public: 
     54        bool operator()(const struct molecule *mol1, const struct molecule *mol2) const 
     55        { 
     56                for (long nelem=LIMELM-1; nelem >= ipHYDROGEN; nelem--)  
     57                { 
     58                        if (mol1->nElem[nelem] > mol2->nElem[nelem]) 
     59                                return false; 
     60                        else if (mol1->nElem[nelem] < mol2->nElem[nelem]) 
     61                                return true;             
     62                }  
     63                return strcmp(mol1->label,mol2->label) < 0; 
     64        } 
     65}; 
    4866 
    4967Mole mole; 
    … …  
    308326 
    309327        /* Create linear list of species and populate it... */ 
    310         mole.list = (struct molecule **)MALLOC((size_t)mole.num_total* 
    311                                                                                                                                                  sizeof(struct molecule *)); 
     328        mole.list.resize(mole.num_total); 
    312329 
    313330        /* ...first active species */ 
    … …  
    322339 
    323340        /* Sort list into a standard ordering */ 
    324         qsort((void *)mole.list,(size_t)mole.num_calc,sizeof(struct molecule *), 
    325                                 mole_cmp); 
     341        sort(&mole.list[0],&mole.list[0]+mole.num_calc,MoleCmp()); 
    326342 
    327343        for (molecule_i p 
    … …  
    359375        } 
    360376 
    361         mole.src = (double *)MALLOC((size_t)mole.num_total*sizeof(double)); 
    362         mole.snk = (double *)MALLOC((size_t)mole.num_total*sizeof(double)); 
     377        mole.zone.resize(mole.num_total); 
    363378 
    364379        make_groups(); 

  • branches/newmole/source/nemala2.cpp

    r1741 r1742  
    2828        { 
    2929                const char *spName = Species[i].chptrSpName; 
    30                 double *g, *ex, *pops, *depart; 
    3130                double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate; 
    32                 double *source, *sink; 
     31                long nlev = Species[i].numLevels_max; 
     32                valarray<double> g(nlev), ex(nlev), pops(nlev), depart(nlev); 
     33                valarray<double> source(nlev), sink(nlev); 
    3334                double cooltl, coolder; 
    3435                double fupsilon; 
    … …  
    3637                int nNegPop; 
    3738                bool lgZeroPop, lgDeBug = false; 
    38                 long nlev = Species[i].numLevels_max; 
    3939 
    4040                /* first find current density (cm-3) of species */ 
    … …  
    7878                } 
    7979 
    80                 g = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 
    81                 ex = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 
    82                 pops = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 
    83                 depart = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 
    84                 source = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 
    85                 sink = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double)); 
    86  
    8780                AulEscp = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));  
    8881                col_str = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));  
    … …  
    248241                        abund,  
    249242                        /* G(nlev) is stat weight of levels */ 
    250                         g,  
     243                        &g[0],  
    251244                        /* EX(nlev) is excitation potential of levels, deg K or wavenumbers 
    252245                         * 0 for lowest level, all are energy rel to ground NOT d(ENER) */ 
    253                         ex,  
     246                        &ex[0],  
    254247                        /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */ 
    255248                        'w', 
    256249                        /* populations [cm-3] of each level as deduced here */ 
    257                         pops,  
     250                        &pops[0],  
    258251                        /* departure coefficient, derived below */ 
    259                         depart, 
     252                        &depart[0], 
    260253                        /* net transition rate, A * esc prob, s-1 */ 
    261254                        &AulEscp,  
    … …  
    272265                        &CollRate, 
    273266                        /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */ 
    274                         source, 
     267                        &source[0], 
    275268                        /* this is an additional destruction rate to continuum, normally zero, units s-1 */ 
    276                         sink, 
     269                        &sink[0], 
    277270                        /* flag saying whether CollRate already done (true), or we need to do it here (false), 
    278271                         * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/ 
    … …  
    340333                                { 
    341334                                        /*Solving for level populations*/ 
    342                                         atom_levelN(nlev,abund,g,ex,'w',pops,depart,&AulEscp,&col_str, 
    343                                                                         &AulDest, &AulPump, &CollRate, source, sink,true,&cooltl,  
     335                                        atom_levelN(nlev,abund,&g[0],&ex[0],'w',&pops[0],&depart[0],&AulEscp,&col_str, 
     336                                                                        &AulDest, &AulPump, &CollRate, &source[0], &sink[0],true,&cooltl,  
    344337                                                                        &coolder, spName, &nNegPop,     &lgZeroPop, lgDeBug ); 
    345338                                         
    … …  
    393386                        } 
    394387                } 
    395  
    396                 /* free vectors */ 
    397                 free( g ); 
    398                 free( ex ); 
    399                 free( pops ); 
    400                 free( depart ); 
    401                 free( source ); 
    402                 free( sink ); 
    403388 
    404389                /* free arrays */ 

  • branches/newmole/source/punch_do.cpp

    r1741 r1742  
    263263                                          timesc.time_H2_Dest_here,  
    264264                                          /* CO destruction timescale */ 
    265                                           1./SDIV((ipCO != -1) ? mole.snk[ipCO] : 0.),  
     265                                          1./SDIV((ipCO != -1) ? mole.zone[ipCO].snk : 0.),  
    266266                                          /* OH destruction timescale */ 
    267                                           1./SDIV((ipOH != -1) ? mole.snk[ipOH] : 0.),  
     267                                          1./SDIV((ipOH != -1) ? mole.zone[ipOH].snk : 0.),  
    268268                                          /* H recombination timescale */ 
    269269                                          1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) ); 

  • branches/newmole/source/radius_increment.cpp

    r1741 r1742  
    253253                if (ipCO != -1) 
    254254                        timesc.BigCOMoleForm = MAX2(  timesc.BigCOMoleForm ,  
    255                                                                                                                                                 1./SDIV(mole.snk[ipCO] )); 
     255                                                                                                                                                1./SDIV(mole.zone[ipCO].snk )); 
    256256        } 
    257257