The main program sets up geometry of solution domain for 2D advective diffusion
with uniform density rho.
Conductivity Gamma is now a function of phi, and the solution is achieved iteratively.
 The boundaries of this solution domain are on finitevolume faces,
but additional nodes are added on the boundaries (Patankar's Practice B, p.69).
 Boundary conditions can be either fixed phi or fixed phi gradient at each face on the boundary.
For example, in bc_e, first column is numerical value of BC,
second column specifies BC type  1 for value, 2 for gradient.
A value of 2 indicates a specified diffusive flux,rather than a phi gradient
(that's because Gamma is changing.)
 Gamma_phi mfile uses 1/phi dependence
 Gamma_phi_bad mfile uses linear dependence
on phi (can get negative Gamma)
 u.m and w.m find velocities on appropriate
volume boundaries

A.m and F_upwind.m implement powerlaw advection scheme
solve.m solves the resulting system of linear equations.