Hi All,
I am trying to solve a system of pdes, I get an error about the use of the
operator *.
Can you please look at my code and see the error I get??
Also, I have in my equations a singularity problem at r=0, is it fine for
fipy to treat this type of problems?
Regards,
Fadoua
from fipy import *
from numpy import memmap
#the constants
nr = 101
lamda = 0.2
alpha = 0.1
xi = 0.9
beta = 0.1
dx = 1./nr
#Diffusion terms
Ds = 0.5
Db = 0.25
steps = 20000
timeStepDuration = 0.9 * dx**2 / (2 * Ds)
#define the grid
mesh = Grid1D(nx = nr, dx = dx)
r = mesh.getCellCenters()
#The variables & the initial conditions
S = CellVariable(name="Sensing chemicals", mesh=mesh, value=0.)
B0 = 0.9
B = CellVariable(name="Bacterial density in the free space", mesh=mesh, value=B0)
valueLeft = 0.
valueRight = 0.
#The boundary conditions
Bcb = (FixedFlux(faces=mesh.getFacesRight(),value=valueRight),
FixedValue(faces=mesh.getFacesLeft(), value=valueLeft))
Bcs = (FixedFlux(faces=mesh.getFacesRight(),value= 0. ),
FixedFlux(faces=mesh.getFacesLeft(), value=5.))
#We define our system pf Differenrial equations
eqB = TransientTerm() == (1/r.getFaceValue()**2) * DiffusionTerm(coeff=Db * r.getFaceValue()**2) - xi * (B.getFaceGrad() * S.getFaceGrad() + B.getFaceValue() * r.getFaceValue() * (r.getFaceValue()**2 * S.getFaceGrad()).getDivergence())
eqS = TransientTerm() == (1/r.getFaceValue()**2) * DiffusionTerm(coeff=Ds * r.getFaceValue()**2) - lamda * S.getFaceValue() + alpha * B.getFaceValue()
viewerB = Viewer(vars=(B,))
viewerS = Viewer(vars=(S,))
#we solve the pde
for step in range(steps):
eqB.solve(var = B, boundaryConditions = Bcb, dt = timeStepDuration)
eqS.solve(var = S, boundaryConditions = Bcs, dt = timeStepDuration)
viewerB.plot()
viewerS.plot()
import time
time.sleep(2)