cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
init_sim_postparse.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 /*InitSimPostparse initialization at start of simulation,
4  * called from cloudy after parser,
5  * called one time per sim in grid
6  * not called for every iteration */
7 #include "cddefines.h"
8 #include "init.h"
9 #include "dense.h"
10 #include "struc.h"
11 #include "elementnames.h"
12 #include "atoms.h"
13 #include "iterations.h"
14 #include "pressure.h"
15 #include "trace.h"
16 #include "radius.h"
17 #include "thermal.h"
18 #include "heavy.h"
19 #include "wind.h"
20 #include "iso.h"
21 #include "h2.h"
22 #include "monitor_results.h"
23 #include "taulines.h"
24 #include "mole.h"
25 #include "rfield.h"
26 #include "continuum.h"
27 
28 /*InitSimPostparse initialization after parser,
29  * called one time for each simulation in a grid,
30  * only before first iteration
31  * sets initial or zeros values before start of a simulation */
32 void InitSimPostparse( void )
33 {
34  DEBUG_ENTRY( "InitSimPostparse()" );
35  static double one=1.0;
36 
38 
39  thermal.thist = 0.;
40  thermal.tlowst = 1e20f;
41  thermal.nUnstable = 0;
42  thermal.lgUnstable = false;
43 
44  colliders.init();
45 
46  mole_global.init();
47 
49 
50  findspecieslocal("e-")->location = &(dense.eden);
52  findspecieslocal("PHOTON")->location = &one;
53  findspecieslocal("CRPHOT")->location = &one;
54  findspecieslocal("CRP")->location = &one;
55 
56  /* >> chng 06 Nov 24 rjrw: Move reaction definitions here so they can be controlled by parsed commands */
58 
59  //mole_cmp_num_in_out_reactions();
60 
61  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
62  (*diatom)->H2_Reset();
63 
64  /* read LTE radiative cooling, if small model in use */
65  if( ! h2.lgEnabled )
67 
68  /* zero out diffuse Lya emission */
69  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
70  for( long ion=0; ion<=nelem; ++ion )
71  Heavy.xLyaHeavy[nelem][ion] = 0;
72 
73  /* convert STOP RADIUS command to STOP THICKNESS command now that inner radius is known */
74  if( iterations.StopRadius[0] > 0. )
75  {
76  for( long j=0; j < iterations.iter_malloc; j++ )
77  {
80  else
81  {
82  fprintf( ioQQQ, " PROBLEM stopping radius is <= inner radius. Bailing out.\n" );
84  }
85  }
86  }
87 
88  /* this term applies centrifugal acceleration if DISK option on wind command */
89  wind.DiskRadius = 0;
90  if( wind.lgDisk )
92 
93  iterations.lgLastIt = false;
94 
96 
97  rfield.lgUSphON = false;
98 
99  /* where cloud is optically thick to brems */
101  rfield.EnergyBremsThin = 0.;
102  rfield.comtot = 0.;
103  rfield.cmcool = 0.;
104  rfield.cinrat = 0.;
109  rfield.EnerGammaRay = 7676.;
110 
111  for( vector<TransitionList>::iterator it = AllTransitions.begin(); it != AllTransitions.end(); ++it )
112  {
113  for( TransitionList::iterator tr = it->begin(); tr != it->end(); ++tr )
114  {
115  tr->Emis().TauTrack().clear();
116  }
117  }
118 
119  /* usually computed in pressure_change, but not called before first zone */
120  if (wind.comass == 0.0)
121  {
122  wind.AccelGravity = 0.f;
123  }
124  else
125  {
126  if (radius.Radius < 0.0)
127  {
128  fprintf(ioQQQ,"PROBLEM: must have positive initial radius for central gravity model\n");
130  }
131  wind.AccelGravity = (realnum)(GRAV_CONST*wind.comass*SOLAR_MASS/POW2(radius.Radius)*
132  (1.-wind.DiskRadius/radius.Radius) );
133  }
134  if( trace.lgTrace && trace.lgWind )
135  {
136  fprintf(ioQQQ," InitSimPostparse sets AccelGravity %.3e lgDisk?%c\n",
137  wind.AccelGravity ,
138  TorF(wind.lgDisk) );
139  }
140 
141  /* pressure related variables */
142  pressure.PresInteg = 0.;
144  pressure.pinzon = 0.;
145 
146  /* update iso sequence level information */
147  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
148  {
149  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
150  {
151  iso_sp[ipISO][nelem].Reset();
152  /* these elements that are turned off */
153  if( nelem>ipHELIUM && !dense.lgElmtOn[nelem] )
154  {
155  iso_sp[ipISO][nelem].numLevels_max = 0;
156  iso_sp[ipISO][nelem].nCollapsed_max = 0;
157  iso_sp[ipISO][nelem].n_HighestResolved_max = 0;
158 
159  iso_sp[ipISO][nelem].numLevels_local = 0;
160  iso_sp[ipISO][nelem].nCollapsed_local = 0;
161  iso_sp[ipISO][nelem].n_HighestResolved_local = 0;
162  }
163  else
164  {
165  iso_update_num_levels( ipISO, nelem );
166  /* must always have at least one collapsed level, unless compiling recomb data file. */
167  ASSERT( iso_sp[ipISO][nelem].nCollapsed_max > 0 || iso_ctrl.lgCompileRecomb[ipISO] );
168  }
169  }
170  }
171 
173  {
174  fprintf(ioQQQ,"\n\n Number of levels in ions treated by iso sequences.\n");
175  fprintf(ioQQQ," ISO Element hi-n(l-resolved) #(l-resolved) n(collapsed)\n");
176  /* option to print number of levels for each element */
177  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
178  {
179  for( long nelem=ipISO; nelem<LIMELM; ++nelem )
180  {
181  /* print number of levels */
182  fprintf(ioQQQ," %s %s %4li %4li %4li \n",
183  iso_ctrl.chISO[ipISO] ,
184  elementnames.chElementSym[nelem],
185  iso_sp[ipISO][nelem].n_HighestResolved_max,
186  iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max,
187  iso_sp[ipISO][nelem].nCollapsed_max);
188  }
189  }
190  }
191  atoms.p2nit = 0.;
192  atoms.d5200r = 0.;
193  atoms.rateMg2 = 0.;
194 
197 
198  /* Initialize pseudo-continua, if requested */
200 
201  /* Initialize species bands, if requested */
203 
204  return;
205 }
void InitSimPostparse(void)
t_mole_global mole_global
Definition: mole.cpp:7
long int iter_malloc
Definition: iterations.h:40
double grain_area
Definition: mole.h:457
double Radius
Definition: radius.h:31
t_thermal thermal
Definition: thermal.cpp:6
void Reset()
Definition: iso.cpp:113
realnum PresInteg
Definition: pressure.h:69
void mole_create_react(void)
double comtot
Definition: rfield.h:275
realnum EnerGammaRay
Definition: rfield.h:448
long int ipEnergyBremsThin
Definition: rfield.h:226
t_struc struc
Definition: struc.cpp:6
t_Heavy Heavy
Definition: heavy.cpp:5
bool lgDisk
Definition: wind.h:74
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
t_monitorresults MonitorResults
realnum EnergyBremsThin
Definition: rfield.h:227
const int NISO
Definition: cddefines.h:311
bool lgCompileRecomb[NISO]
Definition: iso.h:400
realnum PresIntegElecThin
Definition: pressure.h:75
void rfield_opac_zero(long lo, long ihi)
long int nSumErrorCaseMonitor
char TorF(bool l)
Definition: cddefines.h:749
long int nCollapsed_max
Definition: iso.h:518
realnum xLyaHeavy[LIMELM][LIMELM]
Definition: heavy.h:21
void init()
Definition: collision.cpp:49
realnum d5200r
Definition: atoms.h:126
void SpeciesPseudoContCreate()
vector< double > StopThickness
Definition: iterations.h:77
FILE * ioQQQ
Definition: cddefines.cpp:7
molezone * findspecieslocal(const char buf[])
bool lgWind
Definition: trace.h:106
t_dense dense
Definition: global.cpp:15
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int nflux_with_check
Definition: rfield.h:49
Wind wind
Definition: wind.cpp:5
long int n_HighestResolved_local
Definition: iso.h:538
t_trace trace
Definition: trace.cpp:5
double SumErrorCaseMonitor
bool lgUnstable
Definition: thermal.h:65
long int nUnstable
Definition: thermal.h:64
realnum AccelGravity
Definition: wind.h:49
realnum pinzon
Definition: pressure.h:69
realnum rateMg2
Definition: atoms.h:129
void init(void)
Definition: mole.cpp:52
ColliderList colliders(dense)
void iso_update_num_levels(long ipISO, long nelem)
long int n_HighestResolved_max
Definition: iso.h:536
#define POW2
Definition: cddefines.h:979
bool lgTrace
Definition: trace.h:12
bool lgEnabled
Definition: h2_priv.h:352
realnum tlowst
Definition: thermal.h:68
const char * chISO[NISO]
Definition: iso.h:348
t_mole_local mole
Definition: mole.cpp:8
double cmcool
Definition: rfield.h:275
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
t_atoms atoms
Definition: atoms.cpp:5
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
bool lgPrintNumberOfLevels
Definition: iso.h:346
vector< diatomics * > diatoms
Definition: h2.cpp:8
realnum thist
Definition: thermal.h:68
#define cdEXIT(FAIL)
Definition: cddefines.h:482
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
t_iterations iterations
Definition: iterations.cpp:6
t_radius radius
Definition: radius.cpp:5
bool lgElmtOn[LIMELM]
Definition: dense.h:160
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
double extin_mag_V_point
Definition: rfield.h:258
void SpeciesBandsCreate()
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgLastIt
Definition: iterations.h:47
realnum p2nit
Definition: atoms.h:126
double extin_mag_B_point
Definition: rfield.h:258
double extin_mag_V_extended
Definition: rfield.h:262
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
void H2_Read_LTE_cooling_per_H2()
const int ipHELIUM
Definition: cddefines.h:350
double eden
Definition: dense.h:201
double extin_mag_B_extended
Definition: rfield.h:262
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
vector< TransitionList > AllTransitions
Definition: taulines.cpp:9
long int numLevels_max
Definition: iso.h:524
vector< double > StopRadius
Definition: iterations.h:80
double DiskRadius
Definition: wind.h:78
bool lgUSphON
Definition: rfield.h:355
long int nzonePreviousIteration
Definition: struc.h:22
void set_ion_locations()
const int ipHYDROGEN
Definition: cddefines.h:349
double * location
Definition: mole.h:411
long int numLevels_local
Definition: iso.h:529
long int nCollapsed_local
Definition: iso.h:519
realnum comass
Definition: wind.h:14
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
double cinrat
Definition: rfield.h:275