cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
eden_sum.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue
4  *called by ConvBase - ConvEdenIoniz actually updates the electron density
5  * returns 0 if all is ok, 1 if need to abort calc */
6 #include "cddefines.h"
7 #include "hmi.h"
8 #include "trace.h"
9 #include "grainvar.h"
10 #include "rfield.h"
11 #include "mole.h"
12 #include "dense.h"
13 #include "conv.h"
14 
15 /*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue
16  *called by ConvBase - ConvEdenIoniz actually updates the electron density
17  * returns 0 if all is ok, 1 if need to abort calc */
18 int eden_sum(void)
19 {
20  long int i,
21  ion,
22  nelem;
23 
24  realnum sum_all_ions ,
25  sum_metals ,
26  hmole_eden,
27  eden_ions[LIMELM];
28 
29  DEBUG_ENTRY( "eden_sum()" );
30 
31  /* EdenExtra is normally zero, set with EDEN command, to add extra e- */
33 
34  /* sum over all ions */
35  sum_all_ions = 0.;
36  sum_metals = 0.;
37  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
38  {
39  if( nelem==ipLITHIUM )
40  sum_metals = 0.;
41  eden_ions[nelem] = dense.xIonDense[nelem][1];
42  /* >>chng 02 jul 15, upper limit now includes H-like, so all species
43  * in this one loop (except hydrogen and helium ) */
44  /* >>chng 02 apr 26, upper limit had been nelem, not nelem+1, so H-like
45  * metals were missed - ok for H and He since not part of this loop before */
46  for( ion=2; ion <= nelem+1; ion++ )
47  {
48  /* >>chng 96 oct 27, save all contributors to electron density */
49  eden_ions[nelem] += ion*dense.xIonDense[nelem][ion];
50  }
51 
52  sum_all_ions += eden_ions[nelem];
53  sum_metals += eden_ions[nelem];
54  }
55  dense.EdenTrue += sum_all_ions;
56  ASSERT( dense.EdenTrue >= 0. );
57 
58  /* add on molecules */
59  co.comole_eden = 0.;
60  for( i=0; i < mole.num_comole_calc; i++ )
61  {
62  if(COmole[i]->n_nuclei != 1)
63  {
65  }
66  }
67 
68  /* electrons contributed by molecules */
70  ASSERT( dense.EdenTrue >= 0. );
71 
72  /* >>chng 03 nov 28, loop over all H molecules, had just been H- */
73  hmole_eden = 0.;
74  for(i=0;i<N_H_MOLEC;i++)
75  {
76  /* hmi.nElectron is zero for H+ since counted as an ion, -1 for H-, etc */
77  hmole_eden += hmi.Hmolec[i]*hmi.nElectron[i];
78  }
79  /* >>chng 01 jan 08, following logic added to stop H- from forcing
80  * negative electron density during very first search, when H- is badly off */
81  /* >>chng 03 nov 28, add test for lgSearch, had been absent */
82  /* >>chng 04 feb 20, recoil_ion.in shows this important,
83  * update test so prevent neg eden */
84  if( (-hmole_eden) > dense.EdenTrue/4. && conv.lgSearch )
85  {
86  /*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,hmole_eden);*/
87  dense.EdenTrue *= 0.9;
88  }
89  else if( (-hmole_eden) > dense.EdenTrue )
90  {
91  /* >>chng 04 mar 14, add this branch to prevent
92  * negative electron density. occurred in pdr
93  * with large h2, after hmole failure */
94  fprintf(ioQQQ," PROBLEM sum eden from hmole too neg, set to limt. EdenTrue:%.3e hmole_eden:%.3e \n",
95  dense.EdenTrue,
96  hmole_eden);
98  }
99  else
100  {
101  /* account for electrons on all H-bearing molecules */
102  dense.EdenTrue += hmole_eden;
103  }
104  ASSERT(dense.EdenTrue >= 0. );
105 
106  /* this variable is set with the set eden command,
107  * is supposed to override physical electron density */
108  if( dense.EdenSet > 0. )
109  {
111  dense.eden_from_metals = 1.;
112  }
113  else
114  {
115  /* fraction of electrons from si, s, mg, al */
116  /*dense.eden_from_metals = sum_all_ions / SDIV(dense.EdenTrue);*/
117  dense.eden_from_metals = sum_metals / SDIV(dense.EdenTrue);
118  /*fprintf(ioQQQ," debuggg ionfrac %.2f %.3f\n",fnzone,dense.eden_from_metals );*/
119  }
120  /* dense.eden itself is actually changed in ConvEdenIoniz */
121 
122  /* >>chng 00 dec 19, include electrons on grains in total sum */
123  /* >>chng 02 jul 15, only add grain electrons when we have stable soln -
124  * everett.in had very negative grain elec total, so drove total eden negative,
125  * but due to bad inital electron density - do not add grain electrons until we
126  * are close to a solution */
127  /* >>chng 04 sep 26, allow nEdenTrue to become -ve */
129  {
130  /* gv.lgGrainElectrons - should grain electron source/sink be included in overall electron sum?
131  * default is true, set false with no grain electrons command */
133  }
134  else
135  {
136  /* >>chng 05 jan 24, only print if not in search phase */
137  if( !conv.lgSearch )
138  fprintf(ioQQQ,
139  " PROBLEM eden grain neg limt %.2f eden %.4e new eden bef grn %.4e"
140  "grain eden %.4e set edentrue to %.4e Search?%c\n",
141  fnzone,
142  dense.eden,
143  dense.EdenTrue,
144  gv.TotalEden,
145  dense.eden*0.9,
146  TorF(conv.lgSearch));
147 
148  /*dense.EdenTrue = dense.eden*0.9;*/
150  }
151 
152  if( trace.lgTrace || trace.lgESOURCE )
153  {
154  fprintf( ioQQQ,
155  " eden_sum zn:%.2f current:%.4e new true:%.4e ions:%.4e comole:%.4e hmole:%.4e grain:%.4e extra:%.4e sum:%.4e LaOTS:%.4e\n",
156  fnzone ,
157  dense.eden ,
158  dense.EdenTrue ,
159  sum_all_ions ,
160  co.comole_eden ,
161  hmole_eden ,
163  dense.EdenExtra ,
164  sum_all_ions + co.comole_eden + hmole_eden + gv.TotalEden*gv.lgGrainElectrons + dense.EdenExtra,
165  rfield.otslin[2182] );
166 
167  if( trace.lgNeBug )
168  {
169  for(nelem=ipHYDROGEN; nelem < LIMELM; nelem++)
170  {
171  if( nelem==0 )
172  {
173  fprintf( ioQQQ, " eden_sum H -Ne:" );
174  }
175  else if( nelem==10 )
176  {
177  fprintf( ioQQQ, "\n eden_sum Na-Ca:" );
178  }
179  else if( nelem==20 )
180  {
181  fprintf( ioQQQ, "\n eden_sum Sc-Zn:" );
182  }
183  fprintf( ioQQQ, " %.4e", eden_ions[nelem] );
184  if( nelem==29 )
185  {
186  fprintf( ioQQQ, "\n" );
187  }
188  /*if( nelem==9 || nelem==19 || nelem==29 )
189  fprintf( ioQQQ, "\n " );*/
190  }
191  }
192  }
193 
194  /* abort if negative electron density */
195  /*>>chng 04 sep 25, allow negative true elec density so solver can work */
196  if( 0 && dense.EdenTrue < 0. )
197  {
198  fprintf( ioQQQ, "eden_sum finds non-positive electron density.\n" );
199  fprintf( ioQQQ,
200  " eden_sum: EdenTrue to%12.4e fm%12.4e Ne(H):%10.2e Ne(He):%10.2e Ne(C):\n",
201  dense.EdenTrue,
202  dense.eden,
205  ShowMe();
206  cdEXIT(EXIT_FAILURE);
207  }
208  else if( dense.EdenTrue == 0. )
209  {
210  fprintf( ioQQQ, "\nDISASTER PROBLEM eden_sum finds an electron density of zero, this is unphysical.\n" );
211  fprintf( ioQQQ, "There appears to be no source of ionization.\n");
212  fprintf( ioQQQ, "Consider adding background cosmic rays with COSMIC RAY BACKGROUND and BACKGROUND commands.\n\n");
213  lgAbort = true;
214 
215  return 1;
216  }
217 
218  {
219  /*@-redef@*/
220  enum {DEBUG_LOC=false};
221  /*@+redef@*/
222  if( DEBUG_LOC )
223  {
224  fprintf(ioQQQ,"esumdebugg\t%li\t%.2e\t%.2e\n",
225  nzone,
227  }
228  }
229 
230  /* >>chng 05 jan 05, don't let elec den be zero - logs are taken */
232 
233  return 0;
234 }

Generated for cloudy by doxygen 1.8.4