Dear all!
Thank you, it seems the solution is physical (at least the shape). I need to to check it with numbers, however, to be sure it's a proper one... I use fipy revision 2718. The corrections to the Daniel's script I applied to make it work: #x, = mesh.getCellCenters() x = mesh.getCellCenters()[...,0] #convCoeff = FaceVariable(mesh=mesh, value=L(mesh.getFaceCenters()), rank=1) convCoeff = VectorFaceVariable(mesh=mesh, value=L(mesh.getFaceCenters())) fipy doesn't know what "rank=1" means. fipy wanted VectorFaceVariable ... The boundary conditions used in the script are correct. I need to have zero flux at boundaries. By the way, in the solution there are values of order 1e-322 positive as well as negative ones. I guess I can count them as 0, can I? Thanks a lot for your help!
