Re: System of partial, differential and algebraic equations

Mon, 17 Dec 2007 09:47:33 -0800
You already have a (good) answer from Daniel, but here's my take on this: 


On Dec 17, 2007, at 8:14 AM, Etienne Rivard wrote:


I would like to use FiPy for my research. I have set of equations  
to solve that consists mainly of partial differential equations.  
However, it also includes ordinary differential equations and  
algebraic equations. Can FiPy handle that.


Daniel answered this more fully, but the short answer is "yes".


If so, could you point me towards useful documentation?



As Daniel said, our manual includes lots of examples ranging from a  
simple diffusion equation in 1D to coupled sets of PDEs in higher  
dimensions. What exactly do you need to know that isn't covered in  
the examples?



It is a one-dimensional problem. Therefore, I have two independent  
variables: a space coordinate and time. One of the terms contains a  
derivative with respect to position evaluated at a specific  
position. Something like (dy/dr)|R Is that definable in FiPy?


Assuming you have FiPy 1.2, make a 1D mesh:

    >>> from fipy import *
    >>> mesh = Grid1D(nx=100, dx=0.1)

get the coordinates of that mesh:

    >>> r = mesh.getCellCenters()[...,0]

create a y variable with some function of r as its initial value

    >>> y = CellVariable(mesh=mesh, value=r**2)

pick some point R

    >>> R = 4

Get (dy/dr)|R

    >>> y.getGrad()(point=(R,))
    array([ 7.9])


This is actually $\nabla y$, which is a 1D vector. To get the scalar  
value dy/dr, you'd do:


    >>> y.getGrad()(point=(R,))[0]
    7.9


Note that this is not the analytical answer of 2r = 8 because FiPy  
uses cell-centered finite volume. For dx=0.1, the cell centers lie at  
0.05, 0.15, 0.25, etc. The "(point=)" syntax interpolates to find the  
first nearest cell center (and is thus pretty slow, although no big  
deal in 1D), which is at r=3.95 in this case. A much faster way (but  
with no interpolation) is


    >>> dr = 0.1
    >>> y.getGrad()[r == R + dr / 2]
    array([[ 8.1]])


This is actually an array of *all* the points that satisfy r == R +  
dr / 2. Since we're interested in only one (and there is only one)


    >>> y.getGrad()[r == R + dr / 2][0]
    array([ 8.1])
    >>> y.getGrad()[r == R - dr / 2][0]
    array([ 7.9])

give the (dy/dr) at the two points nearest to R = 4.


You can also get the gradient at the face between cells, which  
happens to sit exactly at r = R, with:


    >>> y.getFaceGrad()[rF == 4][0]
    array([ 8.])


Although you can write the equivalent y.getFaceGrad()(point=(R,)),  
that code appears to be buggy.


I don't know enough yet about your equations to know how you'll be using

In all of these cases, the value obtained will not change as FiPy  
evolves y in time. We have code in development that does a better job  
of this sort of thing, but, in the meantime, assuming you wanted to  
solve something like the convective equation $\frac{d y}{d t} = \nabla 
\cdot(y \frac{d y}{d t}|_R)$, you would need to do something like:



    >>> dydrAtR = VectorFaceVariable(mesh=mesh, value=y.getGrad() 
(point=(R,)))

    >>> eq = TransientTerm() == PowerLawConvectionTerm(coeff=dydrAtR)
    >>> for step in range(steps):
    ...     eq.solve(var=y, dt=0.01)
    ...     dydrAtR.setValue(y.getGrad()(point=(R,)))


This is fairly nonlinear, so you'd likely need to "sweep" this  
problem, too. See FAQ 6.5 in our manual for more information on that.

























					Reply via email to