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


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    
45 jpye 1832 /* for the moment, species data are defined in C code, we'll implement something
46     better later on, hopefully. */
47     #include "ammonia.h"
48 jpye 1849 #include "nitrogen.h"
49 jpye 1995 #include "hydrogen.h"
50     #include "water.h"
51 jpye 2110 #include "carbondioxide.h"
52 kchittur 2147 #include "methane.h"
53 jpye 1832
54 jpye 1822 #ifndef ASC_EXPORT
55     # error "Where is ASC_EXPORT?"
56     #endif
57    
58    
59     /*------------------------------------------------------------------------------
60     FORWARD DECLARATIONS
61     */
62    
63 jpye 1825 ExtBBoxInitFunc helmholtz_prepare;
64 jpye 1822 ExtBBoxFunc helmholtz_p_calc;
65 jpye 1825 ExtBBoxFunc helmholtz_u_calc;
66 jpye 1903 ExtBBoxFunc helmholtz_s_calc;
67 jpye 1825 ExtBBoxFunc helmholtz_h_calc;
68 jpye 1903 ExtBBoxFunc helmholtz_a_calc;
69 jpye 2124 ExtBBoxFunc helmholtz_g_calc;
70     ExtBBoxFunc helmholtz_phase_calc;
71 jpye 1822
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 jpye 1825 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 jpye 1903 static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
83 jpye 1825 static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
84 jpye 1903 static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";
85 jpye 2124 static const char *helmholtz_g_help = "Calculate specific Gibbs energy from temperature and density, using Helmholtz fundamental correlation";
86 jpye 1825
87 jpye 2124 static const char *helmholtz_phase_help = "Calculate Maxwell phase criterion residuals using Helmholtz fundamental correlation";
88    
89 jpye 1822 /*------------------------------------------------------------------------------
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 jpye 1920 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
103 jpye 1822
104 jpye 1825 #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 jpye 1822
116 jpye 1909 #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 jpye 1825 CALCFN(helmholtz_u,2,1);
130 jpye 1903 CALCFN(helmholtz_s,2,1);
131 jpye 1825 CALCFN(helmholtz_h,2,1);
132 jpye 1903 CALCFN(helmholtz_a,2,1);
133 jpye 2124 CALCFN(helmholtz_g,2,1);
134     CALCFN(helmholtz_phase,4,3);
135 jpye 1825
136     #undef CALCFN
137    
138 jpye 1822 if(result){
139     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
140     }
141     return result;
142     }
143    
144 jpye 1825 /**
145     'helmholtz_prepare' just gets the data member and checks that it's
146     valid, and stores it in the blackbox data field.
147 jpye 1822 */
148 jpye 1825 int helmholtz_prepare(struct BBoxInterp *bbox,
149 jpye 1822 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 jpye 1849 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 jpye 1995 }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 jpye 2110 }else if(strcmp(comp,"carbondioxide")==0){
185     bbox->user_data = (void*)&helmholtz_data_carbondioxide;
186 kchittur 2147 }else if(strcmp(comp,"methane")==0){
187     bbox->user_data = (void*)&helmholtz_data_methane;
188 jpye 1849 }else{
189 jpye 1995 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
190 jpye 1822 }
191    
192 jpye 1849 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
193 jpye 1822 return 0;
194     }
195    
196 jpye 1825 /*------------------------------------------------------------------------------
197     EVALULATION ROUTINES
198     */
199    
200 jpye 2124 #define CALCPREPARE(NIN,NOUT) \
201 jpye 1825 /* a few checks about the input requirements */ \
202 jpye 2124 if(ninputs != NIN)return -1; \
203     if(noutputs != NOUT)return -2; \
204 jpye 1825 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 jpye 1822 /**
213 jpye 1825 Evaluation function for 'helmholtz_p'.
214 jpye 1909 @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 jpye 1822 @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 jpye 2124 CALCPREPARE(2,1);
228 jpye 1822
229     /* first input is temperature, second is molar density */
230 jpye 1909 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 jpye 1822
238     /* no need to worry about error states etc. */
239 jpye 1825 return 0;
240     }
241 jpye 1822
242 jpye 1825
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 jpye 2124 CALCPREPARE(2,1);
254 jpye 1825
255     /* first input is temperature, second is molar density */
256 jpye 1923 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 jpye 1825
263     /* no need to worry about error states etc. */
264 jpye 1822 return 0;
265     }
266 jpye 1825
267    
268     /**
269 jpye 2124 Evaluation function for 'helmholtz_s'
270 jpye 1825 @param jacobian ignored
271     @return 0 on success
272     */
273 jpye 1903 int helmholtz_s_calc(struct BBoxInterp *bbox,
274     int ninputs, int noutputs,
275     double *inputs, double *outputs,
276     double *jacobian
277     ){
278 jpye 2124 CALCPREPARE(2,1);
279 jpye 1903
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 jpye 1825 int helmholtz_h_calc(struct BBoxInterp *bbox,
294     int ninputs, int noutputs,
295     double *inputs, double *outputs,
296     double *jacobian
297     ){
298 jpye 2124 CALCPREPARE(2,1);
299 jpye 1825
300     /* first input is temperature, second is molar density */
301 jpye 1920 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 jpye 1825
309     /* no need to worry about error states etc. */
310     return 0;
311     }
312    
313    
314 jpye 1903 /**
315 jpye 2124 Evaluation function for 'helmholtz_a'
316 jpye 1903 @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 jpye 2124 CALCPREPARE(2,1);
325 jpye 1825
326 jpye 1903 /* first input is temperature, second is molar density */
327     outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);
328 jpye 1825
329 jpye 1903 /* no need to worry about error states etc. */
330     return 0;
331     }
332    
333    
334 jpye 2124 /**
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 jpye 1903
346 jpye 2124 /* 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