Changeset 1712 for branches/newmole/source

Timestamp:

12/22/07 06:26:02 (11 months ago)

Author:

rjrw

Message:

Mask out grain surface chemistry (was causing convergence failures,
need to investigate further).

Improve diagnostic prints slightly.

Move return of molecular network results to non-network locations out
through the hierarchy (from end of each Newton solver step towards
where it will actually be used).

Location:

branches/newmole/source

Files:

• dense.cpp (modified) (2 diffs)
• mole_eval_balance.cpp (modified) (2 diffs)
• mole_h2_io.cpp (modified) (2 diffs)
• mole_newton_step.cpp (modified) (1 diff)
• mole_solve.cpp (modified) (3 diffs)
• mole_species.cpp (modified) (1 diff)

branches/newmole/source/dense.cpp

@@ -18 +18 @@
                         /* this sum of over the molecules did not include the atom and first
                          * ion that is calculated in the molecular solver */
-                        double sum = dense.xMolecules[nelem];
+                        double sum = 0.;
                         for( ion=0; ion<nelem+2; ++ion )
                         {
                                 sum += dense.xIonDense[nelem][ion];
                         }
-                        if (sum <= SMALLFLOAT && dense.gas_phase[nelem] > SMALLFLOAT)
+                        if (sum+dense.xMolecules[nelem] <= SMALLFLOAT && 
+                                        dense.gas_phase[nelem] > SMALLFLOAT)
                         {
-                                fprintf(ioQQQ,"Problem total density %g of %g nelem %ld(%s)\n",
-                                                                sum,dense.gas_phase[nelem],nelem,
+                                fprintf(ioQQQ,"Problem ions %g moles %g error %g of %g nelem %ld(%s)\n",
+                                                                sum,dense.xMolecules[nelem],
+                                                                sum+dense.xMolecules[nelem]-dense.gas_phase[nelem],
+                                                                dense.gas_phase[nelem],nelem,
                                                                 elementnames.chElementSym[nelem]);
                                 lgOK=false;
@@ -41 +44 @@
                         /* >>chng 04 jul 23, needed to increase to 2e-2 from 1e-2 to get conserve.in to work,
                          * not clear what caused increase in error */
-                        if( fabs(sum-dense.gas_phase[nelem])/sum > 2e-5 )
+                        if( fabs(sum+dense.xMolecules[nelem]-dense.gas_phase[nelem]) > 
+                                        2e-5*dense.gas_phase[nelem] )
                         {
+                                fprintf(ioQQQ,"Problem ions %g moles %g error %g of %g nelem %ld(%s)\n",
+                                                                sum,dense.xMolecules[nelem],
+                                                                sum+dense.xMolecules[nelem]-dense.gas_phase[nelem],
+                                                                dense.gas_phase[nelem],nelem,
+                                                                elementnames.chElementSym[nelem]);
                                 fprintf(ioQQQ,"DEBUG non-conv nelem\t%li\tsum\t%e\ttot gas\t%e\trel err\t%.3e\n",
                                                                 nelem,

branches/newmole/source/mole_eval_balance.cpp

@@ -29 +29 @@
 /*=================================================================*/
 
-#define FABSSWITCH(a) (a)
-
 void mole_eval_balance(long int num_total, double *b, double **c)
 {
@@ -91 +89 @@
                                 if(i!=j)
                                 {
-                                        rate_deriv[i] *= FABSSWITCH(rate->reactants[j]->hevmol);
-                                }
-                        }
-                }
-                
-                rate_tot = rate_deriv[0]*FABSSWITCH(rate->reactants[0]->hevmol);
+                                        rate_deriv[i] *= rate->reactants[j]->hevmol;
+                                }
+                        }
+                }
+                
+                rate_tot = rate_deriv[0]*rate->reactants[0]->hevmol;
                 
                 for(i=0;i<rate->nreactants;i++)

branches/newmole/source/mole_h2_io.cpp

@@ -1818 +1818 @@
                         mole_findrate("H2,OH=>H2O,H") / hmi.H2_rate_destroy / hmi.H2_total,
                         mole_findrate("H2,O=>OH,H") / hmi.H2_rate_destroy / hmi.H2_total,
