Writing ASCEND external relations in C

This HOWTO will guide you through the steps required to access external relation code for ASCEND from the C language. As an example, we will walk through the implementation of an external thermodynamic properties correlation using the modified Benedict-Webb-Rubin correlation, with the added benefit of a useful and reasonably well-documented piece of code. For a little more background on these so-called 'black box' external relations, see the PhD thesis of Kirk Abbott (in Publications).

All the code associated with this HOWTO (and possibly a bit more, over time) is available from models/johnpye/fprops/ in the ModelLibrary,

Decide on the C functions that you want to access

Firstly, you need to provide the C code for the functions that you want to be able to access from ASCEND. In this HOWTO, we're interested in external relations, which means equations or functions that we want to be able to use from within an ASCEND model. The most typical application of external relations is probably thermodynamic property correlations, but other examples include pressure-drop correlations (e.g. pipe friction), and perhaps getting data from lookup tables, etc.

In this case, we want to access the Modified Benedict-Webb-Rubin (MBWR) thermodynamic property correlation for the refrigerant R123. The MBWR correlation data for this refrigerant was published in 1994 by Younglove, and McLinden. We have implemented a function which provides the general form of an MBWR correlation as follows, with the intention that our fluid-specific correlation data would also be passed in to the function as a pointer to a MbwrData struct:

/* from mbwr.h: */
typedef struct MbwrData_struct{

	double rhob_c; /**< critical molar density in mol/L */
	double beta[32]; /**< constants in MBWR for the fluid in question */

} MbwrData;

/* from mbwr.c: */
double mbwr_p(double T, double rhob, MbwrData *data){

    double p;
    /* lots of code here... */
    return p;
}

We can use this code to calculate pressure p as an output, with temperature T and molar density <math>\bar \rho</math> (rhob) as inputs.

The full code for this function is available in models/johnpye/fprops/mbwr.c.

This code can be compiled to a 'shared object file' using a C compiler, or you can write a little 'main' function to drive this function and make some output, just for testing. The following test program includes the data for R123, because we haven't yet implemented a tidy general-purpose system for loading these data from a file or a database etc.

#include "mbwr.h"
#include <stdio.h>
int main(int argc, char **argv){

    if(argc<3)return 1;
    double T = atof(argv[1]); /* K */

    double rhog = atof(argv[2]); /* mol/L? */
    double p;

    /* property data for R123, from J. Phys. Chem. Ref. Data, 23:731-779, 1994. */
    const MbwrData propdata = {
        13.2117771543124 /* rho_c */
        ,{ /* beta[32] */

            -0.657453133659E-02,     0.293479845842E+01,   -0.989140469845E+02
            , 0.201029776013E+05,   -0.383566527886E+07,    0.227587641969E-02

            ,-0.908726819450E+01,    0.434181417995E+04,    0.354116464954E+07
            ,-0.635394849670E-03,    0.320786715274E+01,   -0.131276484299E+04

            ,-0.116360713718E+00,   -0.113354409016E+02,   -0.537543457327E+04
            , 0.258112416120E+01,   -0.106148632128E+00,    0.500026133667E+02

            ,-0.204326706346E+01,   -0.249438345685E+07,   -0.463962781113E+09
            ,-0.284903429588E+06,    0.974392239902E+10,   -0.637314379308E+04

            , 0.314121189813E+06,   -0.145747968225E+03,   -0.843830261449E+07
            ,-0.241138441593E+01,    0.108508031257E+04,   -0.106653193965E-01

            ,-0.121343571084E+02,   -0.257510383240E+03
        }
    };

    p = mbwr_p(T,rhog, propdata);

    fprintf(stderr,"T = %f, rhob = %f --> p = %f", T, rhob, p);

    return 0;
}

Once you're happy that your function is working as expected (and be sure you know what's happening with the units of measurement) then you can get on with the job of wrapping it from ASCEND.

Create an ASCEND wrapper library

