Changeset 1600

Timestamp:

12/01/07 17:09:37 (12 months ago)

Author:

rjrw

Message:

Remove some commented out stuff,

change a conservation test that should no longer be active to an
ASSERT,

instrument convergence of S ion in conv_base which seems a good
bellweather for pdr_co_fully.

Location:

branches/newmole/source

Files:

• conv_base.cpp (modified) (5 diffs)
• ion_solver.cpp (modified) (3 diffs)
• mole_co_solve.cpp (modified) (3 diffs)
• mole_h_step.cpp (modified) (1 diff)

branches/newmole/source/conv_base.cpp

@@ -378 +378 @@
         {
 
-                fprintf(ioQQQ,"Start of loop %g\n",dense.xIonDense[15][0]);
                 check_conserve(__FILE__,__LINE__);
         
@@ -413 +412 @@
                 /*fprintf(ioQQQ,"DEBUG h2\t%.2f\t%.3e\t%.3e", fnzone,hmi.H2_total,findspecies("H2")->hevmol);*/
                 iso_solve( ipH_LIKE );
-                fprintf(ioQQQ,"After H_LIKE %g\n",dense.xIonDense[15][0]);
 
                 /* now do level populations for H2 */
@@ -434 +432 @@
                 /* do all he-like species */
                 iso_solve( ipHE_LIKE );
-                fprintf(ioQQQ,"After He_like %g\n",dense.xIonDense[15][0]);
 
                 /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures
@@ -509 +506 @@
                 IonOxyge();
                 IonNitro();
-                fprintf(ioQQQ,"Before Si %g\n",dense.xIonDense[15][0]);
                 IonSilic();
-                fprintf(ioQQQ,"Before S %g\n",dense.xIonDense[15][0]);
+#define DEBUGRW 1
+                if (DEBUGRW) fprintf(ioQQQ," DEBUG Before S %g\n",dense.xIonDense[15][0]);
                 IonSulph();
-                fprintf(ioQQQ,"After ion %g\n",dense.xIonDense[15][0]);
+                if (DEBUGRW) fprintf(ioQQQ," DEBUG After S %g\n",dense.xIonDense[15][0]);
                 IonChlor();
                 check_conserve(__FILE__,__LINE__);
@@ -519 +516 @@
                 /* do carbon monoxide after oxygen */
                 CO_drive();
-                fprintf(ioQQQ,"After mole %g\n",dense.xIonDense[15][0]);
+                if (DEBUGRW) fprintf(ioQQQ," DEBUG After mole %g\n",dense.xIonDense[15][0]);
+#undef DEBUGRW
 
                 check_conserve(__FILE__,__LINE__);

branches/newmole/source/ion_solver.cpp

@@ -64 +64 @@
         abund_total = SDIV ( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
 
-        /* >>chng 04 nov 22,
-         * if gas nearly all atomic/ionic do not let source - sink terms from 
-         * molecular network force system of balance equations to become 
-         * inhomogeneous
-         * what constitutes a source or sink IS DIFFERENT for hydrogen and the rest 
-         * the H solution must couple with hmole - and its defn of source and sink.  For instance, oxygen charge
-         * transfer goes into source and sink terms for hydrogen.  So we never hose source and sink for H.
-         * for the heavy elements, which couple onto comole, mole.source and sink represent terms that
-         * remove atoms and ions from the ionization ladder.  their presence makes the system of
-         * equations inhomogeneous.  we don't want to do this when comole has a trivial effect on
-         * the ionization balance since the matrix becomes unstable */
-        /* >>chng 04 dec 06, limit from 0.01 to 1e-10 as per NA suggestion */
-#if 0
-        if( nelem>ipHYDROGEN && dense.xMolecules[nelem]/SDIV(dense.gas_phase[nelem]) < 1e-10 )
-        {
-                for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
-                {
-                        mole.source[nelem][ion] = 0.;
-                        mole.sink[nelem][ion] = 0.;
-                }
-        }
-
-        /* protect against case where all gas phase abundances are in molecules, use the
-         * atomic and first ion density from the molecule solver 
-         * >>chng 04 aug 15, NA change from 10 0000 to 10 pre-coef on
-         * FLT_EPSILON for stability in PDR 
-         * the factor 10.*FLT_EPSILON also appears in mole_h_step in total H 
-         * conservation */
-        if( fabs( dense.gas_phase[nelem] -  dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) <
-                10.*FLT_EPSILON )
-        {
-                double renorm;
-                /* >>chng 04 jul 31, add logic to conserve nuclei in fully molecular limit;
-                 * in first calls, when searching for solution, we may be VERY far 
-                 * off, and sum of first ion and atom density may be far larger 
-                 * than difference between total gas and molecular densities,
-                 * since they reflect the previous evaluation of the solution.  Do 
-                 * renorm to cover this case */
-                /* first form sum of all atoms and ions */
-                realnum sum = 0.;
-                for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
-                        sum += dense.xIonDense[nelem][ion];
-                /* now renorm to this sum - this should be unity, and is not if we have
-                 * now conserved particles, due to molecular fraction changing */
-                renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] );
-
-                abund_total = renorm * sum;
-        }
-
-        /* negative residual density */
-        if( abund_total < 0. )
-        {
-                /* >>chng 05 dec 16, do not abort when net atomic/ ionic abundance 
-                 * is negative due to chem net having too much of a species - this 
-                 * happens naturally when ices become well formed, but the code will 
-                 * converge away from it.  Rather, we set the ionization
-                 * convergence flag to "not converge" and soldier on 
-                 * if negative populations do not go away, we will eventually 
-                 * terminate due to convergence failures */
-                /* print error if not search */
-                if( !conv.lgSearch )
-                        fprintf(ioQQQ,
-                                "PROBLEM ion_solver - neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
-                                nelem,
-                                fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
-                                conv.nTotalIoniz );
-                /* fix up is to use half the positive abundance, assuming chem is 
-                 * trying to get too much of this species */
-                abund_total = -abund_total/2.;
-
-                /* say that ionization is not converged - do not abort - but if 
-                 * cannot converge away from negative solution, this will become a 
-                 * convergence failure abort */
-                conv.lgConvIoniz = false;
-                strcpy( conv.chConvIoniz, "neg ion" );
-        }
-#endif
-
         if (dense.gas_phase[nelem] <  dense.xMolecules[nelem]) 
         {
-                fprintf(stdout,"Bad density %li %g %g\n",nelem,dense.gas_phase[nelem],dense.xMolecules[nelem]);
+                fprintf(stdout," DISASTER ion_solver: Bad density %li %g %g\n",nelem,dense.gas_phase[nelem],dense.xMolecules[nelem]);
+                puts( "[Stop in ion_solver]" );
+                cdEXIT(EXIT_FAILURE);
         }
 
