Changeset 685 for trunk/source/ion_solver.cpp

Timestamp:

12/05/06 18:14:43 (2 years ago)

Author:

gary

Message:

dayslow_master.pl - use icl rather than gcc for slow test suite
ion_solver.cpp - undo bug introduced in rev 622. now looks at matrix elements to see if singular
rt_stark.cpp - resolved ambiguous overload
blr_n09_p18_Z20.in - add punch convergence error command
limit_hi_ion.in - renamed func_hi_ion.in to this since it tests code in high ionization limit

Files:

• trunk/source/ion_solver.cpp (modified) (6 diffs)

trunk/source/ion_solver.cpp

@@ -45 +45 @@
         int lgNegPop;
         static int lgMustMalloc=TRUE;
-        static double *sink, *source;
+        static double *sink, *source , *sourcesave;
         int32 nerror;
         static int32 *ipiv;
@@ -221 +221 @@
                         BadMalloc();
                 if( (source=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ))) == NULL )
+                        BadMalloc();
+                if( (sourcesave=(double*)MALLOC( (sizeof(double)*(unsigned)ncell ))) == NULL )
                         BadMalloc();
                 if( (ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)ncell ))) == NULL )
@@ -514 +516 @@
                         MAT( xmatsave, i, j ) = MAT( xmat, i, j );
                 }
+                sourcesave[i] = source[i];
         }
 
@@ -523 +526 @@
         * will not be removed, a possible source of error, but they must not
         * have been significant, given that the molecular fraction is so small */
+        if( !lgHomogeneous && ion_range==2 )
+        {
+                /* solve algebraically */
+                double a = MAT( xmatsave, 0, 0 ), 
+                        b = MAT( xmatsave, 1, 0 ) ,
+                        c = MAT( xmatsave, 0, 1 ) ,
+                        d = MAT( xmatsave, 1, 1 );
+                b = SDIV(b);
+                d = SDIV(d);
+                double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
+                if( fabs(ratio1-ratio2)/MAX2(fratio1,fratio2) <DBL_EPSILON*1e4 )
+                {
+                        abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
+                        lgHomogeneous = TRUE;
+                }
+        }
+#       if 0
         if( nelem==ipHYDROGEN && 
                 fabs(dense.xMolecules[nelem]) / SDIV(dense.gas_phase[ipHYDROGEN]) <DBL_EPSILON*100. )
         {
+                abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
                 lgHomogeneous = TRUE;
         }
+#       endif
 
         /* this is true if no source terms 
@@ -566 +588 @@
         }
 
+        for( i=0; i< ion_range; ++i )
+        {
+                for( j=0; j< ion_range; ++j )
+                {
+                        MAT( xmatsave, i, j ) = MAT( xmat, i, j );
+                }
+                sourcesave[i] = source[i];
+        }
+
         if ( FALSE && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 ) 
         {
@@ -728 +759 @@
                 enum {DEBUG_LOC=FALSE};
                 /*@+redef@*/
-                if( DEBUG_LOC && (nzone >380 ) && nelem == ipHYDROGEN )
-                {
-                        fprintf(ioQQQ,"debuggg\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n", 
+                if( DEBUG_LOC && nelem == ipHYDROGEN )
+                {
+                        fprintf(ioQQQ,"debuggg ion_solver1 \t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n", 
                                 fnzone,
                                 phycon.te,