| 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/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 arguement list does not exist\n"); |
| 202 |
return 1; |
| 203 |
} |
| 204 |
len = gl_length(arglist); |
| 205 |
if (!len) { |
| 206 |
FPRINTF(stderr,"No arguements to external function statement\n"); |
| 207 |
return 1; |
| 208 |
} |
| 209 |
if ((len!=(N_INPUT_ARGS+N_OUTPUT_ARGS))) { |
| 210 |
FPRINTF(stderr,"Number of arguements 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 |
{ |
| 233 |
struct KVALUES_problem *problem; |
| 234 |
if (interp && interp->user_data) { |
| 235 |
problem = (struct KVALUES_problem *)interp->user_data; |
| 236 |
if (problem->components) free((char *)problem->components); |
| 237 |
if (problem->x) free((char *)problem->x); |
| 238 |
if (problem->x) free((char *)problem); |
| 239 |
} |
| 240 |
interp->user_data = NULL; |
| 241 |
} |
| 242 |
|
| 243 |
/** |
| 244 |
This function is one of the registered functions. It operates in |
| 245 |
one modes first_call, and |
| 246 |
creates a KVALUES_problem and calls a number of routines to check |
| 247 |
the arguments (data and arglist) and to cache the information |
| 248 |
processed away in the KVALUES_problem structure. We then attach |
| 249 |
to the hook. |
| 250 |
|
| 251 |
In last_call mode, which should never be seen, |
| 252 |
we delegate to final where we cleanup the problem structure: |
| 253 |
the array of compononent string pointers (we dont own the strings). |
| 254 |
the array of doubles. |
| 255 |
the problem structure itself. |
| 256 |
*/ |
| 257 |
int kvalues_preslv(struct BBoxInterp *interp, |
| 258 |
struct Instance *data, |
| 259 |
struct gl_list_t *arglist |
| 260 |
|
| 261 |
){ |
| 262 |
struct KVALUES_problem *problem; |
| 263 |
int workspace; |
| 264 |
/*struct Instance *arg; |
| 265 |
int nok,c; */ |
| 266 |
|
| 267 |
|
| 268 |
if (interp->task == bb_first_call) { |
| 269 |
if (interp->user_data!=NULL) { /* we have been called before */ |
| 270 |
return 0; /* the problem structure exists */ |
| 271 |
} |
| 272 |
else{ |
| 273 |
problem = CreateProblem(); |
| 274 |
if(CheckArgsOK(interp,data,arglist,problem)) { |
| 275 |
return 1; |
| 276 |
} |
| 277 |
workspace = /* T, x[1..n], P, y[1..n], satp[1..n] */ |
| 278 |
problem->ninputs + problem->noutputs + |
| 279 |
problem->ncomps * 1; |
| 280 |
problem->x = (double *)calloc(workspace, sizeof(double)); |
| 281 |
interp->user_data = (void *)problem; |
| 282 |
return 0; |
| 283 |
} |
| 284 |
} |
| 285 |
if (interp->task == bb_last_call) { |
| 286 |
kvalues_final(interp); |
| 287 |
} |
| 288 |
return 0; |
| 289 |
} |
| 290 |
|
| 291 |
|
| 292 |
/* |
| 293 |
This function provides support to kvalues_fex which is one of the |
| 294 |
registered functions. The input variables are T,x[1..ncomps],P. |
| 295 |
The output variables are y[1..n]. We also load up the x vector |
| 296 |
(our copy) with satP[1..ncomps] as proof of concept that we can |
| 297 |
do (re)initialization. |
| 298 |
*/ |
| 299 |
|
| 300 |
/* |
| 301 |
The name field will be field in vp before the call. |
| 302 |
The rest of the data will be filled the structure |
| 303 |
provided that there are no errors. |
| 304 |
*/ |
| 305 |
static int GetCoefficients(struct Vpdata *vp){ |
| 306 |
if (LookupData(vp)==0) |
| 307 |
return 0; |
| 308 |
else |
| 309 |
return 1; /* error in name lookup */ |
| 310 |
} |
| 311 |
|
| 312 |
|
| 313 |
static int DoCalculation(struct BBoxInterp *interp, |
| 314 |
int ninputs, int noutputs, |
| 315 |
double *inputs, |
| 316 |
double *outputs |
| 317 |
){ |
| 318 |
struct KVALUES_problem *problem; |
| 319 |
int c, offset; |
| 320 |
int ncomps; |
| 321 |
double TotP, T; |
| 322 |
double *liq_frac = NULL; |
| 323 |
double *vap_frac = NULL; |
| 324 |
double *SatP = NULL; |
| 325 |
double *tmp; |
| 326 |
struct Vpdata vp; |
| 327 |
int result = 0; |
| 328 |
|
| 329 |
problem = (struct KVALUES_problem *)interp->user_data; |
| 330 |
ncomps = problem->ncomps; |
| 331 |
assert(ninputs == ncomps+2); |
| 332 |
assert(noutputs == ncomps); |
| 333 |
liq_frac = (double *)calloc(ncomps,sizeof(double)); |
| 334 |
vap_frac = (double *)calloc(ncomps,sizeof(double)); |
| 335 |
SatP = (double *)calloc(ncomps,sizeof(double)); |
| 336 |
|
| 337 |
T = inputs[0]; offset = 1; |
| 338 |
for(c=0;c<ncomps;c++) { |
| 339 |
liq_frac[c] = inputs[c+offset]; |
| 340 |
} |
| 341 |
offset = ncomps+1; |
| 342 |
TotP = inputs[offset]; |
| 343 |
|
| 344 |
for (c=0;c<ncomps;c++) { |
| 345 |
vp.component = problem->components[c]; /* get component name */ |
| 346 |
result = GetCoefficients(&vp); /* get antoines coeffs */ |
| 347 |
SatP[c] = exp(vp.a - vp.b/(T + vp.c)); /* calc satP */ |
| 348 |
vap_frac[c] = SatP[c] * liq_frac[c] / TotP; |
| 349 |
} |
| 350 |
|
| 351 |
#ifdef KVALUES_DEBUG |
| 352 |
FPRINTF(stdout,"Temperature T = %12.8f\n",T); |
| 353 |
FPRINTF(stdout,"TotalPressure P = %12.8f\n",TotP); |
| 354 |
for(c=0;c<ncomps;c++) { |
| 355 |
FPRINTF(stdout,"x[%d] = %12.8g\n",c,liq_frac[c]); |
| 356 |
} |
| 357 |
for (c=0;c<ncomps;c++) { |
| 358 |
FPRINTF(stdout,"y[%d] = %20.8g\n",c,vap_frac[c]); |
| 359 |
} |
| 360 |
#endif /* KVALUES_DEBUG */ |
| 361 |
|
| 362 |
/* |
| 363 |
* Save values, copy the information to the output vector |
| 364 |
* and cleanup. |
| 365 |
*/ |
| 366 |
tmp = problem->x; /* save the input data */ |
| 367 |
*tmp++ = T; |
| 368 |
for (c=0;c<ncomps;c++) { |
| 369 |
*tmp = liq_frac[c]; |
| 370 |
tmp++; |
| 371 |
} |
| 372 |
*tmp++ = TotP; |
| 373 |
|
| 374 |
for (c=0;c<ncomps;c++) { /* save the output data */ |
| 375 |
*tmp = outputs[c] = vap_frac[c]; /* load up the output vector also */ |
| 376 |
tmp++; |
| 377 |
} |
| 378 |
|
| 379 |
for (c=0;c<ncomps;c++) { /* save the internal data */ |
| 380 |
*tmp = SatP[c]; |
| 381 |
tmp++; |
| 382 |
} |
| 383 |
|
| 384 |
free((char *)liq_frac); |
| 385 |
free((char *)vap_frac); |
| 386 |
free((char *)SatP); |
| 387 |
interp->status = calc_all_ok; |
| 388 |
return 0; |
| 389 |
} |
| 390 |
|
| 391 |
int kvalues_fex(struct BBoxInterp *interp, |
| 392 |
int ninputs, int noutputs, |
| 393 |
double *inputs, double *outputs, |
| 394 |
double *jacobian |
| 395 |
){ |
| 396 |
int nok; |
| 397 |
(void)jacobian; |
| 398 |
nok = DoCalculation(interp, ninputs, noutputs, inputs, outputs); |
| 399 |
if (nok) |
| 400 |
return 1; |
| 401 |
else |
| 402 |
return 0; |
| 403 |
} |
| 404 |
|
| 405 |
/** |
| 406 |
Registration function |
| 407 |
*/ |
| 408 |
int kvalues_register(void){ |
| 409 |
int result; |
| 410 |
double epsilon = 1.0e-14; |
| 411 |
|
| 412 |
char kvalues_help[] = "This function does a bubble point calculation\ |
| 413 |
given the mole fractions in the feed, a T and P."; |
| 414 |
|
| 415 |
result = CreateUserFunctionBlackBox("bubblepnt", kvalues_preslv, kvalues_fex, |
| 416 |
NULL,NULL,kvalues_final, |
| 417 |
N_INPUT_ARGS,N_OUTPUT_ARGS, |
| 418 |
kvalues_help,epsilon); |
| 419 |
return result; |
| 420 |
} |
| 421 |
|
| 422 |
#undef N_INPUT_ARGS |
| 423 |
#undef N_OUTPUT_ARGS |
Parent Directory
| 
Revision Log