On Fri, Jan 4, 2013 at 2:51 PM, Salomon Turgman Cohen <[email protected]>wrote:
> Hello developers and users,
> I am interesting in solving a simple 1D diffusion equation using fipy.
> My requirement is that the diffusion coefficient is concentration
> dependent. Following the mesh1D example I know how to do this if I have a
> simple functional form to describe the diffusion's dependence on
> concentration. However, what I have is just 5 data points of diffusion
> versus concentration. What I want to do is a linear interpolation of the
> diffusion in between the data points. I constructed the following function
> as an example:
>
> def diffusion(phi):
> if (phi>=0.0) & (phi<0.902449):
> return phi*(-7.63256E-8) + 1.05240E-7
> elif (phi >= 0.902449) & (phi < 1.807062):
> return phi*(-1.37296E-08) + 4.87503E-08
> elif (phi >= 1.807061587) & (phi < 2.713202143):
> return phi*(-2.23797E-08) + 6.43815E-08
> elif (phi >= 2.713202143) & (phi < 3.570465711):
> return phi*-1.48121E-09 + 7.67966E-09
> elif (phi >= 3.570465711) & (phi < 4.425457264):
> return phi*-4.04495E-10 + 3.83529E-09
> elif (phi >= 4.425457264):
> return phi*-4.62147E-10 + 4.0904E-09
>
> But it cannot handle the arrays obtained from CellVariables command. What
> would be the right way to do this?
>
Hi Salomen,
Sorry for taking so long to reply. The function above won't work at all as
it is not in a vectorized form. There are two issues to think about. The
first it to ensure that a vector operation is used, the second is to ensure
that the dependency of the diffusion coefficient on "phi" is preserved. The
solution to this is a slightly unintuitive, but consists of using ">" and
"*" operators in combination as follows:
def cases(var, values, limits):
assert len(values) == (len(limits) + 1)
newVar = var.copy()
newVar[:] = values[0]
for value, limit in zip(values[1:], limits):
newVar = (var > limit) * (value - newVar) + newVar
return newVar
if __name__ == '__main__':
import fipy as fp
mesh = fp.Grid1D(nx=10, Lx=1.)
var = fp.CellVariable(mesh=mesh)
var[:] = mesh.x
v0 = cases(var, (var / 2, 2., 3.), (0.5, 1.5) )
print v0
var[:] = mesh.x * 2
print v0
var[:] = mesh.x * 3
print v0
This seems to work as you would expect and satisfy the requirements of
variable dependency and vectorization. This is not the most efficient
solution. A more efficient solution would involve subclassing from
CellVariable and using "numpy.where" or fancy indexing to create masks,
rather than using the "*" operators. This would require some more knowledge
of the FiPy internals,.
Cheers
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
