coupling fipy with scipy.integrate.odeint

[Fangohr H.]

Dear all,

we would like to couple fipy with scipy.integrate.odeint.

We attach an example (odeint_example.py) that demonstrates how scipy's odeint 
can be used by providing a function for the right hand side of an equation of 
type

dy
-- = rhs(y,t)
dt

We also attach a diffusion example (filename question.py) which solves a 
diffusion problem with fipy's implicit time integration scheme.

The third file (filename question_odeint.py) contains an attempt to reformulate 
this so that scipy.integrate.odeint is used for the time integration.

We think the current stumbling block is that we don't know how to ask fipy to 
simply compute the dy/dt term, so that we can hand this back to odeint() in the 
function rhs().

Hopefully the question makes sense.

Thank you in advance,

Hans



#Harmonic oscillator example to demonstrate integration of two couple ODEs with odeint



import numpy
from scipy.integrate import odeint


def rhs(y,t):
    x=y[0] #position
    v=y[1] #velocity

    k=1    #spring stiffness
    dxdt = v
    dvdt = -k*x 

    return numpy.array([dxdt,dvdt])


y0 = numpy.array([ 1,0]) #position 0 (equilibrium is at 0), velocity 0

ts = numpy.arange(0,10,1) # time steps for which solution is desired

ys = odeint(rhs,y0,ts)

for t,y in map(None,ts,ys):
    print "time=%7f, pos=%7f, vel=%7f" % (t,y[0],y[1])

#!/usr/bin/env python

import fipy

Ly = 20
Lx = Ly
nx = 20
ny = nx
dx = 1.
dy = dx

mesh = fipy.Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

phi = fipy.CellVariable(name = "solution variable",
                        mesh = mesh,
                        value = 0.0)

#Basic equation
D = 1
eq = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=D) 

#boundary conditions
phi_outside = 0.1
x, y = mesh.getFaceCenters()
facesTopLeft = ((mesh.getFacesLeft() & (y > Ly / 2))
                | (mesh.getFacesTop() & (x < Lx / 2)))
facesBottomRight = ((mesh.getFacesRight() & (y < Ly / 2))
                    | (mesh.getFacesBottom() & (x > Lx / 2)))

BCs = (fipy.FixedValue(faces=facesTopLeft, value=phi_outside),
       fipy.FixedValue(faces=facesBottomRight, value=-phi_outside))

if __name__ == '__main__':
     pass
     #viewer = fipy.Viewer(vars=phi, datamin=-phi_outside*1., datamax=phi_outside*1,colorbar='vertical')
     #viewer.plot()

#and solve the equation by repeatedly looping in time:

explicit_max_timeStepDuration = 0.9 * dx*dy / (2 * D)
timeStepDuration = 10 * explicit_max_timeStepDuration

steps = 20

for step in range(steps):
    eq.solve(var=phi,
             boundaryConditions=BCs,
             dt=timeStepDuration)
    tmp = phi.getValue()
    #print one value to compare with new odeint-timeintegration code (question_odeint.py)
    print step,tmp[42]

    if __name__ == '__main__':
        #viewer.plot()
         pass




#!/usr/bin/env python

import fipy
from scipy.integrate import odeint

Ly = 20
Lx = Ly
nx = 20
ny = nx
dx = 1.
dy = dx

mesh = fipy.Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

phi = fipy.CellVariable(name = "solution variable",
                        mesh = mesh,
                        value = 0.0)

#Basic equation
D = 1
eq = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=D) 


#boundary conditions
phi_outside = 0.1
x, y = mesh.getFaceCenters()
facesTopLeft = ((mesh.getFacesLeft() & (y > Ly / 2))
                | (mesh.getFacesTop() & (x < Lx / 2)))
facesBottomRight = ((mesh.getFacesRight() & (y < Ly / 2))
                    | (mesh.getFacesBottom() & (x > Lx / 2)))

BCs = (fipy.FixedValue(faces=facesTopLeft, value=phi_outside),
       fipy.FixedValue(faces=facesBottomRight, value=-phi_outside))

if __name__ == '__main__':
     #viewer = fipy.Viewer(vars=phi, datamin=-phi_outside*1., datamax=phi_outside*1,colorbar='vertical')
     #viewer.plot()
     pass

#and solve the equation by repeatedly looping in time:

explicit_max_timeStepDuration = 0.9 * dx*dy / (2 * D)
timeStepDuration = 10 * explicit_max_timeStepDuration


def rhs( y, t):
     #set phi to state vector y
     phi.setValue(y)

     #compute time derivative of state vector (HOW??? the code below is a guess but doesn't work)
     eq.cacheRHSvector()
     eq.solve(var=phi,             
              boundaryConditions=BCs)
     dydt = eq.getRHSvector()
     return dydt


steps = 20

#define time steps at which we like odeint() to report back statevector
ts = fipy.arange(0,steps*timeStepDuration,timeStepDuration)


#initial value for system of ODEs
y0 = phi.getValue()

#this does the time integration using scipy.integrate.odeint
ys = odeint(rhs,y0,ts)

#ys is a matrix with 20 rows, and each row as 400 entries
n,tmp= ys.shape
assert n == steps
assert tmp == nx*ny

#simply print one value as function of time to compare against working fipy-time integration code (question.py)
for step in range(steps):
     print step,ys[step,42]