johnpye/fprops/asc_helmholtz.c

/*  ASCEND modelling environment
    Copyright (C) 2008 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//** @file
    Wrapper for helmholtz.c to allow access to the routine from ASCEND.
*/

/* include the external function API from libascend... */
#include <compiler/extfunc.h>

/* include error reporting API as well, so we can send messages to user */
#include <utilities/error.h>

/* for accessing the DATA instance */
#include <compiler/child.h>
#include <general/list.h>
#include <compiler/module.h>
#include <compiler/childinfo.h>
#include <compiler/parentchild.h>
#include <compiler/slist.h>
#include <compiler/type_desc.h>
#include <compiler/packages.h>
#include <compiler/symtab.h>
#include <compiler/instquery.h>
#include <compiler/instmacro.h>
#include <compiler/instance_types.h>

/* the code that we're wrapping... */
#include "helmholtz.h"

/* for the moment, species data are defined in C code, we'll implement something
better later on, hopefully. */
#include "ammonia.h"
#include "nitrogen.h"
#include "hydrogen.h"
#include "water.h"

#ifndef ASC_EXPORT
# error "Where is ASC_EXPORT?"
#endif


/*------------------------------------------------------------------------------
  FORWARD DECLARATIONS
*/

ExtBBoxInitFunc helmholtz_prepare;
ExtBBoxFunc helmholtz_p_calc;
ExtBBoxFunc helmholtz_u_calc;
ExtBBoxFunc helmholtz_s_calc;
ExtBBoxFunc helmholtz_h_calc;
ExtBBoxFunc helmholtz_a_calc;

/*------------------------------------------------------------------------------
  GLOBALS
*/

/* place to store symbols needed for accessing ASCEND's instance tree */
static symchar *helmholtz_symbols[1];
#define COMPONENT_SYM helmholtz_symbols[0]

static const char *helmholtz_p_help = "Calculate pressure from temperature and density, using Helmholtz fundamental correlation";
static const char *helmholtz_u_help = "Calculate specific internal energy from temperature and density, using Helmholtz fundamental correlation";
static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";

/*------------------------------------------------------------------------------
  REGISTRATION FUNCTION
*/

/**
    This is the function called from "IMPORT helmholtz"

    It sets up the functions contained in this external library
*/
extern
ASC_EXPORT int helmholtz_register(){
    int result = 0;

    ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");

#define CALCFN(NAME,INPUTS,OUTPUTS) \
    result += CreateUserFunctionBlackBox(#NAME \
        , helmholtz_prepare \
        , NAME##_calc /* value */ \
        , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
        , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
        , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
        , INPUTS,OUTPUTS /* inputs, outputs */ \
        , NAME##_help /* help text */ \
        , 0.0 \
    ) /* returns 0 on success */

#define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
    result += CreateUserFunctionBlackBox(#NAME \
        , helmholtz_prepare \
        , NAME##_calc /* value */ \
        , NAME##_calc /* derivatives */ \
        , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
        , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
        , INPUTS,OUTPUTS /* inputs, outputs */ \
        , NAME##_help /* help text */ \
        , 0.0 \
    ) /* returns 0 on success */

    CALCFNDERIV(helmholtz_p,2,1);
    CALCFN(helmholtz_u,2,1);
    CALCFN(helmholtz_s,2,1);
    CALCFN(helmholtz_h,2,1);
    CALCFN(helmholtz_a,2,1);

#undef CALCFN

    if(result){
        ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
    }
    return result;
}

