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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3079 - (hide annotations) (download) (as text)
Thu Aug 20 02:17:28 2015 UTC (3 years, 11 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 jacob 2980 /* 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 jacob 3064 by Jacob Shealy, June 25-August 14, 2015
23 jacob 2980
24 jacob 3034 Exports functions for initial model of ideal-solution mixing
25 jacob 2980 */
26    
27 jacob 3064 #include "asc_mixture.h"
28 jacob 2980
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 jacob 2988
39 jacob 2980 #define CALCFN(NAME,INPUTS,OUTPUTS) \
40 jacob 2988 result += CreateUserFunctionBlackBox( #NAME \
41 jacob 3061 , 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 jacob 2980 , INPUTS,OUTPUTS \
47 jacob 3061 , NAME##_help /* help text/documentation */ \
48 jacob 2980 , 0.0 \
49     )
50    
51 jacob 3053 /* CALCFN(mixture_p,2,1); */
52 jacob 2980 CALCFN(mixture_rho,2,1);
53 jacob 3049 CALCFN(mixture_phase_rho,3,1);
54 jacob 3064 CALCFN(mixture_comps_rho,4,1);
55 jacob 3049 CALCFN(mixture_u,2,1);
56     CALCFN(mixture_phase_u,3,1);
57 jacob 3062 CALCFN(mixture_comps_u,4,1);
58 jacob 3049 CALCFN(mixture_h,2,1);
59     CALCFN(mixture_phase_h,3,1);
60 jacob 3062 CALCFN(mixture_comps_h,4,1);
61 jacob 3053 CALCFN(mixture_cp,2,1);
62     CALCFN(mixture_phase_cp,3,1);
63 jacob 3062 CALCFN(mixture_comps_cp,4,1);
64 jacob 3053 CALCFN(mixture_cv,2,1);
65     CALCFN(mixture_phase_cv,3,1);
66 jacob 3062 CALCFN(mixture_comps_cv,4,1);
67 jacob 2980
68 jacob 3053 CALCFN(mixture_s,2,1);
69     CALCFN(mixture_phase_s,3,1);
70 jacob 3064 CALCFN(mixture_comps_s,4,1);
71 jacob 3053 CALCFN(mixture_g,2,1);
72     CALCFN(mixture_phase_g,3,1);
73 jacob 3064 CALCFN(mixture_comps_g,4,1);
74 jacob 3053 CALCFN(mixture_a,2,1);
75 jacob 3064 CALCFN(mixture_phase_a,3,1);
76     CALCFN(mixture_comps_a,4,1);
77 jacob 3053
78 jacob 3064 CALCFN(mixture_count_phases,2,4);
79     CALCFN(mixture_count_components,3,1);
80     CALCFN(mixture_component_frac,4,1);
81 jacob 3053 CALCFN(mixture_component_cnum,4,1);
82 jacob 3061 CALCFN(mixture_dew_p,1,1);
83 jacob 3058 CALCFN(mixture_bubble_p,1,1);
84 jacob 3061 CALCFN(mixture_dew_T,1,1);
85     CALCFN(mixture_bubble_T,1,1);
86 jacob 3069 CALCFN(mixture_state_T_ph,2,1);
87 jacob 3053
88 jacob 2980 if(result){
89     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunctionBlackBox result is %d\n",result);
90     }
91     return result;
92     }
93    
94     /*
95 jacob 3062 Function to prepare persistent data, namely the MixtureSpec struct which
96     will be used in modeling the mixture.
97 jacob 2980 */
98 jacob 2988 int asc_mixture_prepare(struct BBoxInterp *bbox, struct Instance *data, struct gl_list_t *arglist){
99 jacob 2980
100     #define CHECK_TYPE(VAR,TYPE,NAME,TYPENAME) \
101     if(InstanceKind(VAR)!=TYPE){ \
102     ERROR_REPORTER_HERE(ASC_USER_ERROR \
103 jacob 2988 , "DATA member '%s' has type-value %#o, but must have %s (value %#o)" \
104 jacob 2980 , NAME, InstanceKind(VAR), TYPENAME, TYPE); \
105     return 1; \
106 jacob 3034 }else{ \
107 jacob 3062 /* ERROR_REPORTER_HERE(ASC_USER_NOTE, "DATA member %s has correct type %s" */ \
108     /* , NAME, TYPENAME); */ \
109 jacob 2980 }
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 jacob 3034 unsigned i;
119 jacob 2980 struct Instance *npureinst, *compinst, *xinst, *typeinst, *srcinst;
120 jacob 2988 unsigned npure;
121 jacob 2980
122     mix_symbols[0] = AddSymbol("npure");
123     mix_symbols[1] = AddSymbol("components");
124     mix_symbols[2] = AddSymbol("xs");
125 jacob 3079 mix_symbols[3] = AddSymbol("eos");
126 jacob 2980 mix_symbols[4] = AddSymbol("source");
127    
128 jacob 3034 /* 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 jacob 2980 /* number of pure components -- required */
134     npureinst = ChildByChar(data, mix_symbols[NPURE_SYM]);
135 jacob 3034 CHECK_EXIST_TYPE(npureinst, INTEGER_CONSTANT_INST, "npure", "'integer constant'");
136     npure = (int *)(IC_INST(npureinst)->value);
137 jacob 3062 /* ERROR_REPORTER_HERE(ASC_USER_NOTE, "Number of pures is %i", npure); */
138 jacob 2980
139 jacob 2988 const struct gl_list_t *comps_gl;
140 jacob 3034 const struct gl_list_t *xs_gl;
141 jacob 3062 /* const */ char *type = NULL
142     , **comps = ASC_NEW_ARRAY(/* const */ char *, npure)
143     , **srcs = ASC_NEW_ARRAY(/* const */ char *, npure);
144 jacob 3040 double *xs = ASC_NEW_ARRAY(double, npure);
145 jacob 2980
146 jacob 3034 /* Component names -- required */
147 jacob 2980 compinst = ChildByChar(data, mix_symbols[COMP_SYM]);
148 jacob 3034 CHECK_EXIST_TYPE(compinst, ARRAY_INT_INST, "components"
149     , "array indexed with integers");
150 jacob 2980 comps_gl = ARY_INST(compinst)->children;
151    
152 jacob 3034 /*
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 jacob 2980 typeinst = ChildByChar(data, mix_symbols[TYPE_SYM]);
157     if(typeinst){
158 jacob 3034 CHECK_TYPE(typeinst, SYMBOL_CONSTANT_INST, "type", "'symbol constant'");
159     type = SCP(SYMC_INST(typeinst)->value); /* read 'typeinst' into a string */
160 jacob 2980 if(type && strlen(type)==0){
161 jacob 3034 char t[] = "pengrob";
162 jacob 2980 type = t;
163     }
164 jacob 3037 }else{
165     char t[] = "pengrob";
166     type = t;
167 jacob 2980 }
168    
169 jacob 3034 /*
170     Data string representing source -- not required, so check for existence
171     before checking type and assigning to the char array 'source'.
172     */
173 jacob 2980 srcinst = ChildByChar(data, mix_symbols[SOURCE_SYM]);
174     if(srcinst){
175 jacob 3034 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 jacob 3037 }else if(!srcs[0]){
180     srcs[0] = NULL;
181 jacob 2980 }
182 jacob 3037 }else{
183     srcs[0] = NULL;
184 jacob 2980 }
185 jacob 3037 for(i=1;i<npure;i++){
186     srcs[i] = srcs[0];
187     }
188 jacob 2980
189 jacob 3034 /* 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 jacob 2980 /*
196 jacob 3034 Check that the lengths of the arrays 'comps_gl' and 'xs_gl' are equal,
197     and equal to 'npure'
198 jacob 2980 */
199 jacob 3034 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 jacob 2980
230 jacob 3040 /* Read contents of 'comps_gl' and 'xs_gl' into 'comps_child' and 'xs_child' */
231     for(i=0;i<npure;i++){
232 jacob 3064 comps[i] = SCP(SYMC_INST(InstanceChild(compinst, i+1))->value);
233     xs[i] = RC_INST(InstanceChild(xinst, i+1))->value;
234 jacob 3040 }
235 jacob 3034
236     /* Create mixture specification in a MixtureSpec struct */
237     MixtureError merr = MIXTURE_NO_ERROR;
238 jacob 3062 bbox->user_data = (void *) build_MixtureSpec(npure, xs, (void **)comps, type, srcs, &merr);
239 jacob 3079 ERROR_REPORTER_HERE(ASC_USER_NOTE, "The equation of state is %s", type);
240 jacob 3040
241 jacob 3034 return 0;
242 jacob 2980 }
243    
244 jacob 3064 /* ---------------------------------------------------------------------
245 jacob 2980 Functions which calculate results
246     */
247 jacob 3053 #if 0
248 jacob 2980 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 jacob 3053 #endif
260 jacob 2980
261 jacob 3058 /* ---------------------------------------------------------------------
262     Mixture property-calculation functions
263     */
264     /*
265     Find and return the overall mixture density
266     */
267 jacob 2980 int mixture_rho_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
268     double *inputs, double *outputs, double *jacobian){
269 jacob 3053 #define NPHASE PS->phases
270     #define PTYPE PS->ph_type[i]
271 jacob 2980 CALCPREP(2,1);
272 jacob 3053 CALCFLASH;
273 jacob 3062 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS);
274 jacob 2980
275 jacob 3049 double rhos[MS->pures];
276     outputs[0] = mixture_rho(PM, rhos);
277 jacob 3058 #if 0
278 jacob 3062 unsigned i;
279 jacob 3049 ERROR_REPORTER_HERE(ASC_USER_NOTE, "The overall mixture density is %g", outputs[0]);
280 jacob 3053 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 jacob 3049 }
284 jacob 3058 #endif
285 jacob 2980 return 0;
286 jacob 3053 #undef PTYPE
287     #undef NPHASE
288 jacob 2980 }
289 jacob 3049
290 jacob 3058 /*
291     Find and return the mixture density for a single phase
292     */
293 jacob 3049 int mixture_phase_rho_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
294     double *inputs, double *outputs, double *jacobian){
295     CALCPREP(3,1);
296 jacob 3053 CALCFLASH;
297 jacob 3062 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS);
298 jacob 3049
299 jacob 3069 unsigned ph_num = ((unsigned) inputs[2]) - 1; /* number of the phase */
300     double rhos[PS->phases]; /* individual phase densities */
301 jacob 3049
302 jacob 3062 mixture_rho(PM, rhos);
303 jacob 3069 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 jacob 3049 return 0;
310     }
311    
312 jacob 3062 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 jacob 3069 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 jacob 3062 return 0;
334     }
335    
336 jacob 3058 /*
337     Macros used to generate mixture-property functions automatically
338     */
339 jacob 3049 #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 jacob 3053 CALCFLASH; \
344 jacob 3062 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS); \
345 jacob 3049 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 jacob 3053 CALCFLASH; \
355 jacob 3062 PhaseMixState *PM = fill_PhaseMixState(T, p, PS, MS); \
356 jacob 3069 unsigned ph_num = ((unsigned) inputs[2]) - 1; /* number of the phase */ \
357 jacob 3062 double props[PS->phases]; /* internal by-phase property values */ \
358     mixture_##PROP(PM, props, &err); \
359 jacob 3069 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 jacob 3049 return 0; \
366     }
367    
368 jacob 3069 #define RHO PS->PH[ph_num]->rhos[comp_num]
369     #define PPF PS->PH[ph_num]->PF[comp_num]
370    
371 jacob 3062 #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 jacob 3069 unsigned ph_num = ((int) inputs[2]) - 1 /* number of the phase */ \
377     , comp_num = ((int) inputs[3]) - 1; /* number of the component */ \
378 jacob 3062 /* Access the property function from FPROPS for an individual component */ \
379 jacob 3069 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 jacob 3062 return 0; \
393     }
394    
395 jacob 3064 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 jacob 3049
399 jacob 3053 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 jacob 3064 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 jacob 3069 #undef PPF
407     #undef RHOS
408 jacob 3062
409 jacob 3058 /* ---------------------------------------------------------------------
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 jacob 3064 int mixture_count_phases_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
418 jacob 3049 double *inputs, double *outputs, double *jacobian){
419 jacob 3053 CALCPREP(2,4);
420     CALCFLASH;
421 jacob 3049
422 jacob 3053 unsigned i;
423     for(i=1;i<4;i++){
424 jacob 3064 /* Set all outputs to zero initially */
425 jacob 3053 outputs[i] = 0.0;
426     }
427 jacob 3064 outputs[0] = (double) PS->phases;
428 jacob 3049
429 jacob 3053 ERROR_REPORTER_HERE(ASC_USER_NOTE, "There are %u phases", PS->phases);
430     for(i=0;i<PS->phases;i++){
431 jacob 3058 /* For all phases present, set the corresponding phase-fraction output
432     equal to the mass fraction of the phase. */
433 jacob 3053 outputs[i+1] = PS->ph_frac[i];
434 jacob 3064
435 jacob 3053 ERROR_REPORTER_HERE(ASC_USER_NOTE, "\tPhase number %u is %s", i+1,
436     MIX_PHASE_STRING(PS->ph_type[i]));
437     }
438 jacob 3064 return 0;
439     }
440 jacob 3053
441 jacob 3064 /*
442 jacob 3069 Find and return the number of components of some phase within the mixture.
443 jacob 3064 */
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 jacob 3069 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 jacob 3064 }else{
453 jacob 3069 ERROR_REPORTER_HERE(ASC_USER_ERROR, "There is no phase number %u", ph_num+1);
454 jacob 3064 outputs[0] = 0.0;
455     }
456 jacob 3049 return 0;
457     }
458    
459 jacob 3058 /*
460     Find and return the mass fraction of some j_th component within the i_th
461     phase of the mixture.
462     */
463 jacob 3064 int mixture_component_frac_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
464 jacob 3049 double *inputs, double *outputs, double *jacobian){
465 jacob 3053 CALCPREP(4,1);
466     CALCFLASH;
467 jacob 3049
468 jacob 3058 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 jacob 3049 return 0;
483     }
484    
485 jacob 3053 /*
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 jacob 3049 double *inputs, double *outputs, double *jacobian){
491 jacob 3053 CALCPREP(4,1);
492     CALCFLASH;
493 jacob 3049
494 jacob 3058 unsigned i = ((unsigned) inputs[2]) - 1
495     , j = ((unsigned) inputs[3]) - 1;
496 jacob 3053 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 jacob 3049 return 0;
509     }
510 jacob 3058
511     /*
512 jacob 3061 Find and return the mixture dew pressure at some temperature.
513 jacob 3058 */
514 jacob 3061 int mixture_dew_p_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
515 jacob 3058 double *inputs, double *outputs, double *jacobian){
516     CALCPREP(1,1);
517    
518 jacob 3062 int result;
519 jacob 3061 double T = inputs[0] /* the mixture temperature */
520 jacob 3062 , p_d /* the mixture dew pressure (if there is any) */
521     , tol = MIX_XTOL /* tolerance to use in solving for dew pressure */
522     ;
523 jacob 3058
524 jacob 3062 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 jacob 3061 }
552 jacob 3058 }
553    
554     /*
555 jacob 3061 Find and return the mixture bubble pressure at some temperature.
556 jacob 3058 */
557 jacob 3061 int mixture_bubble_p_calc(struct BBoxInterp *bbox, int ninputs, int noutputs,
558 jacob 3058 double *inputs, double *outputs, double *jacobian){
559     CALCPREP(1,1);
560    
561 jacob 3062 int result;
562 jacob 3061 double T = inputs[0] /* the mixture temperature */
563 jacob 3062 , p_b /* the mixture bubble pressure (if there is any) */
564     , tol = MIX_XTOL /* tolerance to use in solving for bubble pressure */
565     ;
566 jacob 3058
567 jacob 3062 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 jacob 3061 }
595 jacob 3058 }
596 jacob 3061
597 jacob 3064 /*
598     Find and return the mixture dew temperature at some pressure.
599     */
600 jacob 3061 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 jacob 3062 int result;
605 jacob 3061 double p = inputs[0] /* the mixture pressure */
606 jacob 3062 , T_d /* the mixture dew temperature (if there is any) */
607     , tol = MIX_XTOL /* tolerance to use in solving for dew temperature */
608     ;
609 jacob 3061
610 jacob 3062 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 jacob 3061 }
638     }
639    
640 jacob 3064 /*
641     Find and return the mixture bubble temperature at some pressure.
642     */
643 jacob 3061 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 jacob 3062 int result;
648 jacob 3061 double p = inputs[0] /* the mixture pressure */
649 jacob 3062 , T_b /* the mixture bubble temperature (if there is any) */
650     , tol = MIX_XTOL /* tolerance to use in solving for dew temperature */
651     ;
652 jacob 3061
653 jacob 3062 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 jacob 3061 }
682     }
683    
684 jacob 3069 /* ---------------------------------------------------------------------
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