Hello Serbulent, thank you for sharing!
Right now I go on with the mesh20x20 example and got a question here:
########################
from fipy import *
nx = 100
ny = nx
dx = 1.
dy = dx
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = CellVariable(name = "solution variable",mesh = mesh,value = 0.5)
D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D)
#print eq
valueTopLeft = 0
valueBottomRight = 1
X, Y = mesh.faceCenters
x,y = mesh.cellCenters
#print X, len(X)
facesTopLeft = ((mesh.facesLeft & (Y > L / 2))
| (mesh.facesTop & (X < L / 2)))
facesBottomRight = ((mesh.facesRight & (Y < L / 2))
| (mesh.facesBottom & (X > L / 2)))
bacteria = (x>9)&(x<22)&(y>9)&(y<22)
bacteria2 = ((x-50)**2+(y-70)**2) < 10
bacteria3 = ((x-60)**2+(y-70)**2) < 10
bacteria4 = ((x-70)**2+(y-70)**2) < 10
print "_____",type(bacteria), type(mesh.facesLeft), type(facesTopLeft)
phi.constrain(valueTopLeft, facesTopLeft)
phi.constrain(valueBottomRight, facesBottomRight)
phi.constrain(0, bacteria)
phi.constrain(10, bacteria2)
phi.constrain(2, bacteria3)
phi.constrain(0.2, bacteria4)
#FixedFlux(mesh.facesLeft,-2)
if __name__ == '__main__':
viewer = Viewer(vars=phi, datamin=0., datamax=2.)
viewer.plot()
timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
steps = 50
for step in range(steps):
if step > 25:
phi.constrain(2.0/(1+step), bacteria2)
eq.solve(var=phi,dt=timeStepDuration)
faceValue = phi.arithmeticFaceValue[bacteria]
print faceValue, len(faceValue)
if __name__ == '__main__':
viewer.plot()
#######################
My questions are:
* Why is my selection "bacteria" of type <type 'numpy.ndarray'>, whereas
the
facesTopLeft is <class 'fipy.variables.binaryOperatorVariable.binOp'> ?
* How can I change the value of phi in the bacteria2-area to the recent
value in another equally sized patch? And how can I manipulate phi at one
single place? It should not be constant afterwards.
Thank you all!
Damian
2013/9/24 Serbulent UNSAL <[email protected]>
> Hi,
>
> If I understand your problem right, it is definitely possible. I apply
> similar scenario for tumor cells and consumption of oxygen and glucose and
> production of H+ ions in tumor micro-environment.
>
> You only need to use a source term (which is a 1D matrix that represents
> your 2D environment). You can set values in a 2D matrix according to
> coordinates and then flatten the matrix. I paste a part of my code below
> that should give you the basic idea.
>
> Serbulent
>
> # This function called in each time step
>
> def calculateO2Concentration(self, o2consumptionMatrix, timeStep):
> chain = itertools.chain(*o2consumptionMatrix) flattenedMatrix =
> list(chain) o2concentrationMatrix =
> self.o2Calculator.calculateO2(flattenedMatrix, timeStep) return
> o2concentrationMatrix;
>
> ---------------------------------------------------------------------------------------------
>
> # and this is the main code which calculates diffusion of nutrient based on
> consumptionfrom fipy import *class o2CalculatorClass: def
> __init__(self,nonDimensionaliser,configObject): self.conf =
> configObject nx = configObject.nx cellSize =
> configObject.cellSize nx = nx/cellSize dx = configObject.dx
> ny = nx dy = dx mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
> boundaryCondition = 0.2 self.phi = CellVariable(name="Oxygen
> Concentration", mesh=mesh, value = boundaryCondition )
> self.nonDimensionalizer = nonDimensionaliser D =
> self.nonDimensionalizer.nonDimensionalO2DiffConst
> self.phi.constrain(boundaryCondition, mesh.facesLeft)
> self.phi.constrain(boundaryCondition, mesh.facesRight)
> self.phi.constrain(boundaryCondition, mesh.facesTop)
> self.phi.constrain(boundaryCondition, mesh.facesBottom)
> proliferationAge = configObject.paStandart self.timeStepDuration =
> 1/proliferationAge self.viewer = Viewer(vars=self.phi)
> self.source = CellVariable(mesh=mesh) #self.sourceCapilary =
> CellVariable(mesh=mesh) self.eq = TransientTerm() ==
> DiffusionTerm(coeff=D) - self.source #+ self.sourceCapilary def
> calculateO2(self, sourceTerm, timeStep):#, sourceCapilary):
> self.source.setValue(sourceTerm)
> #self.sourceCapilary.setValue(sourceCapilary)
> self.eq.solve(var=self.phi, dt=self.timeStepDuration) willImagesSaved
> = self.conf.save_images if willImagesSaved and
> timeStep%self.conf.imageSaveInterval == 0: nameofpic =
> "./images/o2Diffusion_" + str(timeStep) + ".png"
> self.viewer.plot(filename=nameofpic) else:
> self.viewer.plot() #TSVViewer(vars=(phi)).plot(filename="output.txt")
> return self.phi;
>
>
>
>
>
>
> 2013/9/24 Damian Kösters <[email protected]>
>
>> We want to model a petridish of bacteria in a liquid medium. So the cells
>> comsume some substances and produce others. After some time cells in
>> regions with good food supply are going to split. So new cells are set on
>> the dish and we carry on with the calculation using the final values of the
>> last step as initial values of the next.
>>
>> In Comsol we did not find any documentation treating this scenario. Using
>> old values to proceed with a calculation seems not to work when you change
>> the geometry in the mean time. Does this work with FiPy? To me it seems
>> possible (after playing with the mesh20x20 example) to change single points
>> in a cell variable. That would be good news for us.
>>
>> (Furthermore the communication with the Comsol server from Matlab is
>> quite hard to understand and debug. Fipy is much better in this respact.
>> Nonetheless, the Comsol Desktop is better than Comsol's Matlab interface.)
>>
>> Thank you!
>>
>>
>>
>> 2013/9/23 Guyer, Jonathan E. Dr. <[email protected]>
>>
>> On Sep 23, 2013, at 3:47 AM, Damian Kösters <[email protected]>
>>> wrote:
>>>
>>> > I have got a problem I can't solve using COMSOL+Livelink to Matlab.
>>> http://www.comsol.com/community/forums/general/thread/39531/
>>> >
>>> > I just need diffusion calculation and reactions with the substances.
>>> >
>>> > Would you consider this solvable with Fipy?
>>>
>>> Can you describe mathematically what you're trying to do and what's not
>>> working with COMSOL? FiPy can certainly solve reaction-diffusion problems.
>>> I would think COMSOL could, too.
>>>
>>> I don't have any practical experience with COMSOL to follow your problem
>>> definition.
>>>
>>>
>>>
>>> _______________________________________________
>>> fipy mailing list
>>> [email protected]
>>> http://www.ctcms.nist.gov/fipy
>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>>
>>
>>
>> _______________________________________________
>> fipy mailing list
>> [email protected]
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
