I'm using FiPy to model convection/diffusion of a chemical in a given velocity field, and I'm having problems setting this field as ConvectionTerm coefficient.
My velocity field is expressed by two bidimensional arrays, let's say Ugrid for horizontal velocities and Vgrid for vertical ones, representing the values at the faces of each cell. For this reason, I'd say that the proper way to set this field as coeff parameter for the ConvectionTerm would be assigning it to a FaceVariable.
However, I don't get how to pass these two arrays as value for the FaceVariable. Apparently, I have no problems in setting the field as CellVariable by using value = [np.ravel(Ugrid),np.ravel(Vgrid)] and the convection simulation appears to make sense as well, but I don't think this is correct, as I briefly mentioned above.
Any suggestions?
An important consideration with FiPy is that the mesh indexing is not based on grid indexing so mapping from a grid array to cell values should not assume that the mesh is ordered in any particular way. This makes it necessary to use some form of interpolation. Let's first construct a velocity field with grid indexing where the points lie on the grid points (not cell centers).
from fipy import Grid2D, CellVariable, FaceVariable
from fipy import ConvectionTerm, DiffusionTerm, TransientTerm
import numpy as np
from scipy import interpolate
dx = 0.5
nx = 7
dy = 0.5
ny = 7
xy = np.zeros((nx, ny, 2))
xy[:, :, 0] = np.mgrid[0:nx, 0:ny][0] * dx
xy[:, :, 1] = np.mgrid[0:nx, 0:ny][1] * dy
u = xy[..., 0] + xy[..., 1]
v = xy[..., 0] * xy[..., 1]
We now have a (u, v) velocity field each of shape nx, ny and the corresponding coordinates, xy, of shape (nx, ny, 2). Let's say we want to interpolate that to a FiPy mesh with the same domain, but a different grid.
m = Grid2D(nx=3, ny=3, dx=1.0, dy=1.0)
The mesh doesn't necessarily need to align with the grid points for the velocity field. We can then interpolate to that mesh with
xy_interp = np.array(m.faceCenters).swapaxes(0, 1)
u_interp = interpolate.griddata(xy.reshape(-1, 2), u.flatten(), xy_interp, method='cubic')
v_interp = interpolate.griddata(xy.reshape(-1, 2), v.flatten(), xy_interp, method='cubic')
where xy_interp are the mesh's face centers. Note that using griddata requires that the xy_interp are inside the xy otherwise it gives nan values. Once we have the interpolated values, we can set up the FiPy velocity field.
velocity = FaceVariable(mesh=m, rank=1)
velocity[0, :] = u_interp
velocity[1, :] = v_interp
Note that the coefficient for the ConvectionTerm can be a FaceVariable or a CellVariable. Once we have the velocity, we're ready to set up and solve the equation.
var = CellVariable(mesh=m)
eqn = (TransientTerm() + ConvectionTerm(velocity) == DiffusionTerm(1.0))
eqn.solve(var, dt=1.0)
This runs without error for me.