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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2660 - (show annotations) (download) (as text)
Wed Jan 16 05:57:03 2013 UTC (9 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 17098 byte(s)
Working on updating rankine_fprops and associated models to work with new fprops2 code.
Some issue discovered with (p,h) for water (added python/solve_ph1.py to check it).
Next cunit tests to drive some ASCEND models embedding FPROPS (note use of slvreq for this).
1 /* ASCEND modelling environment
2 Copyright (C) 2008 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//** @file
19 Wrapper for FPROPS to allow access from ASCEND.
20 */
21
22 //#define ASC_FPROPS_DEBUG
23
24 /* include the external function API from libascend... */
25 #include <ascend/compiler/extfunc.h>
26
27 /* include error reporting API as well, so we can send messages to user */
28 #include <ascend/utilities/error.h>
29
30 /* for accessing the DATA instance */
31 #include <ascend/compiler/child.h>
32 #include <ascend/general/list.h>
33 #include <ascend/compiler/module.h>
34 #include <ascend/compiler/childinfo.h>
35 #include <ascend/compiler/parentchild.h>
36 #include <ascend/compiler/slist.h>
37 #include <ascend/compiler/type_desc.h>
38 #include <ascend/compiler/packages.h>
39 #include <ascend/compiler/symtab.h>
40 #include <ascend/compiler/instquery.h>
41 #include <ascend/compiler/instmacro.h>
42 #include <ascend/compiler/instance_types.h>
43
44 /* the code that we're wrapping... */
45 #include "fprops.h"
46 #include "sat.h"
47 #include "solve_ph.h"
48
49 /* for the moment, species data are defined in C code, we'll implement something
50 better later on, hopefully. */
51 #include "fluids.h"
52
53 #ifndef ASC_EXPORT
54 # error "Where is ASC_EXPORT?"
55 #endif
56
57
58 /*------------------------------------------------------------------------------
59 FORWARD DECLARATIONS
60 */
61
62 ExtBBoxInitFunc asc_fprops_prepare;
63 ExtBBoxFunc fprops_p_calc;
64 ExtBBoxFunc fprops_u_calc;
65 ExtBBoxFunc fprops_s_calc;
66 ExtBBoxFunc fprops_h_calc;
67 ExtBBoxFunc fprops_a_calc;
68 ExtBBoxFunc fprops_g_calc;
69 ExtBBoxFunc fprops_cp_calc;
70 ExtBBoxFunc fprops_cv_calc;
71 ExtBBoxFunc fprops_w_calc;
72 ExtBBoxFunc fprops_phsx_vT_calc;
73 ExtBBoxFunc fprops_Tvsx_ph_calc;
74
75 #define TCRIT(FLUID) (FLUID->data->T_c)
76 #define TTRIP(FLUID) (FLUID->data->T_t)
77 #define RHOCRIT(FLUID) (FLUID->data->rho_c)
78 #define PCRIT(FLUID) (FLUID->data->p_c)
79
80 /*------------------------------------------------------------------------------
81 GLOBALS
82 */
83
84 /* place to store symbols needed for accessing ASCEND's instance tree */
85 static symchar *fprops_symbols[2];
86 #define COMPONENT_SYM fprops_symbols[0]
87 #define TYPE_SYM fprops_symbols[1]
88
89 static const char *fprops_p_help = "Calculate pressure from temperature and density, using FPROPS";
90 static const char *fprops_u_help = "Calculate specific internal energy from temperature and density, using FPROPS";
91 static const char *fprops_s_help = "Calculate specific entropy from temperature and density, using FPROPS";
92 static const char *fprops_h_help = "Calculate specific enthalpy from temperature and density, using FPROPS";
93 static const char *fprops_a_help = "Calculate specific Helmholtz energy from temperature and density, using FPROPS";
94 static const char *fprops_g_help = "Calculate specific Gibbs energy from temperature and density, using FPROPS";
95 static const char *fprops_cp_help = "Calculate isobaric specific heat from temperature and density, using FPROPS";
96 static const char *fprops_cv_help = "Calculate isochoric specific heat from temperature and density, using FPROPS";
97 static const char *fprops_w_help = "Calculate speed of sound from temperature and density, using FPROPS";
98
99 static const char *fprops_phsx_vT_help = "Calculate p, h, s, x from temperature and density, using FPROPS/Helmholtz eqn";
100
101 static const char *fprops_Tvsx_ph_help = "Calculate T, v, s, x from pressure and enthalpy, using FPROPS/Helmholtz eqn";
102
103 /*------------------------------------------------------------------------------
104 REGISTRATION FUNCTION
105 */
106
107 /**
108 This is the function called from "IMPORT fprops"
109
110 It sets up the functions contained in this external library
111 */
112 extern
113 ASC_EXPORT int fprops_register(){
114 int result = 0;
115
116 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
117
118 #define CALCFN(NAME,INPUTS,OUTPUTS) \
119 result += CreateUserFunctionBlackBox(#NAME \
120 , asc_fprops_prepare \
121 , NAME##_calc /* value */ \
122 , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
123 , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
124 , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
125 , INPUTS,OUTPUTS /* inputs, outputs */ \
126 , NAME##_help /* help text */ \
127 , 0.0 \
128 ) /* returns 0 on success */
129
130 #define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
131 result += CreateUserFunctionBlackBox(#NAME \
132 , asc_fprops_prepare \
133 , NAME##_calc /* value */ \
134 , NAME##_calc /* derivatives */ \
135 , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
136 , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
137 , INPUTS,OUTPUTS /* inputs, outputs */ \
138 , NAME##_help /* help text */ \
139 , 0.0 \
140 ) /* returns 0 on success */
141
142 CALCFNDERIV(fprops_p,2,1);
143 CALCFN(fprops_u,2,1);
144 CALCFN(fprops_s,2,1);
145 CALCFN(fprops_h,2,1);
146 CALCFN(fprops_a,2,1);
147 CALCFN(fprops_g,2,1);
148 CALCFN(fprops_cp,2,1);
149 CALCFN(fprops_cv,2,1);
150 CALCFN(fprops_w,2,1);
151 CALCFN(fprops_phsx_vT,2,4);
152 CALCFN(fprops_Tvsx_ph,2,4);
153
154 #undef CALCFN
155
156 if(result){
157 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
158 }
159 return result;
160 }
161
162 /**
163 'fprops_prepare' just gets the data member and checks that it's
164 valid, and stores it in the blackbox data field.
165 */
166 int asc_fprops_prepare(struct BBoxInterp *bbox,
167 struct Instance *data,
168 struct gl_list_t *arglist
169 ){
170 struct Instance *compinst, *typeinst;
171 const char *comp, *type = NULL;
172
173 fprops_symbols[0] = AddSymbol("component");
174 fprops_symbols[1] = AddSymbol("type");
175
176 /* get the component name */
177 compinst = ChildByChar(data,COMPONENT_SYM);
178 if(!compinst){
179 ERROR_REPORTER_HERE(ASC_USER_ERROR
180 ,"Couldn't locate 'component' in DATA, please check usage of FPROPS."
181 );
182 return 1;
183 }
184 if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
185 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
186 return 1;
187 }
188 comp = SCP(SYMC_INST(compinst)->value);
189 if(comp==NULL || strlen(comp)==0){
190 ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
191 return 1;
192 }
193
194 /* get the component correlation type (FPROPS doesn't mind if none given) */
195 typeinst = ChildByChar(data,TYPE_SYM);
196 if(typeinst){
197 if(InstanceKind(typeinst)!=SYMBOL_CONSTANT_INST){
198 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'type' must be a symbol_constant");
199 return 1;
200 }
201 type = SCP(SYMC_INST(typeinst)->value);
202 CONSOLE_DEBUG("TYPE: %s",type?type:"(null)");
203 if(type && strlen(type)==0)type = NULL;
204 }
205
206 bbox->user_data = (void *)fprops_fluid(comp,type);
207 if(bbox->user_data == NULL){
208 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name/type was not recognised. Check the source-code for for the supported species.");
209 return 1;
210 }
211
212 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s'%s%s%s OK.\n"
213 ,comp
214 ,type?" type '":""
215 ,type?type:""
216 ,type?"'":""
217 );
218 return 0;
219 }
220
221 /*------------------------------------------------------------------------------
222 EVALULATION ROUTINES
223 */
224
225 #define CALCPREPARE(NIN,NOUT) \
226 /* a few checks about the input requirements */ \
227 if(ninputs != NIN)return -1; \
228 if(noutputs != NOUT)return -2; \
229 if(inputs==NULL)return -3; \
230 if(outputs==NULL)return -4; \
231 if(bbox==NULL)return -5; \
232 \
233 /* the 'user_data' in the black box object will contain the */\
234 /* coefficients required for this fluid; cast it to the required form: */\
235 const PureFluid *FLUID = (const PureFluid *)bbox->user_data;\
236 FpropsError err=FPROPS_NO_ERROR;
237
238 /**
239 Evaluation function for 'fprops_p'.
240 @param inputs array with values of inputs T and rho.
241 @param outputs array with just value of p
242 @param jacobian, the partial derivative df/dx, where
243 each row is df[i]/dx[j] over each j for the y_out[i] of
244 matching index. The jacobian array is 1-D, row major, i.e.
245 df[i]/dx[j] -> jacobian[i*ninputs+j].
246 @return 0 on success
247 */
248 int fprops_p_calc(struct BBoxInterp *bbox,
249 int ninputs, int noutputs,
250 double *inputs, double *outputs,
251 double *jacobian
252 ){
253 CALCPREPARE(2,1);
254
255 /* first input is temperature, second is density */
256 if(bbox->task == bb_func_eval){
257 FluidState S = fprops_set_Trho(inputs[0],inputs[1], FLUID, &err);
258 outputs[0] = fprops_p(S, &err);
259 }else{
260 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
261 FluidState S = fprops_set_Trho(inputs[0],inputs[1], FLUID, &err);
262 jacobian[0*1+0] = fprops_dpdT_rho(S, &err);
263 jacobian[0*1+1] = fprops_dpdrho_T(S, &err);
264 }
265
266 /* no need to worry about error states etc. */
267 return 0;
268 }
269
270
271 /**
272 Evaluation function for 'fprops_u'
273 @param jacobian ignored
274 @return 0 on success
275 */
276 int fprops_u_calc(struct BBoxInterp *bbox,
277 int ninputs, int noutputs,
278 double *inputs, double *outputs,
279 double *jacobian
280 ){
281 CALCPREPARE(2,1);
282 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
283
284 /* first input is temperature, second is density */
285 if(bbox->task == bb_func_eval){
286 outputs[0] = fprops_u(S, &err);
287 }else{
288 jacobian[0*1+0] = fprops_dudT_rho(S, &err);
289 jacobian[0*1+1] = fprops_dudrho_T(S, &err);
290 }
291
292 /* no need to worry about error states etc. */
293 return 0;
294 }
295
296
297 /**
298 Evaluation function for 'fprops_s'
299 @param jacobian ignored
300 @return 0 on success
301 */
302 int fprops_s_calc(struct BBoxInterp *bbox,
303 int ninputs, int noutputs,
304 double *inputs, double *outputs,
305 double *jacobian
306 ){
307 CALCPREPARE(2,1);
308 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
309
310 /* first input is temperature, second is density */
311 outputs[0] = fprops_s(S, &err);
312
313 /* no need to worry about error states etc. */
314 return 0;
315 }
316
317
318 /**
319 Evaluation function for 'fprops_h'
320 @param jacobian ignored
321 @return 0 on success
322 */
323 int fprops_h_calc(struct BBoxInterp *bbox,
324 int ninputs, int noutputs,
325 double *inputs, double *outputs,
326 double *jacobian
327 ){
328 CALCPREPARE(2,1);
329 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
330
331 /* first input is temperature, second is density */
332 if(bbox->task == bb_func_eval){
333 outputs[0] = fprops_h(S, &err);
334 }else{
335 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
336 jacobian[0*1+0] = fprops_dhdT_rho(S, &err);
337 jacobian[0*1+1] = fprops_dhdrho_T(S, &err);
338 }
339
340 /* no need to worry about error states etc. */
341 return 0;
342 }
343
344
345 /**
346 Evaluation function for 'fprops_a'
347 @param jacobian ignored
348 @return 0 on success
349 */
350 int fprops_a_calc(struct BBoxInterp *bbox,
351 int ninputs, int noutputs,
352 double *inputs, double *outputs,
353 double *jacobian
354 ){
355 CALCPREPARE(2,1);
356 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
357
358 /* first input is temperature, second is density */
359 outputs[0] = fprops_a(S, &err);
360
361 /* no need to worry about error states etc. */
362 return 0;
363 }
364
365
366 /**
367 Evaluation function for 'fprops_g'
368 @param jacobian ignored
369 @return 0 on success
370 */
371 int fprops_g_calc(struct BBoxInterp *bbox,
372 int ninputs, int noutputs,
373 double *inputs, double *outputs,
374 double *jacobian
375 ){
376 CALCPREPARE(2,1);
377 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
378
379 /* first input is temperature, second is density */
380 outputs[0] = fprops_g(S, &err);
381
382 /* no need to worry about error states etc. */
383 return 0;
384 }
385
386
387 /**
388 Evaluation function for 'fprops_cp'
389 @param jacobian ignored
390 @return 0 on success
391 */
392 int fprops_cp_calc(struct BBoxInterp *bbox,
393 int ninputs, int noutputs,
394 double *inputs, double *outputs,
395 double *jacobian
396 ){
397 CALCPREPARE(2,1);
398 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
399
400 /* first input is temperature, second is density */
401 outputs[0] = fprops_cp(S, &err);
402
403 /* no need to worry about error states etc. */
404 return 0;
405 }
406
407
408 /**
409 Evaluation function for 'fprops_cv'
410 @param jacobian ignored
411 @return 0 on success
412 */
413 int fprops_cv_calc(struct BBoxInterp *bbox,
414 int ninputs, int noutputs,
415 double *inputs, double *outputs,
416 double *jacobian
417 ){
418 CALCPREPARE(2,1);
419 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
420
421 /* first input is temperature, second is density */
422 outputs[0] = fprops_cv(S, &err);
423
424 /* no need to worry about error states etc. */
425 return 0;
426 }
427
428
429 /**
430 Evaluation function for 'fprops_w'
431 @param jacobian ignored
432 @return 0 on success
433 */
434 int fprops_w_calc(struct BBoxInterp *bbox,
435 int ninputs, int noutputs,
436 double *inputs, double *outputs,
437 double *jacobian
438 ){
439 CALCPREPARE(2,1);
440 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
441
442 /* first input is temperature, second is density */
443 outputs[0] = fprops_w(S, &err);
444
445 /* no need to worry about error states etc. */
446 return 0;
447 }
448
449
450 /**
451 Evaluation function for 'fprops_phsx_vT'
452 @return 0 on success
453 */
454 int fprops_phsx_vT_calc(struct BBoxInterp *bbox,
455 int ninputs, int noutputs,
456 double *inputs, double *outputs,
457 double *jacobian
458 ){
459 CALCPREPARE(2,4);
460
461 double rho = 1./inputs[0];
462 double T = inputs[1];
463 double p_sat, rho_f, rho_g;
464
465 if(T < TCRIT(FLUID)){
466 fprops_sat_T(T, &p_sat, &rho_f, &rho_g, FLUID, &err);
467
468 if(rho < rho_f && rho > rho_g){
469 /* saturated */
470 double vf = 1./rho_f;
471 double vg = 1./rho_g;
472 double x = (inputs[0] - vf) /(vg - vf);
473 FluidState Sf = fprops_set_Trho(T, rho_f, FLUID, &err);
474 double sf = fprops_s(Sf, &err);
475 double hf = fprops_h(Sf, &err);
476 FluidState Sg = fprops_set_Trho(T, rho_g, FLUID, &err);
477 double sg = fprops_s(Sg, &err);
478 double hg = fprops_h(Sg, &err);
479 outputs[0] = p_sat;
480 outputs[1] = hf + x * (hg-hf);
481 outputs[2] = sf + x * (sg-sf);
482 outputs[3] = x;
483 /* maybe there was an error solving the saturation state? */
484 return err;
485 }
486 }
487
488 /* non-saturated */
489 FluidState S = fprops_set_Trho(T, rho, FLUID, &err);
490 outputs[0] = fprops_p(S, &err);
491 outputs[1] = fprops_h(S, &err);
492 outputs[2] = fprops_s(S, &err);
493 outputs[3] = rho < RHOCRIT(FLUID) ? 1 : 0;
494 return 0;
495 }
496
497
498 /**
499 Evaluation function for 'fprops_Tvsx_ph'
500 @return 0 on success
501 */
502 int fprops_Tvsx_ph_calc(struct BBoxInterp *bbox,
503 int ninputs, int noutputs,
504 double *inputs, double *outputs,
505 double *jacobian
506 ){
507 CALCPREPARE(2,4);
508
509 static const PureFluid *last = NULL;
510 static double p,h,T,v,s,x;
511 if(last == FLUID && p == inputs[0] && h == inputs[1]){
512 outputs[0] = T;
513 outputs[1] = v;
514 outputs[2] = s;
515 outputs[3] = x;
516 return 0;
517 }
518
519 p = inputs[0];
520 h = inputs[1];
521
522 double hft, pt, rhoft,rhogt;
523 fprops_triple_point(&pt,&rhoft,&rhogt,FLUID,&err);
524 if(err){
525 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve triple point for %s.",FLUID->name);
526 return 5;
527 }
528 FluidState Sft = fprops_set_Trho(TTRIP(FLUID),rhoft,FLUID,&err);
529 hft = fprops_h(Sft, &err);
530 if(h < hft){
531 ERROR_REPORTER_HERE(ASC_PROG_ERR
532 ,"Input enthalpy %f kJ/kg is below triple point liquid enthalpy %f kJ/kg"
533 ,h/1e3,hft/1e3
534 );
535 return 6;
536 }
537
538 if(p < pt){
539 ERROR_REPORTER_HERE(ASC_PROG_ERR
540 ,"Input pressure %f bar is below triple point pressure %f bar"
541 ,p/1e5,pt/1e5
542 );
543 outputs[0] = TTRIP(FLUID);
544 outputs[1] = 1./ rhoft;
545 outputs[2] = FLUID->s_fn(TTRIP(FLUID), rhoft, FLUID->data, &err);
546 outputs[3] = 0;
547 return 7;
548 }
549
550 if(p < PCRIT(FLUID)){
551 double T_sat, rho_f, rho_g;
552
553 fprops_sat_p(p, &T_sat, &rho_f, &rho_g, FLUID, &err);
554 if(err){
555 ERROR_REPORTER_HERE(ASC_PROG_ERR
556 , "Failed to solve saturation state of %s for p = %f bar < pc (= %f bar)"
557 , FLUID->name, p/1e5,PCRIT(FLUID)/1e5
558 );
559 outputs[0] = TTRIP(FLUID);
560 outputs[1] = 1./rhoft;
561 outputs[2] = FLUID->s_fn(TTRIP(FLUID), rhoft, FLUID->data, &err);
562 outputs[3] = 0;
563 return 8;
564 }
565
566 FluidState Sf = fprops_set_Trho(T_sat,rho_f,FLUID,&err);
567 FluidState Sg = fprops_set_Trho(T_sat,rho_g,FLUID,&err);
568 double hf = fprops_h(Sf, &err);
569 double hg = fprops_h(Sg, &err);
570
571 if(hf < h && h < hg){
572 /* saturated */
573 double vf = 1./rho_f;
574 double vg = 1./rho_g;
575 double sf = fprops_s(Sf, &err);
576 double sg = fprops_s(Sg, &err);
577 T = T_sat;
578 x = (h - hf) /(hg - hf);
579 v = vf + x * (vg-vf);
580 s = sf + x * (sg-sf);
581 last = FLUID;
582 outputs[0] = T;
583 outputs[1] = v;
584 outputs[2] = s;
585 outputs[3] = x;
586 #ifdef ASC_FPROPS_DEBUG
587 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Saturated state, p=%f bar, h = %f kJ/kg",p/1e5,h/1e3);
588 #endif
589 return 0;
590 }
591 }
592
593 double rho;
594 fprops_solve_ph(p,h, &T, &rho, 0, FLUID, &err);
595 if(err){
596 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for (p,h): %s",fprops_error(err));
597 return 9;
598 }
599 /* non-saturated */
600 v = 1./rho;
601 FluidState S = fprops_set_Trho(T,rho,FLUID,&err);
602 s = fprops_s(S, &err);
603 x = (v > 1./RHOCRIT(FLUID)) ? 1 : 0;
604 last = FLUID;
605 outputs[0] = T;
606 outputs[1] = v;
607 outputs[2] = s;
608 outputs[3] = x;
609 #ifdef ASC_FPROPS_DEBUG
610 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Non-saturated state, p = %f bar, h = %f kJ/kg",p/1e5,h/1e3);
611 #endif
612 return 0;
613 }
614
615
616

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