/**
   'helmholtz_prepare' just gets the data member and checks that it's
    valid, and stores it in the blackbox data field.
*/
int helmholtz_prepare(struct BBoxInterp *bbox,
       struct Instance *data,
       struct gl_list_t *arglist
){
    struct Instance *compinst;
    const char *comp;

    helmholtz_symbols[0] = AddSymbol("component");

    /* get the data file name (we will look for this file in the ASCENDLIBRARY path) */
    compinst = ChildByChar(data,COMPONENT_SYM);
    if(!compinst){
        ERROR_REPORTER_HERE(ASC_USER_ERROR
            ,"Couldn't locate 'component' in DATA, please check usage of HELMHOLTZ."
        );
        return 1;
    }
    if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
        return 1;
    }
    comp = SCP(SYMC_INST(compinst)->value);
    CONSOLE_DEBUG("COMPONENT: %s",comp);
    if(comp==NULL || strlen(comp)==0){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
        return 1;
    }

    if(strcmp(comp,"ammonia")==0){
        bbox->user_data = (void*)&helmholtz_data_ammonia;
    }else if(strcmp(comp,"nitrogen")==0){
        bbox->user_data = (void*)&helmholtz_data_nitrogen;
    }else if(strcmp(comp,"hydrogen")==0){
        bbox->user_data = (void*)&helmholtz_data_hydrogen;
    }else if(strcmp(comp,"water")==0){
        bbox->user_data = (void*)&helmholtz_data_water;
    }else{
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
    }

    ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
    return 0;
}

/*------------------------------------------------------------------------------
  EVALULATION ROUTINES
*/

#define CALCPREPARE \
    /* a few checks about the input requirements */ \
    if(ninputs != 2)return -1; \
    if(noutputs != 1)return -2; \
    if(inputs==NULL)return -3; \
    if(outputs==NULL)return -4; \
    if(bbox==NULL)return -5; \
    \
    /* the 'user_data' in the black box object will contain the */\
    /* coefficients required for this fluid; cast it to the required form: */\
    HelmholtzData *helmholtz_data = (HelmholtzData *)bbox->user_data

/**
    Evaluation function for 'helmholtz_p'.
    @param inputs array with values of inputs T and rho.
    @param outputs array with just value of p
    @param jacobian, the partial derivative df/dx, where
        each row is df[i]/dx[j] over each j for the y_out[i] of
        matching index. The jacobian array is 1-D, row major, i.e.
        df[i]/dx[j] -> jacobian[i*ninputs+j].
    @return 0 on success
*/
int helmholtz_p_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE;

    /* first input is temperature, second is molar density */
    if(bbox->task == bb_func_eval){
        outputs[0] = helmholtz_p(inputs[0], inputs[1], helmholtz_data);
    }else{
        //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
        jacobian[0*1+0] = helmholtz_dpdT_rho(inputs[0], inputs[1], helmholtz_data);
        jacobian[0*1+1] = helmholtz_dpdrho_T(inputs[0], inputs[1], helmholtz_data);
    }

    /* no need to worry about error states etc. */
    return 0;
}


/**
    Evaluation function for 'helmholtz_u'
    @param jacobian ignored
    @return 0 on success
*/
int helmholtz_u_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE;

    /* first input is temperature, second is molar density */
    if(bbox->task == bb_func_eval){
        outputs[0] = helmholtz_u(inputs[0], inputs[1], helmholtz_data);
    }else{
        jacobian[0*1+0] = helmholtz_dudT_rho(inputs[0], inputs[1], helmholtz_data);
        jacobian[0*1+1] = helmholtz_dudrho_T(inputs[0], inputs[1], helmholtz_data);
    }

    /* no need to worry about error states etc. */
    return 0;
}


/**
    Evaluation function for 'helmholtz_h'
    @param jacobian ignored
    @return 0 on success
*/
int helmholtz_s_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE;

    /* first input is temperature, second is molar density */
    outputs[0] = helmholtz_s(inputs[0], inputs[1], helmholtz_data);

    /* no need to worry about error states etc. */
    return 0;
}


