| Dear FiPy developers, Thanks a lot for sharing the very useful FiPy software and providing reactive support. I’ve got a set of equations for which FiPy should be the right tool to use I think. I have 3 simple convection-reaction-diffusion PDEs (convection-diffusion with reactive source term) in 1-dimenion, which are coupled together and with an additional ODE for the (convection) velocity (which comes from force balance for an active viscous layer). To make it more clear I attach a PDF with the governing equations. I manage to treat convection-reaction-diffusion correctly, but I block on the coupling to the ODE for the convective velocity with FiPy, which includes 1st-order spatial derivatives of the other coupled variables. I tried to treat these spatial derivatives as source terms, but .grad leads to rank-1 terms, while the ODE is written for the velocity defined as a rank-0 CellVariable. In the actual implementation, I try, on the contrary to treat these 1st-order spatial derivatives as a convection term, but I get the error below. I attach my current Python script to make it more clear. Traceback (most recent call last): File "/Users/Turlier/Dropbox/Research/Projects/AdvectionReactionDiffusion/FiPy/coupled-advection-reaction-diffusion-mechanics_1D.py", line 83, in <module> res = Eq.sweep(dt=dt, solver=solver) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/term.py", line 223, in sweep solver = self._prepareLinearSystem(var=var, solver=solver, boundaryConditions=boundaryConditions, dt=dt) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/term.py", line 157, in _prepareLinearSystem buildExplicitIfOther=self._buildExplcitIfOther) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/coupledBinaryTerm.py", line 122, in _buildAndAddMatrices buildExplicitIfOther=buildExplicitIfOther) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/binaryTerm.py", line 68, in _buildAndAddMatrices buildExplicitIfOther=buildExplicitIfOther) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/binaryTerm.py", line 68, in _buildAndAddMatrices buildExplicitIfOther=buildExplicitIfOther) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/binaryTerm.py", line 68, in _buildAndAddMatrices buildExplicitIfOther=buildExplicitIfOther) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/unaryTerm.py", line 99, in _buildAndAddMatrices diffusionGeomCoeff=diffusionGeomCoeff) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/abstractConvectionTerm.py", line 191, in _buildMatrix var, L, b = FaceTerm._buildMatrix(self, var, SparseMatrix, boundaryConditions=boundaryConditions, dt=dt, transientGeomCoeff=transientGeomCoeff, diffusionGeomCoeff=diffusionGeomCoeff) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/faceTerm.py", line 178, in _buildMatrix self._implicitBuildMatrix_(SparseMatrix, L, id1, id2, b, weight['implicit'], var, boundaryConditions, interiorFaces, dt) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/terms/faceTerm.py", line 74, in _implicitBuildMatrix_ L.addAt(numerix.take(coeffMatrix['cell 1 diag'], interiorFaces, axis=-1).ravel(), id1.ravel(), id1.swapaxes(0,1).ravel()) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/matrices/offsetSparseMatrix.py", line 58, in addAt SparseMatrix.addAt(self, vector, id1 + self.mesh.numberOfCells * self.equationIndex, id2 + self.mesh.numberOfCells * self.varIndex) File "/usr/local/lib/python2.7/site-packages/FiPy-unknown-py2.7.egg/fipy/matrices/pysparseMatrix.py", line 246, in addAt self.matrix.update_add_at(vector, id1, id2) IndexError: id1 does not have the same size as b I went through the manual and spent several hours on past discussions on the mailing-list, but I could not find, yet a working solution. Your help would be greatly appreciated. Best, Hervé ---
Hervé Turlier, PhD European Molecular Biology Laboratory, Meyerhofstrasse 1 69117 Heidelberg, Germany [email protected] [email protected] Tel: +49 (0)6221 387-8182 |
1D_AdvectionReactionDiffusionWithMechanics.pdf
Description: Adobe PDF document
#!/usr/bin/env python2.7 # Coupled advection-Reaction-diffusion with mechanics in 1 dimension - matrix form # Typical toy model for actomyosin cortical dynamics
import sys
from fipy import *
# MESH
nx = 100
Lx = 1.
dx = Lx / nx
mesh = Grid1D(nx=nx, dx=dx)
print mesh.cellCenters
x = mesh.cellCenters[0]
# PHYSICAL PARAMETERS
Da = 0.1 # diffusion coefficient for the active myosin pool
Di = 10 # diffusion coefficient for the inactive pool
ka = 1. # myosin activation rate
ki = 1. # myosin inactivation rate
kd = 10. # actin depolymerization rate
eta = 10. # viscosity
zmu = 1. # active stress
gam = 0. # friction term
# INITIATION VALUES
Ca0 = CellVariable(name="initial value", mesh=mesh) # Active myosin pool initial value
sigma = 0.1
Ca0.setValue( 1./(sigma*numerix.sqrt(2.*numerix.pi))*numerix.exp(-0.5*(x/sigma)**2) )
Ci0 = 1. # Inactive myosin pool initial value
A0 = 1. # Actin thickness initial value
V0 = 0. # Initial velocity initial value
# VARIABLES
Ca = CellVariable(name="Active pool", mesh=mesh, hasOld=True, value=Ca0 )
Ci = CellVariable(name="Inactive pool", mesh=mesh, hasOld=True, value=Ci0 )
A = CellVariable(name="Actin thickness", mesh=mesh, hasOld=True, value=A0 )
V = CellVariable(name="Velocity", mesh=mesh, hasOld=True, value=V0 )
# SET UP THE VIEWER
if __name__ == '__main__':
viewer = Matplotlib1DViewer(vars=(Ca,Ci,A), datamin=0., datamax=4.)
viewer.plot()
raw_input("Press <return> to start...")
# BOUNDARY CONDITIONS
# Vanishing velocity at boundaries
V.constrain(0., mesh.facesRight)
V.constrain(0., mesh.facesLeft)
# Zero flux of active myosin (Robin->Neumann because V=0 at boundaries)
Ca.faceGrad.constrain([0.], mesh.facesRight)
Ca.faceGrad.constrain([0.], mesh.facesLeft)
# Zero flux of inactive myosin (Robin->Neumann because V=0 at boundaries)
Ci.faceGrad.constrain([0.], mesh.facesRight)
Ci.faceGrad.constrain([0.], mesh.facesLeft)
# Zero flux of actin
A.faceGrad.constrain([0.], mesh.facesRight)
A.faceGrad.constrain([0.], mesh.facesLeft)
# DEFINE EQUATIONS
R = ka * Ci * A/A0 - ki * Ca # explicit reaction term
EqCa = ( TransientTerm(var=Ca) == - ExponentialConvectionTerm(var=Ca, coeff=(V,)) + DiffusionTerm(var=Ca, coeff=Da) + R )
EqCi = ( TransientTerm(var=Ci) == - ExponentialConvectionTerm(var=Ci, coeff=(V,)) + DiffusionTerm(var=Ci, coeff=Di) - R )
EqA = ( TransientTerm(var=A) == - ConvectionTerm(var=A, coeff=(V,)) - kd*(A-A0) )
EqV = ( ImplicitSourceTerm(var=V, coeff=gam) == DiffusionTerm(var=V, coeff=4*eta*A ) + ConvectionTerm(var=A, coeff=(zmu*Ca,)) )
Eq = EqCa & EqCi & EqA & EqV
# SOLVE EQUATION(S)
elapsed = 0.
dt = 1.e-2
solver = LinearLUSolver(tolerance=1e-10)
while elapsed < 2.:
elapsed += dt
Ca.updateOld()
Ci.updateOld()
A.updateOld()
V.updateOld()
res = 1.
while res > 1.e-3:
res = Eq.sweep(dt=dt, solver=solver)
if __name__ == '__main__':
viewer.plot()
raw_input("Press <return> for next step...")_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
