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
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/ipykernel_23326/474769008.py in <module>
----> 1 get_ipython().run_cell_magic('fortran', '', 'module particles_f90\n\n implicit none\n \n real(8), dimension(:), allocatable :: positions\n real(8), dimension(:), allocatable :: velocities\n \ncontains\n subroutine init_particles( n )\n \n integer, intent(in) :: n\n \n integer :: i\n \n if (.not. allocated(positions)) then\n allocate(positions(n))\n end if\n positions = [(i, i = 1, n, 1)]\n if (.not. allocated(velocities)) then\n allocate(velocities(n))\n end if\n velocities = 1.0\n\n end subroutine init_particles\n \n subroutine push_particles( dt )\n \n real(8), intent(in) :: dt\n \n positions = positions + dt * velocities\n \n end subroutine push_particles\nend module particles_f90\n')
/usr/local/lib/python3.8/site-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
2417 with self.builtin_trap:
2418 args = (magic_arg_s, cell)
-> 2419 result = fn(*args, **kwargs)
2420 return result
2421
/usr/local/lib/python3.8/site-packages/decorator.py in fun(*args, **kw)
230 if not kwsyntax:
231 args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)
233 fun.__name__ = func.__name__
234 fun.__doc__ = func.__doc__
/usr/local/lib/python3.8/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
185 # but it's overkill for just that one bit of state.
186 def magic_deco(arg):
--> 187 call = lambda f, *a, **k: f(*a, **k)
188
189 if callable(arg):
/usr/local/lib/python3.8/site-packages/fortranmagic.py in fortran(self, line, cell)
377 verbosity=args.verbosity)
378 if res != 0:
--> 379 raise RuntimeError("f2py failed, see output")
380
381 self._code_cache[key] = module_name
RuntimeError: f2py failed, see output
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
particles.positions
particles.index = 0
particles.position
particles.index = 1
particles.position
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
for p in particle_array:
print(p.position)
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)
type(geom)
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];
}
}
%%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
import sys
!{sys.executable} -m numpy.f2py --quiet -c cfuncts.c cfuncts.pyf -m cfuncts
import numpy as np
import cfuncts
print(cfuncts.push_particles.__doc__)
n = 10
dt = 0.1
x = np.arange(n, dtype="d")
v = np.ones(n, dtype="d")
cfuncts.push_particles( x, v, dt)
x
References¶
f2py documentation https://docs.scipy.org/doc/numpy/f2py/
Transparents E. Sonnendrucker http://calcul.math.cnrs.fr/Documents/Journees/dec2006/python-fortran.pdf
Documentation Sagemath http://doc.sagemath.org/html/en/thematic_tutorials/numerical_sage/f2py.html
Hans Petter Langtangen : Python Scripting for Computational Science.