gmane.comp.python.fipy

Using a dynamic link library (dll) in python

Etienne Rivard — 2008-09-22T07:12:58

Hello again, I'd like to use a dynamic link library (dll) in a problem I'll solve with FiPy. And not just any DLL. Our chair of research has a license for the excellent REFPROP software by you-know-who. I'm aware of the existence of the ctypes package, and I even read the tutorial. But quite frankly, I didn't understand much. I think there are much prerequisites I'm lacking. The other solution would be a python extension written in C. But I'm not proficient in C and there is no example of using REFPROP in a C program in the documentation. Perhaps this is a very simple problem for one of you? Or would it be better to ask your colleague Mr. Lemmon? Many thanks,

set a viewer

franck kalala — 2008-09-20T22:10:04

Dear All, How do I set a viewer in fipy? Regards 2008/9/12 Daniel Wheeler public.gmane.org>

Face Variable in Equation

Zhiwen Liang — 2008-09-18T21:54:43

Hi all, I wonder if a face variable can be used when defining equations. If I run the following lines (they might not be meaningful, just for example): from fipy import * mesh=Grid1D(nx=2) facev=FaceVariable(mesh=mesh,value=0.) eq=0==DiffusionTerm(0.)+facev I will get "TypeError: The coefficient must be a rank-0 CellVariable or a scalar value.". It seems that I could not use FaceTerm directly. Are there any convenient ways to do this? Thank you in advance for any suggestions. Have a nice day! Zhiwen

convection term with the unknown to a certain power

franck kalala — 2008-09-12T11:34:03

Dear all, I am solving an equation of this type phi_t=div(phi**3*div\theta) in the right hand side of this equation phi is to the power 3. in the manual it is explain how to deal with this term when phi is to the power 1. but here I don't know how to represent it in fipy the second hand side, I try this PowerLowerConvectionTerm((phi**3).getArithmeticFaceValue()*theta.getFaceGrad())) but with tha I don't have expected result for my codes. Is there any other better way to represent such term? Regards 2008/9/3 Daniel Wheeler public.gmane.org>

FiPy: The Virtual Kinetics of Materials Laboratory, version 0.1

R. Edwin Garcia — 2008-09-03T15:30:42

Dear FiPy enthusiasts, It is with great pleasure that we introduce to you the Virtual Kinetics of Materials Laboratory, VKML. VKML is a set of simple FiPy examples that simulate basic aspects of kinetics of materials through an interactive Graphical User Interface. The seminal development includes four examples: Polycrystalline Growth and Coarsening: simulates the growth, impingement, and coarsening of a random distribution of crystallographically oriented nuclei. The user can control every aspect of the model such as the nuclei radius, the size of the simulation cell, and whether the grains are homogeneously dispersed or only on one wall of the simulation. https://www.nanohub.org/tools/vkmlpsgg/ Dendritic Growth: simulates the anisotropic solidification of a single solid seed with an N-fold axis of crystallographic symmetry embedded in an undercooled liquid. The user can specify many material aspects of the solidification process, such as the thermal diffusivity and the strength of the surface tension anisotropy. Default values are physical but arbitrary. This model is based on the phase field method and an example shown in the FiPy manual. https://www.nanohub.org/tools/vkmlggs/ Two-Dimensional Spinodal Decomposition: simulates the time-dependent segregation of two chemical components and its subsequent coarsening, as presented by John Cahn. The default values are physical but arbitrary. https://www.nanohub.org/tools/vkmlsd/ Three-Dimensional Spinodal Decomposition has the same functionality as the 2D version, but has an interactive Three-Dimensional viewer. https://www.nanohub.org/tools/vkmlsd3d/ These modules provide a Graphical User Interface to FiPy, and allow you to perform simulations directly through your web browser. This approach to computing removes the need to install the software on your local machine (unless you really want to), allows you to assess current and potential FiPy applications and instead you only need a web browser to access it and run it. In other words, you can run these simulations (and simulations like this one) from a Windows machine, a Mac, or a Linux box, and you can also run the modules from Michigan, Boston, Japan, or England: from wherever you are. Moreover, if you close your web browser and leave your calculation running, when you come back a few hours later, your calculation will persist. Additionally, if there is something you want to share with a coworker, wherever he or she might be (e.g., the other side of the planet), you can grant him temporary access to your calculation so that the third party can directly see the output (or specify inputs directly into it, without having to travel to where you are). It is a great way to privately (or publicly) collaborate with other people, especially if the users are in different parts of the world. The only requirement to run VKML is to register (registration is 100% free) in the nanoHUB. Please go to: www.nanoHUB.org You can find a link to the working project from: http://www.ctcms.nist.gov/fipy/ Sincerely, R. Edwin Garcia Michael Waters R. Edwin García Assistant Professor of Materials Engineering School of Materials Engineering Purdue University Neil Armstrong Hall of Engineering 701 West Stadium Avenue West Lafayette, IN 47907-2045 phone: (765) 494-0148 email: redwing-olO2ZdjDehc3uPMLIKxrzw< at >public.gmane.org web: http://bicephalous.ecn.purdue.edu/~edwin

