cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mesh.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 
4 #include "cddefines.h"
5 #include "trace.h"
6 #include "iso.h"
7 #include "freebound.h"
8 #include "opacity.h"
9 #include "mesh.h"
10 
12 {
13  DEBUG_ENTRY( "p_ReadResolution()" );
14 
15  if( trace.lgTrace )
16  fprintf( ioQQQ,"p_ReadResolution opening continuum_mesh.ini:");
17 
18  FILE* ioDATA = open_data( "continuum_mesh.ini", "r" );
19 
20  char chLine[INPUT_LINE_LENGTH];
21  /* check that magic number is ok */
22  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
23  {
24  fprintf( ioQQQ, "p_ReadResolution could not read first line of continuum_mesh.ini.\n");
26  }
27 
28  long i = 1;
29  bool lgEOL;
30  /* continuum mesh magic number */
31  long i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
32  long i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
33  long i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
34 
35  bool lgResPower;
36 
37  /* the following is the set of numbers that appear at the start of continuum_mesh.ini */
38  if( i1 == 1 && i2 == 9 && i3 == 29 )
39  // old version of the file (c08 and older), this has pairs: upper limit freq range, resolution
40  // this format is still supported to accomodate users with existing continuum_mesh.ini files.
41  lgResPower = false;
42  else if( i1 == 10 && i2 == 8 && i3 == 8 )
43  // new version of the file (c10 and newer), this has pairs: upper limit freq range, resolving power
44  // resolving power = 1./resolution
45  lgResPower = true;
46  else
47  {
48  fprintf( ioQQQ,
49  "p_ReadResolution: the version of continuum_mesh.ini is not supported.\n" );
50  fprintf( ioQQQ,
51  "I found version number %li %li %li.\n" ,
52  i1 , i2 , i3 );
53  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
55  }
56 
57  size_t n = 0;
58  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
59  {
60  /* skip comment lines */
61  if( chLine[0] != '#')
62  {
63  i = 1;
64  double upper_limit = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
65  double val2 = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
66 
67  // continuum energy could be 0 to indicate high energy bound of the code
68  if( upper_limit == 0. )
69  upper_limit = p_egamry;
70 
71  // all must be positive
72  if( upper_limit <= 0. || val2 <= 0. )
73  {
74  fprintf(ioQQQ, "DISASTER PROBLEM continuum_mesh.ini has a non-positive number.\n");
76  }
77 
78  // convert resolving power into resolution
79  double resolution;
80  if( lgResPower )
81  // resolving power was entered
82  resolution = 1./val2;
83  else
84  // resolution was entered
85  resolution = val2;
86 
87  /* this is option to rescale resolution with set resolution command */
88  resolution *= p_ResolutionScaleFactor;
89 
90  // merge in the major ionization edges
91  while( true )
92  {
93  // if a user-supplied mesh edge agrees within 1/10th of a resolution
94  // element with an edge from p_edges[], we will merge and use the latter
95  double toler = upper_limit*resolution*0.1;
96  if( n < p_edges.size()
97  && fp_equal_tol( p_edges[n].Ryd(), upper_limit, toler ) )
98  {
99  p_RangeUpperLimit.push_back(p_edges[n].Ryd());
100  p_RangeResolution.push_back(resolution);
101  ++n;
102  break;
103  }
104  else if( n < p_edges.size()
105  && p_edges[n].Ryd() < upper_limit )
106  {
107  p_RangeUpperLimit.push_back(p_edges[n].Ryd());
108  p_RangeResolution.push_back(resolution);
109  ++n;
110  }
111  else
112  {
113  p_RangeUpperLimit.push_back(upper_limit);
114  p_RangeResolution.push_back(resolution);
115  break;
116  }
117  }
118  }
119  }
120 
121  ASSERT( n == p_edges.size() );
122 
123  fclose( ioDATA );
124 
125  /* now verify continuum grid is ok - first are all values but the last positive? */
126  for( n=1; n < p_RangeUpperLimit.size(); ++n )
127  {
128  if( p_RangeUpperLimit[n-1] >= p_RangeUpperLimit[n] )
129  {
130  fprintf( ioQQQ,
131  "p_ReadResolution: the range upper limits must be in increasing order.\n" );
133  }
134  }
135  if( !fp_equal(p_RangeUpperLimit.back(), p_egamry) )
136  {
137  fprintf( ioQQQ,
138  "p_ReadResolution: the last range upper limit must be zero.\n" );
140  }
141 
142  return;
143 }
144 
145 void t_mesh::p_SetupMesh(bool lgUnitCell)
146 {
147  DEBUG_ENTRY( "p_SetupMesh()" );
148 
149  if( trace.lgTrace )
150  fprintf( ioQQQ, "p_SetupMesh() setting up the frequency mesh. ResolutionScaleFactor = %g.\n",
152 
153  ASSERT( p_RangeUpperLimit.size() > 0 && p_RangeUpperLimit.size() == p_RangeResolution.size() );
154 
155  // reserve some space for the mesh, no problem if we overrun this though...
156  p_anu_edge.reserve( 8500 );
157 
158  // fill in the edges of the frequency cells, these will have a logarithmic spacing
159  // in each of the frequency ranges with a resolving power set by the user.
160  //
161  // the lower bound of the first cell is exactly equal to emm()
162  // the upper bound of the last cell is exactly equal to egamry()
163  //
164  // the anu and widflx arrays will be derived from the cell boundaries
165  // widflx[i] will simply be the difference between the upper and lower bound
166  // and anu[i] will be exactly in the middle between the two bounds.
167  //
168  // lgUnitCell == false: use highest cell of regular mesh for unit tests
169  // the frequency range available for regular RT will NOT reach egamry().
170  // lgUnitCell == true: create additional cell above regular mesh for unit tests
171  // the frequency range available for regular RT will reach egamry().
172  // this choice currently trips at least one ASSERT.
173  double hi_bound = p_emm;
174  p_anu_edge.push_back( hi_bound );
175  for( size_t n=0; n < p_RangeUpperLimit.size(); ++n )
176  {
177  double lo_bound = hi_bound;
178  hi_bound = p_RangeUpperLimit[n];
179  size_t nbin = size_t(ceil(log(hi_bound/lo_bound)/p_RangeResolution[n]));
180  double bound = lo_bound;
181  double step = exp(log(hi_bound/lo_bound)/nbin);
182  for( size_t i=1; i < nbin; ++i )
183  {
184  bound *= step;
185  p_anu_edge.push_back( bound );
186  }
187  p_anu_edge.push_back( hi_bound );
188  }
189 
190  if( lgUnitCell )
191  {
192  // create an extra cell above the upper limit for the unit test
193  double edge = pow2(p_anu_edge[p_anu_edge.size()-1])/p_anu_edge[p_anu_edge.size()-2];
194  p_anu_edge.push_back( edge );
195  }
196 
197  // anu[i] should be exactly in the middle of p_anu_edge[i] and p_anu_edge[i+1]
198  p_anu.resize( p_anu_edge.size()-1 );
199  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
200  p_anu[i] = (p_anu_edge[i]+p_anu_edge[i+1])/2.;
201 
202  p_widflx.resize( p_anu_edge.size()-1 );
203  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
204  p_widflx[i] = p_anu_edge[i+1] - p_anu_edge[i];
205 
206  p_anu2.resize( p_anu_edge.size()-1 );
207  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
208  p_anu2[i] = pow2(p_anu[i]);
209 
210  p_anu3.resize( p_anu_edge.size()-1 );
211  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
212  p_anu3[i] = pow3(p_anu[i]);
213 
214  p_anusqrt.resize( p_anu_edge.size()-1 );
215  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
216  p_anusqrt[i] = sqrt(p_anu[i]);
217 
218  p_anuln.resize( p_anu_edge.size()-1 );
219  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
220  p_anuln[i] = log(p_anu[i]);
221 
222  p_anulog10.resize( p_anu_edge.size()-1 );
223  for( size_t i=0; i < p_anu_edge.size()-1; ++i )
224  p_anulog10[i] = p_anuln[i]/LN_TEN;
225 
226  return;
227 }
228 
229 // Set up special edges that need to be merged into the frequency mesh
231 {
232  DEBUG_ENTRY( "SetupEdges()" );
233 
234  // a list of major ionization edges that need to be fiddled into the frequency mesh
235  // this list will be sorted at the end, so can be entered in any order
236 
237  // =================================================================================
238  //
239  // NB NB - make sure that all edges below match the values used by the code !!!!
240  // this is checked by ValidateEdges() below.
241  // NB NB - this is especially the case for the H I 1s and He II 2s & 2p edges !!!!
242  // NB NB - the values below are not necessarily the most accurate values possible,
243  // they should match what is used elsewhere in the code.
244  // NB NB - to update the values, print the value used by the code with sufficient
245  // precision and ALWAYS ROUND DOWN to 6 significant digits. It is necessary
246  // to round down to avoid epsilon sensitivity when using ipoint() for the
247  // same edge. The 6 digits are needed to avoid problems with FLT_IS_DBL. If
248  // the value used by the code has less than 6 significant digits, you should
249  // still make it slightly smaller (like for the O III edges).
250  //
251  // =================================================================================
252 
253  // H I 1s ionization potential matching the one used in Cloudy, rel. accuracy ~ 6e-6
254  p_edges.push_back( Energy(0.999460, "Ryd") );
255  // H I 2s ionization potential matching the one used in Cloudy
256  p_edges.push_back( Energy(0.249865, "Ryd") );
257  // He I 1s2 ionization potential matching the one used in Cloudy
258  p_edges.push_back( Energy(1.807139, "Ryd") );
259  // He II 1s ionization potential matching the one used in Cloudy, rel. accuracy ~ 8e-5
260  p_edges.push_back( Energy(3.99931, "Ryd") );
261  // He II 2s & 2p ionization potential, matching the one used in Cloudy
262  // this edge needs to be here to assure that the He II 2s & 2p edges are in a higher
263  // cell than the H I 1s edge. The code in cont_createpointers.cpp relies on this!!
264  p_edges.push_back( Energy(0.9998275, "Ryd") );
265  // ionization potential from O III 2p2 1D
266  p_edges.push_back( Energy(3.85499, "Ryd") );
267  // ionization potential from O III 2p2 1S
268  p_edges.push_back( Energy(3.64599, "Ryd") );
269 
270  sort( p_edges.begin(), p_edges.end() );
271 }
272 
273 // Validate that the p_edges contains values in fair agreement with the rest of the code
275 {
276  DEBUG_ENTRY( "ValidateEdges()" );
277 
278  vector<Energy> cl_edges;
279 
280  // Make sure that the entries here match the ones defined above in SetupEdges()!!
281  // These are the edge energies actually used by Cloudy
282  cl_edges.push_back( Energy( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].xIsoLevNIonRyd ) );
283  cl_edges.push_back( Energy( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].xIsoLevNIonRyd ) );
284  cl_edges.push_back( Energy( iso_sp[ipHE_LIKE][ipHELIUM].fb[ipHe1s1S].xIsoLevNIonRyd ) );
285  cl_edges.push_back( Energy( iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].xIsoLevNIonRyd ) );
286  cl_edges.push_back( Energy( iso_sp[ipH_LIKE][ipHELIUM].fb[ipH2s].xIsoLevNIonRyd ) );
287  cl_edges.push_back( Energy( opac.o3exc ) );
288  cl_edges.push_back( Energy( opac.o3exc3 ) );
289 
290  sort( cl_edges.begin(), cl_edges.end() );
291 
292  if( cl_edges.size() != p_edges.size() )
293  {
294  fprintf( ioQQQ, "DISASTER INTERNAL ERROR: the cl_edges and p_edges arrays have different size.\n" );
295  fprintf( ioQQQ, "In mesh.cpp, make sure the entries in these arrays match one on one.\n\n" );
297  }
298  bool lgErr = false;
299  for( unsigned long i=0; i < p_edges.size(); ++i )
300  {
301  double rel_diff = 1. - p_edges[i].Ryd()/cl_edges[i].Ryd();
302  // the entry in p_edge[i] should be slightly smaller than cl_edge[i]
303  // this avoids problems when calculating ipoint(cl_edge[i])...
304  // it should be at least a factor 1.5*FLT_EPSILON smaller to avoid problems
305  // when using -DFLT_IS_DBL, but is should also not be too large...
306  if( rel_diff < 1.5*double(FLT_EPSILON) || rel_diff > 1.e-5 )
307  {
308  double expected = cl_edges[i].Ryd()*(1. - 5.*double(FLT_EPSILON));
309  fprintf( ioQQQ, "DISASTER INTERNAL ERROR: energy mismatch for p_edges[].\n" );
310  fprintf( ioQQQ, "I expected %.8g but found %.8g.\n", expected, p_edges[i].Ryd() );
311  fprintf( ioQQQ, "Please update p_edges[] to hold %.8g in mesh.cpp.\n\n", expected );
312  lgErr = true;
313  }
314  }
315  if( lgErr )
316  {
318  }
319 }
320 
321 /*CheckMesh perform sanity check confirming that the energy array has been properly filled */
322 void t_mesh::CheckMesh() const
323 {
324  DEBUG_ENTRY( "CheckMesh()" );
325 
326  // is the grid initialized?
327  ASSERT( p_anu.size() > 0 );
328 
329  // check outer bounds of the complete grid
330  ASSERT( fp_equal( anumin(0), emm() ) );
331  ASSERT( fp_equal( anumax(ncells()-1), egamry() ) );
332 
333  bool lgFail = false;
334  double hi_bound = emm();
335  for( size_t n=0; n < p_RangeUpperLimit.size(); ++n )
336  {
337  /* test middle of energy range */
338  double lo_bound = hi_bound;
339  hi_bound = p_RangeUpperLimit[n];
340 
341  double energy = lo_bound*0.5 + hi_bound*0.5;
342  size_t ic = ipointC(energy);
343  if( energy < anumin(ic) )
344  {
345  fprintf( ioQQQ, "CheckMesh middle test low fail\n" );
346  lgFail = true;
347  }
348  else if( energy > anumax(ic) )
349  {
350  fprintf( ioQQQ, "CheckMesh middle test high fail\n" );
351  lgFail = true;
352  }
353 
354  /* test near low bound */
355  energy = lo_bound*0.99 + hi_bound*0.01;
356  ic = ipointC(energy);
357  if( energy < anumin(ic) )
358  {
359  fprintf( ioQQQ, "CheckMesh low test low fail\n" );
360  lgFail = true;
361  }
362  else if( energy > anumax(ic) )
363  {
364  fprintf( ioQQQ, "CheckMesh low test high fail\n" );
365  lgFail = true;
366  }
367 
368  /* test near high bound */
369  energy = lo_bound*0.01 + hi_bound*0.99;
370  ic = ipointC(energy);
371  if( energy < anumin(ic) )
372  {
373  fprintf( ioQQQ, "CheckMesh high test low fail\n" );
374  lgFail = true;
375  }
376  else if( energy > anumax(ic) )
377  {
378  fprintf( ioQQQ, "CheckMesh high test high fail\n" );
379  lgFail = true;
380  }
381  }
382 
383  if( lgFail )
384  {
385  fprintf( ioQQQ, "CheckMesh: sanity check on frequency mesh failed.\n Bailing out\n" );
387  }
388 
389  return;
390 }
double emm() const
Definition: mesh.h:84
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:751
vector< double > p_RangeUpperLimit
Definition: mesh.h:31
const int ipHE_LIKE
Definition: iso.h:65
double o3exc
Definition: opacity.h:295
t_opac opac
Definition: opacity.cpp:5
void CheckMesh() const
Definition: mesh.cpp:322
vector< double > p_anulog10
Definition: mesh.h:47
double p_egamry
Definition: mesh.h:19
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:908
T pow3(T a)
Definition: cddefines.h:992
FILE * ioQQQ
Definition: cddefines.cpp:7
const int ipHe1s1S
Definition: iso.h:43
double o3exc3
Definition: opacity.h:296
void ValidateEdges() const
Definition: mesh.cpp:274
vector< double > p_RangeResolution
Definition: mesh.h:32
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_trace trace
Definition: trace.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
vector< double > p_anu3
Definition: mesh.h:51
size_t ipointC(double anu) const
Definition: mesh.h:152
void p_SetupMesh(bool lgUnitCell)
Definition: mesh.cpp:145
double energy(const genericState &gs)
const int ipH1s
Definition: iso.h:29
vector< double > p_widflx
Definition: mesh.h:44
bool lgTrace
Definition: trace.h:12
vector< double > p_anu
Definition: mesh.h:38
double p_emm
Definition: mesh.h:16
#define EXIT_FAILURE
Definition: cddefines.h:168
vector< Energy > p_edges
Definition: mesh.h:35
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
#define cdEXIT(FAIL)
Definition: cddefines.h:484
vector< double > p_anuln
Definition: mesh.h:48
vector< double > p_anu2
Definition: mesh.h:50
double anumin(size_t i) const
Definition: mesh.h:139
vector< double > p_anu_edge
Definition: mesh.h:41
#define ASSERT(exp)
Definition: cddefines.h:617
const int ipH2s
Definition: iso.h:30
double p_ResolutionScaleFactor
Definition: mesh.h:23
vector< double > p_anusqrt
Definition: mesh.h:49
const int ipH_LIKE
Definition: iso.h:64
T pow2(T a)
Definition: cddefines.h:985
Definition: energy.h:9
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
const int ipHELIUM
Definition: cddefines.h:349
double egamry() const
Definition: mesh.h:88
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
long ncells() const
Definition: mesh.h:80
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:72
void p_SetupEdges()
Definition: mesh.cpp:230
double anumax(size_t i) const
Definition: mesh.h:143
const int ipHYDROGEN
Definition: cddefines.h:348
void p_ReadResolution()
Definition: mesh.cpp:11
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:393