Generating Fortran Code with Python Preprocessor

Hello everyone,

does someone have experience with using Python as a preprocessor in a Fortran project? I thereby not only mean to use it for conditional compilation like with fypp, but to prepare or optimize a certain problem with Python analytically, write down optimized Fortran code and finally compile and run the Fortran code.

Specifically, I have a project where I load an ODE with Python, try to make some pre-optimizations in Sympy and write down the resulting RHS and Jacobian to Fortran-Templates. Afterwards I compile and run them.

Can fypp be used in this scenario too? How should a code be organized in that case and could fpm still be used reasonably? Are there alternatives to consider?

Thanks a lot!

3 Likes

Hi Jannes, welcome to Fortran Discourse!

I’ve used SymPy in a few cases for simplification of common sub-expressions and generating Jacobian matrices analytically. I never had the need to automate it and make it part of the build process. For my purposes I usually had functions that were known in advance, so I only had to generate the code once. I would usually do this in a Jupyter notebook so I could analyze both the symbolic expressions and the generated Fortran output. An example of code generation can be found in the thread: Code Generation Using Sympy (it may require some fixing due to changes in SymPy).

At the moment preprocessor support, other than the built-in compiler ones, is not available in fpm, but @arteevraina, one of the Fortran-lang GSoC students is making steady progress toward it (feat: added basic preprocess table configuration by arteevraina · Pull Request #715 · fortran-lang/fpm · GitHub).

From your post I understand your ODE’s are always defined in Python. In principle you could set up functions to print the Fortran code given a Python function object. You could then import your Python module into fypp at startup. You can accept inputs from the mixed fypp/Fortran source code, although I cannot foresee right now why would that be needed.

Here’s something I could put together in 45 min for evaluating the derivatives of scalar expressions for moisture sorption isotherms:


# Copyright (c) 2022 Ivan Pribec. All rights reserved.
#
# This work is licensed under the terms of the MIT license.  
# For a copy, see <https://opensource.org/licenses/MIT>.

import functools
from jinja2 import Template

import sympy as sp
from sympy.printing.fortran import fcode

#
# Jinja2 template for a Fortran subroutine
#
isotherm_subroutine = '''
{{ doc }}
subroutine {{ name }}( {{ args }} )
{{ argsdeclr }}
{{ body }}
end subroutine

'''
isotherm_tm = Template(isotherm_subroutine)


def apply_cse_and_return_fcode(expr,assign_to,source_format='free',standard=95):
    """Apply common subexpression simplification and return Fortran code
    
    Input:
        expr - a list of expressions we would like to evaluate
        assign_to - the symbols we would like to assign them to (either a tuple of strings or sp.symbols)
    
    Output:
        String containing subroutine or function body
    """
    
    rvar, rexpr = sp.cse(sp.simplify(expr),
        order='none',
        list=False)

    ret = '  real(dp) :: ' + ','.join(str(var) for (var, _) in rvar) + '\n\n'

    for var, var_expr in rvar:
        ret += '  ' + fcode(var_expr,assign_to=var,source_format=source_format,standard=standard) + '\n'
    
    ret += '\n'

    for var, var_expr in zip(assign_to,rexpr):
        ret += '  ' + fcode(var_expr,assign_to=var,source_format=source_format,standard=standard) + '\n'

    return ret



def print_isotherm_subroutine(isotherm, params):
    """Print a Fortran subroutine for a given isotherm expression

    Input:
        isotherm - a decorated isotherm function
        params - a tuple of parameter symbols
    """

    X = sp.Symbol('X')
    T = sp.Symbol('T')

    input_var = (str(X), str(T))
    output_var = ('aw', 'aw_X', 'aw_T')        
    params_var = tuple(str(par) for par in params)

    name = isotherm.__name__

    doc = '\n'.join('!> '+line for line in isotherm.__doc__.split('\n'))

    args = ', '.join(output_var + input_var + params_var)

    argsdeclr = '  real(dp), intent(out) :: '
    argsdeclr += ', '.join(output_var) + '\n'
    
    argsdeclr += '  real(dp), intent(in) :: '
    argsdeclr += ', '.join(input_var) + '\n'
    
    argsdeclr += '  real(dp), intent(in) :: '
    argsdeclr += ', '.join(params_var) + '\n'

    output = isotherm(X,T,*params)
    body = apply_cse_and_return_fcode(output, output_var)


    subroutine = isotherm_tm.render(
        doc=doc,
        name=name,
        args=args,
        argsdeclr=argsdeclr,
        body=body)

    print(subroutine)

def isotherm(func):
    """Decorator for isotherm expressions

    Uses SymPy to find symbolic derivatives of the isotherm
    expression with respect to moisture (X) and temperature (T)
    """

    @functools.wraps(func)
    def wrapper(*args):

        X, T, *params = args

        expr = func(X,T,*params)

        # Generate derivative expressions
        dX = sp.simplify(sp.diff(expr,X))
        dT = sp.simplify(sp.diff(expr,T))

        return expr, dX, dT

    # Extend docstring with return values
    wrapper.__doc__ += """Returns:
        aw - water activity
        diff(aw,X) - derivative of water activity with respect to X
        diff(aw,T) - derivative of water activity with respect to T 
    """

    return wrapper

@isotherm
def isotherm_oswin(X,T,*args):
    """Oswin isotherm (4-parameter version)

    Inputs:
        X - moisture concentration on a dry basis (kg m / kg dry solid)
        T - temperature in degrees Celsius
        a1, a2, b1, b2 - fitting parameters

    """

    a1, a2, b1, b2 = args

    a = a1 + a2*T
    b = b1 + b2*T

    aw = (X/a)**(1/b)/(1 + (X/a)**(1/b))

    return aw

@isotherm
def isotherm_henderson(X,T,*args):
    """Henderson isotherm

    Inputs:
        X - moisture concentration on a dry basis (kg m / kg dry solid)
        T - temperature in degrees Celsius
        A, B, C, D - fitting parameters
    """
    
    A, B, C, D = args
    aw = 1 - sp.exp(-A*(T-B)**C)*X**D

    return aw


if __name__ == '__main__':
    
    print_isotherm_subroutine(
        isotherm=isotherm_oswin,
        params=sp.symbols('a1 a2 b1 b2'))

    print_isotherm_subroutine(
        isotherm=isotherm_henderson,
        params=sp.symbols('A B C D'))

The output of running this Python script is:

!> Oswin isotherm (4-parameter version)
!> 
!>     Inputs:
!>         X - moisture concentration on a dry basis (kg m / kg dry solid)
!>         T - temperature in degrees Celsius
!>         a1, a2, b1, b2 - fitting parameters
!> 
!>     Returns:
!>         aw - water activity
!>         diff(aw,X) - derivative of water activity with respect to X
!>         diff(aw,T) - derivative of water activity with respect to T 
!>     
subroutine isotherm_oswin( aw, aw_X, aw_T, X, T, a1, a2, b1, b2 )
  real(dp), intent(out) :: aw, aw_X, aw_T
  real(dp), intent(in) :: X, T
  real(dp), intent(in) :: a1, a2, b1, b2

  real(dp) :: x0,x1,x2,x3,x4,x5,x6,x7

  x0 = T*a2 + a1
  x1 = 1d0/x0
  x2 = X*x1
  x3 = T*b2 + b1
  x4 = 1d0/x3
  x5 = x2**x4
  x6 = x5 + 1
  x7 = x5/x6**2

  aw = x5/x6
  aw_X = x4*x7/X
  aw_T = -x1*x7*(a2*x3 + b2*x0*log(x2))/x3**2

end subroutine


!> Henderson isotherm
!> 
!>     Inputs:
!>         X - moisture concentration on a dry basis (kg m / kg dry solid)
!>         T - temperature in degrees Celsius
!>         A, B, C, D - fitting parameters
!>     Returns:
!>         aw - water activity
!>         diff(aw,X) - derivative of water activity with respect to X
!>         diff(aw,T) - derivative of water activity with respect to T 
!>     
subroutine isotherm_henderson( aw, aw_X, aw_T, X, T, A, B, C, D )
  real(dp), intent(out) :: aw, aw_X, aw_T
  real(dp), intent(in) :: X, T
  real(dp), intent(in) :: A, B, C, D

  real(dp) :: x0,x1,x2

  x0 = -B + T
  x1 = exp(-A*x0**C)
  x2 = X**D*x1

  aw = 1 - x2
  aw_X = -D*X**(D - 1)*x1
  aw_T = A*C*x0**(C - 1)*x2

end subroutine

The idea of this would be to build a Python library of sorption isotherm expressions as symbolic expressions. These could be either lambdified for use in Python, or printed as optimized Fortran subroutines (shown above). Since the Python definitions are quite terse, we essentially avoid the repetitive work of declaring all the variables in Fortran. As you can see also the docstrings can be kept consistent through-out the Python and Fortran files. The thing I’m not entirely sure about it is the mechanism I used to pass parameters to the functions.

I hope this semi-working example can get the discussion going.

3 Likes