Changeset 1680 for branches/newmole/source

Timestamp:

12/16/07 03:39:51 (11 months ago)

Author:

rjrw

Message:

Fix lines_supra.in (absolute values for abundances in rate matrix).

Enable abundance asserts for dynamics runs -- need to get these more
accurate.

Location:

branches/newmole/source

Files:

• dense.cpp (modified) (1 diff)
• mole_eval_balance.cpp (modified) (1 diff)
• mole_newton_step.cpp (modified) (1 diff)
• pressure_total.cpp (modified) (3 diffs)
• rt_line_one_tauinc.cpp (modified) (1 diff)

branches/newmole/source/dense.cpp

@@ -7 +7 @@
         long int nelem, ion;
         bool lgOK=true;
-
-        if (dynamics.lgAdvection)
-        {
-                // Switch off invariant checks for dynamics -- better to
-                // work out broader invariant which applies here...
-                fixit(); 
-                return true;
-        }
 
         /* this checks that all molecules and ions add up to gas phase abundance

branches/newmole/source/mole_eval_balance.cpp

@@ -90 +90 @@
                                 if(i!=j)
                                 {
-                                        rate_deriv[i] *= rate->reactants[j]->hevmol;
-                                }
-                        }
-                }
-                
-                rate_tot = rate_deriv[0]*rate->reactants[0]->hevmol;
+                                        rate_deriv[i] *= fabs(rate->reactants[j]->hevmol);
+                                }
+                        }
+                }
+                
+                rate_tot = rate_deriv[0]*fabs(rate->reactants[0]->hevmol);
                 
                 for(i=0;i<rate->nreactants;i++)

branches/newmole/source/mole_newton_step.cpp

@@ -194 +194 @@
                 
                 /* add positive ions and neutral atoms: ratios are set by ion_solver,
-                 * we determine sum of two here */
+                 * we determine abundance of the group as a whole here */
                 
                 if (loop == 0) {

branches/newmole/source/pressure_total.cpp

@@ -70 +70 @@
                 pmx;
         realnum den_hmole ,
-                DenseAtomsIons,
-                DensElements;
+                DenseAtomsIons;
 
         DEBUG_ENTRY( "PresTotCurrent()" );
@@ -132 +131 @@
         TempChange(phycon.te , false);
 
-        /* find total number of atoms and ions */
-        DenseAtomsIons = 0.;
-        DensElements = 0.;
-        realnum factor = 1.1f;
-        if( conv.lgSearch )
-                factor = 2.f;
-        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
-        {
-                /* only check element solution if it is on, and ionization 
-                 * has not been set with TABLE ELEMENT command */
-                if( dense.lgElmtOn[nelem]  )
-                {
-                        /* gas_phase is density of all atoms and ions, but not molecules */
-                        DensElements += dense.gas_phase[nelem];
-                        realnum DenseOne = 0;
-                        for( ion=0; ion<=nelem+1; ++ion )
-                                DenseOne += dense.xIonDense[nelem][ion];
-
-                        /* during search phase sums may not add up correctly, and during
-                         * early parts of search phase may be badly off.  This test 
-                         * was introduced as result off fpe due to another
-                         * reason - Te had changed too much during initial search
-                         * for sim in which chemistry was important - fix was to not
-                         * let Te change.  During resulting insanity caused by large
-                         * change, linearization did not work, CO and ionization solvers
-                         * could diverge leading to molecule or ionization population
-                         * larger than 1e38 limit to a realnum, fpe followed.  this is to
-                         * catch that in its early stages */
-                        ASSERT( conv.nTotalIoniz<5 || dense.lgSetIoniz[nelem] ||
-                                DenseOne <= dense.gas_phase[nelem]*factor );
-                        DenseAtomsIons += DenseOne;
-                }
-        }
-        /* DensElements is sum of all gas-phase densities of all elements,
-         * DenseAtomsIons is sum of density of atoms and ions, does not 
-         * include molecules */
-        ASSERT( conv.nTotalIoniz<5 || DenseAtomsIons <= DensElements*factor );
-        ASSERT( DenseAtomsIons > 0. );
+        ASSERT(lgElemsConserved());
 
         /* evaluate the total ionization energy on a scale where neutral atoms
@@ -203 +165 @@
          * never set and only appears in print statements */
         phycon.EnergyBinding = 0.;
+
+        /* find total number of atoms and ions */
+        DenseAtomsIons = 0.;
+        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
+        {
+                /* only check element solution if it is on, and ionization 
+                 * has not been set with TABLE ELEMENT command */
+                if( dense.lgElmtOn[nelem]  )
+                {
+                        /* gas_phase is density of all atoms and ions, but not molecules */
+                        for( ion=0; ion<=nelem+1; ++ion )
+                                DenseAtomsIons += dense.xIonDense[nelem][ion];
+                }
+        }
+        /* DensElements is sum of all gas-phase densities of all elements,
+         * DenseAtomsIons is sum of density of atoms and ions, does not 
+         * include molecules */
+        ASSERT( DenseAtomsIons > 0. );
 
         /* find number of molecules of the heavy elements in the gas phase. 

branches/newmole/source/rt_line_one_tauinc.cpp

@@ -48 +48 @@
 
                         /* sanity check - total fine opacity should be greater than single ine opacity */
-                        ASSERT( fabs(dTau_line) < 0.01 || 
+                        ASSERT( fabs(dTau_line) < 0.01 || wind.windv != 0 ||
                                 rfield.fine_opac_zone[t->Emis->ipFine]*1.1 >= opacity*t->Emis->PopOpc  );
                 }