« back to all changes in this revision
Viewing changes to ascend/packages/kvalues.c
- browse files at revision 1996
- compare with another revision
- download diff
- download tarball
- view history from revision 1996
- Committer: jpye
- Date: 2009-04-28 09:02:48 UTC
- Revision ID: svn-v4:f15cd148-d149-4907-926e-73cf92b3ee03:trunk:2012
Migrate base/generic/packages to ascend/packages.
Move base/generic/test to test/base (awaiting further reorg).
Move base/generic/test to test/base (awaiting further reorg).
- files added:
- ascend/packages
- test/base
- files removed:
- base/generic/packages
- base/generic/test
added
removed
ascend/packages/kvalues.c
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
31
#include <compiler/packages.h>
32
#include <compiler/instance_enum.h>
33
34
35
#include <compiler/expr_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
unsigned long 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 BBoxInterp *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
(void)interp;
199
200
if (!arglist) {
201
FPRINTF(stderr,"External function argument list does not exist\n");
202
return 1;
203
}
204
len = gl_length(arglist);
205
if (!len) {
206
FPRINTF(stderr,"No arguments to external function statement\n");
207
return 1;
208
}
209
if ((len!=(N_INPUT_ARGS+N_OUTPUT_ARGS))) {
210
FPRINTF(stderr,"Number of arguments does not match\n");
211
FPRINTF(stderr,"external function prototype\n");
212
return 1;
213
}
214
215
ninputs = CountNumberOfArgs(arglist,1,N_INPUT_ARGS);
216
noutputs = CountNumberOfArgs(arglist,N_INPUT_ARGS+1,
217
N_INPUT_ARGS+N_OUTPUT_ARGS);
218
n_components = GetNumberOfComponents(arglist);
219
220
problem->ninputs = (int)ninputs;
221
problem->noutputs = (int)noutputs;
222
problem->ncomps = n_components;
223
result = GetComponentData(data,problem); /* get the component names */
224
if (result)
225
return 1;
226
227
return 0;
228
}
229
230
231
static void kvalues_final( struct BBoxInterp *interp){
232
struct KVALUES_problem *problem;
233
if (interp && interp->user_data) {
234
problem = (struct KVALUES_problem *)interp->user_data;
235
if (problem->components) free((char *)problem->components);
236
if (problem->x) free((char *)problem->x);
237
if (problem->x) free((char *)problem);
238
}
239
interp->user_data = NULL;
240
}
241
242
/**
243
This function is one of the registered functions. It operates in
244
one modes first_call, and
245
creates a KVALUES_problem and calls a number of routines to check
246
the arguments (data and arglist) and to cache the information
247
processed away in the KVALUES_problem structure. We then attach
248
to the hook.
249
250
In last_call mode, which should never be seen,
251
we delegate to final where we cleanup the problem structure:
252
the array of compononent string pointers (we dont own the strings).
253
the array of doubles.
254
the problem structure itself.
255
*/
256
int kvalues_preslv(struct BBoxInterp *interp,
257
struct Instance *data,
258
struct gl_list_t *arglist
259
260
){
261
struct KVALUES_problem *problem;
262
int workspace;
263
/*struct Instance *arg;
264
int nok,c; */
265
266
267
if (interp->task == bb_first_call) {
268
if (interp->user_data!=NULL) { /* we have been called before */
269
return 0; /* the problem structure exists */
270
}
271
else{
272
problem = CreateProblem();
273
if(CheckArgsOK(interp,data,arglist,problem)) {
274
return 1;
275
}
276
workspace = /* T, x[1..n], P, y[1..n], satp[1..n] */
277
problem->ninputs + problem->noutputs +
278
problem->ncomps * 1;
279
problem->x = (double *)calloc(workspace, sizeof(double));
280
interp->user_data = (void *)problem;
281
return 0;
282
}
283
}
284
if (interp->task == bb_last_call) {
285
kvalues_final(interp);
286
}
287
return 0;
288
}
289
290
291
/*
292
This function provides support to kvalues_fex which is one of the
293
registered functions. The input variables are T,x[1..ncomps],P.
294
The output variables are y[1..n]. We also load up the x vector
295
(our copy) with satP[1..ncomps] as proof of concept that we can
296
do (re)initialization.
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
if (LookupData(vp)==0)
306
return 0;
307
else
308
return 1; /* error in name lookup */
309
}
310
311
312
static int DoCalculation(struct BBoxInterp *interp,
313
int ninputs, int noutputs,
314
double *inputs,
315
double *outputs
316
){
317
struct KVALUES_problem *problem;
318
int c, offset;
319
int ncomps;
320
double TotP, T;
321
double *liq_frac = NULL;
322
double *vap_frac = NULL;
323
double *SatP = NULL;
324
double *tmp;
325
struct Vpdata vp;
326
int result = 0;
327
328
problem = (struct KVALUES_problem *)interp->user_data;
329
ncomps = problem->ncomps;
330
assert(ninputs == ncomps+2);
331
assert(noutputs == ncomps);
332
liq_frac = (double *)calloc(ncomps,sizeof(double));
333
vap_frac = (double *)calloc(ncomps,sizeof(double));
334
SatP = (double *)calloc(ncomps,sizeof(double));
335
336
T = inputs[0]; offset = 1;
337
for(c=0;c<ncomps;c++) {
338
liq_frac[c] = inputs[c+offset];
339
}
340
offset = ncomps+1;
341
TotP = inputs[offset];
342
343
for (c=0;c<ncomps;c++) {
344
vp.component = problem->components[c]; /* get component name */
345
result = GetCoefficients(&vp); /* get antoines coeffs */
346
SatP[c] = exp(vp.a - vp.b/(T + vp.c)); /* calc satP */
347
vap_frac[c] = SatP[c] * liq_frac[c] / TotP;
348
}
349
350
#ifdef KVALUES_DEBUG
351
FPRINTF(stdout,"Temperature T = %12.8f\n",T);
352
FPRINTF(stdout,"TotalPressure P = %12.8f\n",TotP);
353
for(c=0;c<ncomps;c++) {
354
FPRINTF(stdout,"x[%d] = %12.8g\n",c,liq_frac[c]);
355
}
356
for (c=0;c<ncomps;c++) {
357
FPRINTF(stdout,"y[%d] = %20.8g\n",c,vap_frac[c]);
358
}
359
#endif /* KVALUES_DEBUG */
360
361
/*
362
* Save values, copy the information to the output vector
363
* and cleanup.
364
*/
365
tmp = problem->x; /* save the input data */
366
*tmp++ = T;
367
for (c=0;c<ncomps;c++) {
368
*tmp = liq_frac[c];
369
tmp++;
370
}
371
*tmp++ = TotP;
372
373
for (c=0;c<ncomps;c++) { /* save the output data */
374
*tmp = outputs[c] = vap_frac[c]; /* load up the output vector also */
375
tmp++;
376
}
377
378
for (c=0;c<ncomps;c++) { /* save the internal data */
379
*tmp = SatP[c];
380
tmp++;
381
}
382
383
free((char *)liq_frac);
384
free((char *)vap_frac);
385
free((char *)SatP);
386
interp->status = calc_all_ok;
387
return 0;
388
}
389
390
int kvalues_fex(struct BBoxInterp *interp,
391
int ninputs, int noutputs,
392
double *inputs, double *outputs,
393
double *jacobian
394
){
395
int nok;
396
(void)jacobian;
397
nok = DoCalculation(interp, ninputs, noutputs, inputs, outputs);
398
if (nok)
399
return 1;
400
else
401
return 0;
402
}
403
404
/**
405
Registration function
406
*/
407
int kvalues_register(void){
408
int result;
409
double epsilon = 1.0e-14;
410
411
char kvalues_help[] = "This function does a bubble point calculation\
412
given the mole fractions in the feed, a T and P.";
413
414
result = CreateUserFunctionBlackBox("bubblepnt", kvalues_preslv, kvalues_fex,
415
NULL,NULL,kvalues_final,
416
N_INPUT_ARGS,N_OUTPUT_ARGS,
417
kvalues_help,epsilon);
418
return result;
419
}
420
421
#undef N_INPUT_ARGS
422
#undef N_OUTPUT_ARGS
Loggerhead is a web-based interface for Breezy
Version: 2.0.1