simulkade/fvtool

Finite volume toolbox for Matlab/Octave

FVTool: Finite volume toolbox for Matlab

Tiny Documents 📘.

FVTool in:

Python: PyFVTool
Julia: JFVM.jl

FVTool

This is a finite volume (toy) toolbox for chemical/petroleum engineers. Right now, it can solve a transient convection-diffusion equation with variable velocity field/diffusion coefficients. The discretization schemes include:

• central difference diffusion term
• central difference convection term
• upwind convection term
• TVD convection term with various flux limiters
• transient term
• Dirichlet, Neumann, Robin, and periodic boundary conditions

Which equation do you solve?

You can solve the following PDE (or a subset of it):

with the following boundary conditions:

Believe it or not, the above equations describe the majority of the transport phenomena in chemical and petroleum engineering and similar fields.

How to start

Download the package, start matlab, and run FVToolStartUp

Inspiration

I started writing this tool after playing with FiPy, an amazing python-based finite volume solver. This matlab solver is not a clone, and indeed very limited compared to FiPy. I wrote it to have a very handy tool for testing new ideas (new mathematical models) by solving them in 1D uniform Cartesian grids. Then I extended the code to

• 1D axisymmetric (radial)
• 2D radial (r, theta)
• 2D Cartesian
• 3D Cartesian
• 2D axisymmetric (cylindrical, r, z)
• 3D cylindrical (r, theta, z)

I have overloaded some of the matlab operators to simplify the switch from 1D codes to 2D and 3D.

A simple example

You can solve a diffusion equation, i.e., $ \nabla. (-D \nabla \phi) = 0 $ by running the following code in Matlab:

clc
L = 50;  % domain length
Nx = 20; % number of cells
m = createMesh1D(Nx, 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; % 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); % matrix of coefficients and RHS vector for the BC
M = Mdiff + Mbc; % matrix of coefficients for the PDE
c = solvePDE(m,M, RHSbc); % send M and RHS to the solver
visualizeCells(c); % visualize the results

change the third line to m = createMesh2D(Nx,Nx, L,L); or m = createMesh3D(Nx,Nx,Nx, L,L,L); and see the outcome for yourself.

Examples

There are a few simple examples in the Tutorial folder. You can also find a few more advanced examples (water injection into a heterogeneous oil field, two nonlinear PDEs, coupled fully implicit solution) in the Advanced folder.

Documents

Find some preliminary documents here.

But Matlab is not a free software?

You can use the code in octave. The new (object oriented) version of the code works in Octave 4.0 (with the new classdef function).
I've re-written the code in Julia. It works fine, but the visualization on Windows OS has still some issues.

Questions and bug reports

You can ask your questions by creating a new issue here, or by writing a comment in my blog. You can also ask your question in the Matlab file exchange page of this code. I truly appreciate your feedback and/or contribution.

How to cite:

If you have used the package in your work and you find it usefull, please cite it as:

