From 7538bba4844fbd4913b15bcbf37b0d845421c28c Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Wed, 9 Apr 2014 15:26:01 +0100 Subject: [PATCH 1/3] Added errors for -ve rates, moved macro_pops call out of elements loop --- matom.c | 27 ++++++++++++++++++++++++++- saha.c | 13 ++++++------- 2 files changed, 32 insertions(+), 8 deletions(-) diff --git a/matom.c b/matom.c index bbb82393f..63d4661cb 100755 --- a/matom.c +++ b/matom.c @@ -1546,6 +1546,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 +1576,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 +1611,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 +1628,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 +1669,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); + } } } } @@ -1762,7 +1787,7 @@ macro_pops (xplasma, xne) ion[index_ion].first_nlte_level + ion[index_ion].nlte; index_lvl++) { - this_ion_density += gsl_vector_get (populations, nn); + this_ion_density += gsl_vector_get (populations, conf_to_matrix[index_lvl]); nn++; } diff --git a/saha.c b/saha.c index ee11296ab..ebd223127 100755 --- a/saha.c +++ b/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++) { From 942f16e208cd8e36972f55e8bfdfcda94ba9d943 Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Mon, 14 Apr 2014 13:54:14 +0100 Subject: [PATCH 2/3] Added test to check that populations are solution to rate matrix equation in macro_pops --- matom.c | 87 +++++++++++++++++++++++++++++++++++++++++++++++------ my_linalg.h | 16 ++++++++++ 2 files changed, 93 insertions(+), 10 deletions(-) diff --git a/matom.c b/matom.c index 63d4661cb..f66257499 100755 --- a/matom.c +++ b/matom.c @@ -9,6 +9,7 @@ #include #include #include +//#include #include "my_linalg.h" /***************************************************************************** @@ -1394,7 +1395,7 @@ 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 @@ -1402,8 +1403,10 @@ macro_pops (xplasma, xne) 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 = ¯omain[xplasma->nplasma]; @@ -1692,8 +1695,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 */ @@ -1708,19 +1713,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); @@ -1729,10 +1745,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. */ @@ -1787,7 +1845,14 @@ macro_pops (xplasma, xne) ion[index_ion].first_nlte_level + ion[index_ion].nlte; index_lvl++) { - this_ion_density += gsl_vector_get (populations, conf_to_matrix[index_lvl]); + 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++; } @@ -3140,3 +3205,5 @@ q_recomb (cont_ptr, electron_temperature) return (coeff); } + + diff --git a/my_linalg.h b/my_linalg.h index 5d08ec133..e01d228b6 100755 --- a/my_linalg.h +++ b/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__ */ From d4a410c7c7aa5c0a599b7b11c868084f01d0e534 Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Mon, 14 Apr 2014 13:56:21 +0100 Subject: [PATCH 3/3] Added a comment describing change; --- matom.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matom.c b/matom.c index f66257499..ddef287b9 100755 --- a/matom.c +++ b/matom.c @@ -1375,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