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

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