On May 15, 2012, at 5:15 PM, Yun Tao wrote:

> As expected, only the CentralDifference method fails under a high Peclet 
> value. By the way, is there a way for FiPy to calculate and display that? 
> According to the FiPy doc 
> http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#linear-equations,
>  the advection strength requires the dot product of the advection component 
> and the normal vector at each face, yet the latter value is not obvious to me 
> as an end user not knowing the topology of the solution domain. 

The actual calculation of the Peclet number is performed at

http://matforge.org/fipy/browser/branches/version-2_1/fipy/terms/convectionTerm.py#L143

The calculation of the advection strength is at 

http://matforge.org/fipy/browser/branches/version-2_1/fipy/terms/convectionTerm.py#L123

That code is really old and so now I'd write

  numerix.dot(self.coeff, mesh._getOrientedAreaProjections())



As to visualization, that's going to be tricky. The Peclet number you calculate 
this way is a scalar FaceVariable and we don't have any tools to visualize such 
a thing. I'd be inclined to just visualize the vector field of u * dAP /D (you 
can get dAP with mesh._getCellDistances()). Hmmm... that will tend to 
over-represent advection vectors that lie parallel to a face. Maybe better to 
calculate the scalar Peclet number and then multiply it by the face normals to 
get a vector field that FiPy can plot, e.g.,

>>> conv = -CentralDifferenceConvectionTerm(coeff=b)
>>> diff = DiffusionTerm(coeff=epsilon)
>>> eq = conv + diff == 0

>>> Viewer(vars=conv._getGeomCoeff(mesh) / eq._getDiffusiveGeomCoeff(mesh) * 
>>> mesh._getOrientedFaceNormals())
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Reply via email to