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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 27 - (hide annotations) (download) (as text)
Thu Dec 9 18:22:50 2004 UTC (19 years, 7 months ago) by aw0a
File MIME type: text/x-csrc
File size: 11499 byte(s)
Correcting package checkin
1 aw0a 27 /*********************************************************************\
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