File Exchange

image thumbnail

A Simple Finite Volume Solver for Matlab

version 2.0 (1.77 MB) by

A simple yet general purpose FVM solver for transient convection diffusion PDE

12 Ratings



A simple Finite volume tool
This code is the result of the efforts of a chemical/petroleum engineer to develop a simple tool to solve the general form of convection-diffusion equation:
on simple uniform/nonuniform mesh over 1D, 1D axisymmetric (radial), 2D, 2D axisymmetric (cylindrical), and 3D domains.
The code accepts Dirichlet, Neumann, and Robin boundary conditions (which can be achieved by changing a, b, and c in the following equation) on a whole or part of a boundary:
It also accepts periodic boundary conditions.
The main purpose of this code is to serve as a handy tool for those who try to play with mathematical models, solve the model numerically in 1D, compare it to analytical solutions, and extend their numerical code to 2D and 3D with the minimum number of modifications in the 1D code.
The discretizaion schemes include
* central difference
* upwind scheme for convective terms
* TVD schemes for convective terms with many flux limiters
To get started, go to the `Test` folder and run the test scripts.
A few calculus functions (divergence, gradient, etc) and averaging techniques (arithmetic average, harmonic average, etc) are available, which can be helpful specially for solving nonlinear or coupled equations or implementing explicit schemes.
I have used the code to solve coupled nonlinear systems of PDE. You can find some of them in the Examples/advanced folder.

There are a few functions in the 'PhysicalProperties' folder for the calculation of the physical properties of fluids. Some of them are not mine, which is specified inside the file.

I'll try to update the documents regularly, in the github repository. Please give me your feedback/questions by writing a comment in my weblog: <>
Special thanks: I vastly benefited from the ideas behind Fipy <>, a python-based finite volume solver.

To start the solver, download and extract the zip archive, open and run 'FVToolStartUp' function.
To see the code in action, copy and paste the following in your Matlab command window:

clc; clear;
L = 50; % domain length
Nx = 20; % number of cells
m = createMesh3D(Nx,Nx,Nx, L,L,L);
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=1; % Dirichlet for the left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % Dirichlet for the right boundary
D_val = 1; % value of the diffusion coefficient
D = createCellVariable(m, D_val); % assign the diffusion coefficient to the cells
D_face = harmonicMean(D); % calculate harmonic average of the diffusion coef on the cell faces
Mdiff = diffusionTerm(D_face); % matrix of coefficients for the diffusion term
[Mbc, RHSbc] = boundaryCondition(BC); % matix of coefficients and RHS vector for the BC
M = Mdiff + Mbc; % matrix of cefficients for the PDE
c = solvePDE(m,M, RHSbc); % send M and RHS to the solver
visualizeCells(c); % visualize the results

You can find some animated results of this code in my youtube channel:

Comments and Ratings (22)

Hi all,
Please ask your questions by opening an issue in the github repository ( orby writing a comment in my blog ( I don't have access to this page anymore, since my previous email at TU Delft is not active.
I'm adding a few examples to the github page to answer your questions regarding the Poisson equation and the source terms.
In short, you can have a derivative boundary condition by changing the BC lines in the above code to:
BC.left.a(:) = 1; BC.left.b(:)=0; BC.left.c(:)=1; % Neumann for the left boundary
BC.right.a(:) = 1; BC.right.b(:)=0; BC.right.c(:)=0; % Neumann for the right boundary

A source term can be added by calling the function constantSourceTerm(q) where q is a cell variable. Don't forget to divide the source term by the cell volume.


Just want to say thx

Chu-Yao Chou

Hi I would like to know how to add a source term to the diffusion equation using your solver.


hi, how I can solve 3d poison equation with derivative boundary conditions in pdetool in rectangular channel domain?

Hi Shijie,

I don't visit this page regularly. You can always download the latest version from github:


Shijie liu

Hi Ali,
I don't know why I can't download the attached .zip file. And I have tried many times. Could you send the files to me? My email is Thanks a lot!

Saif Manji

Yuri Feldman

Hi Ali,
Can you please explain the following corrections in your upwind scheme for convection term:
% Also correct for the boundary cells (not the ghost cells)
% Left boundary:
APx(1,:) = APx(1,:)-uw_max(1,:)/(2*DXp(1)); AW(1,:) = AW(1,:)/2;
% Right boundary:
AE(end,:) = AE(end,:)/2; APx(end,:) = APx(end,:) + ue_min(end,:)/(2*DXp(end));
% Bottom boundary:
APy(:,1) = APy(:,1)-vs_max(:,1)/(2*DYp(1)); AS(:,1) = AS(:,1)/2;
% Top boundary:
AN(:,end) = AN(:,end)/2; APy(:,end) = APy(:,end) + vn_min(:,end)/(2*DYp(end));

I will also appreciate if you could refer me to the relevant literature.
Thank you,

Hi Everyone,

Please write your questions or comments preferably in the Github page of the code. My mathworks account is changed and I don't receive notifications anymore.

Sorry Alireza for this late reply. I have a new mathworks account so I don't receive notifications anymore. Yes, you can define a heterogeneous field. For instance, in the example above, replace the mesh creation and D definition lines with:
m=createMesh2D(Nx, Nx, L, L);
D=createCellVariable(m, rand(Nx, Nx));
You can replace the rand(Nx, Nx) with any matrix of Nx x Nx size.




AG (view profile)

Thanks for the magnificent work. Does your code support heterogeneous material properties as well? I am trying to solve a 2D transient heat equation on a domain that has different conductivities and heat capacities and I was hoping your framework could be of help.

Many thanks!

Hi Hongwei, I'm glad you find the code useful. Please let me know if you like to add your application to the example folder. You are always welcome to send a pull request on github.

Hongwei Guo

Thanks a lot !!! Very professional and general code !!! I will try to apply this in electron transport problems !!!

Hongwei Guo

Guoxi HE

Peng Cao

Nikhil Kumar


Martin (view profile)



showdemo function is not available. I will update it later.


added support for 2D radial (r, theta) and 3D cylindrical (r, thetta, z)


updated descriptions


update my weblog address


add youtube channel link to descriptions

MATLAB Release
MATLAB 8.3 (R2014a)

Inspired by: IAPWS_IF97 functional form with no slip

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video