Slow Array Assignment

Zhiwen Liang — 2008-09-01T19:15:04

Dear Fipy developers and users, Thanks again for solving my last problem. This is the most powerful forum I have ever been involved. I do not know if this is an known issue, but I found it very slow to assign fipy variables as elements of the numerix array. Please see below. from fipy import * import fipy.tools.numerix as numerix mesh=Grid3D(nx=10,ny=10,nz=10) Dx=CellVariable(mesh=mesh,value=1.) Dy=CellVariable(mesh=mesh,value=2.) Dz=CellVariable(mesh=mesh,value=3.) zeros=CellVariable(mesh=mesh,value=0.) Ds=numerix.array(((Dx,zeros,zeros),(zeros,Dy,zeros),(zeros,zeros,Dz))) The last line of the above code took 3.5 second while the rest of the code only took 0.002 second, using my laptop. I wonder if it can be improved. Or maybe the way I did it is too clumsy? I define the "Ds" in the code in order to use it as the anisotropic coefficient in the DiffusionTerm. Thank you for any advice. Enjoy the rest of the Labor Day! Regards, Zhiwen

example.diffusion.cricle

franck kalala — 2008-08-26T11:39:41

Hi all I am so sorry to disturb you so much, I am very young in python, in fipy, sure I lent many thing from you since I decided to use fipy, I was this example ************************************************************************************************************** cellSize=0.05 radius=1. lines= ['cellSize = '+str(cellSize) + ';\in', 'radius ='+str(radius) + ';\in', 'point(1)={0, 0, 0, cellSize};\n', 'point(2)={-radius, 0, 0, cellSize};\n', 'point(3)={0, radius, 0, cellSize};\n', 'point(4)={radius, 0, 0, cellSize};\n', 'point(5)={0, -radius, 0, cellSize};\n', 'point(6)={2, 1, 3};\n', 'point(7)={3, 1, 4};\n', 'point(8)={4, 1, 5};\n', 'circle(9)={5, 1, 2};\n', 'Line Loop(10)={6, 7, 8, 9};\n', 'Plane Surface(11)={10};\n', import tempfile (f,geomName)=tempfile.mkstemp('.geo') file=open(geoName, 'w') file.writelines(lines) file.close() import os os.close(f) import sys if sys.platform=='win32': meshName='tmp.msh' else: (f,mesfName)=tempfile.mkstemp('.msh') os.system('gmsh'+geomName+'- -v 0 -format msh -o'+meshName) if sys.platform!='win32': os.close(f) os.remove(geomName) from fipy.meshes.gmshImport import GmshImporter2D mesh=GmshImporter2D(meshName) os.remove(meshName) from fipy.variables.cellVariable import CellVariable phi=CellVariable(name='solution variable', mesh=mesh, value=0) viewer=none if __name__=='__main__' try: from fipy.viewers.gistViewers.gist2DViewer import Gist2DViewer viewer=Gist2DViewer(vars=phi,limits={'datamin':-1, 'datamax':1.}) viewer.plotMesh() raw_input("Irregular circular mesh. Press to proceed...") except: print "Unable to create a Gist2DViewer" D=1. from fipy.terms.transientTerm import TransientTerm from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm eq=TransientTerm()==ImplicitDiffusionTerm(coeff=D) from fipy.tools import numerix exterirorXcoords=numerix.take(mesh.getFaceCenters()[...,0],mesh.getExteriorFaces()) from fipy.boundaryConditions.fixedValue import FixedValue BCs=(FixedValue(faces=mesh.getExteriorFaces(), value=exteriorXcoords),) timeStepDuration=10*0.9*cellSize**2/(2*D) steps=10 for step in range(steps): eq.solve(var=phi,boundaryConditionsStepDuration,dt=timeStepDuration) if viewer is none: viewer.plot() **************************************************************************************** when running this programme I have the error message saying that syntaxe error where it is in red in the above code, then i can not run the code. Thanks again for your help Regards

Surface plot using viewers

franck kalala — 2008-08-26T00:25:52

Hi all How do I plot this surface z=exp(-(x**2+y**2)) using viewers within fipy, for a given 2D mesh. Regards

Solver stagnated error

franck kalala — 2008-08-24T14:24:07

Thanks again for your help, but I still have one error message that I cannot understand, this my code ************************************************************************************************** rom __future__ import division L=1. nx=400 dx=L/nx timeStepDuration=0.012 from fipy.meshes.grid1D import Grid1D mesh=Grid1D(dx=dx, nx=nx) from fipy.variables.cellVariable import CellVariable phi=CellVariable(name="saturation", mesh=mesh) p=CellVariable(name="pressure",mesh=mesh) x=mesh.getCellCenters()[...,0] phi.setValue(.9, where=x<.5) phi.setValue(.3, where=x>=.5) p.setValue(1.) pcs=0 qin=40*((x>=.1)&(x<=.3))+20*((x>=.7)&(x<=.9)) qout=60*((x>=.4)&(x<=.6)) sin=.8 pin=0 q1s=sin*qin+phi*qout q2s=(1-sin)*qin+(1-phi)*qout if __name__=='__main__': import fipy.viewers phiviewer=fipy.viewers.make(vars=(phi), limits={'ymin':.2, 'ymax':1.1}) pviewer=fipy.viewers.make(vars=(p), limits={'ymin':.2, 'ymax':1.1}) phiviewer.plot() pviewer.plot() raw_input("initial condition. press to proceed...") from fipy.terms.transientTerm import TransientTerm from fipy.terms.powerLawConvectionTerm import PowerLawConvectionTerm from fipy.terms.explicitSourceTerm import _ExplicitSourceTerm from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm eq1=ImplicitDiffusionTerm(coeff=1)==_ExplicitSourceTerm(q1s+q2s) eq2=TransientTerm(coeff=2)==_ExplicitSourceTerm(q1s-q2s)+PowerLawConvectionTerm(coeff=2*p.getFaceGrad())-p.getFaceGrad().getDivergence() steps=100 for i in range(steps): eq1.solve(p, dt=timeStepDuration) eq2.solve(phi, dt=timeStepDuration) if __name__ =='__main__': phiviewer.plot() pviewer.plot() ************************************************************************************************************************* when running this code i got the following error message Warning (from warnings module): File "C:/Python25/tsetdipha.py", line 43 eq1.solve(p, dt=timeStepDuration) MaximumIterationWarning: Iterations: 1001. Relative error: 2.96727e+018 Warning (from warnings module): File "C:/Python25/tsetdipha.py", line 43 eq1.solve(p, dt=timeStepDuration) StagnatedSolverWarning: The solver stagnated. Iterations: 20. Relative error: 2.23961e+016 Warning (from warnings module): File "C:/Python25/tsetdipha.py", line 43 eq1.solve(p, dt=timeStepDuration) StagnatedSolverWarning: The solver stagnated. Iterations: 4. Relative error: 4.16036e+016 Warning (from warnings module): File "C:/Python25/tsetdipha.py", line 43 eq1.solve(p, dt=timeStepDuration) StagnatedSolverWarning: The solver stagnated. Iterations: 2. Relative error: 6.11332e+016 This message seems to appear infinitly I don't understand this error, could any one to help me? Regards 2008/8/20 Jonathan Guyer public.gmane.org>

Porting to FiPy-2.0a1. Troubles.

Igor Telezhinsky — 2008-08-24T12:48:06

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: that old fipy failed to create with an error: Traceback (most recent call last): File "espec.py", line 112, in 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")

2D diffusion problem in fipy

franck kalala — 2008-08-22T15:40:08

Thanks. Fipy can solve this diffusdion problem in 2D? u_t=D(u_xx+u_yy) u(x,y,0)=ph(x,y) how to define this problem(the first line) and the structure of the mesh? Regards 2008/8/22 Daniel Wheeler public.gmane.org>

style of line

franck kalala — 2008-08-22T14:30:23

How to set different style of line when plotting several curve? for example if i plot two curves and i want to have different line style for this curve. Regards. Franck 2008/8/20 Jonathan Guyer public.gmane.org>

VectorFaceVariable?

Zhiwen Liang — 2008-08-21T20:45:55

Hi all, I am sorry that I have to trouble you guys again. Thank you in advance for any advice. There used to be a VectorFaceVariable in the old fipy. But after I updated to 2.0a1, it is gone. I wonder if there is a similar type of variable that I can use to replace it. Thanks, Zhiwen

bool source term

franck kalala — 2008-08-20T10:32:51

Thanks to show me how to use characteristic function as a source. But I still have a problem. I am solving these equations 1)eqp=ImplicitDiffusionTerm(coeff=1)==source1, where source1=q1s+q2s(pressure equation) 2)eqphi=TransientTerm(coeff=2)==source2 where source2=q1s-q2s(saturation eqaution) where: qin=40*((phi>=.1)and(phi<=.3))+20*((phi>=.7)and(phi<=.9)) qout=60*((phi>=.4)and(phi<=.6)) sin=.8 q1s=sin*qin+s*qout q2s=(1-sin)*qin+(1-phi)*qout when I run the I program I have the following error: Evaluating diphabis.py Traceback (most recent call last): File "", line 1, in AttributeError: 'bool' object has no attribute 'solve' I don't understand this error, could any one have an idea? see my code below. Regards, Franck ********************************************************************************************************************************************************************************************* nx=400 dx=L/nx timeStepDuration=0.012 from fipy.meshes.grid1D import Grid1D mesh=Grid1D(dx=dx, nx=nx) from fipy.variables.cellVariable import CellVariable phi=CellVariable(name="saturation", mesh=mesh) p=CellVariable(name="pressure",mesh=mesh) x=mesh.getCellCenters()[...,0] phi.setValue(.9, where=x<.5) phi.setValue(.3, where=x>=.5) p.setValue(1.) #pcs=0 #built sources qin=40*((phi>=.1)and(phi<=.3))+20*((phi>=.7)and(phi<=.9)) qout=60*((phi>=.4)and(phi<=.6)) sin=.8 #pin=0 q1s=sin*qin+phi*qout q2s=(1-sin)*qin+(1-phi)*qout from fipy.terms.transientTerm import TransientTerm from fipy.terms.implicitSourceTerm import ImplicitSourceTerm from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm eqp=ImplicitDiffusionTerm(coeff=1)==q1s+q2s eqphi=TransientTerm(coeff=2)==q1s-q2s steps=50 for i in range(steps): eqp.solve(p,dt=timeStepDuration) eqphi.solve(phi,dt=timeStepDuration) if __name__ =='__main__': phiviewer=fipy.viewers.make(vars=(phi), limits={'datamin':0., 'datamax':1.}) pviewer=fipy.viewers.make(vars=(p), limits={'datamin':0., 'datamax':1.}) phiviewer.plot() pviewer.plot() 2008/8/18 Jonathan Guyer public.gmane.org>

