Dear Jonathan & Daniel,
In addition to my earlier question about automatic differentiation on FiPy data
structures, I debugged my code (attached) to realise that the solution in time
is invariant over dx - this just can't be true. I tried saving the fields and
plotting them after simulation to note that V,m,h and n do not change with
cable length (dx). Can this be due to the specific implementation of the
Neumann boundaries? I read in the FAQ section that fixed flux (Neumann)
boundaries sometimes do not function the way they generally should in FiPy? Why
is this the case? If one doesn't know the analytical solution of the PDE, it
may be pretty difficult to find that out, isn't it?
Best,
Biswa
__author__="Biswa Sengupta"
__date__ ="$Jun 1, 2011 10:47:37 AM$"
if __name__ == "__main__":
from fipy import *
from actinact import *
# from FuncDesigner import *
###########################################################################
# squid axon constants
###########################################################################
PI = 3.141528 # dimensionless
cable_length = 5.0 # in cm
cable_radius = 0.025 # in cm
patch_area = (2*PI*cable_radius*cable_length) # in cm^2
C_m_spec = 1.0 # uF/cm^2
e_Na = 50.0 # in mV
gNa_spec = 120.0 # in mS/cm^2
e_K = -77.0 # in mV
gK_spec = 36.0 # in mS/cm^2
e_l = -54.3 # in mV
gL_spec = 1.5 # in mS/cm^2
ra_spec = 30.0 # in Ohm cm
idc = 0.01 # in mA
idc_boun = idc*(ra_spec/(PI*cable_radius*cable_radius)) # in mV/cm (this doesn't have the usual unit of Amperes because idc here is posed as a boundary having the unit of dv/dx)
lamda = sqrt((cable_radius*(10**3))/(2.0*ra_spec*gL_spec)) # cm
taumem = C_m_spec/gL_spec # ms
###########################################################################
nx = 1000
dx = cable_length/nx
mesh = Grid1D(nx = nx, dx = dx)
restvol = -65.0
V = CellVariable(name="Voltage", mesh=mesh, value=restvol,hasOld=1)
m = CellVariable(name="m", mesh=mesh, value=minfy(restvol),hasOld=1)
h = CellVariable(name="h", mesh=mesh, value=hinfy(restvol),hasOld=1)
n = CellVariable(name="n", mesh=mesh, value=ninfy(restvol),hasOld=1)
valueLeft = idc_boun
valueRight = 0.0
# Neumann boundaries
BCs = (FixedFlux(faces=mesh.getFacesLeft(), value=valueLeft), FixedFlux(faces=mesh.getFacesRight(), value=valueRight))
# ionic currents
ina = gNa_spec*m*m*m*h*(V-e_Na)
ik = gK_spec*n*n*n*n*(V-e_K)
il = gL_spec*(V-e_l)
ionic = ina + ik + il
normcurr = -ionic/gL_spec
print "########################################################################### "
print "Cylinder area: ", patch_area, " cm^2"
print "Space constant: ", lamda*(10**4), " um"
print "Time constant: ", taumem, " ms"
print "Current injected: ", 1000*idc/patch_area, " uA/cm^2"
print "Current boundary term: ", idc_boun , " mV/cm"
print "########################################################################### "
# voltage evolution
volEq = TransientTerm(coeff=taumem) == ImplicitDiffusionTerm(coeff=lamda**2) + ImplicitSourceTerm(coeff=normcurr)
# the (in) activation ODEs described as degenerate PDEs
mEq = TransientTerm(coeff=taum(V)) == ImplicitSourceTerm(minfy(V) - m)
hEq = TransientTerm(coeff=tauh(V)) == ImplicitSourceTerm(hinfy(V) - h)
nEq = TransientTerm(coeff=taun(V)) == ImplicitSourceTerm(ninfy(V) - n)
timeStepDuration = .001
steps = 100
numsweep = 5
viewer = Viewer(vars=V, datamin=-80.0, datamax=80.0)
for i in range(steps):
V.updateOld()
m.updateOld()
h.updateOld()
n.updateOld()
for sweep in range(numsweep):
volEq.sweep(var=V,dt=timeStepDuration)
mEq.sweep(var=m,dt=timeStepDuration)
hEq.sweep(var=h,dt=timeStepDuration)
nEq.sweep(var=n,dt=timeStepDuration)
viewer.plot()
filestr = "/Users/biswasengupta/Desktop/finelPDE/python/cable/src/data/file_id_" + str(i+1) + ".txt"
savetxt(filestr, V, fmt="%12.6G")
On 2 Jun 2011, at 15:01, Jonathan Guyer wrote:
>
> On Jun 2, 2011, at 5:41 AM, Biswa Sengupta wrote:
>
>> As per your advice I tried to pose my set of PDE-ODE system as a PDE + 3
>> degenerate PDEs (i.e., the ODEs without the space term). I have tried to
>> work along the examples for time marching but certainly I am not getting
>> something right. I attach herewith the code that I am trying to use. Can you
>> spot some obvious mistake. I get the following error:
>>
>> 4565:src biswasengupta$ python cable.py
>> Traceback (most recent call last):
>> File "cable.py", line 63, in <module>
>> mEq = TransientTerm(coeff=taum(V)) == ImplicitSourceTerm(minfy(V) - m)
>> File "/Users/biswasengupta/Desktop/finelPDE/python/cable/src/actinact.py",
>> line 3, in taum
>> am = 0.1*(V + 40) / (1-exp(-(V+40)/10.0))
>> File
>> "/Library/Frameworks/EPD64.framework/Versions/7.0/lib/python2.7/site-packages/FiPy-2.1.2-py2.7.egg/fipy/variables/variable.py",
>> line 1146, in __float__
>> return float(self.getValue())
>> TypeError: only length-1 arrays can be converted to Python scalars
>
>
> In your actinact.py, replace all instances of
>
> from math import exp
>
> with
>
> from fipy.tools.numerix import exp
>
> The Python math module only understands single floating point numbers and
> should never be used with FiPy. FiPy's numerix, based on NumPy, can deal with
> arrays of floats.
>
> ----
>
> In cable.py, you should remove
>
> from math import *
>
> for the same reason. Also
>
> from numpy import *
>
> is not necessary because all of the same code (with appropriate
> modifications) is obtained when you import fipy. By importing numpy directly,
> you override some changes that FiPy requires. Finally,
>
> from scipy import *
>
> doesn't do anything. You must import the specific submodules you want from
> SciPy.
>
> ----
>
> Once you're past that, you will encounter a few more issues before your
> script will run, some of them quite obscure:
>
> * TypeError: unsupported operand type(s) for *: 'function' and 'float'
>
> This turns out to be caused by
>
> return ninfy
>
> on actinact.py line 41. Change this to
>
> return ninit
>
> and likewise for minfy and hinfy.
>
> * self.matrix.matvec(other, y)
> ValueError: arg 1 must be a 1-dimensional double array of appropriate size.
>
> This is because V holds an array of integers and it needs to be floats.
> Change cable.py
> line 42 to
>
> restvol = -65.
>
> Probably a good idea to change valueLeft and valueRight (but not absolutely
> necessary).
>
> * Finally, the problem runs, but doesn't show anything. This is because
>
> axonlen = 1
> dx = axonlen/nx
>
> gives a grid spacing of zero because of integer math. Change to
>
> axonlen = 1.
>
> The problem now runs. Whether it gives sensible answers is up to you to
> figure out! 8^)
>
> I've attached my modified scripts.
>
>
>
> <cable.py><actinact.py>
