Not for me. If you're seeing something else, then it may come down to either
boundary conditions or parameter values. We probably need to see a minimal
script that exhibits what you're seeing.
I send you attached a minimal script that shows one of the problems
I'm having. Now I'm solving two uncoupled equations. This script
produces the right solution for C (everything is zero) but B starts at
B_0=2 and despite having the boundary conditions fixed at B=B_0 the
solution evolves to B=1 (which violates the boundary conditions).
Am I doing something wrong there?
I also think that the nonlinear terms should be defined as ImplicitSource
terms, but I don't know which variable should I use in the definition, since
one of the terms depends both on B and C.
I would be inclined to make the eqC implicit in B and eqB implicit in C in
order to couple the equations together.
Can you explain me how to do that?
Thanks,
Adrian.
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 9 12:10:04 2012
@author: ajacobo
"""
import fipy as fp
import pylab as pl
#-----Mesh Variables and definition-----#
nx = 21
ny = 21
nz = 21
dx = 1e-7 # [m]
dy = 1e-7 # [m]
dz = 1e-7 # [m]
Lx = dx * nx
Ly = dy * ny
Lz = dz * nz
Mx = Lx/2
My = Ly/2
Mz = Lz/2
nchannels=1
dt = 1e-6
TotTime=1e-5
mesh = fp.Grid3D(dx=dx, dy=dy,dz=dz, nx=nx, ny=ny,nz=nz)
#-----Equation parammeters, variables and definition-----#
D_c = 2e-10 # [m^2 s^-1]
D_b = 2e-11 # [m^2 s^-1]
B_0 = 2 # [mol/m^3]
C = fp.CellVariable(name = "Calcium",
mesh = mesh, value = 0.)
B = fp.CellVariable(name = "Buffer",
mesh = mesh, value = B_0)
eqC = fp.TransientTerm(var=C) == fp.DiffusionTerm(coeff=D_c,var=C)
eqB = fp.TransientTerm(var=B) == fp.DiffusionTerm(coeff=D_b,var=B)
eq = eqC & eqB
#-----Define Boundary Conditions-----#
valueBC_C = 0
valueBC_B = B_0
X, Y, Z = mesh.faceCenters
x=X[0:nx]
facesTop = (mesh.facesTop)
facesBottom = (mesh.facesBottom)
facesFront = (mesh.facesFront)
facesBack = (mesh.facesBack)
facesRight = (mesh.facesRight)
facesLeft = (mesh.facesLeft)
#Apply the boundary conditions
C.constrain(valueBC_C, facesTop)
C.constrain(valueBC_C, facesBottom)
C.constrain(valueBC_C, facesRight)
C.constrain(valueBC_C, facesLeft)
C.constrain(valueBC_C, facesFront)
C.constrain(valueBC_C, facesBack)
B.constrain(valueBC_B, facesTop)
B.constrain(valueBC_B, facesBottom)
B.constrain(valueBC_B, facesRight)
B.constrain(valueBC_B, facesLeft)
B.constrain(valueBC_B, facesFront)
B.constrain(valueBC_B, facesBack)
#-----Solve the equation-----#
steps = int(TotTime/dt)
print steps
time=[]
t=0.
fig1=pl.figure(num=None, figsize=(5.2,4), dpi=100, facecolor='w')
fig1.subplots_adjust(bottom=0.15,left=0.15)
ax = fig1.add_subplot(111)
step=0
for step in range(1,steps):
eq.solve(dt=dt)
ax.plot(x,B([x,[My],[Mz]]))
time.append(t*dt)
t=t+1
pl.show()
#Plot the results
dmin=0.9*min(C)
dmax=1.1*max(C)
viewer = fp.Viewer(vars=C, datamin=dmin, datamax=dmax)
viewer.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
