Added errors for -ve rates, moved macro_pops call out of elements loop

matom.c

@@ -9,6 +9,7 @@
 #include <gsl/gsl_block.h>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_matrix.h>
+//#include <gsl/gsl_blas.h>
 #include "my_linalg.h"
 
 /*****************************************************************************
@@ -1374,7 +1375,8 @@ Returns:   Should compute the fractional level populations for
                        - some treatment of metastables is needed for 
                        doing CNO macro atoms.
 	06may	ksl	57+ -- Modified to reflect plasma structue
-
+    1404 	jm  77a -- Added test to check that populations are solution to rate matrix equation in macro_pops
+                       + more robust erro reporting in general. See pull request #75.
 ************************************************************/
 
 int
@@ -1394,16 +1396,18 @@ macro_pops (xplasma, xne)
   int nn, mm;
   int index_bbu, index_bbd, index_bfu, index_bfd;
   int lower, upper;
-  double this_ion_density;
+  double this_ion_density, level_population;
   double inversion_test;
   double q_ioniz (), q_recomb ();
   int s;			//NEWKSL
   double *a_data, *b_data;
   gsl_permutation *p;		//NEWKSL
   gsl_matrix_view m;		//NEWKSL
   gsl_vector_view b;		//NEWKSL
-  gsl_vector *populations;	//NEWKSL
-  int index_fast_col;
+  gsl_vector *populations, *test_vector;	//NEWKSL
+  gsl_matrix *test_matrix;
+  int index_fast_col, ierr;
+  double test_val;
 
   MacroPtr mplasma;
   mplasma = &macromain[xplasma->nplasma];
@@ -1546,6 +1550,11 @@ macro_pops (xplasma, xne)
 
 		      rate_matrix[lower][lower] += -1. * rate;
 		      rate_matrix[upper][lower] += rate;
+
+		      if (rate < 0.0 || sane_check(rate))
+		      {
+		      	Error("macro_pops: bbu rate is %8.4e in cell/matom %i\n", rate, xplasma->nplasma);
+		      }
 		    }
 
 		  for (index_bbd = 0;
@@ -1571,6 +1580,11 @@ macro_pops (xplasma, xne)
 
 		      rate_matrix[upper][upper] += -1. * rate;
 		      rate_matrix[lower][upper] += rate;
+
+		      if (rate < 0.0 || sane_check(rate))
+		      {
+		      	Error("macro_pops: bbd rate is %8.4e in cell/matom %i\n", rate, xplasma->nplasma);
+		      }
 		    }
 
 
@@ -1601,6 +1615,11 @@ macro_pops (xplasma, xne)
 		      rate_matrix[lower][lower] += -1. * rate;
 		      rate_matrix[upper][lower] += rate;
 
+		      if (rate < 0.0 || sane_check(rate))
+		      {
+		      	Error("macro_pops: bfu rate is %8.4e in cell/matom %i\n", rate, xplasma->nplasma);
+		      }
+
 		      /* Now deal with the stimulated emission. */
 		      /* Lower and upper are the same, but now it contributes in the
 		         other direction. */
@@ -1613,6 +1632,11 @@ macro_pops (xplasma, xne)
 		      rate_matrix[upper][upper] += -1. * rate;
 		      rate_matrix[lower][upper] += rate;
 
+		      if (rate < 0.0 || sane_check(rate))
+		      {
+		      	Error("macro_pops: st. recomb rate is %8.4e in cell/matom %i\n", rate, xplasma->nplasma);
+		      }
+
 		    }
 
 
@@ -1649,6 +1673,11 @@ macro_pops (xplasma, xne)
 
 		      rate_matrix[upper][upper] += -1. * rate;
 		      rate_matrix[lower][upper] += rate;
+
+		      if (rate < 0.0 || sane_check(rate))
+		      {
+		      	Error("macro_pops: bfd rate is %8.4e in cell/matom %i\n", rate, xplasma->nplasma);
+		      }
 		    }
 		}
 	    }
@@ -1667,8 +1696,10 @@ macro_pops (xplasma, xne)
 
 	  /********************************************************************************/
 	  /* The block that follows (down to next line of ***s) is to do the
-	     matix inversion. It uses LU decomposition - the
-	     code for doing this is taken from the GSL manual with very few modifications. */
+	     matrix inversion. It uses LU decomposition - the code for doing this is 
+	     taken from the GSL manual with very few modifications. */
+	  /* here we solve the matrix equation M x = b, where x is our vector containing
+	     level populations as a fraction w.r.t the whole element */
 
 	  /* Replaced inline array allocaation with calloc, which will work with older version of c compilers */
 
@@ -1683,19 +1714,30 @@ macro_pops (xplasma, xne)
 		}
 	    }
 
-	  /* Replaced inline array allocaation with calloc, which will work with older version of c compilers */
+
+	  /* Replaced inline array allocaation with calloc, which will work with older version of c compilers 
+	     calloc also sets the elements to zero, which is required */
 
 	  b_data = (double *) calloc (sizeof (rate), n_macro_lvl);
+
+	  /* replace the first entry with 1.0- this is part of the normalisation constraint */
 	  b_data[0] = 1.0;
 
+      /* create gsl matrix/vector views of the arrays of rates */
 	  m = gsl_matrix_view_array (a_data, n_macro_lvl, n_macro_lvl);	//KSLNEW
 
