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