/[ascend]/branches/jacob/models/johnpye/fprops/asc_mixture.c
ViewVC logotype

Contents of /branches/jacob/models/johnpye/fprops/asc_mixture.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3079 - (show annotations) (download) (as text)
Thu Aug 20 02:17:28 2015 UTC (3 years, 3 months ago) by jacob
File MIME type: text/x-csrc
File size: 22071 byte(s)
Review format of comments used to create doxygen documentaion.  Begin ammending erroneous formatting in previously written doxygen comments in the mixture code.

Added further commentary and documentation to .a4c test files

1 /* ASCEND modelling environment
2 Copyright (C) Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful, but
10 WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation --
17
18 Free Software Foundation, Inc.
19 59 Temple Place - Suite 330
20 Boston, MA 02111-1307, USA.
21 *//*
22 by Jacob Shealy, June 25-August 14, 2015
23
24 Exports functions for initial model of ideal-solution mixing
25 */
26
27 #include "asc_mixture.h"
28
29 /*
30 Register all functions that will be exported
31 */
32 extern ASC_EXPORT int mixture_register(){
33 int result = 0;
34
35 ERROR_REPORTER_HERE(ASC_USER_WARNING,
36 "FPROPS in general, and IN PARTICULAR this mixture module, "
37 "are still EXPERIMENTAL. Use with caution.");
38
39 #define CALCFN(NAME,INPUTS,OUTPUTS) \
40 result += CreateUserFunctionBlackBox( #NAME \
41 , asc_mixture_prepare /* function to initialize */ \
42 , NAME##_calc /* function to find value */ \
43 , (ExtBBoxFunc *)NULL /* derivatives -- none */ \
44 , (ExtBBoxFunc *)NULL /* hessian -- none */ \
45 , (ExtBBoxFinalFunc *)NULL /* finalization -- none */ \
46 , INPUTS,OUTPUTS \
47 , NAME##_help /* help text/documentation */ \
48 , 0.0 \
49 )
50
51 /* CALCFN(mixture_p,2,1); */
52 CALCFN(mixture_rho,2,1);
53 CALCFN(mixture_phase_rho,3,1);
54 CALCFN(mixture_comps_rho,4,1);
55 CALCFN(mixture_u,2,1);
56 CALCFN(mixture_phase_u,3,1);
57 CALCFN(mixture_comps_u,4,1);
58 CALCFN(mixture_h,2,1);
59 CALCFN(mixture_phase_h,3,1);
60 CALCFN(mixture_comps_h,4,1);
61 CALCFN(mixture_cp,2,1);
62 CALCFN(mixture_phase_cp,3,1);
63 CALCFN(mixture_comps_cp,4,1);
64 CALCFN(mixture_cv,2,1);
65 CALCFN(mixture_phase_cv,3,1);
66 CALCFN(mixture_comps_cv,4,1);
67
68 CALCFN(mixture_s,2,1);
69 CALCFN(mixture_phase_s,3,1);
70 CALCFN(mixture_comps_s,4,1);
71 CALCFN(mixture_g,2,1);
72 CALCFN(mixture_phase_g,3,1);
73 CALCFN(mixture_comps_g,4,1);
74 CALCFN(mixture_a,2,1);
75 CALCFN(mixture_phase_a,3,1);
76 CALCFN(mixture_comps_a,4,1);
77
78 CALCFN(mixture_count_phases,2,4);
79 CALCFN(mixture_count_components,3,1);
80 CALCFN(mixture_component_frac,4,1);
81 CALCFN(mixture_component_cnum,4,1);
82 CALCFN(mixture_dew_p,1,1);
83 CALCFN(mixture_bubble_p,1,1);
84 CALCFN(mixture_dew_T,1,1);
85 CALCFN(mixture_bubble_T,1,1);
86 CALCFN(mixture_state_T_ph,2,1);
87
88 if(result){
89 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunctionBlackBox result is %d\n",result);
90 }
91 return result;
92 }
93
94 /*
95 Function to prepare persistent data, namely the MixtureSpec struct which
96 will be used in modeling the mixture.
97 */
98 int asc_mixture_prepare(struct BBoxInterp *bbox, struct Instance *data, struct gl_list_t *arglist){
99
100 #define CHECK_TYPE(VAR,TYPE,NAME,TYPENAME) \
101 if(InstanceKind(VAR)!=TYPE){ \
102 ERROR_REPORTER_HERE(ASC_USER_ERROR \
103 , "DATA member '%s' has type-value %#o, but must have %s (value %#o)" \
104 , NAME, InstanceKind(VAR), TYPENAME, TYPE); \
105 return 1; \
106 }else{ \
107 /* ERROR_REPORTER_HERE(ASC_USER_NOTE, "DATA member %s has correct type %s" */ \
108 /* , NAME, TYPENAME); */ \
109 }
110
111 #define CHECK_EXIST_TYPE(VAR,TYPE,NAME,TYPENAME) \
112 if(! VAR){ \
113 ERROR_REPORTER_HERE(ASC_USER_ERROR \
114 , "Couldn't locate '%s' in DATA, please check usage", NAME); \
115 } \
116 CHECK_TYPE(VAR,TYPE,NAME,TYPENAME)
117
118 unsigned i;
119 struct Instance *npureinst, *compinst, *xinst, *typeinst, *srcinst;
120 unsigned npure;
121
122 mix_symbols[0] = AddSymbol("npure");
123 mix_symbols[1] = AddSymbol("components");
124 mix_symbols[2] = AddSymbol("xs");
125 mix_symbols[3] = AddSymbol("eos");
126 mix_symbols[4] = AddSymbol("source");
127
128 /* check existence of 'data' Instance */
129 if(! data){
130 ERROR_REPORTER_HERE(ASC_USER_ERROR, "Couldn't locate 'data', please check usage");
131 }
132
133 /* number of pure components -- required */
134 npureinst = ChildByChar(data, mix_symbols[NPURE_SYM]);
135 CHECK_EXIST_TYPE(npureinst, INTEGER_CONSTANT_INST, "npure", "'integer constant'");
136 npure = (int *)(IC_INST(npureinst)->value);
137 /* ERROR_REPORTER_HERE(ASC_USER_NOTE, "Number of pures is %i", npure); */
138
139 const struct gl_list_t *comps_gl;
140 const struct gl_list_t *xs_gl;
141 /* const */ char *type = NULL
142 , **comps = ASC_NEW_ARRAY(/* const */ char *, npure)
143 , **srcs = ASC_NEW_ARRAY(/* const */ char *, npure);
144 double *xs = ASC_NEW_ARRAY(double, npure);
145
146 /* Component names -- required */
147 compinst = ChildByChar(data, mix_symbols[COMP_SYM]);
148 CHECK_EXIST_TYPE(compinst, ARRAY_INT_INST, "components"
149 , "array indexed with integers");
150 comps_gl = ARY_INST(compinst)->children;
151
152 /*
153 Component correlation type (equation of state) -- not required, so check
154 for existence before checking type and assigning to the char array 'type'.
155 */
156 typeinst = ChildByChar(data, mix_symbols[TYPE_SYM]);
157 if(typeinst){
158 CHECK_TYPE(typeinst, SYMBOL_CONSTANT_INST, "type", "'symbol constant'");
159 type = SCP(SYMC_INST(typeinst)->value); /* read 'typeinst' into a string */
160 if(type && strlen(type)==0){
161 char t[] = "pengrob";
162 type = t;
163 }
164 }else{
165 char t[] = "pengrob";
166 type = t;
167 }
168
169 /*
170 Data string representing source -- not required, so check for existence
171 before checking type and assigning to the char array 'source'.
172 */
173 srcinst = ChildByChar(data, mix_symbols[SOURCE_SYM]);
174 if(srcinst){
175 CHECK_TYPE(srcinst, SYMBOL_CONSTANT_INST, "source", "'symbol constant'");
176 srcs[0] = SCP(SYMC_INST(srcinst)->value); /* read 'srcinst' into a string */
177 if(srcs[0] && strlen(srcs[0])==0){
178 srcs[0] = NULL;
179 }else if(!srcs[0]){
180 srcs[0] = NULL;
181 }
182 }else{
183 srcs[0] = NULL;
184 }
185 for(i=1;i<npure;i++){
186 srcs[i] = srcs[0];
187 }
188
189 /* Mass fractions -- required */
190 xinst = ChildByChar(data, mix_symbols[X_SYM]);
191 CHECK_EXIST_TYPE(xinst, ARRAY_INT_INST, "xinst"
192 , "array indexed with integers");
193 xs_gl = ARY_INST(xinst)->children;
194
195 /*
196 Check that the lengths of the arrays 'comps_gl' and 'xs_gl' are equal,
197 and equal to 'npure'
198 */
199 if(xs_gl->length!=npure){
200 if(comps_gl->length==xs_gl->length){
201 ERROR_REPORTER_HERE(ASC_USER_ERROR
202 , "The components and mass fractions arrays both differ in length"
203 "\n\tfrom the given number of components 'npure', but are equal in"
204 "\n\tlength to each other. Setting npure = length of the arrays...");
205 npure = xs_gl->length;
206 }else if(comps_gl->length!=npure){
207 ERROR_REPORTER_HERE(ASC_USER_ERROR
208 , "The components and mass fractions arrays both differ in length"
209 "\n\tfrom the given number of components 'npure', and are not equal"
210 "\n\tin length to each other. Setting npure = (smallest length)...");
211 double lens[] = {npure, xs_gl->length, comps_gl->length};
212 npure = min_element(3, lens);
213 }else{
214 ERROR_REPORTER_HERE(ASC_USER_ERROR
215 , "The mass fractions array differs in length from the given number"
216 "\n\tof components 'npure' and the length of the components array."
217 "\n\tSetting npure = (smallest length)...");
218 double lens[] = {npure, xs_gl->length};
219 npure = min_element(2, lens);
220 }
221 }else if(comps_gl->length!=npure){
222 ERROR_REPORTER_HERE(ASC_USER_ERROR
223 , "The components array differs in length from the given number of"
224 "\n\tcomponents 'npure' and the length of the mass fractions array."
225 "\n\tSetting npure = (smallest length)...");
226 double lens[] = {npure, xs_gl->length};
227 npure = min_element(2, lens);
228 }
229
230 /* Read contents of 'comps_gl' and 'xs_gl' into 'comps_child' and 'xs_child' */
231 for(i=0;i<npure;i++){
232 comps[i] = SCP(SYMC_INST(InstanceChild(compinst, i+1))->value);
233 xs[i] = RC_INST(InstanceChild(xinst, i+1))->value;
234 }
235
236 /* Create mixture specification in a MixtureSpec struct */
237 MixtureError merr = MIXTURE_NO_ERROR;
238 bbox->user_data = (void *) build_MixtureSpec(npure, xs, (void **)comps, type, srcs, &merr);
239 ERROR_REPORTER_HERE(ASC_USER_NOTE, "The equation of state is %s", type);
240
241 return 0;
242 }
243
244 /* ---------------------------------------------------------------------
245 Functions which calculate results
246 */
247 #if 0
248 int mixture_p_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
249 double *inputs, double *outputs, double *jacobian){
250 CALCPREP(1,1);
251
252 ERROR_REPORTER_HERE(ASC_USER_NOTE
253 , "Function 'mixture_p_calc' -- this function has no contents as yet.");
254
255 outputs[0] = 0.0;
256
257 return 0;
258 }
259 #endif
260
261 /* ---------------------------------------------------------------------
262 Mixture property-calculation functions
263 */
264 /*
265 Find and return the overall mixture density
266 */
267 int mixture_rho_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
268 double *inputs, double *outputs, double *jacobian){
269 #define NPHASE PS->phases
270 #define PTYPE PS->ph_type[i]
271 CALCPREP(2,1);
272 CALCFLASH;
273 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS);
274
275 double rhos[MS->pures];
276 outputs[0] = mixture_rho(PM, rhos);
277 #if 0
278 unsigned i;
279 ERROR_REPORTER_HERE(ASC_USER_NOTE, "The overall mixture density is %g", outputs[0]);
280 for(i=0;i<NPHASE;i++){
281 ERROR_REPORTER_HERE(ASC_USER_NOTE, "\tThe density of %s phase is %g kg/m3"
282 , MIX_PHASE_STRING(PTYPE), rhos[i]);
283 }
284 #endif
285 return 0;
286 #undef PTYPE
287 #undef NPHASE
288 }
289
290 /*
291 Find and return the mixture density for a single phase
292 */
293 int mixture_phase_rho_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
294 double *inputs, double *outputs, double *jacobian){
295 CALCPREP(3,1);
296 CALCFLASH;
297 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS);
298
299 unsigned ph_num = ((unsigned) inputs[2]) - 1; /* number of the phase */
300 double rhos[PS->phases]; /* individual phase densities */
301
302 mixture_rho(PM, rhos);
303 if(ph_num < PS->phases){
304 outputs[0] = rhos[ph_num]; /* assign density of one phase to output */
305 }else{
306 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u", ph_num+1);
307 outputs[0] = 0.0;
308 }
309 return 0;
310 }
311
312 int mixture_comps_rho_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
313 double *inputs, double *outputs, double *jacobian){
314 CALCPREP(4,1);
315 CALCFLASH;
316
317 unsigned ph_num = ((int) inputs[2]) - 1 /* number of the phase */
318 , comp_num = ((int) inputs[3]) - 1; /* number of the component */
319
320 if(ph_num < PS->phases){
321 if(comp_num < PS->PH[ph_num]->pures){
322 outputs[0] = PS->PH[ph_num]->rhos[comp_num];
323 }else{
324 ERROR_REPORTER_HERE(ASC_USER_ERROR
325 , "There is no component number %u in phase %u."
326 , comp_num+1, ph_num+1);
327 outputs[0] = 0.0;
328 }
329 }else{
330 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u.", ph_num+1);
331 outputs[0] = 0.0;
332 }
333 return 0;
334 }
335
336 /*
337 Macros used to generate mixture-property functions automatically
338 */
339 #define MIX_PROP_EXTFUNC(PROP) \
340 int mixture_##PROP##_calc(struct BBoxInterp *bbox, int ninputs, int noutputs, \
341 double *inputs, double *outputs, double *jacobian){ \
342 CALCPREP(2,1); \
343 CALCFLASH; \
344 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS); \
345 double props[PS->phases]; /* internal by-phase property values */ \
346 outputs[0] = mixture_##PROP(PM, props, &err); \
347 return 0; \
348 }
349
350 #define MIX_PHASE_EXTFUNC(PROP) \
351 int mixture_phase_##PROP##_calc(struct BBoxInterp *bbox, int ninputs, int noutputs, \
352 double *inputs, double *outputs, double *jacobian){ \
353 CALCPREP(3,1); \
354 CALCFLASH; \
355 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS); \
356 unsigned ph_num = ((unsigned) inputs[2]) - 1; /* number of the phase */ \
357 double props[PS->phases]; /* internal by-phase property values */ \
358 mixture_##PROP(PM, props, &err); \
359 if(ph_num < PS->phases){ \
360 outputs[0] = props[ph_num]; /* assign property of one phase to output */ \
361 }else{ \
362 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u.", ph_num+1); \
363 outputs[0] = 0.0; \
364 } \
365 return 0; \
366 }
367
368 #define RHO PS->PH[ph_num]->rhos[comp_num]
369 #define PPF PS->PH[ph_num]->PF[comp_num]
370
371 #define MIX_COMPS_EXTFUNC(PROP) \
372 int mixture_comps_##PROP##_calc(struct BBoxInterp *bbox, int ninputs, int noutputs, \
373 double *inputs, double *outputs, double *jacobian){ \
374 CALCPREP(4,1); \
375 CALCFLASH; \
376 unsigned ph_num = ((int) inputs[2]) - 1 /* number of the phase */ \
377 , comp_num = ((int) inputs[3]) - 1; /* number of the component */ \
378 /* Access the property function from FPROPS for an individual component */ \
379 if(ph_num < PS->phases){ \
380 if(comp_num < PS->PH[ph_num]->pures){ \
381 outputs[0] = fprops_##PROP((FluidState){T, RHO, PPF}, &err); \
382 }else{ \
383 ERROR_REPORTER_HERE(ASC_USER_ERROR \
384 , "There is no component number %u in phase %u." \
385 , comp_num+1, ph_num+1); \
386 outputs[0] = 0.0; \
387 } \
388 }else{ \
389 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u.", ph_num+1); \
390 outputs[0] = 0.0; \
391 } \
392 return 0; \
393 }
394
395 MIX_PROP_EXTFUNC(u); MIX_PROP_EXTFUNC(h); MIX_PROP_EXTFUNC(cp);
396 MIX_PROP_EXTFUNC(cv); MIX_PROP_EXTFUNC(s); MIX_PROP_EXTFUNC(g);
397 MIX_PROP_EXTFUNC(a);
398
399 MIX_PHASE_EXTFUNC(u); MIX_PHASE_EXTFUNC(h); MIX_PHASE_EXTFUNC(cp);
400 MIX_PHASE_EXTFUNC(cv); MIX_PHASE_EXTFUNC(s); MIX_PHASE_EXTFUNC(g);
401 MIX_PHASE_EXTFUNC(a);
402
403 MIX_COMPS_EXTFUNC(u); MIX_COMPS_EXTFUNC(h); MIX_COMPS_EXTFUNC(cp);
404 MIX_COMPS_EXTFUNC(cv); MIX_COMPS_EXTFUNC(s); MIX_COMPS_EXTFUNC(g);
405 MIX_COMPS_EXTFUNC(a);
406 #undef PPF
407 #undef RHOS
408
409 /* ---------------------------------------------------------------------
410 Phase-equilibrium functions
411 */
412 /*
413 Returns the number of phases and mass fraction for each phase
414 (supercritical, vapor, and liquid). Mass fraction is zero if the phase is
415 not present.
416 */
417 int mixture_count_phases_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
418 double *inputs, double *outputs, double *jacobian){
419 CALCPREP(2,4);
420 CALCFLASH;
421
422 unsigned i;
423 for(i=1;i<4;i++){
424 /* Set all outputs to zero initially */
425 outputs[i] = 0.0;
426 }
427 outputs[0] = (double) PS->phases;
428
429 ERROR_REPORTER_HERE(ASC_USER_NOTE, "There are %u phases", PS->phases);
430 for(i=0;i<PS->phases;i++){
431 /* For all phases present, set the corresponding phase-fraction output
432 equal to the mass fraction of the phase. */
433 outputs[i+1] = PS->ph_frac[i];
434
435 ERROR_REPORTER_HERE(ASC_USER_NOTE, "\tPhase number %u is %s", i+1,
436 MIX_PHASE_STRING(PS->ph_type[i]));
437 }
438 return 0;
439 }
440
441 /*
442 Find and return the number of components of some phase within the mixture.
443 */
444 int mixture_count_components_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
445 double *inputs, double *outputs, double *jacobian){
446 CALCPREP(3,1);
447 CALCFLASH;
448
449 unsigned ph_num = ((unsigned) inputs[2]) - 1; /* the number of the phase */
450 if(ph_num < PS->phases){
451 outputs[0] = (double) PS->PH[ph_num]->pures; /* the number of pures */
452 }else{
453 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u", ph_num+1);
454 outputs[0] = 0.0;
455 }
456 return 0;
457 }
458
459 /*
460 Find and return the mass fraction of some j_th component within the i_th
461 phase of the mixture.
462 */
463 int mixture_component_frac_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
464 double *inputs, double *outputs, double *jacobian){
465 CALCPREP(4,1);
466 CALCFLASH;
467
468 unsigned i = ((unsigned) inputs[2]) - 1
469 , j = ((unsigned) inputs[3]) - 1;
470 if(i < PS->phases){
471 if(j < PS->PH[i]->pures){
472 outputs[0] = (double) PS->PH[i]->xs[j];
473 }else{
474 ERROR_REPORTER_HERE(ASC_USER_ERROR
475 , "There is no component number %u in phase %u", j+1, i+1);
476 outputs[0] = 0.0;
477 }
478 }else{
479 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u", i+1);
480 outputs[0] = 0.0;
481 }
482 return 0;
483 }
484
485 /*
486 Find and return the value of the j_th index in member 'c' of some i_th phase 'PH'
487 within the mixture (that is, 'PS->PH[i]->c[j]').
488 */
489 int mixture_component_cnum_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
490 double *inputs, double *outputs, double *jacobian){
491 CALCPREP(4,1);
492 CALCFLASH;
493
494 unsigned i = ((unsigned) inputs[2]) - 1
495 , j = ((unsigned) inputs[3]) - 1;
496 if(i < PS->phases){
497 if(j < PS->PH[i]->pures){
498 outputs[0] = (double) PS->PH[i]->c[j];
499 }else{
500 ERROR_REPORTER_HERE(ASC_USER_ERROR
501 , "There is no component number %u in phase %u", j+1, i+1);
502 outputs[0] = 0.0;
503 }
504 }else{
505 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u", i+1);
506 outputs[0] = 0.0;
507 }
508 return 0;
509 }
510
511 /*
512 Find and return the mixture dew pressure at some temperature.
513 */
514 int mixture_dew_p_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
515 double *inputs, double *outputs, double *jacobian){
516 CALCPREP(1,1);
517
518 int result;
519 double T = inputs[0] /* the mixture temperature */
520 , p_d /* the mixture dew pressure (if there is any) */
521 , tol = MIX_XTOL /* tolerance to use in solving for dew pressure */
522 ;
523
524 result = mixture_dew_pressure(&p_d, MS, T, tol, &err);
525 if(result){
526 switch(result){
527 case 1:
528 ERROR_REPORTER_HERE(ASC_USER_ERROR
529 , "The dew pressure converged on a non-solution point.");
530 break;
531 case 2:
532 ERROR_REPORTER_HERE(ASC_USER_ERROR
533 , "The dew pressure converged on Infinity or NaN.");
534 break;
535 case 3:
536 ERROR_REPORTER_HERE(ASC_USER_ERROR
537 , "The root-finding algorithm that searches for the dew pressure "
538 "\nfailed to converge in the maximum number of iterations.");
539 break;
540 case 4:
541 ERROR_REPORTER_HERE(ASC_USER_ERROR
542 , "There is no dew pressure; all components are supercritical at "
543 "\ntemperature %g K.", T);
544 break;
545 }
546 outputs[0] = MIN_P;
547 return 1;
548 }else{
549 outputs[0] = p_d;
550 return 0;
551 }
552 }
553
554 /*
555 Find and return the mixture bubble pressure at some temperature.
556 */
557 int mixture_bubble_p_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
558 double *inputs, double *outputs, double *jacobian){
559 CALCPREP(1,1);
560
561 int result;
562 double T = inputs[0] /* the mixture temperature */
563 , p_b /* the mixture bubble pressure (if there is any) */
564 , tol = MIX_XTOL /* tolerance to use in solving for bubble pressure */
565 ;
566
567 result = mixture_bubble_pressure(&p_b, MS, T, tol, &err);
568 if(result){
569 switch(result){
570 case 1:
571 ERROR_REPORTER_HERE(ASC_USER_ERROR
572 , "The bubble pressure converged on a non-solution point.");
573 break;
574 case 2:
575 ERROR_REPORTER_HERE(ASC_USER_ERROR
576 , "The bubble pressure converged on Infinity or NaN.");
577 break;
578 case 3:
579 ERROR_REPORTER_HERE(ASC_USER_ERROR
580 , "The root-finding algorithm that searches for the bubble pressure "
581 "\nfailed to converge in the maximum number of iterations.");
582 break;
583 case 4:
584 ERROR_REPORTER_HERE(ASC_USER_ERROR
585 , "There is no bubble pressure; all components are supercritical at "
586 "\ntemperature %g K.", T);
587 break;
588 }
589 outputs[0] = MIN_P;
590 return 1;
591 }else{
592 outputs[0] = p_b;
593 return 0;
594 }
595 }
596
597 /*
598 Find and return the mixture dew temperature at some pressure.
599 */
600 int mixture_dew_T_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
601 double *inputs, double *outputs, double *jacobian){
602 CALCPREP(1,1);
603
604 int result;
605 double p = inputs[0] /* the mixture pressure */
606 , T_d /* the mixture dew temperature (if there is any) */
607 , tol = MIX_XTOL /* tolerance to use in solving for dew temperature */
608 ;
609
610 result = mixture_dew_temperature(&T_d, MS, p, tol, &err);
611 if(result){
612 switch(result){
613 case 1:
614 ERROR_REPORTER_HERE(ASC_USER_ERROR
615 , "The dew temperature converged on a non-solution point");
616 break;
617 case 2:
618 ERROR_REPORTER_HERE(ASC_USER_ERROR
619 , "The dew temperature converged on Infinity or NaN.");
620 break;
621 case 3:
622 ERROR_REPORTER_HERE(ASC_USER_ERROR
623 , "The root-finding algorithm that searches for the dew temperature "
624 "\nfailed to converge in the maximum number of iterations.");
625 break;
626 case 4:
627 ERROR_REPORTER_HERE(ASC_USER_ERROR
628 , "There is no dew temperature; all components are supercritical at "
629 "\npressure %.2f Pa.", p);
630 break;
631 }
632 outputs[0] = MIN_T;
633 return 1;
634 }else{
635 outputs[0] = T_d;
636 return 0;
637 }
638 }
639
640 /*
641 Find and return the mixture bubble temperature at some pressure.
642 */
643 int mixture_bubble_T_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
644 double *inputs, double *outputs, double *jacobian){
645 CALCPREP(1,1);
646
647 int result;
648 double p = inputs[0] /* the mixture pressure */
649 , T_b /* the mixture bubble temperature (if there is any) */
650 , tol = MIX_XTOL /* tolerance to use in solving for dew temperature */
651 ;
652
653 result = mixture_bubble_temperature(&T_b, MS, p, tol, &err);
654 if(result){
655 switch(result){
656 case 1:
657 ERROR_REPORTER_HERE(ASC_USER_ERROR
658 , "The bubble temperature converged on a non-solution point");
659 break;
660 case 2:
661 ERROR_REPORTER_HERE(ASC_USER_ERROR
662 , "The bubble temperature converged on Infinity or NaN.");
663 break;
664 case 3:
665 ERROR_REPORTER_HERE(ASC_USER_ERROR
666 , "The root-finding algorithm that searches for the bubble "
667 "\ntemperature failed to converge in the maximum number of "
668 "\niterations.");
669 break;
670 case 4:
671 ERROR_REPORTER_HERE(ASC_USER_ERROR
672 , "There is no bubble temperature; all components are supercritical "
673 "\nat pressure %.2f Pa.", p);
674 break;
675 }
676 outputs[0] = MIN_T;
677 return 1;
678 }else{
679 outputs[0] = T_b;
680 return 0;
681 }
682 }
683
684 /* ---------------------------------------------------------------------
685 Functions to find Properties from Pressure and Enthalpy (p,h)
686 */
687 int mixture_state_T_ph_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
688 double *inputs, double *outputs, double *jacobian){
689 CALCPREP(2,1);
690
691 int ph_result = 0; /* result of finding temperature from (p,h) conditions */
692 double p = inputs[0] /* pressure */
693 , h = inputs[1] /* enthalpy */
694 , T /* temperature, determined from 'p' and 'h' */
695 , tol = MIX_XTOL /* tolerance to use in solving for dew temperature */
696 ;
697 ph_result = mixture_T_ph(&T, MS, p, h, tol, &err);
698 ROOTSOLVE_SWITCH(ph_result, "system temperature");
699
700 outputs[0] = T;
701 return 0;
702 }

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