-                        mole_findrate("H2*,CH=>CH2,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,
-                        mole_findrate("H2*,O=>OH,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,
-                        mole_findrate("H2*,OH=>H2O,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,
-                        mole_findrate("H2*,C=>CH,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,    
-                        mole_findrate("H2*,C+=>CH+,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,
+                        mole_findrate("H2*,CH=>CH2,H") /  hmi.H2_rate_destroy / hmi.H2_total,
+                        mole_findrate("H2*,O=>OH,H") /  hmi.H2_rate_destroy / hmi.H2_total,
+                        mole_findrate("H2*,OH=>H2O,H") /  hmi.H2_rate_destroy / hmi.H2_total,
+                        mole_findrate("H2*,C=>CH,H") /  hmi.H2_rate_destroy / hmi.H2_total,     
+                        mole_findrate("H2*,C+=>CH+,H") /  hmi.H2_rate_destroy / hmi.H2_total,
                         mole_findrate("H2,CH2=>CH3,H") / hmi.H2_rate_destroy / hmi.H2_total,
                         mole_findrate("H2,CH3=>CH4,H") / hmi.H2_rate_destroy / hmi.H2_total,
                         mole_findrate("H2,CH4+=>CH5+,H") / hmi.H2_rate_destroy / hmi.H2_total,
-                        mole_findrate("H2*,CH2=>CH3,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,
-                        mole_findrate("H2*,CH3=>CH4,H")*hmi.lgLeiden_Keep_ipMH2s /  hmi.H2_rate_destroy / hmi.H2_total,
+                        mole_findrate("H2*,CH2=>CH3,H") /  hmi.H2_rate_destroy / hmi.H2_total,
+                        mole_findrate("H2*,CH3=>CH4,H") /  hmi.H2_rate_destroy / hmi.H2_total,
                         mole_findrate("H2*,O+=>OH+,H") / hmi.H2_rate_destroy / hmi.H2_total,
                         /*heavy elements*/
@@ -1846 +1846 @@
                         mole_findrate("H2,HCN+=>HCNH+,H") / hmi.H2_rate_destroy  / hmi.H2_total
                         );
+                        
+    if (hmi.lgLeiden_Keep_ipMH2s == false)
+                        ASSERT(mole_findrate("H2*,CH=>CH2,H") == 0 &&
+                                                 mole_findrate("H2*,O=>OH,H") == 0 &&
+                                                 mole_findrate("H2*,OH=>H2O,H") == 0 &&
+                                                 mole_findrate("H2*,C=>CH,H") == 0 &&
+                                                 mole_findrate("H2*,C+=>CH+,H") == 0 &&
+                                                 mole_findrate("H2*,CH2=>CH3,H") == 0 &&
+                                                 mole_findrate("H2*,CH3=>CH4,H") == 0 &&
+                                                 mole_findrate("H2*,O+=>OH+,H") == 0    );
                 fprintf(io,"\t%.4e\t%.4e\n",
                         /*H2g/Htot, H2s/Htot chemical network and from big molecule model*/

branches/newmole/source/mole_newton_step.cpp

@@ -491 +491 @@
         *error = (realnum) MIN2(error0,1e30);
 
-        mole_return_cached_species();
-        
         return;
 }

branches/newmole/source/mole_solve.cpp

@@ -91 +91 @@
                 
                 //fprintf(stdout,"Mole zone %ld -- %7.2f loop %d error %15.8g nBad %d\n",nzone,fnzone,i,error,nBad);
-
-                {
-                        /* option to print deduced abundances */
-                        /*@-redef@*/
-                        enum {DEBUG_LOC=false};
-                        /*@+redef@*/
-                        if( DEBUG_LOC /*&& (nzone>68)*/ )
-                        {
-                                fprintf(ioQQQ,"DEBUG mole out\t%i\t%.2f\t%.5e\t%.5e",
-                                        i,
-                                        fnzone,
-                                        phycon.te,
-                                        dense.eden);
-                                fprintf(ioQQQ,
-                                        "\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e", 
-                                        error,
-                                        dense.xIonDense[ipHYDROGEN][0],
-                                        dense.xIonDense[ipHYDROGEN][1],
-                                        mole.sink[ipHYDROGEN][0],
-                                        mole.sink[ipHYDROGEN][1],
-                                        mole.source[ipHYDROGEN][0] , 
-                                        mole.source[ipHYDROGEN][1] ,
-                                        ionbal.RateIonizTot[ipHYDROGEN][0],
-                                        ionbal.RateRecomTot[ipHYDROGEN][0]);
-
-                                /*for( j=0; j<mole.num_calc; j++ )
-                                        fprintf(ioQQQ,"\t%.4e", HMOLEC(j) );*/
-                                fprintf(ioQQQ,"\n" );
-                        }
-                }
-
+        
                 if (error < BIGERROR && nBad == 0 && nPrevBad == 0) /* Stop early if good enough */
                         break;
@@ -144 +114 @@
         }
 
-        ASSERT(lgElemsConserved());
-
+        mole_return_cached_species();
+        
         check_co_ion_converge();
 
@@ -169 +139 @@
         }
 
+        {
+                /* option to print deduced abundances */
+                /*@-redef@*/
+                enum {DEBUG_LOC=false};
+                /*@+redef@*/
+                if( DEBUG_LOC /*&& (nzone>68)*/ )
+                {
+                        fprintf(ioQQQ,"DEBUG mole out\t%i\t%.2f\t%.5e\t%.5e",
+                                                        i,
+                                                        fnzone,
+                                                        phycon.te,
+                                                        dense.eden);
+                        fprintf(ioQQQ,
+                                                        "\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e", 
+                                                        error,
+                                                        dense.xIonDense[ipHYDROGEN][0],
+                                                        dense.xIonDense[ipHYDROGEN][1],
+                                                        mole.sink[ipHYDROGEN][0],
+                                                        mole.sink[ipHYDROGEN][1],
+                                                        mole.source[ipHYDROGEN][0] , 
+                                                        mole.source[ipHYDROGEN][1] ,
+                                                        ionbal.RateIonizTot[ipHYDROGEN][0],
+                                                        ionbal.RateRecomTot[ipHYDROGEN][0]);
+                        
+                        /*for( j=0; j<mole.num_calc; j++ )
+                                fprintf(ioQQQ,"\t%.4e", HMOLEC(j) );*/
+                        fprintf(ioQQQ,"\n" );
+                }
+        }
+        
         return error;
 

branches/newmole/source/mole_species.cpp

@@ -212 +212 @@
                 sp = newspecies("H2Ogrn",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, -238.9f);
                 sp = newspecies("OHgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 38.4f);
-                sp = newspecies("Cgrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 711.2f);
-                sp = newspecies("Ogrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 246.8f);
-                sp = newspecies("Ngrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 470.8f);
-                sp = newspecies("Sgrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 274.7f);
-                sp = newspecies("Hgrn  ",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 216.0f); 
-                sp = newspecies("Sigrn ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 446.0f);
-                sp = newspecies("SiHgrn",MOLECULE,MOLE_ACTIVE,NULL,1.80e-14, 374.9f);
-                sp = newspecies("SiOgrn",MOLECULE,MOLE_ACTIVE,NULL,3.89e-14, -101.6f);
-                sp = newspecies("CHgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.10e-06, 592.5f);
-                sp = newspecies("O2grn ",MOLECULE,MOLE_ACTIVE,NULL,1.91e-07, 0.0f);
-                sp = newspecies("CH2grn",MOLECULE,MOLE_ACTIVE,NULL,8.59e-13, 390.0f);
-                sp = newspecies("CH3grn",MOLECULE,MOLE_ACTIVE,NULL,4.59e-19, 149.0f);
-                sp = newspecies("CH4grn",MOLECULE,MOLE_ACTIVE,NULL,9.12e-28, -66.8f);
-                sp = newspecies("N2grn ",MOLECULE,MOLE_ACTIVE,NULL,2.46e-17, 0.0f);
-                sp = newspecies("NOgrn ",MOLECULE,MOLE_ACTIVE,NULL,5.99e-12, 89.8f);
-                sp = newspecies("S2grn ",MOLECULE,MOLE_ACTIVE,NULL,1.21e-11, 128.3f);
-                sp = newspecies("OCNgrn",MOLECULE,MOLE_ACTIVE,NULL,3.86e-11, 556.0f);
-                sp = newspecies("NHgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.30e-09, 376.5f);
-                sp = newspecies("NH2grn",MOLECULE,MOLE_ACTIVE,NULL,6.78e-20, 191.6f);
-                sp = newspecies("NH3grn",MOLECULE,MOLE_ACTIVE,NULL,1.25e-29, -38.9f);
-                sp = newspecies("CNgrn ",MOLECULE,MOLE_ACTIVE,NULL,3.38e-11, 436.8f);
-                sp = newspecies("HCNgrn",MOLECULE,MOLE_ACTIVE,NULL,1.07e-17, 135.5f);
-                sp = newspecies("HNOgrn",MOLECULE,MOLE_ACTIVE,NULL,9.09e-21, 100.0f);
-                sp = newspecies("HSgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.71e-14, 136.5f);
-                sp = newspecies("CSgrn ",MOLECULE,MOLE_ACTIVE,NULL,6.69e-18, 277.1f);
-                sp = newspecies("NSgrn ",MOLECULE,MOLE_ACTIVE,NULL,3.12e-11, 263.0f);
-                sp = newspecies("SOgrn ",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 5.0f);
-                sp = newspecies("OCSgrn",MOLECULE,MOLE_ACTIVE,NULL,8.83e-24, -142.0f);
-                sp = newspecies("HNCgrn",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 201.0f);
-                sp = newspecies("C2grn ",MOLECULE,MOLE_ACTIVE,NULL,3.71e-09, 817.0f);
-                sp = newspecies("C2Hgrn",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 560.0f);
-                sp = newspecies("C3grn ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 831.0f);
-                sp = newspecies("C3Hgrn",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 602.5f);
-                sp = newspecies("H2grn ",MOLECULE,MOLE_ACTIVE,NULL,4.00E-08, 0.0f);
+                if (0)
+                {
+                        sp = newspecies("Cgrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 711.2f);
+                        sp = newspecies("Ogrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 246.8f);
+                        sp = newspecies("Ngrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 470.8f);
+                        sp = newspecies("Sgrn  ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 274.7f);
+                        sp = newspecies("Hgrn  ",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 216.0f); 
+                        sp = newspecies("Sigrn ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 446.0f);
+                        sp = newspecies("SiHgrn",MOLECULE,MOLE_ACTIVE,NULL,1.80e-14, 374.9f);
+                        sp = newspecies("SiOgrn",MOLECULE,MOLE_ACTIVE,NULL,3.89e-14, -101.6f);
+                        sp = newspecies("CHgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.10e-06, 592.5f);
+                        sp = newspecies("O2grn ",MOLECULE,MOLE_ACTIVE,NULL,1.91e-07, 0.0f);
+                        sp = newspecies("CH2grn",MOLECULE,MOLE_ACTIVE,NULL,8.59e-13, 390.0f);
+                        sp = newspecies("CH3grn",MOLECULE,MOLE_ACTIVE,NULL,4.59e-19, 149.0f);
+                        sp = newspecies("CH4grn",MOLECULE,MOLE_ACTIVE,NULL,9.12e-28, -66.8f);
+                        sp = newspecies("N2grn ",MOLECULE,MOLE_ACTIVE,NULL,2.46e-17, 0.0f);
+                        sp = newspecies("NOgrn ",MOLECULE,MOLE_ACTIVE,NULL,5.99e-12, 89.8f);
+                        sp = newspecies("S2grn ",MOLECULE,MOLE_ACTIVE,NULL,1.21e-11, 128.3f);
+                        sp = newspecies("OCNgrn",MOLECULE,MOLE_ACTIVE,NULL,3.86e-11, 556.0f);
+                        sp = newspecies("NHgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.30e-09, 376.5f);
+                        sp = newspecies("NH2grn",MOLECULE,MOLE_ACTIVE,NULL,6.78e-20, 191.6f);
+                        sp = newspecies("NH3grn",MOLECULE,MOLE_ACTIVE,NULL,1.25e-29, -38.9f);
+                        sp = newspecies("CNgrn ",MOLECULE,MOLE_ACTIVE,NULL,3.38e-11, 436.8f);
+                        sp = newspecies("HCNgrn",MOLECULE,MOLE_ACTIVE,NULL,1.07e-17, 135.5f);
+                        sp = newspecies("HNOgrn",MOLECULE,MOLE_ACTIVE,NULL,9.09e-21, 100.0f);
+                        sp = newspecies("HSgrn ",MOLECULE,MOLE_ACTIVE,NULL,1.71e-14, 136.5f);
+                        sp = newspecies("CSgrn ",MOLECULE,MOLE_ACTIVE,NULL,6.69e-18, 277.1f);
+                        sp = newspecies("NSgrn ",MOLECULE,MOLE_ACTIVE,NULL,3.12e-11, 263.0f);
+                        sp = newspecies("SOgrn ",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 5.0f);
+                        sp = newspecies("OCSgrn",MOLECULE,MOLE_ACTIVE,NULL,8.83e-24, -142.0f);
+                        sp = newspecies("HNCgrn",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 201.0f);
+                        sp = newspecies("C2grn ",MOLECULE,MOLE_ACTIVE,NULL,3.71e-09, 817.0f);
+                        sp = newspecies("C2Hgrn",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 560.0f);
+                        sp = newspecies("C3grn ",MOLECULE,MOLE_ACTIVE,NULL,1.00e-15, 831.0f);
+                        sp = newspecies("C3Hgrn",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 602.5f);
+                        sp = newspecies("H2grn ",MOLECULE,MOLE_ACTIVE,NULL,4.00E-08, 0.0f);
+                }
         }
         sp = newspecies("C2H   ",MOLECULE,MOLE_ACTIVE,NULL,4.00e-15, 560.0f);