/[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 2147 - (show annotations) (download) (as text)
Fri Jan 22 04:15:50 2010 UTC (10 years, 5 months ago) by kchittur
File MIME type: text/x-csrc
File size: 11634 byte(s)


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
45 /* for the moment, species data are defined in C code, we'll implement something
46 better later on, hopefully. */
47 #include "ammonia.h"
48 #include "nitrogen.h"
49 #include "hydrogen.h"
50 #include "water.h"
51 #include "carbondioxide.h"
52 #include "methane.h"
53
54 #ifndef ASC_EXPORT
55 # error "Where is ASC_EXPORT?"
56 #endif
57
58
59 /*------------------------------------------------------------------------------
60 FORWARD DECLARATIONS
61 */
62
63 ExtBBoxInitFunc helmholtz_prepare;
64 ExtBBoxFunc helmholtz_p_calc;
65 ExtBBoxFunc helmholtz_u_calc;
66 ExtBBoxFunc helmholtz_s_calc;
67 ExtBBoxFunc helmholtz_h_calc;
68 ExtBBoxFunc helmholtz_a_calc;
69 ExtBBoxFunc helmholtz_g_calc;
70 ExtBBoxFunc helmholtz_phase_calc;
71
72 /*------------------------------------------------------------------------------
73 GLOBALS
74 */
75
76 /* place to store symbols needed for accessing ASCEND's instance tree */
77 static symchar *helmholtz_symbols[1];
78 #define COMPONENT_SYM helmholtz_symbols[0]
79
80 static const char *helmholtz_p_help = "Calculate pressure from temperature and density, using Helmholtz fundamental correlation";
81 static const char *helmholtz_u_help = "Calculate specific internal energy from temperature and density, using Helmholtz fundamental correlation";
82 static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
83 static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
84 static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";
85 static const char *helmholtz_g_help = "Calculate specific Gibbs energy from temperature and density, using Helmholtz fundamental correlation";
86
87 static const char *helmholtz_phase_help = "Calculate Maxwell phase criterion residuals using Helmholtz fundamental correlation";
88
89 /*------------------------------------------------------------------------------
90 REGISTRATION FUNCTION
91 */
92
93 /**
94 This is the function called from "IMPORT helmholtz"
95
96 It sets up the functions contained in this external library
97 */
98 extern
99 ASC_EXPORT int helmholtz_register(){
100 int result = 0;
101
102 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
103
104 #define CALCFN(NAME,INPUTS,OUTPUTS) \
105 result += CreateUserFunctionBlackBox(#NAME \
106 , helmholtz_prepare \
107 , NAME##_calc /* value */ \
108 , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
109 , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
110 , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
111 , INPUTS,OUTPUTS /* inputs, outputs */ \
112 , NAME##_help /* help text */ \
113 , 0.0 \
114 ) /* returns 0 on success */
115
116 #define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
117 result += CreateUserFunctionBlackBox(#NAME \
118 , helmholtz_prepare \
119 , NAME##_calc /* value */ \
120 , NAME##_calc /* derivatives */ \
121 , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
122 , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
123 , INPUTS,OUTPUTS /* inputs, outputs */ \
124 , NAME##_help /* help text */ \
125 , 0.0 \
126 ) /* returns 0 on success */
127
128 CALCFNDERIV(helmholtz_p,2,1);
129 CALCFN(helmholtz_u,2,1);
130 CALCFN(helmholtz_s,2,1);
131 CALCFN(helmholtz_h,2,1);
132 CALCFN(helmholtz_a,2,1);
133 CALCFN(helmholtz_g,2,1);
134 CALCFN(helmholtz_phase,4,3);
135
136 #undef CALCFN
137
138 if(result){
139 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
140 }
141 return result;
142 }
143
144 /**
145 'helmholtz_prepare' just gets the data member and checks that it's
146 valid, and stores it in the blackbox data field.
147 */
148 int helmholtz_prepare(struct BBoxInterp *bbox,
149 struct Instance *data,
150 struct gl_list_t *arglist
151 ){
152 struct Instance *compinst;
153 const char *comp;
154
155 helmholtz_symbols[0] = AddSymbol("component");
156
157 /* get the data file name (we will look for this file in the ASCENDLIBRARY path) */
158 compinst = ChildByChar(data,COMPONENT_SYM);
159 if(!compinst){
160 ERROR_REPORTER_HERE(ASC_USER_ERROR
161 ,"Couldn't locate 'component' in DATA, please check usage of HELMHOLTZ."
162 );
163 return 1;
164 }
165 if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
166 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
167 return 1;
168 }
169 comp = SCP(SYMC_INST(compinst)->value);
170 CONSOLE_DEBUG("COMPONENT: %s",comp);
171 if(comp==NULL || strlen(comp)==0){
172 ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
173 return 1;
174 }
175
176 if(strcmp(comp,"ammonia")==0){
177 bbox->user_data = (void*)&helmholtz_data_ammonia;
178 }else if(strcmp(comp,"nitrogen")==0){
179 bbox->user_data = (void*)&helmholtz_data_nitrogen;
180 }else if(strcmp(comp,"hydrogen")==0){
181 bbox->user_data = (void*)&helmholtz_data_hydrogen;
182 }else if(strcmp(comp,"water")==0){
183 bbox->user_data = (void*)&helmholtz_data_water;
184 }else if(strcmp(comp,"carbondioxide")==0){
185 bbox->user_data = (void*)&helmholtz_data_carbondioxide;
186 }else if(strcmp(comp,"methane")==0){
187 bbox->user_data = (void*)&helmholtz_data_methane;
188 }else{
189 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
190 }
191
192 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
193 return 0;
194 }
195
196 /*------------------------------------------------------------------------------
197 EVALULATION ROUTINES
198 */
199
200 #define CALCPREPARE(NIN,NOUT) \
201 /* a few checks about the input requirements */ \
202 if(ninputs != NIN)return -1; \
203 if(noutputs != NOUT)return -2; \
204 if(inputs==NULL)return -3; \
205 if(outputs==NULL)return -4; \
206 if(bbox==NULL)return -5; \
207 \
208 /* the 'user_data' in the black box object will contain the */\
209 /* coefficients required for this fluid; cast it to the required form: */\
210 HelmholtzData *helmholtz_data = (HelmholtzData *)bbox->user_data
211
212 /**
213 Evaluation function for 'helmholtz_p'.
214 @param inputs array with values of inputs T and rho.
215 @param outputs array with just value of p
216 @param jacobian, the partial derivative df/dx, where
217 each row is df[i]/dx[j] over each j for the y_out[i] of
218 matching index. The jacobian array is 1-D, row major, i.e.
219 df[i]/dx[j] -> jacobian[i*ninputs+j].
220 @return 0 on success
221 */
222 int helmholtz_p_calc(struct BBoxInterp *bbox,
223 int ninputs, int noutputs,
224 double *inputs, double *outputs,
225 double *jacobian
226 ){
227 CALCPREPARE(2,1);
228
229 /* first input is temperature, second is molar density */
230 if(bbox->task == bb_func_eval){
231 outputs[0] = helmholtz_p(inputs[0], inputs[1], helmholtz_data);
232 }else{
233 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
234 jacobian[0*1+0] = helmholtz_dpdT_rho(inputs[0], inputs[1], helmholtz_data);
235 jacobian[0*1+1] = helmholtz_dpdrho_T(inputs[0], inputs[1], helmholtz_data);
236 }
237
238 /* no need to worry about error states etc. */
239 return 0;
240 }
241
242
243 /**
244 Evaluation function for 'helmholtz_u'
245 @param jacobian ignored
246 @return 0 on success
247 */
248 int helmholtz_u_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 molar density */
256 if(bbox->task == bb_func_eval){
257 outputs[0] = helmholtz_u(inputs[0], inputs[1], helmholtz_data);
258 }else{
259 jacobian[0*1+0] = helmholtz_dudT_rho(inputs[0], inputs[1], helmholtz_data);
260 jacobian[0*1+1] = helmholtz_dudrho_T(inputs[0], inputs[1], helmholtz_data);
261 }
262
263 /* no need to worry about error states etc. */
264 return 0;
265 }
266
267
268 /**
269 Evaluation function for 'helmholtz_s'
270 @param jacobian ignored
271 @return 0 on success
272 */
273 int helmholtz_s_calc(struct BBoxInterp *bbox,
274 int ninputs, int noutputs,
275 double *inputs, double *outputs,
276 double *jacobian
277 ){
278 CALCPREPARE(2,1);
279
280 /* first input is temperature, second is molar density */
281 outputs[0] = helmholtz_s(inputs[0], inputs[1], helmholtz_data);
282
283 /* no need to worry about error states etc. */
284 return 0;
285 }
286
287
288 /**
289 Evaluation function for 'helmholtz_h'
290 @param jacobian ignored
291 @return 0 on success
292 */
293 int helmholtz_h_calc(struct BBoxInterp *bbox,
294 int ninputs, int noutputs,
295 double *inputs, double *outputs,
296 double *jacobian
297 ){
298 CALCPREPARE(2,1);
299
300 /* first input is temperature, second is molar density */
301 if(bbox->task == bb_func_eval){
302 outputs[0] = helmholtz_h(inputs[0], inputs[1], helmholtz_data);
303 }else{
304 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
305 jacobian[0*1+0] = helmholtz_dhdT_rho(inputs[0], inputs[1], helmholtz_data);
306 jacobian[0*1+1] = helmholtz_dhdrho_T(inputs[0], inputs[1], helmholtz_data);
307 }
308
309 /* no need to worry about error states etc. */
310 return 0;
311 }
312
313
314 /**
315 Evaluation function for 'helmholtz_a'
316 @param jacobian ignored
317 @return 0 on success
318 */
319 int helmholtz_a_calc(struct BBoxInterp *bbox,
320 int ninputs, int noutputs,
321 double *inputs, double *outputs,
322 double *jacobian
323 ){
324 CALCPREPARE(2,1);
325
326 /* first input is temperature, second is molar density */
327 outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);
328
329 /* no need to worry about error states etc. */
330 return 0;
331 }
332
333
334 /**
335 Evaluation function for 'helmholtz_g'
336 @param jacobian ignored
337 @return 0 on success
338 */
339 int helmholtz_g_calc(struct BBoxInterp *bbox,
340 int ninputs, int noutputs,
341 double *inputs, double *outputs,
342 double *jacobian
343 ){
344 CALCPREPARE(2,1);
345
346 /* first input is temperature, second is molar density */
347 outputs[0] = helmholtz_g(inputs[0], inputs[1], helmholtz_data);
348
349 /* no need to worry about error states etc. */
350 return 0;
351 }
352
353
354 /**
355 Evaluation function for 'helmholtz_phase'
356 @param jacobian ignored
357 @return 0 on success
358 */
359 int helmholtz_phase_calc(struct BBoxInterp *bbox,
360 int ninputs, int noutputs,
361 double *inputs, double *outputs,
362 double *jacobian
363 ){
364 CALCPREPARE(4,3);
365
366 #define T inputs[0]
367 #define PSAT inputs[1]
368 #define RHOF inputs[2]
369 #define RHOG inputs[3]
370
371 outputs[0] = PSAT - helmholtz_p(T,RHOF,helmholtz_data);
372 outputs[1] = PSAT - helmholtz_p(T,RHOG,helmholtz_data);
373 outputs[2] = helmholtz_g(T,RHOG,helmholtz_data) - helmholtz_g(T,RHOF,helmholtz_data);
374
375 #undef T
376 #undef PSAT
377 #undef RHOF
378 #undef RHOG
379
380 /* we won't worry about error states etc. */
381 return 0;
382 }
383

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