Should be possible to couple scalar and vector equations

[guyer]

As illustrated

If we attempt to represent

$$ \begin{aligned} \frac{\partial C}{\partial t} &= \nabla\cdot(DC \nabla C) + \nabla\cdot(C \vec{v}) \ \frac{\partial \vec{v}}{\partial t} &= \nabla\cdot(D\vec{v} \nabla \vec{v}) \end{aligned} $$

as

import fipy as fp

mesh = fp.Grid2D(nx=10, ny=10)
sca = fp.CellVariable(mesh=mesh, name="sca", hasOld=True)
vec = fp.CellVariable(mesh=mesh, name="vec", rank=1, hasOld=True)

Dsca = 1
Dvec = 2

eq_sca = fp.TransientTerm(var=sca) == fp.DiffusionTerm(coeff=Dsca, var=sca) + fp.ConvectionTerm(coeff=vec, var=sca)
eq_vec = fp.TransientTerm(var=vec) == fp.DiffusionTerm(coeff=[[[Dvec, 0],
                                                               [0, Dvec]]], var=vec)

eq = eq_sca & eq_vec

eq.solve(dt=1)

FiPy raises

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 1 dimension(s) and the array at index 1 has 2 dimension(s)

Solving sequentially is fine:

eq_sca.solve(dt=1)
eq_vec.solve(dt=1)

For stability reasons, it may be necessary to use a sequential approach, anyway (e.g. SIMPLE), but FiPy shouldn't prevent trying to solve coupled.