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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2275 - (show annotations) (download) (as text)
Wed Aug 11 13:58:25 2010 UTC (13 years, 10 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 /* 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 #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 "helmholtz.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 "ammonia.h"
50 #include "nitrogen.h"
51 #include "hydrogen.h"
52 #include "water.h"
53 #include "carbondioxide.h"
54
55 #ifndef ASC_EXPORT
56 # error "Where is ASC_EXPORT?"
57 #endif
58
59
60 /*------------------------------------------------------------------------------
61 FORWARD DECLARATIONS
62 */
63
64 ExtBBoxInitFunc helmholtz_prepare;
65 ExtBBoxFunc helmholtz_p_calc;
66 ExtBBoxFunc helmholtz_u_calc;
67 ExtBBoxFunc helmholtz_s_calc;
68 ExtBBoxFunc helmholtz_h_calc;
69 ExtBBoxFunc helmholtz_a_calc;
70 ExtBBoxFunc helmholtz_g_calc;
71 ExtBBoxFunc helmholtz_phsx_vT_calc;
72 ExtBBoxFunc helmholtz_Tvsx_ph_calc;
73
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 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 static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
85 static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
86 static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";
87 static const char *helmholtz_g_help = "Calculate specific Gibbs energy from temperature and density, using Helmholtz fundamental correlation";
88
89 static const char *helmholtz_phsx_vT_help = "Calculate p, h, s, x from temperature and density, using FPROPS/Helmholtz eqn";
90
91 static const char *helmholtz_Tvsx_ph_help = "Calculate T, v, s, x from pressure and enthalpy, using FPROPS/Helmholtz eqn";
92
93
94 /*------------------------------------------------------------------------------
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 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
108
109 #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
121 #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 CALCFN(helmholtz_u,2,1);
135 CALCFN(helmholtz_s,2,1);
136 CALCFN(helmholtz_h,2,1);
137 CALCFN(helmholtz_a,2,1);
138 CALCFN(helmholtz_g,2,1);
139 CALCFN(helmholtz_phsx_vT,2,4);
140 CALCFN(helmholtz_Tvsx_ph,2,4);
141
142 #undef CALCFN
143
144 if(result){
145 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
146 }
147 return result;
148 }
149
150 /**
151 'helmholtz_prepare' just gets the data member and checks that it's
152 valid, and stores it in the blackbox data field.
153 */
154 int helmholtz_prepare(struct BBoxInterp *bbox,
155 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 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 }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 }else if(strcmp(comp,"carbondioxide")==0){
191 bbox->user_data = (void*)&helmholtz_data_carbondioxide;
192 }else{
193 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
194 }
195
196 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
197 return 0;
198 }
199
200 /*------------------------------------------------------------------------------
201 EVALULATION ROUTINES
202 */
203
204 #define CALCPREPARE(NIN,NOUT) \
205 /* a few checks about the input requirements */ \
206 if(ninputs != NIN)return -1; \
207 if(noutputs != NOUT)return -2; \
208 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 const HelmholtzData *helmholtz_data = (const HelmholtzData *)bbox->user_data
215
216 /**
217 Evaluation function for 'helmholtz_p'.
218 @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 @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 CALCPREPARE(2,1);
232
233 /* first input is temperature, second is density */
234 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
242 /* no need to worry about error states etc. */
243 return 0;
244 }
245
246
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 CALCPREPARE(2,1);
258
259 /* first input is temperature, second is density */
260 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
267 /* no need to worry about error states etc. */
268 return 0;
269 }
270
271
272 /**
273 Evaluation function for 'helmholtz_s'
274 @param jacobian ignored
275 @return 0 on success
276 */
277 int helmholtz_s_calc(struct BBoxInterp *bbox,
278 int ninputs, int noutputs,
279 double *inputs, double *outputs,
280 double *jacobian
281 ){
282 CALCPREPARE(2,1);
283
284 /* first input is temperature, second is density */
285 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 int helmholtz_h_calc(struct BBoxInterp *bbox,
298 int ninputs, int noutputs,
299 double *inputs, double *outputs,
300 double *jacobian
301 ){
302 CALCPREPARE(2,1);
303
304 /* first input is temperature, second is density */
305 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
313 /* no need to worry about error states etc. */
314 return 0;
315 }
316
317
318 /**
319 Evaluation function for 'helmholtz_a'
320 @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 CALCPREPARE(2,1);
329
330 /* first input is temperature, second is density */
331 outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);
332
333 /* no need to worry about error states etc. */
334 return 0;
335 }
336
337
338 /**
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
350 /* first input is temperature, second is density */
351 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
359 /**
360 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 outputs[3] = rho < helmholtz_data->rho_c ? 1 : 0;
400 return 0;
401 }
402
403
404 /**
405 Evaluation function for 'helmholtz_Tvsx_ph'
406 @return 0 on success
407 */
408 int helmholtz_Tvsx_ph_calc(struct BBoxInterp *bbox,
409 int ninputs, int noutputs,
410 double *inputs, double *outputs,
411 double *jacobian
412 ){
413 CALCPREPARE(2,4);
414
415 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
426 p = inputs[0];
427 h = inputs[1];
428
429 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
436 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 return res;
467 }
468
469
470

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