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

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