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

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