/[ascend]/trunk/base/generic/packages/kvalues.c
ViewVC logotype

Contents of /trunk/base/generic/packages/kvalues.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 421 - (show annotations) (download) (as text)
Mon Apr 3 17:26:22 2006 UTC (18 years, 3 months ago) by ben.allan
File MIME type: text/x-csrc
File size: 11856 byte(s)
updated kvalues example to parse and put in model library.
changed import library name to dquote text for more generality.
fixed many, but not all problems in kvalues.c
1 /*********************************************************************\
2 External Distillation Routines.
3 by Kirk Abbott
4 Created: July 4, 1994
5 Version: $Revision: 1.5 $
6 Date last modified: $Date: 1997/07/18 12:20:07 $
7 This file is part of the Ascend Language Interpreter.
8
9 Copyright (C) 1990, 1993, 1994 Thomas Guthrie Epperly, Kirk Abbott.
10
11 The Ascend Language Interpreter is free software; you can redistribute
12 it and/or modify it under the terms of the GNU General Public License as
13 published by the Free Software Foundation; either version 2 of the
14 License, or (at your option) any later version.
15
16 The Ascend Language Interpreter is distributed in hope that it will be
17 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
20
21 You should have received a copy of the GNU General Public License along with
22 the program; if not, write to the Free Software Foundation, Inc., 675
23 Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
24 \*********************************************************************/
25 #include <utilities/ascConfig.h>
26 #include <compiler/compiler.h>
27 #include <compiler/packages.h>
28 #include <compiler/instance_enum.h>
29 #include <compiler/fractions.h>
30 #include <compiler/compiler.h>
31 #include <compiler/dimen.h>
32 #include <compiler/types.h>
33 #include <general/list.h>
34 #include <compiler/sets.h>
35 #include <compiler/atomvalue.h>
36 #include <compiler/parentchild.h>
37 #include <compiler/instquery.h>
38 #include <compiler/instance_name.h>
39 #include <compiler/symtab.h>
40 #include <math.h>
41
42 #define KVALUES_DEBUG 1
43
44 #ifndef EXTERNAL_EPSILON
45 #define EXTERNAL_EPSILON 1.0e-12
46 #endif
47
48 #define N_INPUT_ARGS 3
49 #define N_OUTPUT_ARGS 1
50
51
52 struct Vpdata {
53 char *component;
54 double a, b, c;
55 };
56
57 static struct Vpdata VpTable[] = {
58 "benzene", 15.5381, 2032.73, -33.15,
59 "chloro" , 15.8333, 2477.07, -39.94,
60 "acetone", 15.8737, 2911.32, -56.51,
61 "hexane" , 15.3866, 2697.55, -46.78,
62 "", 0.0, 0.0 , 0.0
63 };
64
65 struct KVALUES_problem {
66 int ncomps;
67 char **components;
68 int ninputs;
69 int noutputs;
70 double *x;
71 };
72
73 static struct KVALUES_problem *CreateProblem(void)
74 {
75 struct KVALUES_problem *result;
76 result = (struct KVALUES_problem *)malloc(sizeof(struct KVALUES_problem));
77 result->ncomps = 0;
78 result->components = NULL;
79 result->ninputs = result->noutputs = 0;
80 result->x = NULL;
81 return result;
82 }
83
84 static int LookupData(struct Vpdata *vp)
85 {
86 char *component;
87 unsigned long c=0L;
88
89 while (strcmp(VpTable[c].component,"")!=0) {
90 if (strcmp(VpTable[c].component,vp->component)==0) {
91 vp->a = VpTable[c].a;
92 vp->b = VpTable[c].b;
93 vp->c = VpTable[c].c;
94 return 0;
95 }
96 c++;
97 }
98 return 1;
99 }
100
101
102 /*
103 * The 2nd. arguement is the array of mole fractions.
104 * the length of this list is then the number of components.
105 */
106 static int GetNumberOfComponents(struct gl_list_t *arglist)
107 {
108 int ncomps;
109 struct gl_list_t *mole_fracs;
110 mole_fracs = (struct gl_list_t *)gl_fetch(arglist,2L);
111 if (mole_fracs) {
112 ncomps = (int)gl_length(mole_fracs);
113 return ncomps;
114 }
115 else{
116 return 0; /* error */
117 }
118 }
119
120 /*
121 *********************************************************************
122 * The following code takes a data instance and tries to decode it
123 * to determine the names of the components used in the model. It will
124 * then create an array of these names and attach it to the problem.
125 * The data is expected to be an instance of a model of the form:
126 *
127 * MODEL kvalues_data;
128 * ncomps IS_A integer;
129 * ncomps := 3;
130 * components[1..ncomps] IS_A symbol;
131 *
132 * components[1] := 'a';
133 * components[2] := 'b';
134 * components[3] := 'c';
135 * END kvalues_data;
136 *
137 *********************************************************************
138 *
139 */
140 static int GetComponentData(struct Instance *data,
141 struct KVALUES_problem *problem)
142 {
143 struct InstanceName name;
144 struct Instance *child;
145 unsigned long pos,c;
146 unsigned long nch;
147 char *component;
148 symchar *s_component = NULL;
149
150 s_component = AddSymbol("components");
151 SetInstanceNameStrPtr(name,s_component);
152 SetInstanceNameType(name,StrName);
153
154 if (!data) {
155 FPRINTF(stderr,"Error: expecting a data instance to be provided\n");
156 return 1;
157 }
158 pos = ChildSearch(data,&name); /* find array head */
159 if (!pos) {
160 FPRINTF(stderr,"Error: cannot find instance \"components\"\n");
161 return 1;
162 }
163 child = InstanceChild(data,pos);
164 if (InstanceKind(child)!= ARRAY_INT_INST) {
165 FPRINTF(stderr,"Error: expecting components to be an array\n");
166 return 1;
167 }
168 nch = NumberChildren(child);
169 if ((nch==0) || (nch!=problem->ncomps)) {
170 FPRINTF(stderr,"Error: Inconsistent component definitions\n");
171 return 1;
172 }
173 if (InstanceKind(InstanceChild(child,1))!=SYMBOL_ATOM_INST) {
174 FPRINTF(stderr,"Error: Expecting an array of symbols\n");
175 return 1;
176 }
177
178 problem->components = (char **)malloc((int)nch*sizeof(char *));
179 for (c=1;c<=nch;c++) {
180 component = (char *)GetSymbolAtomValue(InstanceChild(child,c));
181 problem->components[c-1] = component; /* just peek at it; dont copy it */
182 }
183 return 0;
184 }
185
186
187 static int CheckArgsOK(struct Slv_Interp *slv_interp,
188 struct Instance *data,
189 struct gl_list_t *arglist,
190 struct KVALUES_problem *problem)
191 {
192 unsigned long len,ninputs,noutputs;
193 int n_components;
194 int result;
195
196 if (!arglist) {
197 FPRINTF(stderr,"External function arguement list does not exist\n");
198 return 1;
199 }
200 len = gl_length(arglist);
201 if (!len) {
202 FPRINTF(stderr,"No arguements to external function statement\n");
203 return 1;
204 }
205 if ((len!=(N_INPUT_ARGS+N_OUTPUT_ARGS))) {
206 FPRINTF(stderr,"Number of arguements does not match\n");
207 FPRINTF(stderr,"external function prototype\n");
208 return 1;
209 }
210
211 ninputs = CountNumberOfArgs(arglist,1,N_INPUT_ARGS);
212 noutputs = CountNumberOfArgs(arglist,N_INPUT_ARGS+1,
213 N_INPUT_ARGS+N_OUTPUT_ARGS);
214 n_components = GetNumberOfComponents(arglist);
215
216 problem->ninputs = (int)ninputs;
217 problem->noutputs = (int)noutputs;
218 problem->ncomps = n_components;
219 result = GetComponentData(data,problem); /* get the component names */
220 if (result)
221 return 1;
222
223 return 0;
224 }
225
226
227 /*
228 *********************************************************************
229 * This function is one of the registered functions. It operates in
230 * one of two modes, first_call or last_call. In the former it
231 * creates a KVALUES_problem and calls a number of routines to check
232 * the arguements (data and arglist) and to cache the information
233 * processed away in the KVALUES_problem structure. We then attach
234 * to the hook. *However*, the very first time this function is
235 * called for any given external node, the slv_interp is guaranteed
236 * to be in a clean state. Hence if we find that slv_interp->user_data
237 * is non-NULL, it means that this function has been called before,
238 * regardless of the first_call/last_call flags, and information was
239 * stored there by us. So as not to leak, we either need to deallocate
240 * any storage or just update the information there.
241 *
242 * In last_call mode, the memory associated with the problem needs to
243 * released; this includes:
244 * the array of compononent string pointers (we dont own the strings).
245 * the array of doubles.
246 * the problem structure itself.
247 *
248 *********************************************************************
249 */
250 int kvalues_preslv(struct Slv_Interp *slv_interp,
251 struct Instance *data,
252 struct gl_list_t *arglist)
253 {
254 struct KVALUES_problem *problem;
255 struct Instance *arg;
256 int workspace;
257 int nok,c;
258
259
260 if (slv_interp->first_call) {
261 if (slv_interp->user_data!=NULL) { /* we have been called before */
262 return 0; /* the problem structure exists */
263 }
264 else{
265 problem = CreateProblem();
266 if(CheckArgsOK(slv_interp,data,arglist,problem)) {
267 return 1;
268 }
269 workspace = /* T, x[1..n], P, y[1..n], satp[1..n] */
270 problem->ninputs + problem->noutputs +
271 problem->ncomps * 1;
272 problem->x = (double *)calloc(workspace, sizeof(double));
273 slv_interp->user_data = (void *)problem;
274 return 0;
275 }
276 } else{ /* shutdown */
277 if (slv_interp->user_data) {
278 problem = (struct KVALUES_problem *)slv_interp->user_data;
279 if (problem->components) free((char *)problem->components);
280 if (problem->x) free((char *)problem->x);
281 if (problem->x) free((char *)problem);
282 }
283 slv_interp->user_data = NULL;
284 return 0;
285 }
286 }
287
288
289 /*
290 *********************************************************************
291 * This function provides support to kvalues_fex which is one of the
292 * registered functions. The input variables are T,x[1..ncomps],P.
293 * The output variables are y[1..n]. We also load up the x vector
294 * (our copy) with satP[1..ncomps] as proof of concept that we can
295 * do (re)initialization.
296 *********************************************************************
297 */
298
299 /*
300 * The name field will be field in vp before the call.
301 * The rest of the data will be filled the structure
302 * provided that there are no errors.
303 */
304 static int GetCoefficients(struct Vpdata *vp)
305 {
306 if (LookupData(vp)==0)
307 return 0;
308 else
309 return 1; /* error in name lookup */
310 }
311
312
313 static int DoCalculation(struct Slv_Interp *slv_interp,
314 int ninputs, int noutputs,
315 double *inputs,
316 double *outputs)
317 {
318 struct KVALUES_problem *problem;
319 int c, offset;
320 int ncomps;
321 double TotP, T;
322 double *liq_frac = NULL;
323 double *vap_frac = NULL;
324 double *SatP = NULL;
325 double *tmp;
326 struct Vpdata vp;
327 int result = 0;
328
329 problem = (struct KVALUES_problem *)slv_interp->user_data;
330 ncomps = problem->ncomps;
331 liq_frac = (double *)calloc(ncomps,sizeof(double));
332 vap_frac = (double *)calloc(ncomps,sizeof(double));
333 SatP = (double *)calloc(ncomps,sizeof(double));
334
335 T = inputs[0]; offset = 1;
336 for(c=0;c<ncomps;c++) {
337 liq_frac[c] = inputs[c+offset];
338 }
339 offset = ncomps+1;
340 TotP = inputs[offset];
341
342 for (c=0;c<ncomps;c++) {
343 vp.component = problem->components[c]; /* get component name */
344 result = GetCoefficients(&vp); /* get antoines coeffs */
345 SatP[c] = exp(vp.a - vp.b/(T + vp.c)); /* calc satP */
346 vap_frac[c] = SatP[c] * liq_frac[c] / TotP;
347 }
348
349 #ifdef KVALUES_DEBUG
350 FPRINTF(stdout,"Temperature T = %12.8f\n",T);
351 FPRINTF(stdout,"TotalPressure P = %12.8f\n",TotP);
352 for(c=0;c<ncomps;c++) {
353 FPRINTF(stdout,"x[%d] = %12.8g\n",c,liq_frac[c]);
354 }
355 for (c=0;c<ncomps;c++) {
356 FPRINTF(stdout,"y[%d] = %20.8g\n",c,vap_frac[c]);
357 }
358 #endif /* KVALUES_DEBUG */
359
360 /*
361 * Save values, copy the information to the output vector
362 * and cleanup.
363 */
364 tmp = problem->x; /* save the input data */
365 *tmp++ = T;
366 for (c=0;c<ncomps;c++) {
367 *tmp = liq_frac[c];
368 tmp++;
369 }
370 *tmp++ = TotP;
371
372 for (c=0;c<ncomps;c++) { /* save the output data */
373 *tmp = outputs[c] = vap_frac[c]; /* load up the output vector also */
374 tmp++;
375 }
376
377 for (c=0;c<ncomps;c++) { /* save the internal data */
378 *tmp = SatP[c];
379 tmp++;
380 }
381
382 free((char *)liq_frac);
383 free((char *)vap_frac);
384 free((char *)SatP);
385 slv_interp->status = calc_all_ok;
386 return 0;
387 }
388
389 int kvalues_fex(struct Slv_Interp *slv_interp,
390 int ninputs, int noutputs,
391 double *inputs, double *outputs,
392 double *jacobian)
393 {
394 int nok;
395 nok = DoCalculation(slv_interp, ninputs, noutputs, inputs, outputs);
396 if (nok)
397 return 1;
398 else
399 return 0;
400 }
401
402
403 int KValues_Init (void)
404 {
405 int result;
406
407 char kvalues_help[] = "This function does a bubble point calculation\
408 given the mole fractions in the feed, a T and P.";
409
410 result = CreateUserFunction("bubblepnt",kvalues_preslv, kvalues_fex,
411 NULL,NULL,3,1,kvalues_help);
412 return result;
413 }
414
415 #undef N_INPUT_ARGS
416 #undef N_OUTPUT_ARGS

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