bool source term

franck kalala — 2008-08-20T10:19:03

Thanks to show me how to use characteristic function as a source. But I still have a problem. I am solving these equations 1)eqp=ImplicitDiffusionTerm(coeff=1)==source1, where source1=q1s+q2s 2)eqphi=TransientTerm(coeff=2)==source2 where source2=q1s-q2s where: qin=40*((phi>=.1)and(phi<=.3))+20*((phi>=.7)and(phi<=.9)) qout=60*((phi>=.4)and(phi<=.6)) sin=.8 q1s=sin*qin+s*qout q2s=(1-sin)*qin+(1-phi)*qout when I run the I program I have the following error: Evaluating diphabis.py Traceback (most recent call last): File "", line 1, in AttributeError: 'bool' object has no attribute 'solve' I don't understand this error, could any one have an idea? see my code below. Regards, Franck ********************************************************************************************************************************************************************************************* nx=400 dx=L/nx timeStepDuration=0.012 from fipy.meshes.grid1D import Grid1D mesh=Grid1D(dx=dx, nx=nx) from fipy.variables.cellVariable import CellVariable phi=CellVariable(name="saturation", mesh=mesh) p=CellVariable(name="pressure",mesh=mesh) x=mesh.getCellCenters()[...,0] phi.setValue(.9, where=x<.5) phi.setValue(.3, where=x>=.5) p.setValue(1.) #pcs=0 #built sources qin=40*((phi>=.1)and(phi<=.3))+20*((phi>=.7)and(phi<=.9)) qout=60*((phi>=.4)and(phi<=.6)) sin=.8 #pin=0 q1s=sin*qin+s*qout q2s=(1-sin)*qin+(1-phi)*qout from fipy.terms.transientTerm import TransientTerm from fipy.terms.implicitSourceTerm import ImplicitSourceTerm from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm eqp=ImplicitDiffusionTerm(coeff=1)==q1s+q2s eqs=TransientTerm(coeff=2)==q1s-q2s steps=50 for i in range(steps): eqp.solve(p,dt=timeStepDuration) eqphi.solve(phi,dt=timeStepDuration) if __name__ =='__main__': phiviewer=fipy.viewers.make(vars=(phi), limits={'datamin':0., 'datamax':1.}) pviewer=fipy.viewers.make(vars=(p), limits={'datamin':0., 'datamax':1.}) phiviewer.plot() pviewer.plot() ********************************************************************************************************************************************************************************************** 2008/8/18 Jonathan Guyer public.gmane.org>