/**
    Evaluation function for 'helmholtz_h'
    @param jacobian ignored
    @return 0 on success
*/
int helmholtz_h_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE;

    /* first input is temperature, second is molar density */
    if(bbox->task == bb_func_eval){
        outputs[0] = helmholtz_h(inputs[0], inputs[1], helmholtz_data);
    }else{
        //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
        jacobian[0*1+0] = helmholtz_dhdT_rho(inputs[0], inputs[1], helmholtz_data);
        jacobian[0*1+1] = helmholtz_dhdrho_T(inputs[0], inputs[1], helmholtz_data);
    }

    /* no need to worry about error states etc. */
    return 0;
}


/**
    Evaluation function for 'helmholtz_h'
    @param jacobian ignored
    @return 0 on success
*/
int helmholtz_a_calc(struct BBoxInterp *bbox,
        int ninputs, int noutputs,
        double *inputs, double *outputs,
        double *jacobian
){
    CALCPREPARE;

    /* first input is temperature, second is molar density */
    outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);

    /* no need to worry about error states etc. */
    return 0;
}


1	/* ASCEND modelling environment
2	Copyright (C) 2008 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,
10	but WITHOUT ANY WARRANTY; without even the implied warranty of
11	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12	GNU 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, Inc., 59 Temple Place - Suite 330,
17	Boston, MA 02111-1307, USA.
18	//* @file
19	Wrapper for helmholtz.c to allow access to the routine from ASCEND.
20	*/
21
22	/* include the external function API from libascend... */
23	#include <compiler/extfunc.h>
24
25	/* include error reporting API as well, so we can send messages to user */
26	#include <utilities/error.h>
27
28	/* for accessing the DATA instance */
29	#include <compiler/child.h>
30	#include <general/list.h>
31	#include <compiler/module.h>
32	#include <compiler/childinfo.h>
33	#include <compiler/parentchild.h>
34	#include <compiler/slist.h>
35	#include <compiler/type_desc.h>
36	#include <compiler/packages.h>
37	#include <compiler/symtab.h>
38	#include <compiler/instquery.h>
39	#include <compiler/instmacro.h>
40	#include <compiler/instance_types.h>
41
42	/* the code that we're wrapping... */
43	#include "helmholtz.h"
44
45	/* for the moment, species data are defined in C code, we'll implement something
46	better later on, hopefully. */
47	#include "ammonia.h"
48	#include "nitrogen.h"
49	#include "hydrogen.h"
50	#include "water.h"
51
52	#ifndef ASC_EXPORT
53	# error "Where is ASC_EXPORT?"
54	#endif
55
56
57	/*------------------------------------------------------------------------------
58	FORWARD DECLARATIONS
59	*/
60
61	ExtBBoxInitFunc helmholtz_prepare;
62	ExtBBoxFunc helmholtz_p_calc;
63	ExtBBoxFunc helmholtz_u_calc;
64	ExtBBoxFunc helmholtz_s_calc;
65	ExtBBoxFunc helmholtz_h_calc;
66	ExtBBoxFunc helmholtz_a_calc;
67
68	/*------------------------------------------------------------------------------
69	GLOBALS
70	*/
71
72	/* place to store symbols needed for accessing ASCEND's instance tree */
73	static symchar *helmholtz_symbols[1];
74	#define COMPONENT_SYM helmholtz_symbols[0]
75
76	static const char *helmholtz_p_help = "Calculate pressure from temperature and density, using Helmholtz fundamental correlation";
77	static const char *helmholtz_u_help = "Calculate specific internal energy from temperature and density, using Helmholtz fundamental correlation";
78	static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
79	static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
80	static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";
81
82	/*------------------------------------------------------------------------------
83	REGISTRATION FUNCTION
84	*/
85
86	/**
87	This is the function called from "IMPORT helmholtz"
88
89	It sets up the functions contained in this external library
90	*/
91	extern
92	ASC_EXPORT int helmholtz_register(){
93	int result = 0;
94
95	ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
96
97	#define CALCFN(NAME,INPUTS,OUTPUTS) \
98	result += CreateUserFunctionBlackBox(#NAME \
99	, helmholtz_prepare \
100	, NAME##_calc /* value */ \
101	, (ExtBBoxFunc)NULL / derivatives not provided yet*/ \
102	, (ExtBBoxFunc)NULL / hessian not provided yet */ \
103	, (ExtBBoxFinalFunc)NULL / finalisation not implemented */ \
104	, INPUTS,OUTPUTS /* inputs, outputs */ \
105	, NAME##_help /* help text */ \
106	, 0.0 \
107	) /* returns 0 on success */
108
109	#define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
110	result += CreateUserFunctionBlackBox(#NAME \
111	, helmholtz_prepare \
112	, NAME##_calc /* value */ \
113	, NAME##_calc /* derivatives */ \
114	, (ExtBBoxFunc)NULL / hessian not provided yet */ \
115	, (ExtBBoxFinalFunc)NULL / finalisation not implemented */ \
116	, INPUTS,OUTPUTS /* inputs, outputs */ \
117	, NAME##_help /* help text */ \
118	, 0.0 \
119	) /* returns 0 on success */
120
121	CALCFNDERIV(helmholtz_p,2,1);
122	CALCFN(helmholtz_u,2,1);
123	CALCFN(helmholtz_s,2,1);
124	CALCFN(helmholtz_h,2,1);
125	CALCFN(helmholtz_a,2,1);
126
127	#undef CALCFN
128
129	if(result){
130	ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
131	}
132	return result;
133	}
134
135	/**
136	'helmholtz_prepare' just gets the data member and checks that it's
137	valid, and stores it in the blackbox data field.
138	*/
139	int helmholtz_prepare(struct BBoxInterp *bbox,
140	struct Instance *data,
141	struct gl_list_t *arglist
142	){
143	struct Instance *compinst;
144	const char *comp;
145
146	helmholtz_symbols[0] = AddSymbol("component");
147
148	/* get the data file name (we will look for this file in the ASCENDLIBRARY path) */
149	compinst = ChildByChar(data,COMPONENT_SYM);
150	if(!compinst){
151	ERROR_REPORTER_HERE(ASC_USER_ERROR
152	,"Couldn't locate 'component' in DATA, please check usage of HELMHOLTZ."
153	);
154	return 1;
155	}
156	if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
157	ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
158	return 1;
159	}
160	comp = SCP(SYMC_INST(compinst)->value);
161	CONSOLE_DEBUG("COMPONENT: %s",comp);
162	if(comp==NULL \|\| strlen(comp)==0){
163	ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
164	return 1;
165	}
166
167	if(strcmp(comp,"ammonia")==0){
168	bbox->user_data = (void*)&helmholtz_data_ammonia;
169	}else if(strcmp(comp,"nitrogen")==0){
170	bbox->user_data = (void*)&helmholtz_data_nitrogen;
171	}else if(strcmp(comp,"hydrogen")==0){
172	bbox->user_data = (void*)&helmholtz_data_hydrogen;
173	}else if(strcmp(comp,"water")==0){
174	bbox->user_data = (void*)&helmholtz_data_water;
175	}else{
176	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
177	}
178
179	ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
180	return 0;
181	}
182
183	/*------------------------------------------------------------------------------
184	EVALULATION ROUTINES
185	*/
186
187	#define CALCPREPARE \
188	/* a few checks about the input requirements */ \
189	if(ninputs != 2)return -1; \
190	if(noutputs != 1)return -2; \
191	if(inputs==NULL)return -3; \
192	if(outputs==NULL)return -4; \
193	if(bbox==NULL)return -5; \
194	\
195	/* the 'user_data' in the black box object will contain the */\
196	/* coefficients required for this fluid; cast it to the required form: */\
197	HelmholtzData helmholtz_data = (HelmholtzData )bbox->user_data
198
199	/**
200	Evaluation function for 'helmholtz_p'.
201	@param inputs array with values of inputs T and rho.
202	@param outputs array with just value of p
203	@param jacobian, the partial derivative df/dx, where
204	each row is df[i]/dx[j] over each j for the y_out[i] of
205	matching index. The jacobian array is 1-D, row major, i.e.
206	df[i]/dx[j] -> jacobian[i*ninputs+j].
207	@return 0 on success
208	*/
209	int helmholtz_p_calc(struct BBoxInterp *bbox,
210	int ninputs, int noutputs,
211	double inputs, double outputs,
212	double *jacobian
213	){
214	CALCPREPARE;
215
216	/* first input is temperature, second is molar density */
217	if(bbox->task == bb_func_eval){
218	outputs[0] = helmholtz_p(inputs[0], inputs[1], helmholtz_data);
219	}else{
220	//ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
221	jacobian[0*1+0] = helmholtz_dpdT_rho(inputs[0], inputs[1], helmholtz_data);
222	jacobian[0*1+1] = helmholtz_dpdrho_T(inputs[0], inputs[1], helmholtz_data);
223	}
224
225	/* no need to worry about error states etc. */
226	return 0;
227	}
228
229
230	/**
231	Evaluation function for 'helmholtz_u'
232	@param jacobian ignored
233	@return 0 on success
234	*/
235	int helmholtz_u_calc(struct BBoxInterp *bbox,
236	int ninputs, int noutputs,
237	double inputs, double outputs,
238	double *jacobian
239	){
240	CALCPREPARE;
241
242	/* first input is temperature, second is molar density */
243	if(bbox->task == bb_func_eval){
244	outputs[0] = helmholtz_u(inputs[0], inputs[1], helmholtz_data);
245	}else{
246	jacobian[0*1+0] = helmholtz_dudT_rho(inputs[0], inputs[1], helmholtz_data);
247	jacobian[0*1+1] = helmholtz_dudrho_T(inputs[0], inputs[1], helmholtz_data);
248	}
249
250	/* no need to worry about error states etc. */
251	return 0;
252	}
253
254
255	/**
256	Evaluation function for 'helmholtz_h'
257	@param jacobian ignored
258	@return 0 on success
259	*/
260	int helmholtz_s_calc(struct BBoxInterp *bbox,
261	int ninputs, int noutputs,
262	double inputs, double outputs,
263	double *jacobian
264	){
265	CALCPREPARE;
266
267	/* first input is temperature, second is molar density */
268	outputs[0] = helmholtz_s(inputs[0], inputs[1], helmholtz_data);
269
270	/* no need to worry about error states etc. */
271	return 0;
272	}
273
274
275	/**
276	Evaluation function for 'helmholtz_h'
277	@param jacobian ignored
278	@return 0 on success
279	*/
280	int helmholtz_h_calc(struct BBoxInterp *bbox,
281	int ninputs, int noutputs,
282	double inputs, double outputs,
283	double *jacobian
284	){
285	CALCPREPARE;
286
287	/* first input is temperature, second is molar density */
288	if(bbox->task == bb_func_eval){
289	outputs[0] = helmholtz_h(inputs[0], inputs[1], helmholtz_data);
290	}else{
291	//ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
292	jacobian[0*1+0] = helmholtz_dhdT_rho(inputs[0], inputs[1], helmholtz_data);
293	jacobian[0*1+1] = helmholtz_dhdrho_T(inputs[0], inputs[1], helmholtz_data);
294	}
295
296	/* no need to worry about error states etc. */
297	return 0;
298	}
299
300
301	/**
302	Evaluation function for 'helmholtz_h'
303	@param jacobian ignored
304	@return 0 on success
305	*/
306	int helmholtz_a_calc(struct BBoxInterp *bbox,
307	int ninputs, int noutputs,
308	double inputs, double outputs,
309	double *jacobian
310	){
311	CALCPREPARE;
312
313	/* first input is temperature, second is molar density */
314	outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);
315
316	/* no need to worry about error states etc. */
317	return 0;
318	}
319
320
321