Python class linked with fortran module

Python class linked with fortran module#

%load_ext fortranmagic

%%fortran
module particles_f90

    implicit none
    
    real(8), dimension(:), allocatable :: positions
    real(8), dimension(:), allocatable :: velocities
            
contains
    subroutine init_particles( n )
    
        integer, intent(in) :: n
                
        integer :: i
        
        if (.not. allocated(positions)) then
            allocate(positions(n))
        end if
        positions = [(i, i = 1, n, 1)]
        if (.not. allocated(velocities)) then
            allocate(velocities(n))
        end if
        velocities = 1.0

    end subroutine init_particles
 
    subroutine push_particles( dt )
        
        real(8), intent(in) :: dt
    
        positions = positions + dt * velocities
        
    end subroutine push_particles
end module particles_f90

The Python class#

class Particles(object):
    
    def __init__(self, n):
        self.index       = 0
        self.numberof    = n
        particles_f90.init_particles( n)
        self.positions  = particles_f90.positions
        self.velocities = particles_f90.velocities
        
    @property 
    def position(self):      
        return self.positions[self.index]
    
    @property 
    def velocity(self):      
        return self.velocities[self.index]

Access to Fortran data from Python#

particles = Particles(10)
particles.velocities 

array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

particles.positions

array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.])

particles.index = 0
particles.position

np.float64(1.0)

particles.index = 1
particles.position

np.float64(2.0)

Create an Iterator class#

class ParticleArray(object):
    
    def __init__(self, particles):
        self.particles = particles
        self.numberof = particles.numberof
        
    def __getitem__(self, index): 
        self.particles.index = index 
        return self.particles
    
    def __len__(self): 
        return self.numberof
    
    def __iter__(self): 
        for i in range(self.numberof):
            self.particles.index = i
            yield self.particles

particle_array = ParticleArray(particles)
particle_array[0].position

np.float64(1.0)

for p in particle_array:
    print(p.position)

1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0

Fortran derived type#

%%fortran
module mesh

implicit none
type :: geometry
    real(8) :: xmin, xmax, dx            ! coordinates of origin and grid size
    integer :: nx                        ! number of grid points
    real(8), dimension(:), pointer :: x  ! coordinates of points
end type geometry

contains

subroutine create(geom, xmin, xmax, nx)

    !f2py integer(8), intent(out) :: geom
    type(geometry), pointer :: geom
    real(8), intent(in) :: xmin, xmax
    integer, intent(in) :: nx
            
    real(8) :: dx
            
    integer :: i
            
    allocate(geom)
    geom%xmin = xmin
    geom%xmax = xmax
    geom%dx = ( xmax - xmin ) / (nx-1) 
    geom%nx = nx
    allocate(geom%x(nx))
    do i=1,nx
        geom%x(i)=geom%xmin+(i-1)*geom%dx
    end do

end subroutine create

subroutine view(geom)
    !f2py integer(8), intent(in) :: geom
    type(geometry), pointer :: geom
    print*, 'nx = ', geom%nx
    print*, geom%xmin, geom%xmax
    print*, geom%x(:)
end subroutine view

subroutine get_size(geom, nx)

    !f2py integer(8), intent(in) :: geom
    type(geometry), pointer :: geom
    integer, intent(out) :: nx
    
    nx = geom%nx
    
end subroutine get_size


end module mesh

geom = mesh.create(0.0, 1.0, 10)

mesh.get_size(geom)

10

type(geom)

int

f2py with C code#

  • Signature file is mandatory

  • intent(c) must be used for all variables and can be set globally.

  • Function name is declared with intent(c)

%rm -rf cfuncts*

%%file cfuncts.c

void push_particles(double* positions, double* velocities, double dt, int n){
    for (int i=0; i<n; i++){
       positions[i] = positions[i] + dt * velocities[i];
        
    }
} 

Writing cfuncts.c

%%file cfuncts.pyf

python module cfuncts 
    interface
        subroutine push_particles(positions, velocities, dt, n) 
            intent(c):: push_particles
            intent(c)
            integer, optional, depend(velocities) :: n = len(velocities)
            real(8), dimension(n),  intent(inplace)  :: positions 
            real(8), dimension(n),  intent(in) :: velocities
            real(8), intent(in) :: dt
        end subroutine push_particles
    end interface
end python module cfuncts

Writing cfuncts.pyf

import sys
!{sys.executable} -m numpy.f2py --quiet -c cfuncts.c cfuncts.pyf -m cfuncts

/usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/f2py/f2py2e.py:723: VisibleDeprecationWarning: 
distutils has been deprecated since NumPy 1.26.x
Use the Meson backend instead, or generate wrappers without -c and use a custom build script
  builder = build_backend(

/usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/cmd.py:66: SetuptoolsDeprecationWarning: setup.py install is deprecated.
!!

        ********************************************************************************
        Please avoid running ``setup.py`` directly.
        Instead, use pypa/build, pypa/installer or other
        standards-based tools.

        See https://blog.ganssle.io/articles/2021/10/setup-py-deprecated.html for details.
        ********************************************************************************

!!
  self.initialize_options()
WARN: Could not locate executable armflang

Removing build directory /tmp/tmpg4mkb5lb

import numpy as np
import cfuncts

print(cfuncts.push_particles.__doc__)

push_particles(positions,velocities,dt,[n])

Wrapper for ``push_particles``.

Parameters
----------
positions :  rank-1 array('d') with bounds (n)
velocities : input rank-1 array('d') with bounds (n)
dt : input float

Other Parameters
----------------
n : input int, optional
    Default: len(velocities)

n = 10
dt = 0.1
x = np.arange(n, dtype="d")
v = np.ones(n, dtype="d")
cfuncts.push_particles( x, v, dt)

x

array([0.1, 1.1, 2.1, 3.1, 4.1, 5.1, 6.1, 7.1, 8.1, 9.1])

References#