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

Annotation of /trunk/models/johnpye/fprops/asc_helmholtz.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2275 - (hide annotations) (download) (as text)
Wed Aug 11 13:58:25 2010 UTC (13 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 13771 byte(s)
Fixing error return from asc_helmholtz function.
The bug with 'defaultall' functionality is still here... see boiler_simple_test mode in models/johnpye/fprops/rankine_fprops.a4c.
1 jpye 1822 /* 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 helmholtz.c to allow access to the routine from ASCEND.
20     */
21    
22     /* include the external function API from libascend... */
23 jpye 2018 #include <ascend/compiler/extfunc.h>
24 jpye 1822
25     /* include error reporting API as well, so we can send messages to user */
26 jpye 2018 #include <ascend/utilities/error.h>
27 jpye 1822
28     /* for accessing the DATA instance */
29 jpye 2018 #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 jpye 1822
42     /* the code that we're wrapping... */
43     #include "helmholtz.h"
44 jpye 2262 #include "sat.h"
45 jpye 2272 #include "solve_ph.h"
46 jpye 1822
47 jpye 1832 /* for the moment, species data are defined in C code, we'll implement something
48     better later on, hopefully. */
49     #include "ammonia.h"
50 jpye 1849 #include "nitrogen.h"
51 jpye 1995 #include "hydrogen.h"
52     #include "water.h"
53 jpye 2110 #include "carbondioxide.h"
54 jpye 1832
55 jpye 1822 #ifndef ASC_EXPORT
56     # error "Where is ASC_EXPORT?"
57     #endif
58    
59    
60     /*------------------------------------------------------------------------------
61     FORWARD DECLARATIONS
62     */
63    
64 jpye 1825 ExtBBoxInitFunc helmholtz_prepare;
65 jpye 1822 ExtBBoxFunc helmholtz_p_calc;
66 jpye 1825 ExtBBoxFunc helmholtz_u_calc;
67 jpye 1903 ExtBBoxFunc helmholtz_s_calc;
68 jpye 1825 ExtBBoxFunc helmholtz_h_calc;
69 jpye 1903 ExtBBoxFunc helmholtz_a_calc;
70 jpye 2124 ExtBBoxFunc helmholtz_g_calc;
71 jpye 2236 ExtBBoxFunc helmholtz_phsx_vT_calc;
72 jpye 2272 ExtBBoxFunc helmholtz_Tvsx_ph_calc;
73 jpye 1822
74     /*------------------------------------------------------------------------------
75     GLOBALS
76     */
77    
78     /* place to store symbols needed for accessing ASCEND's instance tree */
79     static symchar *helmholtz_symbols[1];
80     #define COMPONENT_SYM helmholtz_symbols[0]
81    
82 jpye 1825 static const char *helmholtz_p_help = "Calculate pressure from temperature and density, using Helmholtz fundamental correlation";
83     static const char *helmholtz_u_help = "Calculate specific internal energy from temperature and density, using Helmholtz fundamental correlation";
84 jpye 1903 static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
85 jpye 1825 static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
86 jpye 1903 static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";
87 jpye 2124 static const char *helmholtz_g_help = "Calculate specific Gibbs energy from temperature and density, using Helmholtz fundamental correlation";
88 jpye 1825
89 jpye 2236 static const char *helmholtz_phsx_vT_help = "Calculate p, h, s, x from temperature and density, using FPROPS/Helmholtz eqn";
90    
91 jpye 2272 static const char *helmholtz_Tvsx_ph_help = "Calculate T, v, s, x from pressure and enthalpy, using FPROPS/Helmholtz eqn";
92 jpye 2124
93 jpye 2272
94 jpye 1822 /*------------------------------------------------------------------------------
95     REGISTRATION FUNCTION
96     */
97    
98     /**
99     This is the function called from "IMPORT helmholtz"
100    
101     It sets up the functions contained in this external library
102     */
103     extern
104     ASC_EXPORT int helmholtz_register(){
105     int result = 0;
106    
107 jpye 1920 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
108 jpye 1822
109 jpye 1825 #define CALCFN(NAME,INPUTS,OUTPUTS) \
110     result += CreateUserFunctionBlackBox(#NAME \
111     , helmholtz_prepare \
112     , NAME##_calc /* value */ \
113     , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
114     , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
115     , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
116     , INPUTS,OUTPUTS /* inputs, outputs */ \
117     , NAME##_help /* help text */ \
118     , 0.0 \
119     ) /* returns 0 on success */
120 jpye 1822
121 jpye 1909 #define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
122     result += CreateUserFunctionBlackBox(#NAME \
123     , helmholtz_prepare \
124     , NAME##_calc /* value */ \
125     , NAME##_calc /* derivatives */ \
126     , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
127     , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
128     , INPUTS,OUTPUTS /* inputs, outputs */ \
129     , NAME##_help /* help text */ \
130     , 0.0 \
131     ) /* returns 0 on success */
132    
133     CALCFNDERIV(helmholtz_p,2,1);
134 jpye 1825 CALCFN(helmholtz_u,2,1);
135 jpye 1903 CALCFN(helmholtz_s,2,1);
136 jpye 1825 CALCFN(helmholtz_h,2,1);
137 jpye 1903 CALCFN(helmholtz_a,2,1);
138 jpye 2124 CALCFN(helmholtz_g,2,1);
139 jpye 2236 CALCFN(helmholtz_phsx_vT,2,4);
140 jpye 2272 CALCFN(helmholtz_Tvsx_ph,2,4);
141 jpye 1825
142     #undef CALCFN
143    
144 jpye 1822 if(result){
145     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
146     }
147     return result;
148     }
149    
150 jpye 1825 /**
151     'helmholtz_prepare' just gets the data member and checks that it's
152     valid, and stores it in the blackbox data field.
153 jpye 1822 */
154 jpye 1825 int helmholtz_prepare(struct BBoxInterp *bbox,
155 jpye 1822 struct Instance *data,
156     struct gl_list_t *arglist
157     ){
158     struct Instance *compinst;
159     const char *comp;
160    
161     helmholtz_symbols[0] = AddSymbol("component");
162    
163     /* get the data file name (we will look for this file in the ASCENDLIBRARY path) */
164     compinst = ChildByChar(data,COMPONENT_SYM);
165     if(!compinst){
166     ERROR_REPORTER_HERE(ASC_USER_ERROR
167     ,"Couldn't locate 'component' in DATA, please check usage of HELMHOLTZ."
168     );
169     return 1;
170     }
171     if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
172     ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
173     return 1;
174     }
175     comp = SCP(SYMC_INST(compinst)->value);
176     CONSOLE_DEBUG("COMPONENT: %s",comp);
177     if(comp==NULL || strlen(comp)==0){
178     ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
179     return 1;
180     }
181    
182 jpye 1849 if(strcmp(comp,"ammonia")==0){
183     bbox->user_data = (void*)&helmholtz_data_ammonia;
184     }else if(strcmp(comp,"nitrogen")==0){
185     bbox->user_data = (void*)&helmholtz_data_nitrogen;
186 jpye 1995 }else if(strcmp(comp,"hydrogen")==0){
187     bbox->user_data = (void*)&helmholtz_data_hydrogen;
188     }else if(strcmp(comp,"water")==0){
189     bbox->user_data = (void*)&helmholtz_data_water;
190 jpye 2110 }else if(strcmp(comp,"carbondioxide")==0){
191     bbox->user_data = (void*)&helmholtz_data_carbondioxide;
192 jpye 1849 }else{
193 jpye 1995 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
194 jpye 1822 }
195    
196 jpye 1849 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
197 jpye 1822 return 0;
198     }
199    
200 jpye 1825 /*------------------------------------------------------------------------------
201     EVALULATION ROUTINES
202     */
203    
204 jpye 2124 #define CALCPREPARE(NIN,NOUT) \
205 jpye 1825 /* a few checks about the input requirements */ \
206 jpye 2124 if(ninputs != NIN)return -1; \
207     if(noutputs != NOUT)return -2; \
208 jpye 1825 if(inputs==NULL)return -3; \
209     if(outputs==NULL)return -4; \
210     if(bbox==NULL)return -5; \
211     \
212     /* the 'user_data' in the black box object will contain the */\
213     /* coefficients required for this fluid; cast it to the required form: */\
214 jpye 2272 const HelmholtzData *helmholtz_data = (const HelmholtzData *)bbox->user_data
215 jpye 1825
216 jpye 1822 /**
217 jpye 1825 Evaluation function for 'helmholtz_p'.
218 jpye 1909 @param inputs array with values of inputs T and rho.
219     @param outputs array with just value of p
220     @param jacobian, the partial derivative df/dx, where
221     each row is df[i]/dx[j] over each j for the y_out[i] of
222     matching index. The jacobian array is 1-D, row major, i.e.
223     df[i]/dx[j] -> jacobian[i*ninputs+j].
224 jpye 1822 @return 0 on success
225     */
226     int helmholtz_p_calc(struct BBoxInterp *bbox,
227     int ninputs, int noutputs,
228     double *inputs, double *outputs,
229     double *jacobian
230     ){
231 jpye 2124 CALCPREPARE(2,1);
232 jpye 1822
233 jpye 2236 /* first input is temperature, second is density */
234 jpye 1909 if(bbox->task == bb_func_eval){
235     outputs[0] = helmholtz_p(inputs[0], inputs[1], helmholtz_data);
236     }else{
237     //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
238     jacobian[0*1+0] = helmholtz_dpdT_rho(inputs[0], inputs[1], helmholtz_data);
239     jacobian[0*1+1] = helmholtz_dpdrho_T(inputs[0], inputs[1], helmholtz_data);
240     }
241 jpye 1822
242     /* no need to worry about error states etc. */
243 jpye 1825 return 0;
244     }
245 jpye 1822
246 jpye 1825
247     /**
248     Evaluation function for 'helmholtz_u'
249     @param jacobian ignored
250     @return 0 on success
251     */
252     int helmholtz_u_calc(struct BBoxInterp *bbox,
253     int ninputs, int noutputs,
254     double *inputs, double *outputs,
255     double *jacobian
256     ){
257 jpye 2124 CALCPREPARE(2,1);
258 jpye 1825
259 jpye 2236 /* first input is temperature, second is density */
260 jpye 1923 if(bbox->task == bb_func_eval){
261     outputs[0] = helmholtz_u(inputs[0], inputs[1], helmholtz_data);
262     }else{
263     jacobian[0*1+0] = helmholtz_dudT_rho(inputs[0], inputs[1], helmholtz_data);
264     jacobian[0*1+1] = helmholtz_dudrho_T(inputs[0], inputs[1], helmholtz_data);
265     }
266 jpye 1825
267     /* no need to worry about error states etc. */
268 jpye 1822 return 0;
269     }
270 jpye 1825
271    
272     /**
273 jpye 2124 Evaluation function for 'helmholtz_s'
274 jpye 1825 @param jacobian ignored
275     @return 0 on success
276     */
277 jpye 1903 int helmholtz_s_calc(struct BBoxInterp *bbox,
278     int ninputs, int noutputs,
279     double *inputs, double *outputs,
280     double *jacobian
281     ){
282 jpye 2124 CALCPREPARE(2,1);
283 jpye 1903
284 jpye 2236 /* first input is temperature, second is density */
285 jpye 1903 outputs[0] = helmholtz_s(inputs[0], inputs[1], helmholtz_data);
286    
287     /* no need to worry about error states etc. */
288     return 0;
289     }
290    
291    
292     /**
293     Evaluation function for 'helmholtz_h'
294     @param jacobian ignored
295     @return 0 on success
296     */
297 jpye 1825 int helmholtz_h_calc(struct BBoxInterp *bbox,
298     int ninputs, int noutputs,
299     double *inputs, double *outputs,
300     double *jacobian
301     ){
302 jpye 2124 CALCPREPARE(2,1);
303 jpye 1825
304 jpye 2236 /* first input is temperature, second is density */
305 jpye 1920 if(bbox->task == bb_func_eval){
306     outputs[0] = helmholtz_h(inputs[0], inputs[1], helmholtz_data);
307     }else{
308     //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
309     jacobian[0*1+0] = helmholtz_dhdT_rho(inputs[0], inputs[1], helmholtz_data);
310     jacobian[0*1+1] = helmholtz_dhdrho_T(inputs[0], inputs[1], helmholtz_data);
311     }
312 jpye 1825
313     /* no need to worry about error states etc. */
314     return 0;
315     }
316    
317    
318 jpye 1903 /**
319 jpye 2124 Evaluation function for 'helmholtz_a'
320 jpye 1903 @param jacobian ignored
321     @return 0 on success
322     */
323     int helmholtz_a_calc(struct BBoxInterp *bbox,
324     int ninputs, int noutputs,
325     double *inputs, double *outputs,
326     double *jacobian
327     ){
328 jpye 2124 CALCPREPARE(2,1);
329 jpye 1825
330 jpye 2236 /* first input is temperature, second is density */
331 jpye 1903 outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);
332 jpye 1825
333 jpye 1903 /* no need to worry about error states etc. */
334     return 0;
335     }
336    
337    
338 jpye 2124 /**
339     Evaluation function for 'helmholtz_g'
340     @param jacobian ignored
341     @return 0 on success
342     */
343     int helmholtz_g_calc(struct BBoxInterp *bbox,
344     int ninputs, int noutputs,
345     double *inputs, double *outputs,
346     double *jacobian
347     ){
348     CALCPREPARE(2,1);
349 jpye 1903
350 jpye 2236 /* first input is temperature, second is density */
351 jpye 2124 outputs[0] = helmholtz_g(inputs[0], inputs[1], helmholtz_data);
352    
353     /* no need to worry about error states etc. */
354     return 0;
355     }
356    
357    
358 jpye 2236
359 jpye 2124 /**
360 jpye 2236 Evaluation function for 'helmholtz_phsx_vT'
361     @return 0 on success
362     */
363     int helmholtz_phsx_vT_calc(struct BBoxInterp *bbox,
364     int ninputs, int noutputs,
365     double *inputs, double *outputs,
366     double *jacobian
367     ){
368     CALCPREPARE(2,4);
369    
370     double rho = 1./inputs[0];
371     double T = inputs[1];
372     double p_sat, rho_f, rho_g;
373    
374     if(T < helmholtz_data->T_c){
375     int res = fprops_sat_T(T, &p_sat, &rho_f, &rho_g, helmholtz_data);
376    
377     if(rho < rho_f && rho > rho_g){
378     /* saturated */
379     double vf = 1./rho_f;
380     double vg = 1./rho_g;
381     double x = (inputs[0] - vf) /(vg - vf);
382     double sf = helmholtz_s(T,rho_f, helmholtz_data);
383     double hf = helmholtz_h(T,rho_f, helmholtz_data);
384     double sg = helmholtz_s(T,rho_g, helmholtz_data);
385     double hg = helmholtz_h(T,rho_g, helmholtz_data);
386     outputs[0] = p_sat;
387     outputs[1] = hf + x * (hg-hf);
388     outputs[2] = sf + x * (sg-sf);
389     outputs[3] = x;
390     /* maybe there was an error solving the saturation state? */
391     return res;
392     }
393     }
394    
395     /* non-saturated */
396     outputs[0] = helmholtz_p(T,rho, helmholtz_data);
397     outputs[1] = helmholtz_h(T,rho, helmholtz_data);
398     outputs[2] = helmholtz_s(T,rho, helmholtz_data);
399 jpye 2274 outputs[3] = rho < helmholtz_data->rho_c ? 1 : 0;
400 jpye 2236 return 0;
401     }
402    
403    
404     /**
405 jpye 2272 Evaluation function for 'helmholtz_Tvsx_ph'
406 jpye 2124 @return 0 on success
407     */
408 jpye 2272 int helmholtz_Tvsx_ph_calc(struct BBoxInterp *bbox,
409 jpye 2124 int ninputs, int noutputs,
410     double *inputs, double *outputs,
411     double *jacobian
412     ){
413 jpye 2272 CALCPREPARE(2,4);
414 jpye 2124
415 jpye 2272 static const HelmholtzData *last = NULL;
416     static double p,h,T,v,s,x;
417     int res;
418     if(last == helmholtz_data && p == inputs[0] && h == inputs[1]){
419     outputs[0] = T;
420     outputs[1] = v;
421     outputs[2] = s;
422     outputs[3] = x;
423     return 0;
424     }
425 jpye 2124
426 jpye 2272 p = inputs[0];
427     h = inputs[1];
428 jpye 2124
429 jpye 2272 if(p < fprops_pc(helmholtz_data)){
430     double T_sat, rho_f, rho_g;
431     res = fprops_sat_p(p, &T_sat, &rho_f, &rho_g, helmholtz_data);
432     if(res)return 1;
433     double hf = helmholtz_h(T_sat,rho_f, helmholtz_data);
434     double hg = helmholtz_h(T_sat,rho_g, helmholtz_data);
435 jpye 2124
436 jpye 2272 if(hf < h && h < hg){
437     /* saturated */
438     double vf = 1./rho_f;
439     double vg = 1./rho_g;
440     double sf = helmholtz_s(T_sat,rho_f, helmholtz_data);
441     double sg = helmholtz_s(T_sat,rho_g, helmholtz_data);
442     T = T_sat;
443     x = (h - hf) /(hg - hf);
444     v = vf + x * (vg-vf);
445     s = sf + x * (sg-sf);
446     last = helmholtz_data;
447     outputs[0] = T;
448     outputs[1] = v;
449     outputs[2] = s;
450     outputs[3] = x;
451     return 0;
452     }
453     }
454    
455     double rho;
456     res = fprops_solve_ph(p,h, &T, &rho, 0, helmholtz_data);
457     /* non-saturated */
458     v = 1./rho;
459     s = helmholtz_s(T,rho, helmholtz_data);
460     x = (v > 1./helmholtz_data->rho_c) ? 1 : 0;
461     last = helmholtz_data;
462     outputs[0] = T;
463     outputs[1] = v;
464     outputs[2] = s;
465     outputs[3] = x;
466 jpye 2275 return res;
467 jpye 2124 }
468    
469 jpye 2272
470    

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