# 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 a `Use data block` call in the Aria Material block.

> - `uq.getUserRealRegionData(arr, var_name)` - Retrieve a data array from a data block bound with a `Use data block` call in the Aria Region block.

> - `uq.getUserRealInstanceData(arr, var_name)` - Retrieve a data array from a data block bound with a `Use data block` call 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 a `Use data block` call, 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