Right after "from fipy import *" put in "from fipy.terms.diffusionTerm import DiffusionTermNoCorrection as DiffusionTerm" and replace the solver "LinearBicgstabSolver(tolerance=1e-9, iterations=500)" with "LinearLUSolver(tolerance=1e-9, iterations=500)" and things should be cool. You should be able to change the solver to what you want, but LinearBicgstabSolve threw an error when I tried it.
I noticed a few things, one of the cells (74) only has one neighbor, this is one of the cells that keeps growing in value and diverging. The diffusion term I suggested doesn't correct for non-orthogonality so there will be some error due to that, but it also won't diverge. The non-orthogonality correction could probably do with more faces so it can get a better read on the tangential component of the gradient. BTW This is entirely speculation as to the reason for the issue at this point. Also the non-orthogonality correction (as well as the anisotropy) has an explicit contribution (the tangential contribution to the flux) so time steps could have some upper stability condition as well. This may have contributed to the divergence. Anyhow, the changes I suggested above should (finger's crossed) get you some results. It may actually be worth changing the default DiffusionTerm to the uncorrected version as the default for fipy. Let me think about that. On Thu, Feb 5, 2009 at 4:06 PM, A.S.Reeve <[email protected]> wrote: > Daniel, > > If the importer is messing things up, it only does it for gmsh generated > grids. I've attached my rectilinear grid that follows the gmsh file > conventions. It imports fine and runs well with heterogeneous grids. > > Andy > > Andrew Reeve > Associate Prof. > Dept. of Earth Sciences > University of Maine > 207-581-2353 > > On Thu, 5 Feb 2009, Daniel Wheeler wrote: > >> >> Thanks for that. We actually need to update fipy to use mayavi 2. >> >> I removed the extrude function and ran the problem as a plain 2D >> problem and it diverges just the same (I lost my wager). For this >> problem the matrix should be symmetric so I'll check that first. I >> doubt that the non-orthogonality would cause the solution to diverge, >> it makes me think that the gmsh importer is messing something up. >> >> On Thu, Feb 5, 2009 at 11:35 AM, A.S.Reeve <[email protected]> wrote: >>> >>> Daniel, >>> >>> Here's the bit of code (below) I'm using to sort the cell vertices so >>> hexahedrons can be plotted in mayavi2. I don't think this is very >>> general, >>> and it might not work properly if the hexahedrons are far from >>> rectilinear. >>> I'm not using the fipy viewer, and instead create vtk files to be viewed >>> externally. >>> >>> The rectilinear grid I created that seems to be working was also >>> extruded, >>> which is why I assumed the extrusion routine was O.K. >>> >>> I'm a bit under the gun to get this working and turn in some results. >>> I'll >>> look into making a 3-D grid with gmsh, once I clear this project off my >>> plate. >>> >>> Andy >>> >>> cellVerts=mesh._getCellVertexIDs().data >>> vertCoords=mesh.getVertexCoords() >>> >>> #order=[0,1,2,3,4,5] >>> order=[0,1,3,2,4,5,7,6] >>> >>> cellVerts2=array(cellVerts[:]) >>> #for i in range(8): >>> # cellVerts2[i,:]=cellVerts[order[i],:] >>> >>> cellVerts2=transpose(cellVerts2) >>> vertCoords=transpose(vertCoords) >>> >>> dtype=[('x',float),('y',float),('z',float),('id',float)] >>> ###FiPy mixes up cell verts...need to sort them to be able to plot >>> ###sensibly in mayavi >>> for i in range(mesh.getNumberOfCells()): >>> IDs=array([[j] for j in cellVerts2[i]]) >>> cell=[vertCoords[j] for j in cellVerts2[i]] >>> cell=hstack((cell,IDs)) >>> cell=array([tuple(cell[j]) for j in range(8)],dtype=dtype) >>> cell.sort(order=['z','x','y']) >>> cellVerts2[i]=array([cell[j][3] for j in order]) >>> >>> VisGrid=VtkData( >>> UnstructuredGrid(vertCoords,hexahedron=cellVerts2) >>> ,CellData(Scalars(head,'head')) >>> ) >>> VisGrid.tofile('head_data','ascii') >>> >>> >>> Andrew Reeve >>> Associate Prof. >>> Dept. of Earth Sciences >>> University of Maine >>> 207-581-2353 >>> >>> On Thu, 5 Feb 2009, Daniel Wheeler wrote: >>> >>>> >>>> On Thu, Feb 5, 2009 at 9:05 AM, A.S.Reeve <[email protected]> wrote: >>>>> >>>>> I re-meshed my problem using a rectilinear grid, used the same script, >>>>> and >>>>> everything works fine now. Now that I've worked through this, I recall >>>>> having this issue before. It seems like the gmsh generated grids only >>>>> work >>>>> for homogeneous or nearly homogeneous problem domains. If this is true, >>>>> perhaps there should be some sort of a warning issued when importing >>>>> gmsh >>>>> grids? >>>> >>>> Hi Andrew, Hope all is well. I'd wager on the extrude function rather >>>> than gmsh import as being the source of the problem. I wrote that when >>>> you were here that time and I haven't looked at it since. One thing to >>>> do is to create a 3D mesh in gmsh and use that without using the >>>> extrude function. Have you already tried that? I'm going to look into >>>> it, but anything you have will be helpful. >>>> >>>>> I also recall there being problems plotting hexahedrons in mayavi. Is >>>>> this >>>>> true, or am I imagining this? If its true, I could send you the bit of >>>>> code >>>>> I'm using to plot up a rectilinear hexahedral grid for inclusion in >>>>> FiPy. >>>> >>>> Yeah, we had problems with cell vertex ordering when we wrote the >>>> mayavi viewer. I'm not sure we ever sorted that one out so any code >>>> you have that deals with that issue would be very helpful. Thanks! >>>> >>>> -- >>>> Daniel Wheeler >>>> >>> >>> >> >> >> >> -- >> Daniel Wheeler > -- Daniel Wheeler
