| 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 |
Parent Directory
| 
Revision Log