Dear FiPy developers and users!
Thank you for the answers to all my previous questions!
I have NEW troubles coming with NEW version. I will be very grateful for your
help with solution of these troubles.
I'm using fipy for the solution of convection problem. I was glad with the
result of old fipy, however I wanted to make the solution more efficient. This
could be done by introducing logarithmic mesh to my problem that would result
in less steps in space and faster response from fipy. I was adviced to
introduce the logarithmic mesh with:
>>> mesh=Grid1D(dx=10**numerix.arange(-3., 3., .1))
that old fipy failed to create with an error:
Traceback (most recent call last):
File "espec.py", line 112, in <module>
espec(500)
File "espec.py", line 56, in espec
mesh = Grid1D(dx = dx)
File "/usr/lib/python2.5/site-packages/fipy/meshes/grid1D.py", line 53, in
Grid1D
return grid1D.Grid1D(dx = dx, nx = nx)
File "/usr/lib/python2.5/site-packages/fipy/meshes/numMesh/grid1D.py", line
87, in __init__
Mesh1D.__init__(self, vertices, faces, cells)
File "/usr/lib/python2.5/site-packages/fipy/meshes/numMesh/mesh.py", line 71,
in __init__
_CommonMesh.__init__(self)
File "/usr/lib/python2.5/site-packages/fipy/meshes/common/mesh.py", line 64,
in __init__
self._calcTopology()
File "/usr/lib/python2.5/site-packages/fipy/meshes/numMesh/mesh.py", line
303, in _calcTopology
self._calcFaceCellIDs()
File "/usr/lib/python2.5/site-packages/fipy/meshes/numMesh/mesh.py", line
327, in _calcFaceCellIDs
numerix.put(firstRow, cellFaceIDsFlat[::-1], array[::-1])
TypeError: 'FlatIter' object is unsubscriptable
I decided to switch to the new version (FiPy-2.0a1 install via svn from the
trunk) and first tested if it solves my equation similar to the old version. I
was REALLY surprised to find completely different solution that is no way
correct. I met as well other problems, but let it solve step by step. I attach
here the code of my script with comments showing the changes I made for the
FiPy-2.0a1.
Is it possible to have the old and the new version on the same PC? How to tell
the script which version to use?
After installation FiPy-2.0a1 didn't want to run because README.txt was not in
/usr/lib/python2.5/site-packages/FiPy-2.0a1-py2.5.egg/. I copied it there
manually.
Why does FiPy-2.0a1 solve the equation differently from FiPy-1.2.1?
Kind regards,
Igor.
The code:
(NB:Pygist shows negative values as positive in log scale!)
-----------------------------------------------
from fipy import *
from numpy import *
from scipy import log10
from scipy import exp
from scipy import sqrt
nsteps = 1000
Emax = 1000.
Ee0=5.11e-7
n=122.5
B=0.0001
Ucmb=0.25
A=5.075e-5
Ke=1
alpha=-2
Ecut = 10
yr=3.156E+7
def L(E):
Loss = (1e-12)*(0.65e-3)*B*B*(1/(Ee0*Ee0))*E*E+ \
(1e-12)*(2.66e-14)*Ucmb*(1/(Ee0*Ee0))*E*E+ \
(1.370e-16)*n*E*log10(E/Ee0+0.36)+ \
(8.978e-13)*n*Ee0
return Loss*yr*1/1.
def Jin(E):
return Ke*A*E**alpha*exp(-E/Ecut)
def espec(t):
dx = Emax / nsteps
mesh = Grid1D(dx = dx, nx=nsteps)
x, = mesh.getCellCenters()
#x = mesh.getCellCenters()[...,0] #worked in old version
convCoeff = FaceVariable(mesh=mesh, value=L(mesh.getFaceCenters()), rank=1)
#convCoeff = VectorFaceVariable(mesh=mesh, value=L(mesh.getFaceCenters()))
#worked in old version
#initial conditions
phi = CellVariable(mesh=mesh, name = "variable", value=Jin(x))
phiInitial = CellVariable(mesh=mesh, name="variable initial", value=Jin(x))
viewer = make(vars=(phi,phiInitial),limits={'datamin': -1e-35,
'datamax':1e0})
viewer.xlog=1
viewer.ylog=1
valueLeft = 0
valueRight = 0
boundaryConditions = (FixedFlux(faces=mesh.getFacesLeft(),
value=valueLeft), FixedFlux(faces=mesh.getFacesRight(), value=valueRight))
eq = TransientTerm() == PowerLawConvectionTerm(convCoeff)
timeStepDuration = 0.5
i=0
time=0
while time <= t:
timeStep = timeStepDuration
if time>= 50:
timeStep = timeStepDuration*10
if time>=500:
timeStep = timeStepDuration*100
if time>=5000:
timeStep = timeStepDuration*1000
eq.solve(var=phi, boundaryConditions=boundaryConditions, dt=timeStep,
solver=LinearLUSolver(tolerance=1.0e-15))
i=i+1
time=time+timeStep
viewer.plot()
TSVViewer(vars=(phi)).plot(filename="espectrum.tsv")
espec(10000)
raw_input("Press enter to close")