matrix diffusive coefficient and characteristic source function

franck kalala — 2008-08-18T12:56:42

Dear all, How can solve this problem 1) \nabla(S.\nabla phi)=source where S is a 2 by 2 matrice? 1) \nabla(\nabla phi)=source=q(phi) where source=q(phi)=charactristic function of the intervall [0.3,.7] how do I set a characteristic function in fipy ie my given source. Franck Regards.

Detalization of the mesh

Igor Telezhinsky — 2008-08-17T23:22:03

Dear all, Thank you for previous answers! Now I have a few tricky questions: Lets say I have a mesh mesh = Grid1D(dx = dx, nx = nsteps) x = mesh.getCellCenters()[...,0] x is now a space variable of my problem. How many steps in space is it essential to use for the convection problem? (see PDE solution problem thread for details). In what way the number of steps in space changes the solution precision? I would like to have small steps in space at least near 0 (i.e. 0, 0.001, 0.002 ... ). This is needed for my further research. I experimented a bit and found that fipy slows down dramatically when number of steps in space increases (as N^2)! That's why I thought of logarithmic scale so that all space dependent functions would now depend on 10^x. I met a problem, however: ARE THERE NEGATIVE NUMBERS ON THE MESH? For example, I would like to have a range of x (-3,3) what would respond to the range (0.001,1000) with 1000000 steps if had to do this detalization in linear scale. This number of steps is impossible for fipy ... Could anybody help in solving this issue? Thanks in advance for all answers. Cheers, Igor.

