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

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