@misc{ali_akbar_eftekhari_2015_32745,
  author       = {Ali Akbar Eftekhari and Kai Schüller and Ferran Brosa Planella and Martinus Werts and Behzad Hosseinzadeh},
  title        = {FVTool: a finite volume toolbox for Matlab},
  month        = oct,
  year         = 2015,
  doi          = {10.5281/zenodo.32745},
  url          = {https://doi.org/10.5281/zenodo.32745}
}

I will also appreciate it if you write me a couple of lines about how you have used it in your research. It encourages me to maintain the code.

fvtool's Issues

A bug in 3D Cartesian TVD scheme

There is a bug in the 3D TVD scheme. I had never tested the function. Now it is solved.

Spherical coordinate

The derivations are ready to implement.

convection diffusion with time-space dependent coefficients

Hello!
I would like to solve a convection diffusion equation with a point source and time-and-space-dependent coefficients $D=D(x,t)$, u=(x,t), $x\in R^d$, $d=1,2,3.$

Is it possible to use ode45 or similar solver instead of doing loop in time if have time-dependent coefficients? Or would you recommend to use time-loop?

I tried to modify example diffusionODEtutorial but unfortunately it does not work.
The example convectionODEexample does not work either...

Best regards,
Anna

Boundary conditions

Dear Dr. Eftekhari,

Thank you for provided us with a flexible quite simple software. I just do not get how to use the boundary conditions or which kind of approach It uses.
In a 1D (so far i am interested in that), I see that the matrix will be:
| a1 a2 0 .....0 |
A =| .....O ........... | where a1 = -(BC.left.b/2 + BC.left.a/dx_1)
| 0 .....0 a3 a4 | a2 = -(BC.left.b/2 + BC.left.a/dx_1)
a3 = BC.right.b/2 + BC.right.a/dx_end
a4 = BC.right.b/2 - BC.right.a/dx_end

and the vector of boundary conditions:
|q1 | q1 = -(BC.left.c)
q = | 0 |
|q2 | q2 = BC.right.c

So I must implement the values of a, b, and c for each side:

If I want a dirichlet boundary such as c(0, t) =C0; (left side constant), how do it do it? (equal for the right side?)

and

If I want a flux boundary (cauchy boundary) c(0, t) = C0 + (D/v) dC(0, t)/dx ; (left side constant), how do it do it? (equal for the right side?)

and a no flux boundary such as v = 0 and dC(0, t)/dx =0;(left side constant), how do it do it? (equal for the right side?)

Thank a lot

Best Regards,

Daniel

Explicit integration transport equation 1D

Dear Dr. Eftekhari,
I have another question which is probably kind of silly, although it could be quite helpful to me.
If I want to solve a classical transport advection equation (without sources), I could end up with the following discrete system:

E*(C_n1 -C_n) / DT = L*( (1-thetha)*Cn_1+ (thetha)*C_n) + q

where E is a diagonal matrix with the values of the porosity
C_n1 is the vector of dependent variable at time n+1
C_n is the vector of dependent variable at time n
Dt is the time step
L is a (probably tridiagonal matrix) matrix related to the spatial discretization for advection and diffusion.
thetha is a factor from 1 to 0, when 1 forward euler and when 0 implicit euler.
q is a vector of boundary conditions

In the propose discretization the spatial discretization order is not clear but the temporal is obviously first order (although if thetha is 1/2 it will be second order). I would like to work with a time explicit scheme, namely forward euler and thetha 1. First, for a 1D problem

In the tutorial folder, most of the problems in 1D are solved by an implicit approach, except (I think) ' diffusiontutorialExplicit .m'. Where the solution is given by the line:

% step 1: calculate divergence term
RHS = divergenceTerm(Dave._gradientTerm(c_old));
% step 2: calculate the new value for internal cells
c_old_int=internalCells(c_old);
c_val_int = dt_excludeGhostRHS(meshstruct, RHS+constantSourceTerm(c_old))+c_old_int(:);

I have tried to add advection with the following modification:

% step 1: calculate divergence term
RHSd = divergenceTerm(Dave._gradientTerm(c_old));
RHSv = divergenceTerm(u._arithmeticMean(c_old));
% step 2: calculate the new value for internal cells
c_old_int=internalCells(c_old);
c_val_int = subdt*excludeGhostRHS(m, -RHSv+RHSd+constantSourceTerm(c_old))+c_old_int(:);

Although the results are not satisfactory, I wonder if it is because the time step is big, the implementation that I use is wrong or maybe something to do with boundary conditions?
Which ways there are to use a simple forward euler once the spatial discretization has been done with FVT tool?

Thanks again

Best,
Daniel

solving reaction-diffusion(-advection) equations?

Hello,
Many thanks for this great piece of code. I am just wondering if it could be used also for reaction diffusion problems? like for the Fisher-KPP equation u_t = Laplacian(phi(u)) + u(1-u) where phi(u)=u^m. While I can solve u_t = Laplacian(phi(u)) with FVTool, I don't know how to treat the nonlinear reaction term. How about systems of nonlinear PDEs?
Thanks a lot!
Cheers,
Jan

diffusion equation

Hello! I was desperate to find this tool.

It will be very helpful to me. Thank you!

I have to solve the PDE systems and the thing I want to know is that FVTool can solve this kind of
problem. (d means gradient and Dc means the diffusion coefficient(constant) of concentration C
B means the concentration of B)
-> dCdt = d/dx((Dc-C) * heaviside(Dc-C) * dC/dx) + k1 * B * C

That is, can FVTool solve not only "dCdt=d^2U/dx^2" but also
"dCdt=d/dx(A(C)*dC/dx)" ?

• A(C) is a certain function of concentration C.

Please share your wisdom with me.

Add mesh visualization functionality

Not really important, but may be helpful for preparing presentations and documents.

advection-diffusion-reaction

Hi, I need to solve the 1D advection diffusion equation for a species concentration, with a sink term. At the left BC there is a constant concentration for all time steps, at the right BC, the advective flux equals the diffusive flux initially. I want to add an optional sink term "k(phi-phi_max)", which should be zero for cells when phi<phi_max, and which should have a value that sets the cells concentration equal to phi_max when concentration reaches phi_max, not allowing the concentration to rise above phi_max. Sorry for the long description, but is this possible with your code?

Spherical coordinates not working

I tried running the diffusiontutorial_spherical code, but I commented out line 20, and uncommented line 19 (switching BCs from left to right), and the computation no longer performs as expected. As I am just making a symmetric change in the BCs, is there a reason it doesn't work?

Thanks!

1D heat transfer through the composite wall

Ozel - 2014 - Effect of insulation location on dynamic heat-tran.pdf

Optimization of insulation thickness using implicit finite difference method under steady periodic condition.

Implicit vs explicit Euler

Hi Ali,

I was wondering why, while using backward Euler in time, the convection term obtained by TVD usees the pervious iterate (phi_old)... Would it be more natural to use explicit Euler?

[See e.g. ConvectionTVDBurgers and diffconv_variable_diff in Examples]

Thanks,
Anna

Size o linear system of equations

I noticed that when specifying 2-D problem with NxNy nodes, the coefficient matrix for the linear system is NxN, with N=(Nx+2)(Ny+2). I would expect something like NxN with N=(Nx+2)*(Ny+2). I understand the two extra points are probably related by implementation of boundary conditions, but I don't quite understand why the coefficient matrix is as large as it is. Please let me know if I am missing out on something fundamental.

Non-uniform grid support

My Julia implementation supports non-uniform grids. I need one full day to add the same feature to FVTool.

Radial (R, Theta) and 3D cylindrical (R, Theta, Z) coordinates

Support for radial and 3D cylindrical coordinates. The visualization routines must be updated. A few example to clarify the boundary conditions for axi-symmetric coordinates.

single phase flow and transport

Hi, very interesting suite of codes! I am looking for an example where both single flow and the ADE are solved but I haven't been able to find any. Is there any example that I could look at? Thanks in advance!

Changing bounds of mesh

Hi

Is it possible to change the bounds of a mesh? For instance in:
m = createMeshCylindrical2D(Nr, Nz, Lr, Lz)

Setting the values for Lr and Lz cause the mesh to extend from "0 to Lr" and "0 to Lz". Is it possible to change this so that the mesh extends from, say, "x to Lr" and "y to Lz" ? Thanks in advance!

Diffusion equation

Hello,

First of all, thank you for this tool! It is being very helpful.

I wanted to ask if there is any way I can input a tensor in D in the diffusion term, and if so how to do it.

Thanks!

Convenience functions for spatially variable cell and face variables

John pointed out this, and I added two new functions to return the coordinate of the cell centers and face centers. An script that describes these functions is here: https://github.com/simulkade/FVTool/blob/master/Examples/Tutorial/variables_in_space.m

What's the difference between Cylindrical2D and Radial2D?

Hello,
Thank you for sharing FVTool, I've benefited a lot from it.

After reading some of the codes, I was a little confused by 2-D Cylindrical coordinate system. I know 3-D Cylindrical coordinate system(r, θ, z), but don't know what 2-D Cylindrical means.

If the 2-D Cylindrical coordinate system is based on the r-z plane, what is the difference between it and the 2-D Cartesian coordinate system(x-y)?
If the 2-D Cylindrical coordinate system is based on the r-θ plane, what is the difference between it and the 2-D Radial coordinate system(r-θ)?

In addition, can this tool solve equations in 3-D Spherical grids?

I'll appreciate it if you can make me clear about these questions.
Thanks!