OK. Here's a reduced form of the code that still produces errors from
time to time. Sometimes it finishes without problems, sometimes the
variable values become 1.#QNAN (when I print out C_1b), and sometimes I
get an error at line:
Res.append(eq_Ab.sweep(var=Ab, boundaryConditions=Ab_BCs, dt=dt))
and the error is: IndexError: id1 is not the same size as b (line 322,
sparseMatrix.py)
--------- code -------------
from fipy import *
import sys,os,re
from math import exp
# sphere radius (um),
R = 300
nx = 60
dx = R * 1. / nx
D = 14
init_R_1a = 0
init_R_1b = 540
Ab0 = 100
k_on_a = 0
k_off_a = 0
k_on_b = 1e-4
k_off_b = 1e-5
k_x_a = 0
k_negx_a = 0
k_x_b = 0
k_negx_b = 0
epsilon = 0.12
tao_alpha = 3 * 3600
tao_beta = 50 * 3600
alpha = 0.27
beta = 0.73
######## start of variables
mesh = Grid1D(nx=nx, dx=dx)
Ab = CellVariable(name="unbound antibody", mesh=mesh, value=0., hasOld=1)
R_1a = CellVariable(name="unbound receptor a on cell 1", mesh=mesh,
value=init_R_1a, hasOld=1)
R_1b = CellVariable(name="unbound receptor b on cell 1", mesh=mesh,
value=init_R_1b, hasOld=1)
C_1a = CellVariable(name="bounded receptor a-antibody complex on cell
1", mesh=mesh, value=0., hasOld=1)
C_1b = CellVariable(name="bounded receptor b-antibody complex on cell
1", mesh=mesh, value=0., hasOld=1)
C_1ab = CellVariable(name="bounded receptor a+b-antibody complex on cell
1", mesh=mesh, value=0., hasOld=1)
r = mesh.getCellCenters()[...,0]
r2 = CellVariable(mesh=mesh, value=r**2)
time = Variable()
######## start of equations
eq_Ab = TransientTerm(coeff = r2) == ImplicitDiffusionTerm(coeff = D*r2) \
- ImplicitSourceTerm(coeff = r2 * (k_on_a/epsilon) *
R_1a) \
- ImplicitSourceTerm(coeff = r2 * (k_on_b/epsilon) *
R_1b) \
+ (r2 * k_off_a * C_1a) + (r2 * k_off_b * C_1b)
eq_C_1a = TransientTerm() == ((k_on_a/epsilon) * Ab * R_1a) -
ImplicitSourceTerm(coeff = k_off_a) \
- ImplicitSourceTerm(coeff =
k_x_b * R_1b) \
+ k_negx_b * C_1ab
eq_C_1b = TransientTerm() == ((k_on_b/epsilon) * Ab * R_1b) -
ImplicitSourceTerm(coeff = k_off_b) \
- ImplicitSourceTerm(coeff =
k_x_a * R_1a) \
+ k_negx_a * C_1ab
eq_C_1ab = TransientTerm() == k_x_a * C_1b * R_1a + k_x_b * C_1a * R_1b \
-
ImplicitSourceTerm(coeff = k_negx_a) \
-
ImplicitSourceTerm(coeff = k_negx_b)
eq_R_1a = TransientTerm() == - ImplicitSourceTerm(coeff =
(k_on_a/epsilon) * Ab) + (k_off_a * C_1a) - ImplicitSourceTerm(coeff =
(k_x_a) * C_1b) + (k_negx_a * C_1ab)
eq_R_1b = TransientTerm() == - ImplicitSourceTerm(coeff =
(k_on_b/epsilon) * Ab) + (k_off_b * C_1b) - ImplicitSourceTerm(coeff =
(k_x_b) * C_1a) + (k_negx_b * C_1ab)
####### resolution for sweeping variables
desiredRes = 1e-2
endTime = 21 * 3600
# timestep size
dt = 600
for step in range(1):
time.setValue(0.)
while time() < endTime:
print 'time: ' + str(time())
# remember to updateOld at beginning of each timestep
Ab.updateOld()
R_1a.updateOld()
R_1b.updateOld()
C_1a.updateOld()
C_1b.updateOld()
Ab_p = Ab0 * (alpha * 2**(-time() / tao_alpha) + beta *
2**(-time() / tao_beta))
Ab_r = epsilon * Ab_p
print "Ab_r: %s" % Ab_r
# boundary condition for [Ab]
Ab_BCs = (FixedValue(faces=mesh.getFacesRight(),value=Ab_r), \
FixedFlux(faces=mesh.getFacesLeft(),value=0.))
# keep sweeping all variables until residuals are all smaller
than desiredRes
while 1:
Res = []
Res.append(eq_Ab.sweep(var=Ab, boundaryConditions=Ab_BCs,
dt=dt))
Res.append(eq_C_1a.sweep(var=C_1a, dt=dt))
Res.append(eq_C_1b.sweep(var=C_1b, dt=dt))
Res.append(eq_R_1a.sweep(var=R_1a, dt=dt))
Res.append(eq_R_1b.sweep(var=R_1b, dt=dt))
if len(filter(lambda x: x > desiredRes, Res))==0:
break
print "max(R_1b): %s, max(C_1ab): %s" % (max(R_1b), max(C_1ab))
print "C_1b: %s" % C_1b
time.setValue(time() + dt)
-------- code -------------------
Daniel Wheeler wrote:
Hi Elizabeth,
On 8/15/07, Elizabeth Tseng <[EMAIL PROTECTED]> wrote:
Hi,
I was using the PDE solver with about 10 variables and a system of
equations. I noticed that sometimes, after sweeping, the variable values
become 1.#QNAN. What confuses me even more is that it shows up randomly.
If I re-run the time steps, sometimes it will go all the way to the end
without giving me NANs, and sometimes, after a few time steps, 1.#QNAN
pops up.
Very strange. If you reduce the input file to the simplest form possible that
still produces the bug and post it, we could take a look and see if there isn't
something that is obviously wrong.
One of the equations I'm solving is like this:
eq_C_1b = TransientTerm() == ((k_on_b/epsilon) * Ab * R_1b) -
ImplicitSourceTerm(coeff = k_off_b) \
-
ImplicitSourceTerm(coeff = k_x_a * R_1a) \
+ k_negx_a * C_1ab
eq_C_1ab = xxx.....
Because I'm actually testing different control sets, sometimes, some
variables (ex: C_1ab) and constants (ex: k_x_a) will be purposely set to
0, but the equations (ex: eq_C1_ab) still gets sweeped and the variable
updated. Could this be the cause?
You can set the variable c_1ab to zero and then solve eq_C_1ab for C_1ab.
That shouldn't be a problem at all. Likewise setting any coefficient to zero
shouldn't be a problem.
Cheers
Thanks,
Elizabeth
