4.5.7. User Subroutines
Warning
The use of user subroutines should be considered a last resort, especially for production simulations - most custom capabilities can be implemented using GPU-compatible, well-tested capabilities like string functions, user expressions, and multi-dimensional tables.
While most customization can be accomplished using String Functions, Aria still supports the inclusion of several types of user subroutines. It is worth noting that simulations that use user subroutines will not run fully on the GPU on platforms where they are available and there may be a significant performance penalty relative to using string functions in those cases.
Aria offers support of C and Fortran user subroutines through fixed-interface signatures. Callback data query requests from the subroutines are also supported, again through fixed-interface function signatures. It is important to note that C subroutines provide some level of variable-type safety and bounds checking that is not supported with Fortran subroutines.
Usage of both C and Fortran subroutines within Aria are described in the sections which follow.
4.5.7.1. Building Usersubs
The process for building user subroutines is the same for C and Fortran subroutines. When you invoke the sierra --make command, the specified input file is scanned for imported shared object (.so) files and they are compiled from source files of the same name.
sierra --make aria -i input.i
Note
The shared object generally does not have to be re-compiled between Aria runs unless you change the contents of the C or F file, however it is a good practice to re-compile them when you change what version of Aria you are running.
Subroutines can also be compiled using the buildsub command, which should produce an equivalent output
buildsub -e aria -i input.i
4.5.7.1.1. Static Linking
Fortran subroutines also have the option of being statically linked instead using a shared object file. To use this option, include the user subroutine file = myfile.F line in the Aria input file (in the same scope as the load user plugin file command). When this option is used, a new Aria executable is created in a UserSubsProject/bin directory created at the path where the sierra --make command was invoked.
4.5.7.2. C Subroutines
A C user subroutine must be a text file of suffix .C that will be compiled to a shared object file of the same name with the .so extension. For example, a file named Aucal_density.C would be compiled to a file named Aucal_density.so and referenced from the Aria input file using the line command
load user plugin file ./Aucal_Density.so
Note that, since Aucal_density.so is a filename, its input file name is case-sensitive, and considerations must be made concerning its location path.
An example C file to provide a subroutine for density is shown below.
#include "Aria_Calore_User_Sub_Support.h"
int
aucal_density(
UserQuery & user_query,
ElementIds elemid,
int n,
int nint,
int spatial_dimension,
CoordinateElementIntegrationPoints coords,
TemperatureElementIntegrationPoints temp,
ScalarElementIntegrationPoints rho)
{
// iterate elements and integration points
for (int j = 0; j < nint; ++j) {
for (int i = 0; i < n; ++i) {
rho(i,j) = 0.1;
}
}
return 0;
}
RegisterFunction(ElemUserSubC, aucal_density, "bulk_density");
In the code above, the actual function (named aucal_density) is declared with the required signature for an element scalar. The function is then registered for use in Aria with the following command (which must be after the function declaration)
RegisterFunction(ElemUserSubC, aucal_density, "bulk_density");
The registration function takes three arguments:
ElemUserSubC- The signature of the function being registered.
aucal_density- The name of the function being registered.
"bulk_density"- A name to refer to the subroutine by in the Aria input file.
Once registered, the subroutine can be referred to in Aria using the registered name and type:
density = calore_user_sub name = bulk_density type=element
Subroutine registration can also be done in a single function in the user C file, using
extern "C" void bulk_prop_subs()
{
RegisterFunction(ElemUserSubC, aucal_density, "bulk_density");
}
When using this method, the command to load the plugin in the Aria input file must refer to this registration function (unless the registration function uses the default name of dl_register)
load user plugin file ./Aucal_Density.so using function bulk_prop_subs
4.5.7.2.1. Signatures
There are a limited number of recognized subroutine signatures supported for Aria C user functions, described in the following sections. Scaled signatures are used to combine the basic user subroutine functionality with field variable scaling. With the exception of the node signature, these signatures are used for subroutines that operate on quadrature point values. The primary differences between the different subroutine types are the signature (the type and order of arguments that the function takes) and the signature type (the first argument passed to the RegisterFunction call).
Node Signature
The node subroutine signature is generally used for Dirichlet boundary conditions.
/**
* Prototype for the nodal signature.
*
* user_query: query function
* node_id: array of global ids of the nodes
* num_nodes: number of nodes in the array
* spatial_dimension: number of nodes in the array
* coordinate_array: array of coordinates of the current nodes
* node_data_array: array of value to be calculated
*/
int
AcalNodeUserSubC(
UserQuery & user_query,
int node_id,
int num_nodes,
int spatial_dimension,
CoordinateNodes coordinate_array,
DataNodes node_data_array);
It would be registered in the C file using the NodeUserSubC signature argument, as shown in the following example subroutine
#include "Aria_Calore_User_Sub_Support.h"
int
dirichsub(
UserQuery& user_query,
int nodeid,
int nnode,
int spatial_dimension,
CoordinateNodes coords,
DataNodes temperature)
{
for (int i = 0; i < nnode; ++i) {
double ynode = coords(1,i);
if ( ynode < 0.5 ) {
temperature(i) = 1.5;
} else {
temperature(i) = 0.1;
}
}
return 0;
}
RegisterFunction(NodeUserSubC, dirichsub, "dirich_bc");
A node subroutine would be invoked in the Aria input file by setting type=node
BC Dirichlet for Temperature on surface_4 = Calore_User_Sub name=dirich_bc type=node
Element Signature
The element signature is used to define material properties, flux boundary conditions, and other boundary properties.
/**
* Prototype for the scalar element signature.
*
* user_query: query function
* elem_ids: pointer to the global id of the current element
* num_elements: number of elements in the workset
* nint: number of gauss points
* coordinate_array: coordinates array of the current node (3, nelem, nint)
* temperature_array: temperature array (nelem, nint)
* elem_data_array: array of values to be calculated (nelem, nint)
*/
int
AcalElemUserSubC(
UserQuery & user_query,
ElementIds elem_ids,
int num_elements,
int nint,
int spatial_dimension,
CoordinateElementIntegrationPoints coordinate_array,
TemperatureElementIntegrationPoints temperature_array,
ScalarElementIntegrationPoints element_data_array)
It would be registered in the C file using the ElemUserSubC signature argument, as shown in the following example subroutine
#include "Aria_Calore_User_Sub_Support.h"
int
aucal_density(
UserQuery & user_query,
ElementIds elemid,
int nelem,
int nint,
int spatial_dimension,
CoordinateElementIntegrationPoints coords,
TemperatureElementIntegrationPoints temp,
ScalarElementIntegrationPoints rho)
{
for (int j = 0; j < nint; ++j) {
for (int i = 0; i < nelem; ++i) {
rho(i,j) = 0.1 + temp(i,j);
}
}
return 0;
}
RegisterFunction(ElemUserSubC, aucal_density, "aucal_density");
Scaled Element Signature
This variant of the element subroutine employs an additional field variable that is supplied through the function signature. Here the field variable is used to perform some logical decision on how the result is assigned and the name of the Field variable to be used is specified from the input file. The call signature for scaled scalar element-based user subroutines is
/**
* Prototype for the scaled scalar element signature.
*
* user_query: query function
* elem_ids: pointer to the global id of the current element
* num_elements: number of elements in the workset
* nint: number of gauss points
* coordinate_array: coordinates array of the current node (3, nelem, nint)
* temperature_array: temperature array (nelem, nint)
* scaling_array: scaling array (nelem, nint)
* elem_data_array: array of values to be calculated (nelem, nint)
*/
int
AcalScaledElemUserSubC(
UserQuery & user_query,
ElementIds elem_ids,
int num_elements,
int nint,
int spatial_dimension,
CoordinateElementIntegrationPoints coordinate_array,
TemperatureElementIntegrationPoints temperature_array,
ScalarElementIntegrationPoints scaling_array,
ScalarElementIntegrationPoints element_data_array)
and it is registered using the ScaledElemUserSubC signature argument.
Scaled user subroutines are typically used in Calore block-style commands
begin convective flux boundary condition convect
add surface surface_2
reference temperature = 180.0
scaled convective coefficient subroutine is aucal_htcoef fieldname = c0
end convective flux boundary condition convect
One should note that the number of arguments for the “Scaled Element Signature” and “Reference Temperature Element Signature” are identical. Attempts to use a user subroutine in the wrong context will be detected while parsing the input file.
Reference Temperature Element Signature
In some problems the surface heat transfer coefficient is determined to be a function of reference temperature rather than just the surface temperature. One example of this would be when modeling problems of natural convection where one computes a Rayleigh number based upon the difference between surface temperature and reference temperature to determine the local Nusselt number and a local heat transfer coefficient. In these cases it is sometimes more convenient to utilize an existing model for reference temperature and use those values to evaluate the heat transfer coefficient based upon additional criteria within a user subroutine. Here the interface for the Element Signature is modified to also supply interpolated values of the heat transfer coefficient at the quadrature points with the following call signature
/**
* Prototype for the scalar reference temperature dependent
* heat transfer coefficient signature.
*
* user_query: query function
* elem_ids: pointer to the global id of the current element
* num_elements: number of elements in the workset
* nint: number of gauss points
* coordinate_array: coordinates array of the current node (3, nelem, nint)
* ref_temperature_array: reference temperature array (nelem, nint)
* temperature_array: temperature array (nelem, nint)
* elem_data_array: array of values to be calculated (nelem, nint)
*/
int
AcalElemUserSubHtC(
UserQuery & user_query,
ElementIds elem_ids,
int num_elements,
int nint,
int spatial_dimension,
CoordinateElementIntegrationPoints coordinate_array,
ScalarElementIntegrationPoints ref_temperature_array,
TemperatureElementIntegrationPoints temperature_array,
ScalarElementIntegrationPoints element_data_array)
and it is registered using the ElemUserSubHtcC signature argument. This would typically be used in a Calore-style convective BC
begin convective flux boundary condition convect
add surface surface_2
reference temperature = 180.0
ref temp convective coefficient subroutine is aucal_htcoef
end convective flux boundary condition convect
Scaled Reference Temperature Element Signature
The scaled version of the reference temperature element signature has an additional argument for the scale field
/**
* Prototype for the scaled scalar reference temperature dependent
* heat transfer coefficient signature.
*
* user_query: query function
* elem_ids: pointer to the global id of the current element
* num_elements: number of elements in the workset
* nint: number of gauss points
* coordinate_array: coordinates array of the current node (3, nelem, nint)
* ref_temperature_array: reference temperature array (nelem, nint)
* temperature_array: temperature array (nelem, nint)
* scale_elem_data_array: scaling value array (nelem, nint)
* elem_data_array: array of values to be calculated (nelem, nint)
*/
int
AcalScaledElemUserSubHtcc(
UserQuery & user_query,
ElementIds element_ids,
int num_elements,
int num_integration_points,
int spatial_dimension,
CoordinateElementIntegrationPoints coordinate_array,
ScalarElementIntegrationPoints ref_temperature_array,
TemperatureElementIntegrationPoints temperature_array,
ScalarElementIntegrationPoints scaling_array,
ScalarElementIntegrationPoints flux_array)
and is registered using the ScaledElemUserSubHtcC signature type. This would typically be used in a Calore-style convective BC
begin convective flux boundary condition convect
add surface surface_2
reference temperature = 180.0
scaled ref_temp convective coefficient subroutine = aucal_htcoef fieldname = c0
end convective flux boundary condition convect
4.5.7.2.2. User Data
All the subroutine signatures take a UserQuery object, which can be used to retrieve a variety of data from the Aria input file and simulation state. If you add a data block in the input file (at the sierra level)
Begin data block mydata
real foo = (293. 433. 533. 623. 773. 1273. )
real bar = 12.0
End
Then you can retrieve these values in the user subroutine using the UserQuery object. The data block can be “bound” to different scopes to allow access without using it’s name. For example, to bind it to the region scope you would use the following command in the Aria Region.
Use data block mydata
Similarly, to bind it to an Aria material you would use that command in the Aria Material block.
A list of the common methods that can be called on a UserQuery object are
uq.currentTime()- Return the current time
uq.currentTimeStepSize()- Return the current time step size
uq.getUserRealMaterialData(arr, var_name)- Retrieve a data array from a data block bound with aUse data blockcall in the Aria Material block.
uq.getUserRealRegionData(arr, var_name)- Retrieve a data array from a data block bound with aUse data blockcall in the Aria Region block.
uq.getUserRealInstanceData(arr, var_name)- Retrieve a data array from a data block bound with aUse data blockcall in the relevant IC/BC/Source block.
uq.getUserRealData(arr, var_name, data_block_name)- Retrieve a data array from a data block not bound with aUse data blockcall, instead using the data block name directly in the function call.
uq.getGlobalVariable(var, name)- Retrieve the value of a global variable by name
uq.getGlobalVariableValue<double>(name)- Return a global variable value by name (specify the variable type as the template parameter)
For example, to retrieve data from a block bound to the material scope, you would use
RealArray1d foo(6);
RealArray1d bar(1);
user_query.getUserRealMaterialData(foo, "foo");
user_query.getUserRealMaterialData(bar, "bar");
double t = user_query.currentTime();
double p = user_query.getGlobalVariableValue("enc_pressure");
4.5.7.3. Fortran Subroutines
An Aria Fortran user subroutine must be a text file of suffix .F or .f, in FORTRAN77 format. The syntax for dynamically linked Fortran subroutines is the same as for C ones:
load user plugin file Aucal_ConvFlux.so using function conv_flux_register
where conv_flux_register is the name of the Fortran subroutine used for user subroutine registration.
For statically linked subroutines, the additional user subroutine file = myfile.F command should be included in the input file.
4.5.7.3.1. Signatures
There are only two valid Fortran subroutine signatures supported in Aria - node and element routines. A Fortran subroutine requires a registration function to be defined to register functions for use in Aria. Like the C subroutine registration, the call to fmwkusersub takes three arguments, the signature, the function, and a string to use to refer to it in the Aria input file.
c=============================================================
c Must have a registration subroutine
c
subroutine material_register
implicit none
external aucal_fort_elem_sub_f
external my_rho
call fmwkusersub(aucal_fort_elem_sub_f, my_rho,
& 'my_rho')
end
c=============================================================
c Actual density calculation
c
subroutine my_rho(objid,nelem,nint,coords,temperature,
& rho,time,rdat,idat,ierror)
implicit none
integer nelem,nint
integer objid(nelem)
double precision coords(3,nelem,nint)
double precision temperature(nelem,nint)
double precision rho(nelem,nint)
double precision time
double precision rdat
double precision idat
integer ierror
integer i, j
ierror = 0
do i=1,nelem
do j=1,nint
rho(i,j) = 0.1 + 0.001*temperature(i,j)
end do
end do
return
end
This registration function is then referred to in the Aria input file using
load user plugin file mysub.so using function material_register
and the subroutine can be referenced in a material property as
density = Fortran sub_name = my_rho
Node Signature
An example nodal Fortran routine is shown below - for applying a Dirichlet BC.
c
c Subroutine to handle coordinate and time-dependent
c Temperature Dirichlet BC
c
subroutine dirich_test(objid,nnodes,coords,temperature,
& time,rdat,idat,ierror)
implicit none
integer nnodes
integer objid
double precision coords(2,nnodes)
double precision temperature(nnodes)
double precision time
integer ierror
double precision rdat(1)
integer idat(1)
c
c INPUT ARGUMENTS:
c objid(nelem): array of mesh object IDs.
c nelem: number of elements in the workset
c nint: Number of integration points per element
c coords(2,n,nint): array of the coordinates (x,y,z)
c of the Gauss points
c rdat() real parameter data
c idat() integer data
c
c OUTPUT ARGUMENTS:
c temperature: temperature at the node
c ierror: user defined error code,
c zero means success.
c
integer i, j
ierror = 0
do i=1,nnodes
if( coords(2,i) .lt. 0.5 ) then
temperature(i) = 1.5
else
temperature(i) = 0.1
end if
end do
return
end
The command to use this in the Aria input file would be
Begin Temperature Boundary Condition bcblock1
Temperature fortran subroutine = dirich_test
Add surface surface_2
End Temperature Boundary Condition bcblock1
Element Signature
A Fortran element subroutine can be used to apply a flux, define a source, or define material properties. The example below shows how to use it to define a heat flux boundary condition
c
c Subroutine to handle coordinate, time and temperature
c dependent flux
c
subroutine flux_test(objid,nelem,nint,coords,temperature,
& flux,time,rdat,idat,ierror)
implicit none
integer nelem, nint
integer objid(nelem)
double precision coords(3,nelem,nint)
double precision temperature(nelem,nint)
double precision flux(nelem,nint)
double precision time
integer ierror
double precision rdat(2)
integer idat(3)
c
c INPUT ARGUMENTS:
c objid(nelem): array of mesh object ID's.
c nelem: number of elements in the workset
c nint: Number of integration points per element in the workset
c coords(3,n,nint): array of the coordinates (x,y,z) of the Gauss points
c temperature(n,nint): temperature at the Gauss points
c rdat() real parameter data
c idat() integer data
c
c OUTPUT ARGUMENTS:
c flux: array of length ( n, nint) containing
c the heat flux coefficient
c at integration points
c ierror: user defined error code for calore to test.
c zero means success.
c
integer i, j
ierror = 0
do i=1,nelem
do j=1,nint
if( coords(3,i,j) .gt. 4.0 ) then
flux(i,j) = 4000.0
else
flux(i,j) = 0.0
end if
end do
end do
return
end
This is referenced in the Aria input file using the following commands
Begin Heat Flux Boundary Condition C1
Add Surface surface_1
Flux Fortran Subroutine = ht_flux_test
integer data 8, 12, 14
real data 1.0 2.0
end