Changeset 684

Timestamp:

12/05/06 13:55:14 (2 years ago)

Author:

rjrw

Message:

Patch back to deal with SEGVs on previous step.

Still need to a) reorganize startup b) improve treatment of
decoupling within the network.

Location:

branches/newmole/source

Files:

• mole_co_etc.cpp (modified) (1 diff)
• mole_co_reactions.cpp (modified) (1 diff)
• mole_co_solve.cpp (modified) (1 diff)

branches/newmole/source/mole_co_etc.cpp

@@ -189 +189 @@
         sp = newspecies("H2CCl+",MOLECULE,ACTIVE,NULL,0.00e+00);
         sp = newspecies("ClO+  ",MOLECULE,ACTIVE,NULL,4.00e-15);
-        if( gv.lgDustOn && mole.lgGrain_mole_deplete )
+        if(1 || (gv.lgDustOn && mole.lgGrain_mole_deplete ))
         {
                 sp = newspecies("COgrn ",MOLECULE,ACTIVE,NULL,1.00e-15);

branches/newmole/source/mole_co_reactions.cpp

@@ -1686 +1686 @@
 
         /* >>chng 06 feb 28, turn off this rate when no grain molecules */
-        if (mole.lgGrain_mole_deplete)
+        /* >>chng 06 dec 05 rjrw: do this in newreact rather than rate */
+        for( nd=0; nd < gv.nBin; nd++ )
         {
-                for( nd=0; nd < gv.nBin; nd++ )
-                {
-                        /* >>chng 06 mar 04, update expression for projected grain surface area, PvH */
-                        den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
-                }
+                /* >>chng 06 mar 04, update expression for projected grain surface area, PvH */
+                den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
         }
 

branches/newmole/source/mole_co_solve.cpp

@@ -50 +50 @@
         float abundan;
         struct chem_element_s *element;
+        double sum;
 
         CO_step();                  /* Calculate the matrix elements */
+
+        /* Ugly hack to deal with species which have become uncoupled */
+        for (i=0;i<mole.num_comole_calc;i++)
+        {
+                sum = 0.;
+                for (j=0;j<mole.num_comole_calc;j++)
+                {
+                        sum = sum+fabs(mole.c[i][j]);
+                }
+                if (sum < SMALLFLOAT && fabs(mole.b[i]) < SMALLFLOAT)
+                {
+                        mole.c[i][i] = 1.;
+                }
+        }
 
         cartot_mol = dense.xMolecules[ipCARBON] + findspecies("C")->hevmol + findspecies("C+")->hevmol;