Mesh Problem

Zhiwen Liang — 2008-08-15T17:54:39

Dear Fipy Developers, Thanks again for solving my problem last time. I ran into some error when I tried to translate a mesh to another position. Please see the following lines: Traceback (most recent call last): File "", line 1, in File "/Library/Python/2.5/site-packages/FiPy-2.0a1-py2.5.egg/fipy/meshes/numMesh/uniformGrid2D.py", line 358, in getCellCenters return centers.reshape((2, self.numberOfCells), order="FORTRAN") + self.origin ValueError: shape mismatch: objects cannot be broadcast to a single shape I wonder if anybody has seen this. You might have known it from the error that the fipy version is 2.0a1 and python version is 2.5. I am running this on a Mac OS X 10.5.4. And the Scipy version is 0.6.0, numpy version is 1.0.1. Thank you in advance for any advice. Regards, Zhiwen

the contaminant transport problem

franck kalala — 2008-08-14T00:08:24

Dear all, I am trying to solve with fipy the contaminant transport problem which looks like: Latex: \frac{partial phi}{partial t}=\nabla(S\nabla phi)-\nabla(phi.v) where S is the identity matrix times delta ie S=delta x I , where delta is a given parameter. and v=(v1,v2) a given vector the analytical solution is: phi(x,y,t)=(1/F)exp((-50*(x-x0-v1*t)**2+(y-y0-v2*t)**2) /F) where F=(200*nu*t+1) and nu is a given parameter and x=(x0,y0) a given point. this analytical solution is Gaussian peak starting at the point x=(x0,y0), being transported by the convective field v=(v1,v2), and diffusing. we consider this domain: (0,3)x(0,3) and T=2 I write a program and I got some plots, but these plots don't give what I am expecting. here is my program: I get some problem when the diffusion coefficient is a matrix and not a scalar and the analytical solution is a function of x,y and t ************************************************************************************************************************************** gamma=0.1 nu=.1 convCoeff=v=(0.8, 0.4) x=[0.5,1.35] L=3 nx=30 ny=nx dx=3./nx dy=dx from fipy.meshes.grid2D import Grid2D from fipy import numerix from scipy import * from scipy.linalg import * diffcoeff=S=gamma*numerix.array([[1, 0],[0 ,1]]) #my coefficient here is a matrix, not a scalar as in all example in the documentation I don't know if this a good way to set this matrix in fipy mesh=Grid2D(dx=dx,dy=dy,nx=nx,ny=ny) from fipy.variables.cellVariable import CellVariable phi=CellVariable(name="solution variable", mesh=mesh) #initializing the solution variable to zero #here I am building my equation in fipy, is this correct? from fipy.terms.transientTerm import TransientTerm from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm from fipy.terms.exponentialConvectionTerm import ExponentialConvectionTerm diffTerm= ImplicitDiffusionTerm(coeff=diffcoeff) eqX=TransientTerm()==diffTerm-ExponentialConvectionTerm(coeff=convCoeff, diffusionTerm=diffTerm) xx=mesh.getCellCenters() D=norm(S)# because S is a matrix I take it matrix norm and affect it to D, I don't know if this is correct, to take the norm of S, or maybe there another way to deal with this matrix coefficient xcoord=xx[...,0] ycoord=xx[...,1] phi.setValue(numerix.exp(-50*(xcoord-x[0])**2+(ycoord-x[1])**2))# here I am setting the initial condition at t=0 # is this time step duration correct? where D is the matrix norm of S, should I use the matrix norm of S? timeStepDuration=0.9*dx**2/(2*D) steps=100 t=timeStepDuration*steps phiAnalytical=CellVariable(name="analytical value", mesh=mesh) try: phiAnalytical.setValue((1/(200*nu*t+1))*numerix.exp((-50*((xcoord-x[0]-v[0]*t)**2+(ycoord-x[1]-v[1]*t)**2))/(200*nu*t+1))) except importError: print "The Scipy library is not available to test the solution to the transient diffusion equation" if __name__ =='__main__': from fipy import viewers viewer=viewers.make(vars=(phi, phiAnalytical), limits={'datamin':0, 'datamax':1}) viewer.plot() for step in range(steps): eqX.solve(var=phi,dt=timeStepDuration) if __name__ == '__main__': viewer.plot() print phi.allclose(phiAnalytical, atol = 7e-4) Moreover, I need to use the final time T=2 in my code, Thanks for your help in giving me some hint to deal with this problem Regards Franck

