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
#-------------------------------------------------------------------------------
# Name: module1
# Purpose:
#
# Author: msaudreau
#
# Created: 15/05/2011
# Copyright: (c) msaudreau 2011
# Licence: <your licence>
#-------------------------------------------------------------------------------
#!/usr/bin/env python
def main():
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 = Lx/dx
ny = 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
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()
