cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
heat_save.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 /*SaveHeat save contributors to local heating, with save heat command, called by punch_do */
4 #include "cddefines.h"
5 #include "thermal.h"
6 #include "radius.h"
7 #include "conv.h"
8 #include "dense.h"
9 #include "taulines.h"
10 #include "phycon.h"
11 #include "elementnames.h"
12 #include "dynamics.h"
13 #include "save.h"
14 #include "lines_service.h"
15 
16 /*SaveHeat save contributors to local heating, with save heat command,
17  * called by punch_do */
18 void SaveHeat(FILE* io)
19 {
20  DEBUG_ENTRY( "SaveHeat()" );
21 
22  vector<realnum>SaveVal(LIMELM*LIMELM, -FLT_MAX);
23  vector<long>ipsave(LIMELM*LIMELM, INT_MIN);
24  vector<long>jpsave(LIMELM*LIMELM, INT_MIN);
25  vector<long>ipOrdered(LIMELM*LIMELM);
26  vector<string> chLabel(LIMELM*LIMELM);
27 
28  double cool_total = thermal.ctot;
29  double heat_total = thermal.htot;
30 
31  /* >>chng 06 mar 17, comment out following block and replace with this
32  * removing dynamics heating & cooling and report only physical
33  * heating and cooling
34  * NB the heating and cooling as punched no longer need be
35  * equal for a converged model */
36  cool_total -= dynamics.Cool();
37  heat_total -= dynamics.Heat();
38  if( false )
39  {
40  if(dynamics.Cool() > dynamics.Heat())
41  {
42  cool_total -= dynamics.Heat();
43  heat_total -= dynamics.Heat();
44  }
45  else
46  {
47  cool_total -= dynamics.Cool();
48  heat_total -= dynamics.Cool();
49  }
50  }
51 
52  long ipnt = 0;
53 
54  /* heat sources are saved in a 2d square array */
55  /* WeakHeatCool set with 'set weakheatcool' command
56  * default is 0.05 */
57  for( long i=0; i < LIMELM; i++ )
58  {
59  for( long j=0; j < LIMELM; j++ )
60  {
61  if( safe_div( thermal.heating(i,j), heat_total, 0. ) > SMALLFLOAT )
62  {
63  ipsave[ipnt] = i;
64  jpsave[ipnt] = j;
65  SaveVal[ipnt] = (realnum)( safe_div( thermal.heating(i,j), heat_total, 0. ) );
66  ipnt++;
67  }
68  }
69  }
70 
71  /* now check for possible line heating not counted in 1,23
72  * there should not be any significant heating source here
73  * since they would not be counted in derivative correctly */
74  for( long i=0; i < thermal.ncltot; i++ )
75  {
76  if( safe_div( thermal.heatnt[i], heat_total, 0. ) > save.WeakHeatCool )
77  {
78  realnum awl;
79  awl = thermal.collam[i];
80  /* force to save wavelength convention as printout */
81  if( awl > 100000 )
82  awl /= 10000;
83  fprintf( io, " Negative coolant was %s %.2f %.2e\n",
84  thermal.chClntLab[i], awl, safe_div( thermal.heatnt[i], heat_total, 0. ) );
85  }
86  }
87 
88  if( !conv.lgConvTemp )
89  {
90  fprintf( io, "#>>>> Temperature not converged.\n" );
91  }
92  else if( !conv.lgConvEden )
93  {
94  fprintf( io, "#>>>> Electron density not converged.\n" );
95  }
96  else if( !conv.lgConvIoniz() )
97  {
98  fprintf( io, "#>>>> Ionization not converged.\n" );
99  }
100  else if( !conv.lgConvPres )
101  {
102  fprintf( io, "#>>>> Pressure not converged.\n" );
103  }
104 
105  /* this is mainly to keep the compiler from flagging possible paths
106  * with j not being set */
107  long i = INT_MIN;
108  long j = INT_MIN;
109  /* following loop tries to identify strongest agents and turn to labels */
110  for( long k=0; k < ipnt; k++ )
111  {
112  /* generate labels that make sense in printout
113  * if not identified with a specific agent, will print indices as [i][j],
114  * heating is thermal.heating(i,j) */
115  i = ipsave[k];
116  j = jpsave[k];
117  /* i >= j indicates agent is one of the first LIMELM elements */
118  if( i >= j )
119  {
120  if( dense.xIonDense[i][j] == 0. && thermal.heating(i,j)>SMALLFLOAT )
121  fprintf(ioQQQ,"DISASTER assert about to be thrown - search for hit it\n");
122  /*fprintf(ioQQQ,"DEBUG %li %li %.2e %.2e\n", i , j ,
123  dense.xIonDense[i][j],
124  thermal.heating(i,j));*/
125  ASSERT( dense.xIonDense[i][j] > 0. || thermal.heating(i,j)<SMALLFLOAT );
126  /* this is case of photoionization of atom or ion */
127  chLabel[k] = elementnames.chElementSym[i];
128  chLabel[k] += elementnames.chIonStage[j];
129  }
130  /* notice that in test i and j are swapped from order in heating array */
131  else if( i == 0 && j == 1 )
132  {
133  /* photoionization from all excited states of Hydrogenic species,
134  * heating(0,1) */
135  chLabel[k] = "Hn=2";
136  }
137  else if( i == 0 && j == 3 )
138  {
139  /* collisional ionization of all iso-seq from all levels,
140  * heating(0,3) */
141  chLabel[k] = "Hion";
142  }
143  else if( i == 0 && j == 7 )
144  {
145  /* UTA ionization heating(0,7) */
146  chLabel[k] = " UTA";
147  }
148  else if( i == 0 && j == 8 )
149  {
150  /* thermal.heating(0,8) is heating due to collisions within
151  * X of H2, code var hmi.HeatH2Dexc_used */
152  chLabel[k] = "H2vH";
153  }
154  else if( i == 0 && j == 17 )
155  {
156  /* thermal.heating(0,17) is heating due to photodissociation
157  * heating of X within H2,
158  * code var hmi.HeatH2Dish_used */
159  chLabel[k] = "H2dH";
160  }
161  else if( i == 0 && j == 9 )
162  {
163  /* CO dissociation, co.CODissHeat, heating(0,9) */
164  chLabel[k] = "COds";
165  }
166  else if( i == 0 && j == 20 )
167  {
168  /* extra heat thermal.heating(0,20)*/
169  chLabel[k] = "extH";
170  }
171  else if( i == 0 && j == 21 )
172  {
173  /* pair heating thermal.heating(0,21)*/
174  chLabel[k] = "pair";
175  }
176  else if( i == 0 && j == 11 )
177  {
178  /* free free heating, heating(0,11) */
179  chLabel[k] = "H FF";
180  }
181  else if( i == 0 && j == 12 )
182  {
183  /* heating coolant (not line), physical cooling process, often a bug, heating(0,12) */
184  chLabel[k] = "Hcol";
185  }
186  else if( i == 0 && j == 13 )
187  {
188  /* grain photoelectric effect, heating(0,13) */
189  chLabel[k] = "GrnP";
190  }
191  else if( i == 0 && j == 14 )
192  {
193  /* grain collisions, heating(0,14) */
194  chLabel[k] = "GrnC";
195  }
196  else if( i == 0 && j == 15 )
197  {
198  /* H- heating, heating(0,15) */
199  chLabel[k] = "H- ";
200  }
201  else if( i == 0 && j == 16 )
202  {
203  /* H2+ heating, heating(0,16) */
204  chLabel[k] = "H2+ ";
205  }
206  else if( i == 0 && j == 18 )
207  {
208  /* H2 photoionization heating, heating(0,18) */
209  chLabel[k] = "H2ph";
210  }
211  else if( i == 0 && j == 19 )
212  {
213  /* Compton heating, heating(0,19) */
214  chLabel[k] = "Comp";
215  }
216  else if( i == 0 && j == 22 )
217  {
218  /* line heating, heating(0,22) */
219  chLabel[k] = "line";
220  }
221  else if( i == 0 && j == 23 )
222  {
223  /* iso-sequence line heating - all elements together,
224  * heating(0,23) */
225  chLabel[k] = "Hlin";
226  }
227  else if( i == 0 && j == 24 )
228  {
229  /* charge transfer heating, heating(0,24) */
230  chLabel[k] = "ChaT";
231  }
232  else if( i == 1 && j == 3 )
233  {
234  /* helium triplet line heating, heating(1,3) */
235  chLabel[k] = "He3l";
236  }
237  else if( i == 1 && j == 5 )
238  {
239  /* advective heating, heating(1,5) */
240  chLabel[k] = "adve";
241  }
242  else if( i == 1 && j == 6 )
243  {
244  /* cosmic ray heating thermal.heating(1,6)*/
245  chLabel[k] = "CR H";
246  }
247  else if( i == 25 && j == 27 )
248  {
249  /* Fe 2 line heating, heating(25,27) */
250  chLabel[k] = "Fe 2";
251  }
252  else
253  {
254  ostringstream oss;
255  oss << "[" << i << "][" << j << "]";
256  chLabel[k] = oss.str();
257  }
258  }
259 
260  /* now sort by decreasing importance */
261  /*spsort netlib routine to sort array returning sorted indices */
262  UNUSED int nFail;
263  spsort(/* input array to be sorted */
264  &SaveVal[0],
265  /* number of values in x */
266  ipnt,
267  /* permutation output array */
268  &ipOrdered[0],
269  /* flag saying what to do - 1 sorts into increasing order, not changing
270  * the original routine */
271  -1,
272  /* error condition, should be 0 */
273  &nFail);
274 
275  /*>>chng 06 jun 06, change start of save to give same info as cooling
276  * as per comment by Yumihiko Tsuzuki */
277  /* begin the print out with zone number, total heating and cooling */
278  fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e",
280  phycon.te,
281  heat_total,
282  cool_total );
283 
284  for( long k=0; k < ipnt; k++ )
285  {
286  int ip = ipOrdered[k];
287  i = ipsave[ip];
288  j = jpsave[ip];
289  ASSERT( i<LIMELM && j<LIMELM );
290  if(k > 4 && safe_div( thermal.heating(i,j), heat_total, 0. ) < save.WeakHeatCool )
291  break;
292  fprintf( io, "\t%s\t%.7f ",
293  chLabel[ip].c_str(), SaveVal[ip] );
294  }
295  fprintf( io, " \n" );
296 
297  /* a negative pointer in the heating array is probably a problem,
298  * indicating that some line has become a heat source */
299  bool lgHeatLine = false;
300 
301  /* check if any lines were major heat sources */
302  for( i=0; i < ipnt; i++ )
303  {
304  /* heating(22,0) is line heating - identify line if important */
305  if( ipsave[ipOrdered[i]] == 0 && jpsave[ipOrdered[i]] == 22 )
306  lgHeatLine = true;
307  }
308 
309  if( lgHeatLine )
310  {
311  long level = -1;
312  /* a line was a major heat source - identify it */
313  TransitionProxy t = FndLineHt(&level);
314  if( safe_div( t.Coll().heat(), heat_total, 0. ) > 0.005 )
315  {
316  ASSERT( t.associated() );
317  double TauIn = t.Emis().TauIn();
318  double Pump = t.Emis().pump();
319  double EscP = t.Emis().Pesc();
320  double CS = t.Coll().col_str();
321  /* ratio of line to total heating */
322  double ColHeat = safe_div( t.Coll().heat(), heat_total, 0. );
323 
324  fprintf( io, " LHeat lv%2ld %s TIn%10.2e Pmp%9.1e EscP%9.1e CS%9.1e Hlin/tot%10.2e\n",
325  level, chLineLbl(t).c_str(), TauIn, Pump, EscP, CS, ColHeat );
326  }
327  }
328  return;
329 }
bool lgConvEden
Definition: conv.h:195
double htot
Definition: thermal.h:169
t_thermal thermal
Definition: thermal.cpp:6
void SaveHeat(FILE *io)
Definition: heat_save.cpp:18
double Cool()
Definition: dynamics.cpp:2207
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
char chIonStage[LIMELM+1][CHARS_ION_STAGE]
Definition: elementnames.h:29
const realnum SMALLFLOAT
Definition: cpu.h:246
bool lgConvPres
Definition: conv.h:192
t_conv conv
Definition: conv.cpp:5
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
bool lgConvIoniz() const
Definition: conv.h:108
t_dynamics dynamics
Definition: dynamics.cpp:42
t_dense dense
Definition: global.cpp:15
bool associated() const
Definition: transition.h:54
t_elementnames elementnames
Definition: elementnames.cpp:5
double Heat()
Definition: dynamics.cpp:2193
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
realnum WeakHeatCool
Definition: save.h:488
double & heat() const
Definition: collision.h:224
char chClntLab[NCOLNT][NCOLNT_LAB_LEN+1]
Definition: thermal.h:108
EmissionList::reference Emis() const
Definition: transition.h:447
const TransitionProxy FndLineHt(long int *level)
float realnum
Definition: cddefines.h:124
realnum & Pesc() const
Definition: emission.h:580
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
double depth_mid_zone
Definition: radius.h:31
t_radius radius
Definition: radius.cpp:5
double heating(long nelem, long ion)
Definition: thermal.h:186
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
#define ASSERT(exp)
Definition: cddefines.h:613
long int ncltot
Definition: thermal.h:106
realnum collam[NCOLNT]
Definition: thermal.h:103
bool lgConvTemp
Definition: conv.h:189
const int LIMELM
Definition: cddefines.h:308
realnum & col_str() const
Definition: collision.h:191
CollisionProxy Coll() const
Definition: transition.h:463
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
double heatnt[NCOLNT]
Definition: thermal.h:104
realnum & TauIn() const
Definition: emission.h:470
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1222
double ctot
Definition: thermal.h:130
#define UNUSED
Definition: cpu.h:14
double & pump() const
Definition: emission.h:530