To create a wrapper for your function, you will essentially need to create a shared library (.DLL on Windows, or .SO file on Linux) of a specific kind that can be 'dynamically loaded' at runtime by ASCEND, when ASCEND sees an IMPORT statement in your model code. Using tools like SCons isn't too hard, and if you do it carefully, you will be able to use the same technique to build your shared library on Windows, Linux and Mac.

The first important feature of the shared library that you will be creating is that it will have a 'registration function'. This function is the entry point for the shared library you are creating, and this function will do the job of telling ASCEND what external relations you are going to be providing to it. Note that it is possible to embed external methods, external relations, and external anything else you fancy in a single shared library in ASCEND; there is no requirement that you have a single shared library for each single external relation you are creating.

The registration function will look something like this:

/* forward declarations */
ExtBBoxInitFunc mbwr_p_prepare;
ExtBBoxFunc mbwr_p_calc;

/* registration function for the MBWR shared library */
extern
ASC_EXPORT int mbwr_register(){
	const char *mbwr_help = "Modified Benedict-Webb-Rubin correlation for thermodynamic properties";

	int result = 0;

	ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Initialising MBWR...\n");

	result += CreateUserFunctionBlackBox("mbwr_p"
		, mbwr_p_prepare
		, mbwr_p_calc /* value */

		, (ExtBBoxFunc*)NULL /* derivatives not provided yet*/
		, (ExtBBoxFunc*)NULL /* hessian not provided yet */
		, (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */

		, 2,1 /* inputs, outputs */
		, mbwr_help
		, 0.0

	); /* returns 0 on success */

	if(result){
		ERROR_REPORTER_HERE(ASC_PROG_ERR,"CreateUserFunction result = %d\n",result);

	}
	return result;
}

The key part of that function is that it contains a call to the function CreateUserFunctionBlackBox, and this function tells ASCEND about two functions, mbwr_p_prepare and mbwr_p_calc that we are yet to write. It tells ASCEND that we are providing a 'black box' user function, which is another way of saying 'external relation'. It says that this external relation will have 2 inputs and 1 output. It gives some 'help text' to ASCEND to allow the user to identify what it is that they have loaded. It also gives the name of the provided external relation as 'mbwr_p'; this is the name that you will be able to use to access this function when writing an ASCEND model from within the ASCEND modelling language.

You note that we add the output of CreateUserFunctionBlackBox to the variable 'result'. This is done because typically we would have multiple functions being registered in a single external library, and each registration could potentially return a non-zero error code. So we add all the result codes together and if we get a non-zero total, then we output a message.

Note also that we are using the ERROR_REPORTER_HERE function (from compiler/error.h in libascend) to convey output back to ASCEND. This is much better than outputting to the console, as output via this function will be visible in the GUI, and it has 'status' associated with it too (e.g. ASC_PROG_NOTE is a 'note' for the user, and ASC_PROG_ERR is an error being reported, which would normally be shown in red).

Preparation and Calculation wrapper functions

In the previous section, we had forward declarations for two functions that we haven't written yet:

/* forward declarations */
ExtBBoxInitFunc mbwr_p_prepare;

ExtBBoxFunc mbwr_p_calc;

These functions are where we really get into passing information to and from our external code.

'Prepare' function

This function is the first of the two that we'll implement. We should first highlight the difference between the 'registration' and the 'preparation' function. The registration function is called when ASCEND sees the IMPORT statement for our shared library, for example IMPORT "mbwr";. We can, if we want, do some preparation of specific stuff at that time, but we don't yet know anything about which functions are going to be used, and how often, and with what values.

On the other hand, the preparation function is called when the external relation is first seen in your code, for example

myprops: mbwr_p(
    T, rhob : INPUT;

    p : OUTPUT;
    mydata : DATA
);

At this point, ASCEND will know that the 'mbwr_p' function is going to be used in this model, so it's appropriate that if there is any data that needs to be preloaded or precalculated, that we should do it now, inside the 'prepare' function.

The 'prepare' function has the following structure:

int mbwr_p_prepare(struct BBoxInterp *bbox,
	   struct Instance *data,

	   struct gl_list_t *arglist
){
   /* implementation here... */
   return statuscode;

}

In particular, we are passed a pointer to a BBoxInterp structure, as well as a 'data' instance (more later), and an 'arglist'. The arglist is a list with n_inputs + n_outputs elements, each of which itself is a list of instances (allowing either single values or arrays to be passed to your external relation). You can use the arglist to check that your input arguments have the correct data type, units, etc, although in common practise we don't both to do that.

The 'data' instance is important, as it allows you to pass data to your external relation that specifies how the relation should work, but in the form of constants or parameters, as opposed to variables that will actually vary during the modelling calculations. This is typically used to pass a filename for a data file, or a name of a substance, so that the external relation will load that file or substance data before the calculation commences.

In our particular example, we plan to implement a general-purpose property calculation system here based on the MBWR correlation, so we will make use of the data instance as a place where the name of the substance can be passed into ASCEND.

Our 'mbwr_p_prepare' function therefore contains the following code, which looks for a member called 'component' inside the DATA instance, and asserts that its value be equal to "R123", because we've only implemented support for this one fluid so far.

int mbwr_p_prepare(struct BBoxInterp *bbox,
	   struct Instance *data,

	   struct gl_list_t *arglist
){
	struct Instance *compinst;
	const char *comp;

	mbwr_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 MBWR."
		);
		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,"R123")!=0){
		ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component must be 'R123' at this stage (only one component supported)");

	}

	ERROR_REPORTER_HERE(ASC_PROG_NOTE,"PREPARING MBWR_P...\n");

	bbox->user_data = (void*)&mbwr_r123;

	return 0;
}

