cloudy trunk
Loading...
Searching...
No Matches
cool_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/*CoolSave save coolants */
4#include "cddefines.h"
5#include "thermal.h"
6#include "dynamics.h"
7#include "radius.h"
8#include "conv.h"
9#include "phycon.h"
10#include "save.h"
11#include "grainvar.h"
12#include "hmi.h"
13#include "coolheavy.h"
14
15
16/* this is limit to number of coolants to print out */
17static const int IPRINT = 100;
18
19/*CoolSave save coolants */
20void CoolSave(FILE * io, char chJob[])
21{
22 long int i,
23 ip,
24 is;
25
26 int nFail;
27
28 double cset,
29 cool_total,
30 heat_total;
31
33 *csav,
34 *sgnsav;
35 long int *index;
36
37 DEBUG_ENTRY( "CoolSave()" );
38
39 /* cannot do one-time init since thermal.ncltot can change */
40 index = (long int *)CALLOC((size_t)thermal.ncltot,sizeof(long int));
41 csav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum));
42 sgnsav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum));
43
44 cool_total = thermal.ctot;
45 heat_total = thermal.htot;
46
47 /* >>chng 06 mar 17, comment out following block and replace with this
48 * removing dynamics heating & cooling and report only physical
49 * heating and cooling
50 * NB the heating and cooling as punched no longer need be
51 * equal for a converged model */
52 cool_total -= dynamics.Cool();
53 heat_total -= dynamics.Heat();
54# if 0
55 if(dynamics.Cool > dynamics.Heat())
56 {
57 cool_total -= dynamics.Heat();
58 heat_total -= dynamics.Heat();
59 }
60 else
61 {
62 cool_total -= dynamics.Cool;
63 heat_total -= dynamics.Cool;
64 }
65# endif
66
67 /* cset will be weakest cooling to consider
68 * WeakHeatCool set with 'set weakheatcool' command
69 * default is 0.05 */
70 cset = cool_total*save.WeakHeatCool;
71
72 /* first find all strong lines, both + and - sign */
73 ip = thermal.ncltot;
74
75 for( i=0; i < ip; i++ )
76 {
77 csav[i] = (realnum)(MAX2(thermal.cooling[i],thermal.heatnt[i])/
78 SDIV(cool_total));
79
80 /* save sign to remember if heating or cooling line */
81 if( thermal.heatnt[i] == 0. )
82 {
83 sgnsav[i] = 1.;
84 }
85 else
86 {
87 sgnsav[i] = -1.;
88 }
89 }
90
91 /* order strongest to weakest */
92 /* now sort by decreasing importance */
93 /*spsort netlib routine to sort array returning sorted indices */
94 spsort(
95 /* input array to be sorted */
96 csav,
97 /* number of values in x */
98 ip,
99 /* permutation output array */
100 index,
101 /* flag saying what to do - 1 sorts into increasing order, not changing
102 * the original routine */
103 -1,
104 /* error condition, should be 0 */
105 &nFail);
106
107 /* warn if tcovergence failure occurred */
108 if( !conv.lgConvTemp )
109 {
110 fprintf( io, "#>>>> Temperature not converged.\n" );
111 }
112 else if( !conv.lgConvEden )
113 {
114 fprintf( io, "#>>>> Electron density not converged.\n" );
115 }
116 else if( !conv.lgConvIoniz() )
117 {
118 fprintf( io, "#>>>> Ionization not converged.\n" );
119 }
120 else if( !conv.lgConvPres )
121 {
122 fprintf( io, "#>>>> Pressure not converged.\n" );
123 }
124
125 if( strcmp(chJob,"EACH") == 0 )
126 {
127 /* begin the print out with zone number, total heating and cooling */
128 fprintf( io, "%.5e\t%.4e\t%.4e",
129 radius.depth_mid_zone,
130 phycon.te,
131 cool_total );
132 double debug_ctot = 0.;
133
134 for( int i=0 ; i <= LIMELM ; i++)
135 {
136 fprintf( io, "\t%.4e", thermal.elementcool[i] );
137 debug_ctot += thermal.elementcool[i];
138 }
139
140 /*fprintf( io, "\t%.4e", thermal.dust );
141 fprintf( io, "\t%.4e", thermal.H2cX );
142 fprintf( io, "\t%.4e", thermal.CT_C );
143 fprintf( io, "\t%.4e", thermal.H_fb );
144 fprintf( io, "\t%.4e", thermal.H2ln );
145 fprintf( io, "\t%.4e", thermal.HDro );
146 fprintf( io, "\t%.4e", thermal.H2p );
147 fprintf( io, "\t%.4e", thermal.FF_c );
148 fprintf( io, "\t%.4e", thermal.hvFB );
149 fprintf( io, "\t%.4e", thermal.eeff );
150 fprintf( io, "\t%.4e", thermal.adve );
151 fprintf( io, "\t%.4e", thermal.Comp );
152 fprintf( io, "\t%.4e", thermal.Extr );
153 fprintf( io, "\t%.4e", thermal.Expn );
154 fprintf( io, "\t%.4e", thermal.Cycl );
155 fprintf( io, "\t%.4e", thermal.Hvin );
156 fprintf( io, " \n" );*/
157
158 fprintf( io, "\t%.4e", MAX2(0.,gv.GasCoolColl) );
159 debug_ctot += MAX2(0.,gv.GasCoolColl);
160
161 fprintf( io, "\t%.4e", MAX2(0.,-hmi.HeatH2Dexc_used) );
162 debug_ctot += MAX2(0.,-hmi.HeatH2Dexc_used);
163
164 fprintf( io, "\t%.4e", thermal.char_tran_cool );
165 debug_ctot += thermal.char_tran_cool;
166
167 fprintf( io, "\t%.4e", hmi.hmicol );
168 debug_ctot += hmi.hmicol;
169
170 fprintf( io, "\t%.4e", CoolHeavy.h2line );
171 debug_ctot += CoolHeavy.h2line;
172
173 fprintf( io, "\t%.4e", CoolHeavy.HD );
174 debug_ctot += CoolHeavy.HD;
175
176 fprintf( io, "\t%.4e", CoolHeavy.H2PlsCool );
177 debug_ctot += CoolHeavy.H2PlsCool;
178
179 fprintf( io, "\t%.4e", MAX2(0.,CoolHeavy.brems_cool_net) );
180 debug_ctot += MAX2(0.,CoolHeavy.brems_cool_net);
181
182 fprintf( io, "\t%.4e", CoolHeavy.heavfb );
183 debug_ctot += CoolHeavy.heavfb;
184
185 fprintf( io, "\t%.4e", CoolHeavy.eebrm );
186 debug_ctot += CoolHeavy.eebrm;
187
188 fprintf( io, "\t%.4e", CoolHeavy.tccool );
189 debug_ctot += CoolHeavy.tccool;
190
191 fprintf( io, "\t%.4e", CoolHeavy.cextxx );
192 debug_ctot += CoolHeavy.cextxx;
193
194 fprintf( io, "\t%.4e", CoolHeavy.expans );
195 debug_ctot += CoolHeavy.expans;
196
197 fprintf( io, "\t%.4e", CoolHeavy.cyntrn );
198 debug_ctot += CoolHeavy.cyntrn;
199
200 fprintf( io, "\t%.4e", CoolHeavy.colmet );
201 debug_ctot += CoolHeavy.colmet;
202
203 fprintf( io, "\t%.4e", thermal.dima );
204 debug_ctot += thermal.dima;
205 fprintf( io, " \n" );
206 if( fabs( (debug_ctot - cool_total)/cool_total ) > 1e-10 )
207 {
208 fprintf( ioQQQ , "PROBLEM with the SAVE EACH COOLING output\n" );
209 fprintf( ioQQQ , "PROBLEM One or more coolants have been lost, the sum of the reported cooling is %.4e\n", debug_ctot );
210 fprintf( ioQQQ , "PROBLEM The total cooling is %.4ee\n", cool_total );
211 fprintf( ioQQQ , "PROBLEM The difference is %.4e\n", cool_total - debug_ctot );
213 }
214 }
215 else if( strcmp(chJob,"COOL") == 0 )
216 {
217 /*>>chng 06 jun 06, change start of save to give same info as heating
218 * as per comment by Yumihiko Tsuzuki */
219 /* begin the print out with zone number, total heating and cooling */
220 fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e",
221 radius.depth_mid_zone,
222 phycon.te,
223 heat_total,
224 cool_total );
225
226 /* print only up to IPRINT, which is defined above */
227 ip = MIN2( ip , IPRINT );
228
229 /* now print the coolants
230 * keep sign of coolant, for strong negative cooling
231 * order is ion, wavelength, fraction of total */
232 for( is=0; is < ip; is++ )
233 {
234 if(is > 4 && (thermal.cooling[index[is]] < cset && thermal.heatnt[index[is]] < cset))
235 break;
236 fprintf( io, "\t%s %.1f\t%.7f",
237 thermal.chClntLab[index[is]],
238 thermal.collam[index[is]],
239 sign(csav[index[is]],sgnsav[index[is]]) );
240 }
241 fprintf( io, " \n" );
242 }
243 else
245
246 /* finished, now free space */
247 free(sgnsav);
248 free(csav);
249 free(index);
250
251 return;
252}
253
FILE * ioQQQ
Definition cddefines.cpp:7
#define MIN2
Definition cddefines.h:761
#define CALLOC
Definition cddefines.h:510
const int LIMELM
Definition cddefines.h:258
T sign(T x, T y)
Definition cddefines.h:800
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition service.cpp:1100
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_conv conv
Definition conv.cpp:5
void CoolSave(FILE *io, char chJob[])
Definition cool_save.cpp:20
static const int IPRINT
Definition cool_save.cpp:17
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
t_dynamics dynamics
Definition dynamics.cpp:44
GrainVar gv
Definition grainvar.cpp:5
t_hmi hmi
Definition hmi.cpp:5
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