/[ascend]/trunk/models/johnpye/fprops/solve_ph.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/solve_ph.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2680 - (show annotations) (download) (as text)
Mon Jan 28 06:30:25 2013 UTC (9 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 10136 byte(s)
Working on problem with solve_ph. Could be that one of the deriv routines is wrong in the saturation region?
1 /*
2 ASCEND modelling environment
3 Copyright (C) 2004-2010 John Pye
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include "solve_ph.h"
20 #include "fprops.h"
21 #include "sat.h"
22 #include "derivs.h"
23 #include "rundata.h"
24
25 #include <stdio.h>
26 #include <math.h>
27 #include <stdlib.h>
28
29 #define SOLVE_PH_DEBUG
30 #define SOLVE_PH_ERRORS
31
32 //#define FPE_DEBUG
33 //#define ASSERT_DEBUG
34
35 #ifdef ASSERT_DEBUG
36 # include <assert.h>
37 #else
38 # define assert(ARGS...)
39 #endif
40
41 #include <setjmp.h>
42 #include <signal.h>
43
44 #ifdef FPE_DEBUG
45 #define _GNU_SOURCE
46 #include <fenv.h>
47 int feenableexcept(int excepts);
48 int fedisableexcept(int excepts);
49 int fegetexcept(void);
50 #endif
51
52 #ifdef SOLVE_PH_DEBUG
53 # include "color.h"
54 # define MSG(FMT, ...) \
55 color_on(stderr,ASC_FG_BRIGHTRED);\
56 fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
57 color_on(stderr,ASC_FG_BRIGHTBLUE);\
58 fprintf(stderr,"%s: ",__func__);\
59 color_off(stderr);\
60 fprintf(stderr,FMT "\n",##__VA_ARGS__)
61 #else
62 # define MSG(ARGS...) ((void)0)
63 #endif
64
65 #ifdef SOLVE_PH_ERRORS
66 # include "color.h"
67 # define ERRMSG(STR,...) \
68 color_on(stderr,ASC_FG_BRIGHTRED);\
69 fprintf(stderr,"ERROR:");\
70 color_off(stderr);\
71 fprintf(stderr," %s:%d:" STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
72 #else
73 # define ERRMSG(ARGS...) ((void)0)
74 #endif
75
76 int fprops_region_ph(double p, double h, const PureFluid *fluid, FpropsError *err){
77 double Tsat, rhof, rhog;
78 double p_c = fluid->data->p_c;
79
80 if(p >= p_c)return FPROPS_NON;
81
82 fprops_sat_p(p, &Tsat, &rhof, &rhog, fluid, err);
83 if(*err){
84 *err = FPROPS_SAT_CVGC_ERROR;
85 return FPROPS_ERROR;
86 }
87
88 double hf = fluid->h_fn(Tsat, rhof, fluid->data, err);
89 if(h <= hf)return FPROPS_NON;
90
91 double hg = fluid->h_fn(Tsat,rhog, fluid->data, err);
92 if(h >= hg)return FPROPS_NON;
93
94 return FPROPS_SAT;
95 }
96
97 typedef void SignalHandler(int);
98
99 #ifdef FPE_DEBUG
100 jmp_buf mark;
101 void fprops_fpe(int sig){
102 MSG("Catching signal %d!",sig);
103 feclearexcept(FE_DIVBYZERO|FE_INVALID|FE_OVERFLOW);
104 longjmp(mark, -1);
105 }
106 #endif
107
108 void fprops_solve_ph(double p, double h, double *T, double *rho, int use_guess
109 , const PureFluid *fluid, FpropsError *err
110 ){
111 double Tsat, rhof, rhog, hf, hg;
112 double T1, rho1;
113 int subcrit_pressure = 0;
114 int liquid_iteration = 0;
115 double rhof_t;
116 double p_c = fluid->data->p_c;
117
118 MSG("Solving for p=%f bar, h=%f kJ/kgK (EOS type %d, '%s')",p/1e5,h/1e3,fluid->type,fluid->name);
119
120 #ifdef FPE_DEBUG
121 feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
122 SignalHandler *old = signal(SIGFPE,&fprops_fpe);
123 (void)old;/* not used for anything at this stage */
124 int jmpret = setjmp(mark);
125 if(jmpret==0){
126 #endif
127 MSG("p_c = %f bar",p_c/1e5);
128 if(p < p_c){
129 /* TODO what about testing for p >= p_t? */
130 MSG("Calculate saturation Tsat(p < p_c) with p = %f",p);
131 fprops_sat_p(p, &Tsat, &rhof, &rhog, fluid, err);
132 if(*err){
133 ERRMSG("Unable to solve saturation state");
134 *err = FPROPS_SAT_CVGC_ERROR;
135 return;
136 }
137 hf = fluid->h_fn(Tsat, rhof, fluid->data, err);
138 hg = fluid->h_fn(Tsat, rhog, fluid->data, err);
139 MSG("at p = %f bar, T_sat = %f, rhof = %f, hf = %f kJ/kg, hg = %f",p/1e5,Tsat,rhof, hf/1e3,hg/1e3);
140
141 if(hf <= h && h <= hg){
142 MSG("SATURATION REGION");
143 /* saturation region... easy */
144 double x = (h - hf)/(hg - hf);
145 *rho = 1./(x/rhog + (1.-x)/rhof);
146 *T = Tsat;
147 return;
148 }
149
150 subcrit_pressure = 1;
151 if(h < hf){
152 liquid_iteration = 1;
153 if(!use_guess){
154 *T = Tsat;
155 *rho = rhof;
156 MSG("h < hf; LIQUID GUESS: T = %f, rho = %f",*T, *rho);
157 }
158 }else if(!use_guess){
159 *T = 1.1 * Tsat;
160 *rho = rhog * 0.5;
161 MSG("GAS GUESS: T = %f, rho = %f",*T, *rho);
162 }
163 }else{ /* p >= p_c */
164 /* FIXME: still some problems here at very high pressures */
165 if(!use_guess){
166 /* FIXME we should cache/precalculate hc, store in rundata. */
167 double hc = fluid->h_fn(fluid->data->T_c, fluid->data->rho_c, fluid->data, err);
168 assert(!isnan(hc));
169 MSG("hc = %f kJ/kgK",hc/1e3);
170 if(h < 0.8*hc){ /* FIXME should make use of h_ref here in some way? Otherwise arbitrary... */
171 #if 0
172 double Tsat1, psat1, rhof1, rhog1;
173 int res = fprops_sat_p(p, &Tsat1, &rhof1, &rhog1, fluid);
174 MSG("Unable to solve saturation at p = %f",p);
175 *T = Tsat1;
176 *rho = rhof1;
177 #else
178 MSG("h < 0.9 hc... using saturation Tsat(hf) for starting guess");
179 double Tsat1, psat1, rhof1, rhog1;
180 fprops_sat_hf(h, &Tsat1, &psat1, &rhof1, &rhog1, fluid, err);
181 if(*err){
182 MSG("Unable to solve Tsat(hf)");
183 /* accuracy of the estimate of Tsat(hf) doesn't matter
184 very much so we can ignore this error. */
185
186 /* try initialising to T_t, rho_f(T_t) */
187 Tsat1 = fluid->data->T_t;
188 MSG("T_t = %f",Tsat1);
189 fprops_sat_T(Tsat1, &psat1, &rhof1, &rhog1, fluid, err);
190 if(*err)MSG("Unable to solve rhof(Tt)");
191 }
192 *T = Tsat1;
193 *rho = rhof1;
194 #endif
195 }else{
196 *T = fluid->data->T_c * 1.01;
197 *rho = fluid->data->rho_c * 1.05;
198 }
199 MSG("SUPERCRITICAL GUESS: T = %f, rho = %f", *T, *rho);
200 }
201 }
202
203 MSG("STARTING NON-SAT ITERATION");
204 //*rho = 976.82687191126922;
205 //*T = 344.80371310850518;
206
207 T1 = *T;
208 rho1 = *rho;
209 assert(!isnan(T1));
210 assert(!isnan(rho1));
211
212 if(liquid_iteration){
213 double pt,rhogt;
214 fprops_triple_point(&pt, &rhof_t, &rhogt, fluid, err);
215 if(*err){
216 ERRMSG("Unable to solve triple point liquid density.");
217 *err = FPROPS_SAT_CVGC_ERROR;
218 return;
219 }
220 }
221
222 /* try our own home-baked newton iteration */
223 int i = 0;
224 *err = FPROPS_NO_ERROR;
225 double delta_T = 0;
226 double delta_rho = 0;
227 MSG("STARTING ITERATION");
228 MSG("rhof_t = %f",rhof_t);
229 while(i++ < 200){
230 double p1 = fluid->p_fn(T1,rho1, fluid->data, err);
231 if(*err){
232 MSG("Got an error ('%s') in p_fn calculation",fprops_error(*err));
233 }
234 double h1 = fluid->h_fn(T1,rho1, fluid->data, err);
235 assert(!isnan(h1));
236
237 if(i >= 2){
238 int nred = 10;
239 while(p1 <= 0 && nred){
240 rho1 = rho1 - delta_rho;
241 T1 = T1 - delta_T;
242 delta_rho *= 0.5;
243 delta_T *= 0.5;
244 rho1 = rho1 + delta_rho;
245 T1 = T1 + delta_T;
246 p1 = fluid->p_fn(T1,rho1,fluid->data, err);
247 MSG("Set smaller step as p < 0. T1 = %f, rho1 = %f --> p1 = %f",T1, rho1, p1);
248 nred--;
249 }
250 }
251
252 MSG(" %d: T = %f, rho = %f\tp = %f bar, h = %f kJ/kg", i, T1, rho1, p1/1e5, h1/1e3);
253 //MSG(" p error = %f bar",(p1 - p)/1e5);
254 //MSG(" h error = %f kJ/kg",(h1 - h)/1e3);
255
256 if(p1 < 0){
257 MSG("p1 < 0, reducing dT, drho");
258 T1 -= (delta_T *= 0.5);
259 rho1 -= (delta_rho *= 0.5);
260 continue;
261 }
262
263 if(fabs(p1 - p) < 1e-4 && fabs(h1 - h) < 1e-8){
264 MSG("Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);
265 *T = T1;
266 *rho = rho1;
267 return;
268 }
269 /* calculate step, we're solving log(p1) in this code... */
270 double f = log(p1) - log(p);
271 double g = h1 - h;
272 assert(!isnan(f));
273 assert(!isnan(g));
274 assert(rho1 != 0);
275 assert(p1 != 0);
276 assert(!isnan(fprops_non_dZdT_v('p', T1,rho1, fluid,err)));
277 double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, fluid,err);
278 double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, fluid,err);
279 double g_T = fprops_non_dZdT_v('h', T1,rho1, fluid,err);
280 double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, fluid,err);
281 assert(!isnan(f_T));
282
283 if(isnan(f_rho)){
284 ERRMSG(" rho1 = %f, T1 = %f",rho1, T1);
285 }
286 assert(!isnan(f_rho));
287
288 assert(!isnan(g_T));
289 assert(!isnan(g_rho));
290
291 double det = g_rho * f_T - f_rho * g_T;
292 assert(det!=0);
293 assert(!isnan(det));
294 MSG(" df/dT = %e\t\tdf/drho = %e",f_T, f_rho);
295 MSG(" dg/dT = %e\t\tdg/drho = %e",g_T, g_rho);
296
297 delta_T = -1./det * (g_rho * f - f_rho * g);
298 delta_rho = -1./det * (f_T * g - g_T * f);
299 assert(!isnan(delta_T));
300 assert(!isnan(delta_rho));
301 MSG(" dT = %f", delta_T);
302 MSG(" drho = %f", delta_rho);
303
304 if(subcrit_pressure){
305 if(h > hg){
306 /* vapour */
307 int i = 0;
308 while(rho1 + delta_rho > rhog && i++ < 20){
309 delta_rho *= 0.5;
310 }
311 }else{
312 /* liquid */
313 while(rho1 + delta_rho < rhof && i++ < 20){
314 delta_rho *= 0.5;
315 }
316 if(T1 + delta_T < fluid->data->T_t) delta_T = 0.5 * (T1 + fluid->data->T_t);
317 }
318 }else{
319 #if 0
320 /* supercritical... stay above critical temperature of density > rho_crit */
321 if(rho1 + delta_rho > fluid->data->rho_c && T1 + delta_T < fluid->data->T_c){
322 delta_T = 0.5 *(T1 + fluid->data->T_c);
323 }
324 #endif
325 }
326 /* don't go too dense */
327 #if 1
328 if(liquid_iteration){
329 MSG("rho1 = %f, delta_rho = %f, rho1+delta_rho = %f", rho1, delta_rho, rho1+delta_rho);
330 if(rho1 + delta_rho > rhof_t){
331 MSG("Limit rho to be less than rhof_t");
332 delta_rho = rhof_t - rho1;
333 }
334 }
335 #endif
336
337 /* don't go too hot */
338 if(T1 + delta_T > 5000) delta_T = 5000 - T1;
339
340 /* avoid huge step */
341 while(fabs(delta_T / T1) > 0.7){
342 MSG("Reduce delta_T");
343 delta_T *= 0.5;
344 }
345
346 /* NOTE if step limit is less than 0.5, we're getting errors */
347 while(fabs(delta_rho / rho1) > 0.7){
348 MSG("Reduce delta_rho");
349 delta_rho *= 0.5;
350 }
351
352 T1 = T1 + delta_T;
353 if(T1 < fluid->data->T_t)T1 = fluid->data->T_t;
354 rho1 = rho1 + delta_rho;
355 }
356 #ifdef FPE_DEBUG
357 }else{
358 /* an FPE occurred */
359 MSG("An FPE occurred");
360 abort();
361 }
362 #endif
363
364 *T = T1;
365 *rho = rho1;
366 ERRMSG("Iteration failed for '%s' with p = %.12e, h = %.12e",fluid->name, p,h);
367 *err = FPROPS_NUMERIC_ERROR;
368 return;
369
370 #if 0
371 int res = fprops_nonsolver('p','h',p,h,T,rho,fluid);
372 ERRMSG("Iteration failed in nonsolver");
373 return res;
374 #endif
375 }
376

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22