Sorry,
Here is the code. In this case, if we just set velo=epsilon.grad, then
the code works well. But if we use velo = DOT(D,epsilon.grad), it doesn't.
(Actually, D is the unity Matrix ((1,0),(0,1)).)
Zebo
2014-06-06 9:28 GMT-05:00 Daniel Wheeler <[email protected]>:
> On Thu, Jun 5, 2014 at 12:23 PM, Zebo LI <[email protected]> wrote:
> > Hi,
> >
> > I currently working on a diffusion problem using Fipy. And I need to
> > calculate the convection velocity via the product of a matrix and a
> vector.
> > I used numeric.array and numeric.NUMERIX.dot, but it does not work well,
> as
> > shown in the attached file.
>
> You probably just need to reshape the arrays properly.
>
> Any chance you could cut and paste the whole code or subset that runs
> and produces the error so I can reproduce the error on my computer
> (and not as an image preferably)?
>
> Cheers,
>
> Daniel
>
> --
> Daniel Wheeler
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
#!/usr/bin/env python
r"""Solve a two-dimensional diffusion problem in a square domain.
.. index:: Grid2D
>>> from fipy import *
>>> nx = 100
>>> ny = nx
>>> dx = 0.1
>>> dy = dx
>>> L = dx * nx
>>> mesh = Grid2D(dx=dx, nx=nx , dy = dy, ny = ny)
We create a :class:`~fipy.variables.cellVariable.CellVariable` and initialize it to zero:
>>> phi = CellVariable(name = "solution variable",
... mesh = mesh,
... value = 0.)
and then create a diffusion equation. This is solved by default with an
iterative conjugate gradient solver.
>>> X= mesh.cellCenters[0]
>>> Xface = mesh.faceCenters[0]
>>> Y= mesh.cellCenters[1]
>>> Yface = mesh.faceCenters[1]
>>> mu0 = 0
>>> epsilon = 3.14 /((X-5)**2+(Y-5)**2+1)
>>> epsilonface = 3.14 /((Xface-5)**2+(Yface-5)**2+1)
>>> D0 = 4.5
>>> D = D0*numerix.array(((1,0),(0,1)))
>>> DOT = numerix.NUMERIX.dot
velo = DOT(D,epsilon.grad)
>>> velo = epsilon.grad
>>> eq = TransientTerm(coeff=1.0,var=phi) == DiffusionTerm(coeff=D,var=phi) + ConvectionTerm(coeff=velo,var=phi)
>>> valueRight = numpy.exp(mu0-epsilonface)
>>> valueTop = numpy.exp(mu0-epsilonface)
>>> valueLeft = numpy.exp(mu0-epsilonface)
>>> valueBottom = numpy.exp(mu0-epsilonface)
>>> phi.constrain(valueRight, mesh.facesRight)
>>> phi.constrain(valueLeft, mesh.facesLeft)
>>> phi.constrain(valueTop,mesh.facesTop)
>>> phi.constrain(valueBottom,mesh.facesBottom)
>>> phi.setValue(numpy.exp(mu0-epsilon))
>>> viewer = Viewer(vars=phi, datamin=0.000, datamax=1.00)
and solve the equation by repeatedly looping in time:
>>> timeStepDuration = 10 * 0.9 * dx**2 / (2 * D0)*10
>>> steps = 100
>>> for step in range(steps):
... eq.solve(var=phi,
... dt=timeStepDuration)
... viewer.plot()
>>> viewer._promptForOpinion()
"""
__docformat__ = 'restructuredtext'
##from fipy.tools.profiler.profiler import Profiler
##from fipy.tools.profiler.profiler import calibrate_profiler
if __name__ == '__main__':
import fipy.tests.doctestPlus
import numpy
exec(fipy.tests.doctestPlus._getScript())
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
