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

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