cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cddefines.h"
4 #include "mole.h"
5 #include "parser.h"
6 
9 
11 {
12  DEBUG_ENTRY( "t_mole_global::zero()" );
13  /* flag to turn off molecular network */
14  lgNoMole = false;
15  lgNoHeavyMole = false;
16  /* capture of molecules onto grain surfaces - formation of ices
17  * flag says to include this process - turned off with the
18  * command NO GRAIN MOLECULES */
19  lgGrain_mole_deplete = true;
20  /* flag saying that H2O water destruction rate went to zero */
21  lgH2Ozer = false;
22  /* option to turn on the UMIST rates, naturally this will be 1, set to zero
23  with the set UMIST rates command */
24  lgLeidenHack = false;
25  /* option to use diffuse cloud chemical rates from Table 8 of
26  * >> refer Federman, S. R. & Zsargo, J. 2003, ApJ, 589, 319
27  * By default, this is false - changed with set chemistry command */
28  lgFederman = true;
29  /* option to use effective temperature as defined in
30  * >> refer Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
31  * By default, this is false - changed with set chemistry command */
32  lgNonEquilChem = false;
36  lgProtElim = true;
40  lgNeutrals = true;
41  /* option to use H2 continuum dissociation cross sections computed by P.C. Stancil
42  * By default, this is true - changed with "set H2 continuum dissociation xxx" command
43  * options are "Stancil" or "AD69" */
44  lgStancil = false;
45  // all isotopes are currently disabled by default
46  lgTreatIsotopes.resize( LIMELM );
47  fill( lgTreatIsotopes.begin(), lgTreatIsotopes.end(), false );
48 
49 }
50 /*=================================================================*/
51 /*mole_Init called from cdInit to initialize CO routines */
53 {
54  DEBUG_ENTRY( "Mole::init()" );
55 
56  /* prevent memory leaks */
57  /* \todo this is a temporary fix for PR14. We should improve the overall design
58  * of this code to prevent valid pointers being overwritten in a second call to mole_Init */
59  static bool lgmole_Init_called=false;
60  static long int num_total_MALLOC=-1;
61  if(! lgmole_Init_called )
62  {
63  /* say that we have been called */
64  lgmole_Init_called = true;
65 
66  make_species();
69  mole.species.resize( num_total );
70  num_total_MALLOC = num_total;
71  }
72 
73  if( num_total>num_total_MALLOC )
74  {
75  /* number of species has increased since last time - this can't happen
76  * tsuite / programs / comp4 has 95 first time, 98 second time */
77  fprintf(ioQQQ,"DISASTER - the number of species in the chemistry network has increased. This is not allowed.\n");
78  fprintf(ioQQQ,"This could happen if an element was initially turned off or grains not included, then the element or grains was included. There are not allowed.\n");
79  fprintf(ioQQQ,"Sorry.\n");
81  }
82 
83  for( long i=0; i<num_total; ++i )
84  {
85  mole.species[i].zero();
86  mole.species[i].index = i;
87  }
88  mole.elec = 0.;
89 }
90 
92 {
93  /* these are source and sink terms for heavy element ionization balance from the
94  * chemistry */
95  source =
96  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
97  sink =
98  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
99  xMoleChTrRate =
100  (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
101 
102  for( long nelem=0; nelem<LIMELM; ++nelem )
103  {
104  /* chemistry source and sink terms for ionization ladders */
105  source[nelem] =
106  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
107  sink[nelem] =
108  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
109  xMoleChTrRate[nelem] =
110  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
111  for( long ion=0; ion<nelem+2; ++ion )
112  {
113  xMoleChTrRate[nelem][ion] =
114  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
115  }
116  }
117 }
118 
120 {
121  for( long nelem=0; nelem< LIMELM; ++nelem )
122  {
123  /* these have one more ion than above */
124  for( long ion=0; ion<nelem+2; ++ion )
125  {
126  /* zero out the source and sink arrays */
127  source[nelem][ion] = 0.;
128  sink[nelem][ion] = 0.;
129  for( long ion2=0; ion2<nelem+2; ++ion2 )
130  {
131  xMoleChTrRate[nelem][ion][ion2] = 0.;
132  }
133  }
134  }
135  for( long i=0; i < mole_global.num_calc; i++ )
136  {
137  species[i].column = 0.;
138  }
139 
140  fill( reaction_rks.begin(), reaction_rks.end(), 0. );
141  fill( old_reaction_rks.begin(), old_reaction_rks.end(), 0. );
142 }
143 
145 {
146  DEBUG_ENTRY( "ParseChemistry()" );
147 
148  Symbol s = p.getSymbol();
149  if (s.toktype == Symbol::EOSTAT ||
150  s.toktype == Symbol::ERROR)
151  {
152  fprintf(ioQQQ,"Error, CHEMISTRY requires option to be specified\n");
154  }
155 
156  for(;;)
157  {
158  if (s.toktype == Symbol::NAME)
159  {
160  if (!strncmp("REAC",s.value.c_str(),4))
161  {
162  s = p.getSymbol();
163  vector<string> reactions;
164  if (s.toktype == Symbol::STRING)
165  {
166  reactions.push_back(s.value);
167  }
168  else if (s.toktype == Symbol::OPERATOR && s.value == "(")
169  {
170  for (;;)
171  {
172  s = p.getSymbol();
173  if (s.toktype != Symbol::STRING)
174  {
175  fprintf(ioQQQ,"No reactions found for CHEMistry REACtion command\n");
176  fprintf(ioQQQ,"Reactions needs to be included in quotes \"\"\n");
178  }
179  reactions.push_back(s.value);
180  s = p.getSymbol();
181  if (s.toktype == Symbol::OPERATOR)
182  {
183  if (s.value == ")")
184  {
185  break;
186  }
187  else if (s.value != ",")
188  {
189  fprintf(ioQQQ,"Syntax error for CHEMistry REACtion command\n");
191  }
192  }
193  else
194  {
195  fprintf(ioQQQ,"Syntax error for CHEMistry REACtion command\n");
197  }
198  }
199  }
200  else
201  {
202  fprintf(ioQQQ,"No reaction found for CHEMistry REACtion command\n");
203  fprintf(ioQQQ,"Reaction needs to be included in quotes \"\"\n");
205  }
206  s = p.getSymbol();
207  if (s.toktype == Symbol::NAME)
208  {
209  if (s.value == "OFF")
210  {
211  for (vector<string>::iterator it = reactions.begin();
212  it != reactions.end(); ++it)
213  {
214  mole_global.offReactions[*it] = true;
215  }
216  }
217  else
218  {
219  fprintf(ioQQQ,"Error, CHEMISTRY REACTION option %s "
220  "not known\n",s.value.c_str());
222  }
223  }
224  else
225  {
226  fprintf(ioQQQ,"Error, CHEMISTRY REACTION needs action\n");
228  }
229  }
230  }
231  else
232  {
233  fprintf(ioQQQ,"Error, unknown option to CHEMISTRY: %s\n",
234  s.value.c_str());
236  }
237  s = p.getSymbol();
238  if (s.toktype == Symbol::EOSTAT ||
239  s.toktype == Symbol::ERROR)
240  {
241  break;
242  }
243  if (s.toktype != Symbol::OPERATOR || s.value != ",")
244  {
245  fprintf(ioQQQ,"Syntax error, CHEMISTRY option [, option]\n");
247  }
248  s = p.getSymbol();
249  }
250 }
251 
t_mole_global mole_global
Definition: mole.cpp:7
bool lgStancil
Definition: mole.h:337
vector< double > reaction_rks
Definition: mole.h:470
vector< bool > lgTreatIsotopes
Definition: mole.h:359
int num_total
Definition: mole.h:362
int num_calc
Definition: mole.h:362
void make_species(void)
void zero()
Definition: mole.cpp:119
bool lgH2Ozer
Definition: mole.h:331
void mole_make_list(void)
bool lgProtElim
Definition: mole.h:347
FILE * ioQQQ
Definition: cddefines.cpp:7
Definition: parser.h:43
Symbol getSymbol()
Definition: parser.cpp:824
bool lgNonEquilChem
Definition: mole.h:342
map< string, bool > offReactions
Definition: mole.h:367
#define MALLOC(exp)
Definition: cddefines.h:554
double ** sink
Definition: mole.h:464
void init(void)
Definition: mole.cpp:52
vector< double > old_reaction_rks
Definition: mole.h:471
bool lgLeidenHack
Definition: mole.h:334
void zero()
Definition: mole.cpp:10
double ** source
Definition: mole.h:464
t_mole_local mole
Definition: mole.cpp:8
realnum *** xMoleChTrRate
Definition: mole.h:466
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
bool lgNeutrals
Definition: mole.h:352
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
bool lgGrain_mole_deplete
Definition: mole.h:356
void alloc()
Definition: mole.cpp:91
string value
Definition: parser.h:38
double elec
Definition: mole.h:460
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
enum tokens toktype
Definition: parser.h:37
bool lgNoHeavyMole
Definition: mole.h:328
bool lgFederman
Definition: mole.h:336
bool lgNoMole
Definition: mole.h:325
Definition: parser.h:34
void mole_make_groups(void)
void ParseChemistry(Parser &p)
Definition: mole.cpp:144