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.