'Calculate' function

The 'calculate' function is where ASCEND instructs our external relation to actually be evaluated. This is the part that gets called over and over again during ASCEND's solver iterations, so it's written for efficient access to the values that need to be used. The structure of the 'calculate' function is as follows:

int mbwr_p_calc(struct BBoxInterp *bbox,
		int ninputs, int noutputs,

		double *inputs, double *outputs,
		double *jacobian

){

In actual fact, the external relation API in ASCEND has been designed in such a way that the same function can be used both for function evaluation as well as calculation of partial derivatives (Jacobian matrix). If you want to use the API in this way, they your 'calculate' function needs to look in the 'bbox' struct to see the value of bbox->task as noted in the declaration of BBoxInterp. In our case we're only suing our function for simple evaluation, so we don't need to check bbox->task, and we will also ignore the jacobian parameter.

Our resulting wrapper function contains a few simple checks then then makes the call to our mbwr_p function (the one that we're wrapping) using the appropriate values from the inputs array, and places the result in the correct position in the outputs array. Note that this function doesn't include any arguments for the DATA parameter; it is assumed that you've already dojne what you needed with the DATA parameters, and they won't have changed since then. Also note that we make use of the bbox->user_data general-purpose pointer to indicate where our fluid-specific property data is located.

int mbwr_p_calc(struct BBoxInterp *bbox,
		int ninputs, int noutputs,

		double *inputs, double *outputs,
		double *jacobian

){
	/* 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: */
	MbwrData *mbwr_data = (MbwrData *)bbox->user_data;

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

	/* no need to worry about error states etc, so always return 'success': */
	return 0;
}

Note that in your 'calculate' function you will typically need to convert from the units of measurement used in external code to the 'base units' used by ASCEND. In most cases this means giving values in SI base units, such as kg, m, s, K or A, or combinations of those, such as m³/s.

If it is possible that your external function may return error codes, you have the opportunity of passing those errors back to ASCEND, or alternatively you can choose to ignore the errors. In many cases Newton iteration will work best if you can evaluate smooth curves even outside their 'acceptable range', and not return errors when values are outside those ranges.

Building and linking the code

At this point, we have all the main parts of the wrapper code for accessing our mbwr_p function from ASCEND. You can see the resulting file models/johnpye/fprops/asc_mbwr.c. It now remains to get the code compiled and linked to ASCEND. Because we have been using some code from within libascend, we need to link our dynamically library back to libascend, so this is an opportunity to make use of a standardised approach that ASCEND provides for creating derivative software.

At this point we have two files, mbwr.c (containing the code-to-be-wrapped) and asc_mbwr.c (the ASCEND wrapper code for our function). To build the shared library, we add the following SConstruct file that will make use of SCons to generate a shared library for ASCEND on our system. The only requirement is that the program 'ascend-config' should already be on our PATH, which might require some tweaking in the case of Windows.

Invalid language.

You need to specify a language like this: <source lang="html">...</source>

Supported languages for syntax highlighting:

a4c, abap, abc, abnf, actionscript, ada, agda, alan, algol, ampl, amtrix, applescript, arc, arm, as400cl, ascend, asciidoc, asp, aspect, assembler, ats, autohotkey, autoit, avenue, awk, ballerina, bat, bbcode, bcpl, bibtex, biferno, bison, blitzbasic, bms, bnf, boo, c, carbon, ceylon, charmm, chill, chpl, clean, clearbasic, clipper, clojure, clp, cmake, cobol, coffeescript, coldfusion, conf, cpp2, critic, crk, crystal, cs_block_regex, csharp, css, d, dart, delphi, diff, dockerfile, dts, dylan, ebnf, ebnf2, eiffel, elixir, elm, email, erb, erlang, euphoria, exapunks, excel, express, factor, fame, fasm, felix, fish, fortran77, fortran90, frink, fsharp, fstab, fx, gambas, gdb, gdscript, go, graphviz, haml, hare, haskell, haxe, hcl, html, httpd, hugo, icon, idl, idlang, inc_luatex, informix, ini, innosetup, interlis, io, jam, jasmin, java, javascript, js_regex, json, jsp, jsx, julia, kotlin, ldif, less, lhs, lilypond, limbo, lindenscript, lisp, logtalk, lotos, lotus, lua, luban, makefile, maple, markdown, matlab, maya, mercury, meson, miranda, mod2, mod3, modelica, moon, ms, msl, mssql, mxml, n3, nasal, nbc, nemerle, netrexx, nginx, nice, nim, nix, nsis, nxc, oberon, objc, ocaml, octave, oorexx, org, os, oz, paradox, pas, pdf, perl, php, pike, pl1, plperl, plpython, pltcl, po, polygen, pony, pov, powershell, pro, progress, ps, psl, pure, purebasic, purescript, pyrex, python, q, qmake, qml, qu, r, rebol, rego, rexx, rnc, rpg, rpl, rst, ruby, rust, s, sam, sas, scad, scala, scilab, scss, shellscript, slim, small, smalltalk, sml, snmp, snobol, solidity, spec, spn, sql, squirrel, styl, svg, swift, sybase, tcl, tcsh, terraform, tex, toml, tsql, tsx, ttcn3, txt, typescript, upc, vala, vb, verilog, vhd, vimscript, vue, wat, whiley, wren, xml, xpp, yaiff, yaml, yaml_ansible, yang, zig, znn

Providing ASCEND has been installed on your system (Note: including 'devel' files, if you're installing from an RPM package or Windows installer), you should now be able to just run SCons to build the shared library, like so:

Invalid language.

You need to specify a language like this: <source lang="html">...</source>

Supported languages for syntax highlighting:

a4c, abap, abc, abnf, actionscript, ada, agda, alan, algol, ampl, amtrix, applescript, arc, arm, as400cl, ascend, asciidoc, asp, aspect, assembler, ats, autohotkey, autoit, avenue, awk, ballerina, bat, bbcode, bcpl, bibtex, biferno, bison, blitzbasic, bms, bnf, boo, c, carbon, ceylon, charmm, chill, chpl, clean, clearbasic, clipper, clojure, clp, cmake, cobol, coffeescript, coldfusion, conf, cpp2, critic, crk, crystal, cs_block_regex, csharp, css, d, dart, delphi, diff, dockerfile, dts, dylan, ebnf, ebnf2, eiffel, elixir, elm, email, erb, erlang, euphoria, exapunks, excel, express, factor, fame, fasm, felix, fish, fortran77, fortran90, frink, fsharp, fstab, fx, gambas, gdb, gdscript, go, graphviz, haml, hare, haskell, haxe, hcl, html, httpd, hugo, icon, idl, idlang, inc_luatex, informix, ini, innosetup, interlis, io, jam, jasmin, java, javascript, js_regex, json, jsp, jsx, julia, kotlin, ldif, less, lhs, lilypond, limbo, lindenscript, lisp, logtalk, lotos, lotus, lua, luban, makefile, maple, markdown, matlab, maya, mercury, meson, miranda, mod2, mod3, modelica, moon, ms, msl, mssql, mxml, n3, nasal, nbc, nemerle, netrexx, nginx, nice, nim, nix, nsis, nxc, oberon, objc, ocaml, octave, oorexx, org, os, oz, paradox, pas, pdf, perl, php, pike, pl1, plperl, plpython, pltcl, po, polygen, pony, pov, powershell, pro, progress, ps, psl, pure, purebasic, purescript, pyrex, python, q, qmake, qml, qu, r, rebol, rego, rexx, rnc, rpg, rpl, rst, ruby, rust, s, sam, sas, scad, scala, scilab, scss, shellscript, slim, small, smalltalk, sml, snmp, snobol, solidity, spec, spn, sql, squirrel, styl, svg, swift, sybase, tcl, tcsh, terraform, tex, toml, tsql, tsx, ttcn3, txt, typescript, upc, vala, vb, verilog, vhd, vimscript, vue, wat, whiley, wren, xml, xpp, yaiff, yaml, yaml_ansible, yang, zig, znn

The resulting file, 'libmbwr_ascend.so' (on Linux) is the file that will be loaded when ASCEND sees an IMPORT "mbwr"; statement in your model file.

Obviously if you want, you can just type these compiler and linker commands manually, without using a build script. Or you can use GNU Make or CMake as you wish.

Using it

To make use of the external relation that we've now created, we can create a simple model file that makes a single evaluation of the external relation:

(* test of the fluid properties routines currently implemented in asc_mbwr.c *)
REQUIRE "atoms.a4l";
IMPORT "johnpye/fprops/mbwr";

MODEL mbwrconf;
	component IS_A symbol_constant;
	component :== 'R123';

END mbwrconf;

MODEL fprops_test;

	p IS_A pressure;

	T IS_A temperature;
	rhob IS_A molar_density;

	propsdata IS_A mbwrconf;

	props: mbwr_p(
		T, rhob : INPUT;

		p : OUTPUT;
		propsdata : DATA
	);

METHODS
METHOD on_load;
	FIX rhob, T;
	rhob := 1 {mol/m^3};

	T := 400 {K};
END on_load;
END fprops_test;

Note in particular here the way in which we've set up the 'mbwrconf' model containing the details of the fluid that we're using, and way that the 'DATA' parameter is set in the external relation call.

Note that the syntax for external relations is somewhat limited; you can't include any other terms in these 'equations': the whole inputs-to-outputs relation is defined 'external' to the ASCEND model.

The other issue here is the IMPORT statement. This statement uses the ASCENDLIBRARY environment variable to hunt for your shared library before it attempts to import it. In this case, the 'mbwr' model has been created inside the standard ASCEND ModelLibrary, so I didn't need to change the value of ASCENDLIBRARY. But sometimes you will need to do that before running ASCEND, eg:

Invalid language.

You need to specify a language like this: <source lang="html">...</source>

Supported languages for syntax highlighting:

a4c, abap, abc, abnf, actionscript, ada, agda, alan, algol, ampl, amtrix, applescript, arc, arm, as400cl, ascend, asciidoc, asp, aspect, assembler, ats, autohotkey, autoit, avenue, awk, ballerina, bat, bbcode, bcpl, bibtex, biferno, bison, blitzbasic, bms, bnf, boo, c, carbon, ceylon, charmm, chill, chpl, clean, clearbasic, clipper, clojure, clp, cmake, cobol, coffeescript, coldfusion, conf, cpp2, critic, crk, crystal, cs_block_regex, csharp, css, d, dart, delphi, diff, dockerfile, dts, dylan, ebnf, ebnf2, eiffel, elixir, elm, email, erb, erlang, euphoria, exapunks, excel, express, factor, fame, fasm, felix, fish, fortran77, fortran90, frink, fsharp, fstab, fx, gambas, gdb, gdscript, go, graphviz, haml, hare, haskell, haxe, hcl, html, httpd, hugo, icon, idl, idlang, inc_luatex, informix, ini, innosetup, interlis, io, jam, jasmin, java, javascript, js_regex, json, jsp, jsx, julia, kotlin, ldif, less, lhs, lilypond, limbo, lindenscript, lisp, logtalk, lotos, lotus, lua, luban, makefile, maple, markdown, matlab, maya, mercury, meson, miranda, mod2, mod3, modelica, moon, ms, msl, mssql, mxml, n3, nasal, nbc, nemerle, netrexx, nginx, nice, nim, nix, nsis, nxc, oberon, objc, ocaml, octave, oorexx, org, os, oz, paradox, pas, pdf, perl, php, pike, pl1, plperl, plpython, pltcl, po, polygen, pony, pov, powershell, pro, progress, ps, psl, pure, purebasic, purescript, pyrex, python, q, qmake, qml, qu, r, rebol, rego, rexx, rnc, rpg, rpl, rst, ruby, rust, s, sam, sas, scad, scala, scilab, scss, shellscript, slim, small, smalltalk, sml, snmp, snobol, solidity, spec, spn, sql, squirrel, styl, svg, swift, sybase, tcl, tcsh, terraform, tex, toml, tsql, tsx, ttcn3, txt, typescript, upc, vala, vb, verilog, vhd, vimscript, vue, wat, whiley, wren, xml, xpp, yaiff, yaml, yaml_ansible, yang, zig, znn

In the ASCEND window, you should see some debug output from the 'prepare' function, so you will know that your library was successfully loaded. Then you can solve the model and check that your model is giving results that actually make sense.

It is worth thinking about setting up Model Self-testing when you create new external libraries, so that you have a quick way of checking that subsequent changes to your code or other ASCEND code haven't broken your output or changed any values from what is expected.

Adding derivatives

At this point it needs to be explained that the Newton-based solver algorithm that is being used here by ASCEND to solve more complex models will often fail with an implementation like this. This is a result of the Newton algorithm, which requires partial derivatives of the residual for each equation in the system of equations with respect to each variable in the system. When an external relation is provided, ASCEND's usual system of automatic symbolic differentiation of equations can't be used. Instead, ASCEND calculates a rather crude finite-difference derivative, and this derivative can often be sufficiently inaccurate that it affects ASCEND's ability to converge to a solution.

For this reason, when you declare an external relation to ASCEND, you have the opportunity to provide a function which can calculate exact values of these partial derivatives (which end up as elements in the system's Jacobian matrix).

To set an element within the jacobian (derivative) matrix, use the following 'row-major' syntax (see ascend/compiler/extfunc.h):

for(i = 0; i< noutputs; ++i){

  for(j = 0; j< ninputs; ++j){
    jacobian[i*ninputs + j] = /* insert value of ∂F_i/∂x_j here */

  }
}

See also

• FPROPS (documentation of the 'completed' library)