How to write equation terms with independent variable in FiPy?

Question

I am looking to solve a diffusion equation using FiPy and have read some of their documentation, but can't seem to find anything that relates to writing a diffusion term that includes additional terms that are functions of the independent variable (i.e. space). The closest thing that I found was on the FAQ, where they suggest rewriting additional terms as a ConvectionTerm. However, I believe this only applies to the case where the additional terms are functions of the solution variable rather than the independent variable. For example, I am trying to solve a 1D diffusion equation with the following diffusion term (where derivatives are w.r.t. to the independent variable x, and y is the solution variable):

D * sin(x) * Div_x {sin(x) * Grad_x {y}}

I feel that this is a pretty simple expression, but I can't find how to express it in FiPy notation. Any help would be hugely appreciated!

Exact Code:

from fipy import Variable,FaceVariable,CellVariable,Grid1D,ImplicitSourceTerm,TransientTerm,DiffusionTerm,Viewer,ConvectionTerm
from fipy.tools import numerix

D = 1
c0 = 1
ka = 1
r0 = 1

nx = 100
dx = 2*math.pi/100
mesh = Grid1D(nx=nx, dx=dx)
conc = CellVariable(name="concentration", mesh=mesh, value=0.) # This is the "phi" in the docs
valueLeft = c0
valueRight = 0
conc.constrain(valueRight, mesh.facesRight)
conc.constrain(valueLeft, mesh.facesLeft)
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100
show_per_steps = 50

A = 1 / (r0**2 * numerix.sin(mesh.x)[0])
dA = -(numerix.cos(mesh.x)[0])/(r0**2 * numerix.sin(mesh.x)[0]**2)
dsindA = (numerix.cos(mesh.x)[0])**3/(numerix.sin(mesh.x)[0])**2
eqX = TransientTerm() + ImplicitSourceTerm(ka) == DiffusionTerm(D*A*numerix.sin(mesh.x)[0]) - ConvectionTerm(D*dA*numerix.cos(mesh.x)[0])+ D*conc*dsindA

from builtins import range
for step in range(steps):
    eqX.solve(var=conc, dt=timeStepDuration)
    if __name__ == '__main__' and step % show_per_steps == 0:
        viewer = Viewer(vars=(conc), datamin=0., datamax=c0)
        viewer.plot()

wd15 · Accepted Answer · 2019-09-20T14:14:33.807

2

FiPy allows the terms' coefficients to be functions of space. For example, the following works in FiPy,

from fipy import Grid1D, CellVariable, Viewer
from fipy import TransientTerm, numerix, DiffusionTerm
from fipy import LinearLUSolver

m = Grid1D(nx=100, Lx=numerix.pi / 4.)

v = CellVariable(mesh=m)
v[:] = m.x**2

eqn = TransientTerm() == DiffusionTerm(numerix.sin(m.x))
vi = Viewer(v, colorbar=None)
vi.plot()
solver = LinearLUSolver()
for i in range(10):
    eqn.solve(v, dt=0.1, solver=solver)
    vi.plot()
    print('step', i)

input('stopped')

In the above, the diffusion coefficient is a function of space. The m.x is a CellVariable that holds the cell center positions. numerix is used which enables operations on FiPy variables in the same way as Numpy does for Numpy arrays.

Now, in the above question there is a sin(x) outside of the derivative which is not allowed in FiPy. Everything needs to fit inside of the derivative to work with FiPy. So, we need to rewrite the term so that all the coefficients are inside of the derivative. For any general case, we can write

$f \partial_i \left( g \partial_i \phi \right ) = \partial_i \left(g f \partial_i \phi \right ) - \partial_i \left(g \phi \partial_i f \right ) + \phi \partial_i \left( g \partial_i f)$

which allows us to use a diffusion, convection and source to represent the term in FiPy. If f=sin(x) and g=sin(x) then the FiPy code would be

s = numerix.sin(m.cellCenters)
c = numerix.cos(m.cellCenters)
eqn = ... + DiffusionTerm(D * s[0] * s[0]) - ConvectionTerm(D * s * c) + D * y * (c[0] * c[0] - s[0] * s[0])

The ... are included as I don't know the full equation.

edited Sep 20 '19 at 14:14

answered Sep 17 '19 at 14:57

wd15

1,068
5
8

Thanks for the answer! This was exactly what I was looking for. However, when I use numerix.sin(mesh.x) as the coefficient in ConvectionTerm, I am getting the error: "VectorCoeffError: The coefficient must be a vector value." It works fine for DiffusionTerm. – tooty44 Sep 19 '19 at 03:28
Sorry about that. FiPy is a little confusing in 1D as the convection term expects a vector rather than a scalar for its coefficient, which is irrelevant in 1D, but does require the coefficient to have shape (1, nx). I've edited the answer to reflect this. The `s` and `c` variables are now vector quantities, (1, nx), rather than scalars, (nx,). Hence, the diffusion and source terms use `s[0]` and `c[0]` to get scalar quantities. – wd15 Sep 20 '19 at 14:19
Thanks, I think adding the indexing helped clear up some of the errors, but for some reason I still get `VectorCoeffError: The coefficient must be a vector value` when solving the equation. I've added my code to the question in case that it might help clarify the issue. I am trying to solve diffusion on a spherical surface from a point source and so have everything is terms of the altitude angle, theta. – tooty44 Sep 21 '19 at 03:38
Here's a [working version of your code](https://gist.github.com/wd15/54f6d412a6f5e9f5ad5abd2f73549842). – wd15 Sep 23 '19 at 19:54

How to write equation terms with independent variable in FiPy?

1 Answers1