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)

Reply via email to