@@ -455 +379 @@
 
                         /* these are the external source and sink terms */
-                        /* source first */
                         /* need negative sign to get positive pops */
                         source[i] -= mole.source[nelem][ion];
-
                         totsrc += mole.source[nelem][ion];
-                        /* sink next */
-                        /*MAT( xmat, i, i ) += mole.sink[nelem][ion];*/
                         MAT( xmat, i, i ) -= mole.sink[nelem][ion];
                 }
@@ -478 +398 @@
                 {                        
                         ion = i+ion_low;
-                        /*MAT( xmat, i, i ) += dynamics.Rate;*/
                         MAT( xmat, i, i ) -= dynamics.Rate;
-                        /*src[i] += dynamics.Source[nelem][ion];*/
                         source[i] -= dynamics.Source[nelem][ion];
                         /* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(xmat,i,i),

branches/newmole/source/mole_co_solve.cpp

@@ -73 +73 @@
                                 mole.source[nelem][ion] = mole.src[element_list[nelem]->ipMl[ion]];
                                 mole.sink[nelem][ion] = mole.snk[element_list[nelem]->ipMl[ion]];
-                                // dense.xIonDense[nelem][ion] = (realnum) COmole[element_list[nelem]->ipMl[ion]]->hevmol;
+                                if (dense.IonLow[nelem] > ion)
+                                        dense.IonLow[nelem] = ion;
+                                if (dense.IonHigh[nelem] < ion)
+                                        dense.IonHigh[nelem] = ion;                                     
                         }
                         else
@@ -81 +84 @@
                         }
                 }
-                if (element_list[nelem]->ipMl[0] != -1)
-                {
-                        dense.IonLow[nelem] = 0;
-                        dense.IonHigh[nelem] = MAX2( dense.IonHigh[nelem] , 1 );
-                }
         }
 
         /* if populations not conserved then not converged */
-#       define EPS_MOLE 0.1
-        if( dense.lgElmtOn[ipHYDROGEN] &&
-                dense.xMolecules[ipHYDROGEN] > dense.gas_phase[ipHYDROGEN]*(1.+EPS_MOLE) )
-        {
-                conv.lgConvIoniz = false;
-                sprintf( conv.chConvIoniz, "COcon%2i",ipHYDROGEN );
-                conv.BadConvIoniz[0] = dense.xMolecules[ipHYDROGEN];
-                conv.BadConvIoniz[1] = dense.gas_phase[ipHYDROGEN];
-        }
-        else
-        {
-                for (nelem=0;nelem<LIMELM; ++nelem)
-                {
-                        element = element_list[nelem];
-                        if( dense.lgElmtOn[nelem] &&
-                                        dense.xMolecules[nelem] > dense.gas_phase[nelem]*(1.+EPS_MOLE) )
-                        {
-                                conv.lgConvIoniz = false;
-                                sprintf( conv.chConvIoniz, "COcon%2li",nelem );
-                                conv.BadConvIoniz[0] = dense.xMolecules[nelem];
-                                conv.BadConvIoniz[1] = dense.gas_phase[nelem];
-                                break;
-                        }
-                }
-        }
+#       define EPS_MOLE 0.01
+        for (nelem=0;nelem<LIMELM; ++nelem)
+        {
+                element = element_list[nelem];
+                if( dense.lgElmtOn[nelem] )
+                        ASSERT(dense.xMolecules[nelem] <= dense.gas_phase[nelem]*(1.+EPS_MOLE) );
+        }
+#       undef EPS_MOLE
 
         DEBUG_EXIT( "CO_solve()" );
@@ -122 +103 @@
 }
 
-#       undef EPS_MOLE
 
 // #define SVD

branches/newmole/source/mole_h_step.cpp

@@ -613 +613 @@
 
         CO_return_cached_species();
+
         hmi.H2_total = 0.;
         /* this is where the transition struc expects to find the H2 abundance */