Marc, Attached is a fixed script. I did two things. One was to make
sure that nx and ny are set as integers. I think they were floats,
which was causing some subtle issue or other. I thought we auto cast
them to ints internally, but who knows. That fixed the IndexError. I
also added a pylab.figure() statement as you already had a plot. The
second plot didn't seem to show up and adding the statement opened up
a new window. Seems to at least do something now. No idea if it's what
you want though. Good luck.
On Thu, May 19, 2011 at 10:34 AM, Marc Saudreau
<[email protected]> wrote:
> I've got float64 in both cases !!
> May my *.py file help you ?
>
> Thanks
>
> Marc
>
> Daniel Wheeler <[email protected]> a écrit :
>
>>
>> What do you get if you print numerix.array(S).dtype?
>>
>> On Thu, May 19, 2011 at 9:13 AM, Marc Saudreau
>> <[email protected]> wrote:
>>>
>>> Hi Daniel,
>>>
>>> thanks a lot. Your example works fine when I use it alone. But when I
>>> want
>>> to use it in
>>> my code, an error occurs with the instruction v([S, S], order=1) :
>>> __call__ C:\FiPy-2.1.2\fipy\variables\cellVariable.py 205
>>> IndexError: arrays used as indices must be of integer (or boolean) type
>>>
>>> I tried to get the solution by myself but to be honest this one is a
>>> little
>>> bit disturbing.
>>> Any solution ?
>>>
>>> Thanks again for your help and time
>>>
>>> Marc
>>>
>>>
>>> Daniel Wheeler <[email protected]> a écrit :
>>>
>>>>
>>>> Hi Marc,
>>>>
>>>> CellVariable's __call__ method extracts data. So you can do the
>>>> following:
>>>>
>>>> from fipy import *
>>>> import pylab
>>>>
>>>> L = 1.
>>>> n = 10
>>>> m = Grid2D(nx=n, ny=n, dx=L / n, dy=L / n)
>>>> x, y = m.cellCenters
>>>> v = CellVariable(mesh=m, value = x * y)
>>>> N = 50
>>>> S = numerix.arange(N) / float(N - 1) * L
>>>> pylab.plot(S, v([S, S], order=1))
>>>> pylab.show()
>>>> raw_input('stopped')
>>>>
>>>> Basically, you can extract any data from "v" with "v([x, y])" where
>>>> "x" and "y" define a set of coordinates.
>>>>
>>>> Cheers
>>>>
>>>> On Tue, May 17, 2011 at 4:47 AM, Marc Saudreau
>>>> <[email protected]> wrote:
>>>>>
>>>>> Hi everyone,
>>>>>
>>>>> I'm trying to use Fipy to solve a 2D heat transfer problem, and I need
>>>>> to
>>>>> extract from my 2D grid, some 1D profiles. For instance I need to plot
>>>>> temperature versus y axis for x given.
>>>>> I had a look at Fipy examples but I did not find any way to do it.
>>>>> Thanks a lot for any help.
>>>>>
>>>>> Best regards,
>>>>>
>>>>> Marc
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Daniel Wheeler
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>>
>>
>>
>>
>> --
>> Daniel Wheeler
>>
>>
>>
>
>
--
Daniel Wheeler
#-------------------------------------------------------------------------------
# Name: module1
# Purpose:
#
# Author: msaudreau
#
# Created: 15/05/2011
# Copyright: (c) msaudreau 2011
# Licence: <your licence>
#-------------------------------------------------------------------------------
#!/usr/bin/env python
def main():
import fipy
from fipy import *
import pylab
import numpy
#Define plotting function
dataMin = 24.
dataMax = 26.
#Define geometrical attributes
#Grid dimension (m)
Lx = 1
Ly = 0.5
#Grid spatial step
dx = 0.02
dy = 0.01
#Grid points
nx = int(Lx/dx)
ny = int(Ly/dy)
mesh = fipy.Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
#Define physical attributes
#Air
Tair = 20. #Temperature (Celsius)
Cpair = 1.2 #heat capacity
rhoair = 1.2 #density
#Material
Tmat = 25. #Temperature (Celsius)
Cpleaf = 2700 #heat capacity (J/kgK)
rholeaf = 1000 #density (kg/m3) The product rhoCp equals 2.7MJ/Km3 for a leaf (see Bajon et al., 2005, Infrared Physics & Techn.)
kleaf = 0.0762 #heat conductivity (W/mK) (see Casada et al., 1989, Ame. Soc. Agri. Eng.)
aleaf = kleaf/(Cpleaf*rholeaf) #heat diffusivity
Cpmat = 1400 #heat capacity (J/kgK)
rhomat = 1000 #density (kg/m3)
kmat = 20. #heat conductivity (W/mK)
amat = kmat/(Cpmat*rhomat) #heat diffusivity
print "Material heat diffusivity:",amat,"m2/s"
print "Material diffusion time L*L/diffusivity coeff. y axis", Ly*Ly/amat,' seconds'
print "Material diffusion time L*L/diffusivity coeff. x axis", Lx*Lx/amat,' seconds'
## raw_input("Press <return> to proceed...")
#Test the diffusion process
#Heat fluxes are set to zero (adiabatic conditions). Initial temperature field To is
#constant in x direction and linear in the y direction.
#The steady-state temperature value should be the average of To
#Define the initial temperature field
x, y = mesh.getCellCenters()
T1 = CellVariable(name ="Temperature", mesh=mesh,hasOld=True, value = Tmat + 2*(y-Ly/2)/Ly)
print T1.getGlobalValue()
viewer1 = MatplotlibViewer(vars=T1, name ="Temperature", datamin=dataMin, datamax=dataMax)
viewer1.plot()
raw_input("Press <return> to proceed...")
#Plot 1D profile
#Variable to plot
v = CellVariable(mesh=mesh, value = T1 )
N = 100
S = numpy.arange(N) / float(N - 1)*Ly
print numpy.array(S).dtype
print v
print S
pylab.figure()
pylab.plot(S, v([S, S], order=1))
pylab.show()
#Define heat fluxes / Boundary conditions
#y = 0 | Adiabatic wall
BottomFlux = 0.
#y = Ly | Adiabatic wall
TopFluxValue = 0 #W/m2
#x = 0 | Adiabatic wall
LeftFlux = 0. #W/m2
#x = Lx | Adaibatic wall
RightFlux = 0. #W/m2
LeftBC = FixedFlux(faces=mesh.getFacesLeft(), value=LeftFlux)
RightBC = FixedFlux(faces=mesh.getFacesRight(), value=RightFlux)
BottomBC = FixedFlux(faces=mesh.getFacesBottom(), value=BottomFlux)
TopBC = FixedFlux(faces=mesh.getFacesTop(), value=TopFluxValue)
BCs = (LeftBC, RightBC, BottomBC,TopBC)
#Find steady state solution by solving the transient equation
eqn = ImplicitDiffusionTerm(coeff=kmat) == TransientTerm(coeff = rhomat*Cpmat)
#Define time step
timestep = 100*dy*dy/amat
timemax = Ly*Ly/amat
steps = int32(floor(timemax/timestep))+1
print "time step", timestep
print "time max", timemax
print "step number", steps
raw_input("Time step defined. Press <return> to proceed...")
for step in range(steps):
T1.updateOld()
res = 1
while res > 1e-6:
res = eqn.sweep(T1, boundaryConditions=BCs, dt=timestep)
print "step, residual :", step, res
viewer1.plot()
print "min(T), max(T):", min(T1), max(T1)
print "T =",T1
print "Tmat", Tmat
raw_input("Transient equation solved. Press <return> to proceed...")
if __name__ == '__main__':
main()
