Start up a Maxima session (I run Fedora 9, but I'm sure most other distributions have it, I think they even have a windows installer). It's easy to define our matrix operator, invert it, and apply it:

L : matrix(

[1,0,0, -dt, 0, 0],

[0,1,0, 0, -dt, 0],

[0,0,1, 0, 0, -dt],

[0,0,0, 1, a*B[3]/c,-a*B[2]/c],

[0,0,0,-a*B[3]/c, 1, a*B[1]/c],

[0,0,0, a*B[2]/c,-a*B[1]/c, 1] );

Linv : invert(L);

X: [x,y,z,vx,vy,vz];

P : a*[0,0,0,E[1],E[2],E[3]];

RHS : X - P;

LinvRHS : factor( expand( Linv.RHS ) );

Simple as that we have a symbolic representation of our inverted operator acting on our state vector. Now we want Fortran90 expressions.

load("f90");

fname : "back_euler_inv_op.f90";

f : openw(fname);

printf( f, "! These fortran expressions generated automatically from

! the maxima batch file magnetron.mac ~%");

for i : 1 thru length(LinvRHS) do

( printf( f, "x_n1(~d)=",i),

with_stdout( f, f90(LinvRHS[i,1]) ) );

close(f);

There is a small bug in the "f90" module (as of version 5.15) that causes it to use the Fortran77 style continuation when writing out a matrix, so that is why each element must be written out with a separate call to f90().

Well, that's neat, but I need to appease the quiche eaters, so wouldn't it be better if these expressions, and our inner loops (the time integration loop for this case) in Fortran were callable from something flexible and trendy like Python? That's where F2Py comes in. Generating a module for Python is one simple call:

f2py -c back_euler_sub.f90 -m back_euler

Where my subroutine using the Maxima generated expressions is in the file "back_euler_sub.f90" and the module will be named "back_euler".

The python script to use our new module is quite short:

from numpy import *

from back_euler import *

B = zeros( 3 )

E = zeros( 3 )

x_n = zeros( 6 )

x_n1 = zeros( 6 )

q = -4.8032e-10

m = 9.1094e-28

c = 2.9979e10

dt = 3e-10

nt = 50

E[1] = -1.0/3.0

B[2] = 33.6412

X = zeros( (6,nt), dtype='double', order='F' )

X = back_euler(nt,dt,q,m,c,E,B)

Now our low-level number crunching and inner loop is in nice, fast compiled F90 and we can call it from a high-level, object-oriented language to do things like zero finding (as in the magnetron post), or optimization or something else entirely.

## No comments:

## Post a Comment