Plotting flux

Rashmi Ranjan Mohanty — 2008-08-13T15:36:50

Hello Everyone, How to plot mass flux across the mesh in one dimension? flux = M * chemPotential.getFaceGrad() It may be a really trivial question,but I could not plot it. I could print it to a file using TSVViewer, but even then it prints flux value at mesh nodes rather than cell centers. Thanks. Rashmi

Source plus Convection Equation

Igor Telezhinsky — 2008-08-12T09:13:24

Dear all! Thank to all of you for the answers to my previous questions! I'm solving the convection equation (see PDE solution problem thread for details). My equation looks like eq = TransientTerm() == PowerLawConvectionTerm(convCoeff) Now I want to add a VARIABLE SOURCE to it. A constant source is described in manual and works fine. To add a variable source I do the following: #define a function that describes my source def S(x): E=x if E>10: return 5e-7*E**alpha*exp(-E/Ecut) else: return 0 #fipy said it wants a CellVariable for the source term Source = CellVariable(mesh=mesh) Source.setValue(S(x)) #change the equation: eq = TransientTerm() == PowerLawConvectionTerm(convCoeff) + Source Problems/Questions: 1. Is it correct way to introduce the space variable dependent Source to my problem? 2. I have if/else condition in S(x) function. The script works without if/else with no errors, but with condition the following error appears: File "conv_s_1208.py", line 59, in Source.setValue(S(x)) File "conv_s_1208.py", line 48, in S if E>10: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() It seems that S(x) gets the whole array of x values of the mesh and in this case can't compare. Originally I wouldn't like that S(x) function was designed to take an array argument. Shouldn't the setValue method take care of it? (i.e. if the function is written for the scalar argument and an array x[] is passed, setValue should make a new array where y[i]=S(x[i]). This new array will be then the Source term.) If setValue method can't take care of it, is there a way to overcome this difficulty and do not change S(x) function? Is there a way in python to check what object (array, float ...) is a passed to the function as an argument? Thanks to everybody for the answers!