cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sanity_check.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*SanityCheck, check that various parts of the code still work, called by Cloudy after continuum
4  * and optical depth arrays are set up, but before initial temperature and ionization */
5 #include "cddefines.h"
6 #include "dense.h"
7 #include "elementnames.h"
8 #include "continuum.h"
9 #include "rfield.h"
10 #include "iso.h"
11 #include "opacity.h"
12 #include "hydro_bauman.h"
13 #include "heavy.h"
14 #include "trace.h"
15 #include "freebound.h"
16 #include "integrate.h"
17 #include "atmdat_gaunt.h"
18 #include "cloudy.h"
19 
20 #ifndef NDEBUG
21 namespace {
22  class my_Integrand: public std::unary_function<double, double>
23  {
24  public:
25  double operator() (double x) const
26  {
27  return sin(x);
28  }
29  // the constructor is needed to avoid warnings by the Sun Studio compiler
30  my_Integrand() {}
31  };
32 
33  // create a version of sin with C++ linkage
34  // this is done because we need a C++ pointer to this routine below
35  double mySin(double x)
36  {
37  return sin(x);
38  }
39 }
40 #endif
41 
42 /* NB - this routine must not change any global variables - any that are changed as part of
43  * a test must be reset, so that the code retains state */
44 
47 STATIC void SanityCheckGaunt(long Z, double loggam2, double logu, double refval, double relerr);
48 
49 /* chJob is either "begin" or "final"
50  * "begin is before code starts up
51  * "final is after model is complete */
52 void SanityCheck( const char *chJob )
53 {
54  DEBUG_ENTRY( "SanityCheck()" );
55 
56  if( strcmp(chJob,"begin") == 0 )
57  {
59  }
60  else if( strcmp(chJob,"final") == 0 )
61  {
63  }
64  else
65  {
66  fprintf(ioQQQ,"SanityCheck called with insane argument.\n");
68  }
69 }
70 
72 {
73  /* PrtComment also has some ending checks on sanity */
74 }
75 
77 {
78  bool lgOK=true;
79  int lgFlag;// error return for spsort, 0 success, >=1 for errors
80  int32 ner, ipiv[3];
81  long i ,
82  j ,
83  nelem ,
84  ion ,
85  nshells;
86  double *A;
87 
88  /* this will be charge to the 4th power */
89  double Aul ,
90  error,
91  Z4;
92 
93  long n;
94 
95  const int NDIM = 10;
96  double xMatrix[NDIM][NDIM] , yVector[NDIM] ,
97  rcond;
98  realnum *fvector;
99  long int *ipvector;
100 
101  DEBUG_ENTRY( "SanityCheck()" );
102 
103  /*********************************************************
104  * *
105  * confirm that various part of cloudy still work *
106  * *
107  *********************************************************/
108 
109  /* if this is no longer true at end, we have a problem */
110  lgOK = true;
111 
112  /*********************************************************
113  * *
114  * check that all the Lyas As are ok *
115  * *
116  *********************************************************/
117  for( nelem=0; nelem<LIMELM; ++nelem )
118  {
119  /* this element may be turned off */
120  if( dense.lgElmtOn[nelem] )
121  {
122  /* H_Einstein_A( n, l, np, lp, iz ) - all are on physics scale */
123  Aul = H_Einstein_A( 2, 1, 1, 0, nelem+1 );
124  /*fprintf(ioQQQ,"%li\t%.4e\n", nelem+1, iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul() );*/
125  if( fabs(Aul - iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul() ) /Aul > 0.01 )
126  {
127  fprintf(ioQQQ," SanityCheck found insane H-like As.\n");
128  lgOK = false;
129  }
130  }
131  }
132 
133  /*********************************************************
134  * *
135  * check that gaunt factors are good *
136  * *
137  *********************************************************/
138  // our Gaunt factors are merged relativistic and non-relativistic data
139  // for high gamma^2 we have non-relativistic data, these are compared to Sutherland (1998)
140  // for low gamma^2 we have relativistic data, these are compared to Nozawa+ (1998)
141  // more elaborate checks are done in the unit test suite
142 
143  // these data were taken from
144  // >>refer HI gff Sutherland R.S., 1998, MNRAS, 300, 321
145  SanityCheckGaunt( 1, 4., -4., 2.7008, 0.0004 );
146  SanityCheckGaunt( 1, 4., 4., 1.1040, 0.0004 );
147  SanityCheckGaunt( 1, 4., 0., 1.0237, 0.0004 );
148  SanityCheckGaunt( 2, 3., -3., 2.2126, 0.0004 );
149  SanityCheckGaunt( 2, 1., -1., 1.5088, 0.0004 );
150  // these data were taken from
151  // >>refer gff Nozawa S., Itoh N., Kohyama Y., 1998, ApJ, 507, 530
152  SanityCheckGaunt( 1, -4., -4., 6.120, 0.005 );
153  SanityCheckGaunt( 2, -4., -4., 8.934, 0.005 );
154  SanityCheckGaunt( 1, -3., -1., 1.852, 0.005 );
155  SanityCheckGaunt( 2, -3., -1., 1.921, 0.005 );
156  SanityCheckGaunt( 1, 0., 3., 0.2001, 0.005 );
157  SanityCheckGaunt( 2, 0., 3., 0.2041, 0.005 );
158 
159  /*********************************************************
160  * *
161  * check some transition probabililties for he-like ions *
162  * *
163  *********************************************************/
164  for( nelem=1; nelem<LIMELM; ++nelem )
165  {
166  /* the helike 9-1 transition, A(3^3P to 2^3S) */
167  double as[]={
168  /* updated with Johnson values */
169  0. ,9.47e+006 ,3.44e+008 ,1.74e+009 ,5.51e+009 ,1.34e+010 ,
170  2.79e+010 ,5.32E+010 ,8.81e+010 ,1.46E+011 ,2.15e+011 ,3.15e+011 ,
171  4.46e+011 ,6.39E+011 ,8.26e+011 ,1.09e+012 ,1.41e+012 ,1.86E+012 ,
172  2.26e+012 ,2.80e+012 ,3.44e+012 ,4.18e+012 ,5.04e+012 ,6.02e+012 ,
173  7.14e+012 ,8.40e+012 ,9.83e+012 ,1.14e+013 ,1.32e+013 ,1.52e+013
174  };
175 
176  if( iso_sp[ipHE_LIKE][nelem].numLevels_max > 8 && dense.lgElmtOn[nelem])
177  {
178  /* following used to print current values of As */
179  if( fabs( as[nelem] - iso_sp[ipHE_LIKE][nelem].trans(9,1).Emis().Aul() ) /as[nelem] > 0.025 )
180  {
181  fprintf(ioQQQ,
182  " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
183  nelem,
184  as[nelem] ,
185  iso_sp[ipHE_LIKE][nelem].trans(9,1).Emis().Aul() );
186  lgOK = false;
187  }
188  }
189  }
190 
191  /* only do this test if case b is not in effect */
192  if( !opac.lgCaseB )
193  {
194 
195  for( i = 0; i <=110; i++ )
196  {
197  double DrakeTotalAuls[111] = {
198  -1.0000E+00, -1.0000E+00, -1.0000E+00, 1.02160E+07,
199  1.02160E+07, 1.02160E+07, 1.80090E+09, 2.78530E+07,
200  1.82990E+07, 1.05480E+07, 7.07210E+07, 6.37210E+07,
201  5.79960E+08, 1.60330E+07, 1.13640E+07, 7.21900E+06,
202  3.11920E+07, 2.69830E+07, 1.38380E+07, 1.38330E+07,
203  2.52270E+08, 9.20720E+06, 6.82220E+06, 4.56010E+06,
204  1.64120E+07, 1.39290E+07, 7.16030E+06, 7.15560E+06,
205  4.25840E+06, 4.25830E+06, 1.31150E+08, 5.62960E+06,
206  4.29430E+06, 2.95570E+06, 9.66980E+06, 8.12340E+06,
207  4.19010E+06, 4.18650E+06, 2.48120E+06, 2.48120E+06,
208  1.64590E+06, 1.64590E+06, 7.65750E+07, 3.65330E+06,
209  2.84420E+06, 1.99470E+06, 6.16640E+06, 5.14950E+06,
210  2.66460E+06, 2.66200E+06, 1.57560E+06, 1.57560E+06,
211  1.04170E+06, 1.04170E+06, 7.41210E+05, 7.41210E+05,
212  4.84990E+07, 2.49130E+06, 1.96890E+06, 1.39900E+06,
213  4.16900E+06, 3.46850E+06, 1.79980E+06, 1.79790E+06,
214  1.06410E+06, 1.06410E+06, 7.02480E+05, 7.02480E+05,
215  4.98460E+05, 4.98460E+05, -1.0000E+00, -1.0000E+00,
216  3.26190E+07, 1.76920E+06, 1.41440E+06, 1.01460E+06,
217  2.94830E+06, 2.44680E+06, 1.27280E+06, 1.27140E+06,
218  7.52800E+05, 7.52790E+05, 4.96740E+05, 4.96740E+05,
219  3.51970E+05, 3.51970E+05, -1.0000E+00, -1.0000E+00,
220  -1.0000E+00, -1.0000E+00, 2.29740E+07, 1.29900E+06,
221  1.04800E+06, 7.57160E+05, 2.16090E+06, 1.79030E+06,
222  9.33210E+05, 9.32120E+05, 5.52310E+05, 5.52310E+05,
223  3.64460E+05, 3.64460E+05, 2.58070E+05, 2.58070E+05,
224  -1.0000E+00, -1.0000E+00, -1.0000E+00, -1.0000E+00,
225  -1.0000E+00, -1.0000E+00, 1.67840E+07};
226 
227  if( DrakeTotalAuls[i] > 0. &&
228  i < iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max)
229  {
230  if( fabs( DrakeTotalAuls[i] - (1./iso_sp[ipHE_LIKE][ipHELIUM].st[i].lifetime()) ) /DrakeTotalAuls[i] > 0.001 )
231  {
232  fprintf(ioQQQ,
233  " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
234  i,
235  DrakeTotalAuls[i],
236  (1./iso_sp[ipHE_LIKE][ipHELIUM].st[i].lifetime()) );
237  lgOK = false;
238  }
239  }
240  }
241  }
242 
243  /*********************************************************
244  * *
245  * check the threshold photoionization cs for He I *
246  * *
247  *********************************************************/
248  if( dense.lgElmtOn[ipHELIUM] )
249  {
250  /* HeI photoionization cross sections at threshold for lowest 20 levels */
251  const int NHE1CS = 20;
252  double he1cs[NHE1CS] =
253  {
254  5.480E-18 , 9.253E-18 , 1.598E-17 , 1.598E-17 , 1.598E-17 , 1.348E-17 ,
255  8.025E-18 , 1.449E-17 , 2.852E-17 , 1.848E-17 , 1.813E-17 , 2.699E-17 ,
256  1.077E-17 , 2.038E-17 , 4.159E-17 , 3.670E-17 , 3.575E-17 , 1.900E-17 ,
257  1.900E-17 , 4.175E-17
258  };
259 
260  /* loop over levels and check on photo cross section */
261  j = MIN2( NHE1CS+1 , iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max -iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max );
262  for( n=1; n<j; ++n )
263  {
264  /* above list of levels does not include the ground */
265  i = iso_sp[ipHE_LIKE][ipHELIUM].fb[n].ipOpac;
266  ASSERT( i>0 );
267  /*fprintf(ioQQQ,"%li\t%lin", n , i );*/
268  /* >>chng 02 apr 10, from 0.01 to 0.02, values stored
269  * where taken from calc at low contin resolution, when continuum
270  * resolution changed this changes too */
271  /*fprintf(ioQQQ,"%li %.2e\n", n,( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] );*/
272  /* >>chng 02 jul 16, limt from 0.02 to 0.04, so that "set resolution 4" will work */
273  /* >>chng 04 may 18, levels 10 and 11 are about 12% off - because of energy binning, chng from 0.08 to 0.15 */
274  if( fabs( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
275  {
276  fprintf(ioQQQ,
277  " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
278  n,
279  he1cs[n-1] ,
280  opac.OpacStack[i - 1]);
281  fprintf(ioQQQ,
282  " n=%li, l=%li, s=%li\n",
283  iso_sp[ipHE_LIKE][ipHELIUM].st[n].n() ,
284  iso_sp[ipHE_LIKE][ipHELIUM].st[n].l() ,
285  iso_sp[ipHE_LIKE][ipHELIUM].st[n].S());
286  lgOK = false;
287  }
288  }
289  }
290 
291  for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
292  {
293  long nelem = ipISO;
294  /* Check for agreement between on-the-fly and interpolation calculations
295  * of recombination, but only if interpolation is turned on. */
296  if( !iso_ctrl.lgNoRecombInterp[ipISO] )
297  {
298  /* check the recombination coefficients for ground state */
299  error = fabs( iso_recomb_check( ipISO, nelem , 0 , 7500. ) );
300  if( error > 0.01 )
301  {
302  fprintf(ioQQQ,
303  " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n",
304  iso_ctrl.chISO[ipISO],
305  elementnames.chElementSym[nelem],
306  0,
307  error );
308  lgOK = false;
309  }
310 
311  /* check the recombination coefficients for ground state of the root of each iso sequence */
312  error = fabs( iso_recomb_check( ipISO, nelem , 1 , 12500. ) );
313  if( error > 0.01 )
314  {
315  fprintf(ioQQQ,
316  " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n",
317  iso_ctrl.chISO[ipISO],
318  elementnames.chElementSym[nelem],
319  1,
320  error );
321  lgOK = false;
322  }
323  }
324  }
325 
326  /*********************************************************
327  * *
328  * check out the sorting routine *
329  * *
330  *********************************************************/
331 
332  const int NSORT = 100 ;
333 
334  fvector = (realnum *)MALLOC((NSORT)*sizeof(realnum) );
335 
336  ipvector = (long *)MALLOC((NSORT)*sizeof(long int) );
337 
338  nelem = 1;
339  /* make up some unsorted values */
340  for( i=0; i<NSORT; ++i )
341  {
342  nelem *= -1;
343  fvector[i] = (realnum)nelem * ((realnum)NSORT-i);
344  }
345 
346  /*spsort netlib routine to sort array returning sorted indices */
347  spsort(fvector,
348  NSORT,
349  ipvector,
350  /* flag saying what to do - 1 sorts into increasing order, not changing
351  * the original routine */
352  1,
353  &lgFlag);
354 
355  if( lgFlag ) lgOK = false;
356 
357  for( i=1; i<NSORT; ++i )
358  {
359  /*fprintf(ioQQQ," %li %li %.0f\n",
360  i, ipvector[i],fvector[ipvector[i]] );*/
361  if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
362  {
363  fprintf(ioQQQ," SanityCheck found insane sort\n");
364  lgOK = false;
365  }
366  }
367 
368  free( fvector );
369  free( ipvector);
370 
371 # if 0
372  ttemp = (realnum)sqrt(phycon.te);
373  /* check that the temperatures make sense */
374  if( fabs(ttemp - phycon.sqrte )/ttemp > 1e-5 )
375  {
376  fprintf(ioQQQ , "SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
377  phycon.te ,
378  sqrt(phycon.te) ,
379  phycon.sqrte ,
380  fabs(sqrt(phycon.te) - phycon.sqrte ) );
381  lgOK = false;
382  }
383 # endif
384 
385  /*********************************************************
386  * *
387  * confirm the validity of the mesh *
388  * *
389  *********************************************************/
390 
392  rfield.CheckMesh();
393 
394  /*********************************************************
395  * *
396  * confirm that hydrogen einstein As are still valid *
397  * *
398  *********************************************************/
399  for( nelem=0; nelem<2; ++nelem )
400  {
401  /* this element may be turned off */
402  if( dense.lgElmtOn[nelem] )
403  {
404  /*Z4 = (double)(POW2(nelem+1)*POW2(nelem+1));*/
405  /* form charge to the 4th power */
406  Z4 = (double)(nelem+1);
407  Z4 *= Z4;
408  Z4 *= Z4;
409  /* H Lya */
410  double ans1 = iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul();
411  double ans2 = 6.265e8*Z4;
412  if( fabs(ans1-ans2)/ans2 > 1e-3 )
413  {
414  fprintf(ioQQQ , "SanityCheck finds insane A for H Lya %g %g nelem=%li\n",
415  ans1 , ans2 , nelem );
416  lgOK = false;
417  }
418  }
419  }
420 
421  /* check that hydrogenic branching ratios add up to unity */
422  for( nelem=0; nelem<LIMELM; ++nelem )
423  {
424  if( dense.lgElmtOn[nelem] )
425  {
426  int ipHi, ipLo;
427  for( ipHi=4; ipHi< iso_sp[ipH_LIKE][nelem].numLevels_max-iso_sp[ipH_LIKE][nelem].nCollapsed_max; ++ipHi )
428  {
429  double sum = 0.;
430  for( ipLo=0; ipLo<ipHi; ++ipLo )
431  {
432  sum += iso_sp[ipH_LIKE][nelem].BranchRatio[ipHi][ipLo];
433  }
434  if( fabs(sum-1.)>0.01 )
435  {
436  fprintf(ioQQQ ,
437  "SanityCheck H branching ratio sum not unity for nelem=%li upper level=%i sum=%.3e\n",
438  nelem, ipHi, sum );
439  lgOK = false;
440  }
441  }
442  }
443  }
444 
445  /* check that hydrogenic lifetimes are sane (compare inverse sum of As with closed form of lifetime) */
446  for( nelem=0; nelem<LIMELM; ++nelem )
447  {
448  if( dense.lgElmtOn[nelem] )
449  {
450  int ipHi, ipLo;
451  for( ipHi=1; ipHi< iso_sp[ipH_LIKE][nelem].numLevels_max-iso_sp[ipH_LIKE][nelem].nCollapsed_max; ++ipHi )
452  {
453  double inverse_sum = 0.;
454  double sum = 0.;
455  long ipISO = ipH_LIKE;
456 
457  /* we do not have an accurate closed form for l=0 lifetimes. Everything else should be very accurate. */
458  if( L_(ipHi)==0 )
459  continue;
460 
461  for( ipLo=0; ipLo<ipHi; ++ipLo )
462  {
463  sum += iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul();
464  }
465 
466  inverse_sum = 1./sum;
467 
468  double lifetime = iso_state_lifetime( ipH_LIKE, nelem, N_(ipHi), L_(ipHi) );
469  double percent_error = (1. - inverse_sum/lifetime)*100;
470  /* The closed form seems to consistently yield lifetimes roughly 0.2% smaller than the inverse sum
471  * Transition probabilities go like energy cubed. The cube of the finite mass rydberg is about 0.16% less than unity.
472  * is that the difference? */
473  /* for now enforce to better than 0.5% */
474  if( fabs(percent_error) > 0.5 )
475  {
476  fprintf(ioQQQ ,
477  "SanityCheck hydrogenic lifetime not consistent with closed form for nelem=%2li upper level=%2i inverse-sum= %.3e lifetime= %.3e error= %.2f%%\n",
478  nelem, ipHi, inverse_sum, lifetime, percent_error );
479  lgOK = false;
480  }
481  }
482  }
483  }
484 
485 
486  /* check photo cross sections for H */
487  long ipISO = ipH_LIKE;
488  nelem = 0;
489  for( long n=1; n <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++n )
490  {
491  // This sanity check is of little value, as the cross-section is calculated
492  // with the same routine here as when populating the opacity stack.
493  // However, the cell mid-point energy used in populating OpacStack can be
494  // greater than the energy implied by the relative energy below.
495  // For very Rydberg, very yrare levels, the cross-section can fall off fast
496  // enough over that energy difference to fail this check.
497  // Will reported a sim that failed the check for some levels with n >= 71.
498  // So just don't bother checking above 60.
499  if( n > 60 )
500  break;
501 
502  double rel_photon_energy = 1. + FLT_EPSILON*2.;
503  for( long l=0; l < n; l++ )
504  {
505  double cs;
506  long index = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][2];
507 
508  cs = H_photo_cs( rel_photon_energy, n, l, nelem+1 );
509 
510  error = fabs(cs - opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1] )/
511  ( (cs + opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1] )/2.);
512 
513  if( error > 0.05 )
514  {
515  fprintf(ioQQQ , "SanityCheck finds insane H photo cs, n,l = %li, %li, expected, found = %e, %e, error = %e\n",
516  n, l, opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1], cs, error );
517  lgOK = false;
518  }
519  }
520  }
521 
522  /*********************************************************
523  * *
524  * confirm that Gaussian integration routines still work *
525  * *
526  *********************************************************/
527  ASSERT( fp_equal( qg32( 0., PI, mySin ) , 2. ) );
528 
529 #ifndef NDEBUG
530  /* And again with the new structure */
531  my_Integrand func;
533 #endif
534  ASSERT( fp_equal( mySine.sum( 0., PI ), 2. ) );
535 
536  /*********************************************************
537  * *
538  * confirm that exponential integral routines still work *
539  * *
540  *********************************************************/
541 
542  /* check that first and second exponential integrals are ok,
543  * step through range of values, beginning with following */
544  double x = 1e-3;
545  do
546  {
547  {
548  /* check that fast e1 routine is ok */
549  double ans1 = e1(x);
550  double ans2 = expn( 1 , x );
551  error = safe_div(ans1-ans2,ans1+ans2,0.);
552  if( fabs(error) > 3e-7 )
553  {
554  fprintf(ioQQQ , "SanityCheck finds insane E1 %g %g %g error %g\n",
555  x , ans1 , ans2, error );
556  lgOK = false;
557  }
558  }
559  {
560  /* check that e2 is ok */
561  double ans1 = e2(x);
562  double ans2 = expn( 2 , x );
563  error = safe_div(ans1-ans2,ans1+ans2,0.);
564  if( fabs(error) > 1e-8 )
565  {
566  fprintf(ioQQQ , "SanityCheck finds insane E2 %g %g %g error %g\n",
567  x , ans1 , ans2, error );
568  lgOK = false;
569  }
570  }
571  /* now increment x */
572  x *= 2.;
573  /* following limit set by values not underflowing to zero */
574  } while( x < 700. );
575 
576  /*********************************************************
577  * *
578  * confirm that matrix inversion routine still works *
579  * *
580  *********************************************************/
581 
582  /* these are the answer, chosen to get xvec 1,2,3 */
583  yVector[0] = 1.;
584  yVector[1] = 3.;
585  yVector[2] = 3.;
586 
587  /* zero out the main matrix */
588  for(i=0;i<3;++i)
589  {
590  for( j=0;j<3;++j )
591  {
592  xMatrix[i][j] = 0.;
593  }
594  }
595 
596  /* remember that order is column, row, alphabetical order, rc */
597  xMatrix[0][0] = 1.;
598  xMatrix[0][1] = 1.;
599  xMatrix[1][1] = 1.;
600  xMatrix[2][2] = 1.;
601 
602  /* this is the default matrix solver */
603  /* this test is the 1-d matrix with 2-d macro simulation */
604  /* LDA is right dimension of matrix */
605 
606  /* MALLOC space for the 1-d array */
607  A = (double*)MALLOC( sizeof(double)*NDIM*NDIM );
608 
609  /* copy over the main matrix */
610  for( i=0; i < 3; ++i )
611  {
612  for( j=0; j < 3; ++j )
613  {
614  A[i*NDIM+j] = xMatrix[i][j];
615  }
616  }
617 
618  ner = 0;
619 
620  /*void DGETRF(long,long,double*,long,long[],long*);*/
621  /*void DGETRF(int,int,double*,int,int[],int*);*/
622  getrf_wrapper(3, 3, A, NDIM, ipiv, &ner);
623  if( ner != 0 )
624  {
625  fprintf( ioQQQ, " SanityCheck DGETRF error\n" );
627  }
628 
629  /* usage DGETRS, 'N' = no transpose
630  * order of matrix,
631  * number of cols in bvec, =1
632  * array
633  * leading dim of array */
634  /*void DGETRS(char,int,int,double*,int,int[],double*,int,int*);*/
635  getrs_wrapper('N', 3, 1, A, NDIM, ipiv, yVector, 3, &ner);
636 
637  if( ner != 0 )
638  {
639  fprintf( ioQQQ, " SanityCheck DGETRS error\n" );
641  }
642  /* release the vector */
643  free( A );
644 
645  /* now check on validity of the solution, demand that this
646  * simple problem have come within a few epsilons of the
647  * correct answer */
648 
649  /* find largest deviation */
650  rcond = 0.;
651  for(i=0;i<3;++i)
652  {
653  x = fabs( yVector[i]-i-1.);
654  rcond = MAX2( rcond, x );
655  /*printf(" %g ", yVector[i]);*/
656  }
657 
658  if( rcond>DBL_EPSILON)
659  {
660  fprintf(ioQQQ,
661  "SanityCheck found too large a deviation in matrix solver = %g \n",
662  rcond);
663  /* set flag saying that things are not ok */
664  lgOK = false;
665  }
666  /* end matrix inversion check */
667 
668  /* these pointers were set to INT_MIN in ContCreatePointers,
669  * then set to valid numbers in ipShells and OpacityCreate1Element
670  * this checks that all values have been properly filled */
671  for( nelem=0; nelem<LIMELM; ++nelem )
672  {
673  /* must reset state of code after tests performed, remember state here */
674  realnum xIonF[NISO][LIMELM];
675  double hbn[NISO][LIMELM], hn[NISO][LIMELM];
676 
677  if( dense.lgElmtOn[nelem] )
678  {
679  /* set these abundances so that opacities can be checked below */
680  hbn[ipH_LIKE][nelem] = iso_sp[ipH_LIKE][nelem].st[0].DepartCoef();
681  hn[ipH_LIKE][nelem] = iso_sp[ipH_LIKE][nelem].st[0].Pop();
682  xIonF[ipH_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipH_LIKE];
683 
684  iso_sp[ipH_LIKE][nelem].st[0].DepartCoef() = 0.;
685  iso_sp[ipH_LIKE][nelem].st[0].Pop() = 1.;
686  dense.xIonDense[nelem][nelem+1-ipH_LIKE] = 1.;
687 
688  if( nelem > ipHYDROGEN )
689  {
690 
691  hbn[ipHE_LIKE][nelem] = iso_sp[ipHE_LIKE][nelem].st[0].DepartCoef();
692  hn[ipHE_LIKE][nelem] = iso_sp[ipHE_LIKE][nelem].st[0].Pop();
693  xIonF[ipHE_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipHE_LIKE];
694 
695  /* this does not exist for hydrogen itself */
696  iso_sp[ipHE_LIKE][nelem].st[0].DepartCoef() = 0.;
697  iso_sp[ipHE_LIKE][nelem].st[0].Pop() = 1.;
698  dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = 1.;
699  }
700 
701  for( ion=0; ion<=nelem; ++ion )
702  {
703  /* loop over all shells that are defined */
704  for( nshells=0; nshells<Heavy.nsShells[nelem][ion]; ++nshells )
705  {
706  for( j=0; j<3; ++j )
707  {
708  /* >>chng 00 apr 05, array index is on fortran scale so must be
709  * >= 1. This test had been <0, correct for C. Caught by Peter van Hoof */
710  if( opac.ipElement[nelem][ion][nshells][j] <=0 )
711  {
712  /* this is not possible */
713  fprintf(ioQQQ,
714  "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n",
715  nelem , ion , nshells, j );
716  fprintf(ioQQQ,
717  "value was %li \n", opac.ipElement[nelem][ion][nshells][j] );
718  /* set flag saying that things are not ok */
719  lgOK = false;
720  }
721  }
722  }
723 
724  if( nelem > 1 )
725  {
726  vector<realnum> saveion;
727  saveion.resize(nelem+2);
728  /* check that photoionization cross sections are ok */
729  for( j=0; j <= (nelem + 1); j++ )
730  {
731  saveion[j] = dense.xIonDense[nelem][j];
732  dense.xIonDense[nelem][j] = 0.;
733  }
734 
735  dense.xIonDense[nelem][ion] = 1.;
736 
737  OpacityZero();
738  opac.lgRedoStatic = true;
739 
740  /* generate opacity with standard routine - this is the one
741  * called in OpacityAddTotal to make opacities in usual calculations */
742  OpacityAdd1Element(nelem);
743 
744  /* this starts one beyond energy of threshold since cs may be zero there */
745  for( j=Heavy.ipHeavy[nelem][ion]; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
746  {
747  if( opac.opacity_abs[j]+opac.OpacStatic[j] < FLT_MIN )
748  {
749  /* this is not possible */
750  fprintf(ioQQQ,
751  "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n",
752  nelem , ion );
753  fprintf(ioQQQ,
754  "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n",
755  opac.opacity_abs[j] ,
756  opac.OpacStatic[j] ,
757  nelem ,
758  ion ,
759  rfield.anu(j));
760  /* set flag saying that things are not ok */
761  lgOK = false;
762  break;
763  }
764  }
765  /* reset the ionization distribution */
766  for( j=0; j <= (nelem + 1); j++ )
767  {
768  dense.xIonDense[nelem][j] = saveion[j];
769  }
770 
771  }
772  }
773  iso_sp[ipH_LIKE][nelem].st[ipH1s].DepartCoef() = hbn[ipH_LIKE][nelem];
774  iso_sp[ipH_LIKE][nelem].st[ipH1s].Pop() = hn[ipH_LIKE][nelem];
775  dense.xIonDense[nelem][nelem+1-ipH_LIKE] = xIonF[ipH_LIKE][nelem];
776 
777  if( nelem > ipHYDROGEN )
778  {
779  iso_sp[ipHE_LIKE][nelem].st[ipHe1s1S].DepartCoef() = hbn[ipHE_LIKE][nelem];
780  iso_sp[ipHE_LIKE][nelem].st[ipHe1s1S].Pop() = hn[ipHE_LIKE][nelem];
781  dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = xIonF[ipHE_LIKE][nelem];
782  }
783  }
784  }
785 
786 
787  /*********************************************************
788  * *
789  * everything is done, all checks make, did we pass them?*
790  * *
791  *********************************************************/
792 
793  if( lgOK )
794  {
795  /*return if ok */
796  if( trace.lgTrace )
797  {
798  fprintf( ioQQQ, " SanityCheck returns OK\n");
799  }
800  return;
801  }
802 
803  else
804  {
805  /* stop since problem encountered, lgEOF set false */
806  fprintf(ioQQQ , "SanityCheck finds insanity so exiting\n");
807  ShowMe();
809  }
810 }
811 
812 STATIC void SanityCheckGaunt(long Z, double loggam2, double logu, double refval, double relerr)
813 {
814  DEBUG_ENTRY( "SanityCheckGaunt()" );
815  double gam2 = exp10(loggam2);
816  double u = exp10(logu);
817  double ZZ = double(Z);
818  double Te = TE1RYD*pow2(ZZ)/gam2;
819  double anu = pow2(ZZ)*u/gam2;
820  double val = t_gaunt::Inst().gauntff( Z, Te, anu );
821  if( fabs(val/refval-1.) > relerr )
822  {
823  fprintf( ioQQQ, "Gaunt sanity check failed: log(gam2) = %e, log(u) = %e, Z = %ld, "
824  "expected gff = %e, got gff = %e\n", loggam2, logu, Z, refval, val );
826  }
827 }
void OpacityAdd1Element(long int ipZ)
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:222
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
double * OpacStack
Definition: opacity.h:163
double * opacity_abs
Definition: opacity.h:103
qList st
Definition: iso.h:482
double exp10(double x)
Definition: cddefines.h:1383
const int ipHE_LIKE
Definition: iso.h:65
t_opac opac
Definition: opacity.cpp:5
t_Heavy Heavy
Definition: heavy.cpp:5
void CheckMesh() const
Definition: mesh.cpp:322
double * OpacStatic
Definition: opacity.h:122
double iso_recomb_check(long ipISO, long nelem, long level, double temperature)
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:310
double e1(double x)
long int KshellLimit
Definition: continuum.h:98
long int nCollapsed_max
Definition: iso.h:518
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
void SanityCheck(const char *chJob)
const int ipHe1s1S
Definition: iso.h:43
double expn(int n, double x)
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:807
double anu(size_t i) const
Definition: mesh.h:111
void ValidateEdges() const
Definition: mesh.cpp:274
t_dense dense
Definition: global.cpp:15
static t_gaunt & Inst()
Definition: cddefines.h:209
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_trace trace
Definition: trace.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
#define MALLOC(exp)
Definition: cddefines.h:556
multi_arr< double, 2 > BranchRatio
Definition: iso.h:480
long int n_HighestResolved_max
Definition: iso.h:536
#define L_(A_)
Definition: iso.h:23
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
const int ipH1s
Definition: iso.h:29
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
EmissionList::reference Emis() const
Definition: transition.h:447
#define N_(A_)
Definition: iso.h:22
t_continuum continuum
Definition: continuum.cpp:6
const char * chISO[NISO]
Definition: iso.h:348
t_rfield rfield
Definition: rfield.cpp:9
bool lgCaseB
Definition: opacity.h:173
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
STATIC void SanityCheckFinal()
#define cdEXIT(FAIL)
Definition: cddefines.h:484
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1015
double sum(double min, double max) const
Definition: integrate.h:44
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
bool lgElmtOn[LIMELM]
Definition: dense.h:160
double gauntff(long Z, double Te, double anu)
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
const int ipH2p
Definition: iso.h:31
STATIC void SanityCheckBegin()
#define ASSERT(exp)
Definition: cddefines.h:617
double qg32(double, double, double(*)(double))
Definition: service.cpp:1271
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
T pow2(T a)
Definition: cddefines.h:985
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
const int ipHELIUM
Definition: cddefines.h:349
bool lgRedoStatic
Definition: opacity.h:159
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
double e2(double x)
#define S(I_, J_)
long int numLevels_max
Definition: iso.h:524
double sqrte
Definition: phycon.h:58
void ShowMe(void)
Definition: service.cpp:205
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:348
long int nflux
Definition: rfield.h:48
realnum & Aul() const
Definition: emission.h:668
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1318
void OpacityZero(void)
Definition: opacity_zero.cpp:8
STATIC void SanityCheckGaunt(long Z, double loggam2, double logu, double refval, double relerr)
bool lgNoRecombInterp[NISO]
Definition: iso.h:405