+	  /* these are used for testing the solution below */
+	  test_matrix = gsl_matrix_alloc(n_macro_lvl, n_macro_lvl);
+	  test_vector = gsl_vector_alloc (n_macro_lvl);
+
+	  gsl_matrix_memcpy(test_matrix, &m.matrix); 	// create copy for testing 
+
 	  b = gsl_vector_view_array (b_data, n_macro_lvl);	//KSLNEW
 
+	  /* the populations vector will be a gsl vector which stores populations */
 	  populations = gsl_vector_alloc (n_macro_lvl);	//KSLNEW
 
 
-
 	  p = gsl_permutation_alloc (n_macro_lvl);	//NEWKSL
 
 	  gsl_linalg_LU_decomp (&m.matrix, p, &s);
@@ -1704,10 +1746,52 @@ macro_pops (xplasma, xne)
 
 	  gsl_permutation_free (p);
 
-	  free (a_data);	//NEWKSL
-	  free (b_data);	//NEWKSL
 
 	  /**********************************************************/
+
+      /* JM 140414 -- before we clean, we should check that the populations vector we 
+         have just created really is a solution to the matrix equation */
+
+      /* declaration contained in my_linalg.h, taken from gsl library */
+      ierr = gsl_blas_dgemv(CblasNoTrans, 1.0, test_matrix, populations, 0.0, test_vector);
+
+      if (ierr != 0)
+        {
+      	  Error("macro_pops: bad return when testing matrix solution to rate equations.\n");
+        }
+
+      /* now cycle through and check the solution to y = m * populations really is
+         (1, 0, 0 ... 0) */
+
+      for (mm = 0; mm < n_macro_lvl; mm++)
+        {
+
+          /* get the element of the vector we want to check */
+          test_val = gsl_vector_get (test_vector, mm);
+
+          if (mm == 0)		// 0th element should be 1
+          {
+      	    if ( (fabs(test_val - 1.0)) > EPSILON)
+      	      Error("macro_pops: test vector has first element = %8.4e, should be 1.0\n",
+      	      	     test_val);    
+          }     
+
+          else 				// all other elements should be zero
+          {
+          	if ( fabs(test_val)  > EPSILON)
+      	      Error("macro_pops: test vector has element %i = %8.4e, should be 0.0\n",
+      	      	     mm, test_val);    
+          }
+        }
+
+        /* free memory */
+        free (a_data);	
+	    free (b_data);	
+	    free (test_vector);
+	    free (test_matrix);
+
+
+
 	  /* MC noise can cause population inversions (particularly amongst highly excited states)
 	     which are never a good thing and most likely unphysical.
 	     Therefor let's follow Leon's procedure (Lucy 2003) and remove inversions. */
@@ -1762,7 +1846,14 @@ macro_pops (xplasma, xne)
 		   ion[index_ion].first_nlte_level + ion[index_ion].nlte;
 		   index_lvl++)
 		{
-		  this_ion_density += gsl_vector_get (populations, nn);
+		  level_population = gsl_vector_get (populations, conf_to_matrix[index_lvl]);	
+		  this_ion_density += level_population;
+
+		  /* JM 140409 -- add error check for negative populations */
+		  if (level_population < 0.0 || sane_check(level_population))
+		    Error("macro_pops: level %i has frac. pop. %8.4e in cell %i\n", 
+		    	   index_lvl, level_population, xplasma->nplasma);
+
 		  nn++;
 		}
 
@@ -3115,3 +3206,5 @@ q_recomb (cont_ptr, electron_temperature)
 
   return (coeff);
 }
+
+

my_linalg.h

@@ -400,5 +400,21 @@ int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
 int gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D);
 
 
+/* JM 140414 -- I've had to include a few definitions from gsl/gsl_cblas.h and gsl/gsl_blas_types.h here
+   instead of including the files, as some of the variables conflicted with python */
+
+typedef  enum CBLAS_TRANSPOSE   CBLAS_TRANSPOSE_t;
+enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
+
+/* we use this function to multiply matrices together and check if our populations matrix
+   is a solution to the rate equations */
+int  gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA,
+                     double alpha,
+                     const gsl_matrix * A,
+                     const gsl_vector * X,
+                     double beta,
+                     gsl_vector * Y);
+
+
 __END_DECLS
 #endif /* __GSL_LINALG_H__ */

saha.c

@@ -653,19 +653,18 @@ lucy (xplasma)
 
 	  lucy_mazzali1 (nh, t_r, t_e, www, nelem, xplasma->ne,
 			 xplasma->density, xne, newden);
+  }
 
-	  /* Re solve for the macro atom populations with the current guess for ne */
-          /* JM1308 -- note that unlike lucy mazzali above, here we actually modify the xplasma
-	     structure for those ions which are being treated as macro ions. This means that the
-	     the newden array will contain wrong values for these particular macro ions, but due
-             to the if loop at the end of this subroutine they are never passed to xplasma */
-
+    /* Re solve for the macro atom populations with the current guess for ne */
+    /* JM1308 -- note that unlike lucy mazzali above, here we actually modify the xplasma
+       structure for those ions which are being treated as macro ions. This means that the
+       the newden array will contain wrong values for these particular macro ions, but due
+       to the if loop at the end of this subroutine they are never passed to xplasma */
 	  if (geo.macro_ioniz_mode == 1)
 	    {
 	      macro_pops (xplasma, xne);
 	    }
 
-	}
 
       for (nion = 0; nion < nions; nion++)
 	{