These are some notes that may be helpful to anyone who wants to play with this tool. (FNYAQ = Frequently Not-Yet-Asked Questions)
This is the general form of the equation that can be solved using JFVM code: $$ \alpha \frac {\partial \phi} {\partial t} + \nabla.(\mathbf{u} \phi) + \nabla.(-D \nabla \phi) + \beta \phi = \phi_0$$ The boundary condition is written in the general form of $$a \nabla \phi . \mathbf{n} + b \phi = c$$ Simple and compact as the above equations seem, they represent a diverse number of physical systems. Yes. you can guess that I'm a chemical engineer.
This code is not a exactly a PDE solver. It is a discretization tool, which uses finite volume method, to convert each term in the mentioned equations to a linear system of equations: $$M_t(\phi-\phi_{old})+M_{conv} \phi-M_{diff}\phi+M_{ls}\phi+M_{bc}\phi=RHS_{0}+RHS_{bc}$$ This allows the user to easily handle any equation with transient-convection-diffusion-source terms or even a system of coupled PDE's.
Because it's mass conservative, and I like it (remember the ChemEng thing).
The main purpose is to have a toy tool that allows the user to play with 1D models. Then if the mentioned user is happy with 1D results, easily switch it to 2D and 3D or switch to other coordinate systems:
The discretization schemes include
FiPy. I've learned a lot from using it. It's a fantastic Python package, which was a bit too big for what I wanted to do. It has still some issues with boundary conditions in cylindrical coordinate. I use more or less the same names that FiPy uses in describing FVM terms.
As a first example, I try to solve the diffusion equation, i.e., $$ \frac {\partial c} {\partial t} + \nabla.(-D \nabla c)=0$$ The boundary condition is written in the general form of $$a \nabla \phi . \mathbf{n} + b \phi = c$$ We use Dirichlet boundary condition for the left and right boundaries, i.e., $c=1.0$ at the left boundary and $c=1.0$ at the right boundary. The initial condition is $c=0.0$ at $t=0.$
CellValue of the initial conditionCellValue of the diffusion coefficient [yes, can be hetrogeneous]averaging scheme to calculate the average diffusion coefficient on the cell faces; assign it to a FaceValueCellValue in the new time step
In [1]:
using JFVM
It loads all the functions, and also loads PyPlot and PyCall, which is used to load mayavi.mlab as a variable called mayavis (short for mayavi visualization tools). one and two dimensional visualizations are done in matplotlib and 3D visualization is done in mayavi, so don't forget to install both packages. In ubuntu, you can simply type
sudo apt-get install python-matplotlib mayavi2before using the code.
The created mesh can be uniform (same cell sizes) or nonuniform. The uniform mesh is defined based on the size of the domain and the number of cells in each direction, and the nonuniform mesh is defined based on the position of cell faces in each direction. For example, you can create a uniform 1D mesh with 10 cells on a 5 m domain by writing:
In [2]:
Nx = 10 # number of cells in x direction
Lx = 5.0 # [m] domain length in x direction
m_uni = createMesh1D(Nx, Lx)
Out[2]:
You can generate a nonuniform mesh by explicitly specifying the location of each cell face in the x direction. Let's see the location of cell faces in the uniform mesh we created above:
In [3]:
m_uni.facecenters.x
Out[3]:
Now I can define a different cell face vector to define a new nonuniform mesh:
In [4]:
x = [0.0, 0.1, 0.3, 0.35, 1.0, 3.0, 3.4, 4.5, 5.0]
m_nonuni = createMesh1D(x)
Out[4]:
Other commands for creating mesh structures are:
In [5]:
Nx = 3 # number of cells in x-direction
Ny = 4 # number of cells in y-direction
Nz = 5 # number of cells in z-direction
Nr = 6 # number of cells in r-direction [radial]
Ntheta = 7 # number of cells in \theta-direction [2d radial, 3d cylindrical and spherical]
Lx = 1.0 # domain length in x-direction
Ly = 2.0 # domain length in y-direction
Lz = 3.0 # domain length in z-direction
Lr = 4.0 # domain length in r-direction
Ltheta = 2*pi # domain length in \theta-direction [in radian]
x=[0.0 , 0.4, 0.6, 1.0] # cell face positions in the x-direction
y=[0.0 , 0.4, 0.6, 1.0, 2.0] # cell face positions in the y-direction
z=[0.0 , 0.4, 0.6, 1.0, 2.0, 3.0] # cell face positions in the z-direction
r=[0.0 , 0.4, 0.6, 1.0, 2.0, 3.0, 4.0] # cell face positions in the r-direction
theta=5.0./[0.0 , 0.4, 0.6, 1.0, 2.0, 3.0, 4.0, 5.0]*Ltheta # cell face positions in the \theta-direction
# 1D cartesian uniform mesh:
m_1d_cartesian_uniform = createMesh1D(Nx, Lx)
# 1D cartesian nonuniform mesh
m_1d_cartesian_nonuniform = createMesh1D(x)
# 1D radial uniform mesh
m_1d_radial_uniform = createMeshCylindrical1D(Nr, Lr)
# 1D radial nonuniform mesh
m_1d_radial_nonuniform = createMeshCylindrical1D(r)
# 2D cartesian uniform mesh
m_2d_cartesian_uniform = createMesh2D(Nx, Ny, Lx, Ly)
# 2D cartesian nonuniform mesh
m_2d_cartesian_nonuniform = createMesh2D(x, y)
# 2D cylindrical uniform mesh
m_2d_cylindrical_uniform = createMeshCylindrical2D(Nr, Nz, Lr, Lz)
# 2D cylindrical nonuniform mesh
m_2d_cylindrical_nonuniform = createMeshCylindrical2D(r, z)
# 2D radial uniform mesh
m_2d_radial_uniform = createMeshRadial2D(Nr, Ntheta, Lr, Ltheta)
# 2D radial nonuniform mesh
m_2d_radial_nonuniform = createMeshRadial2D(r, theta)
# 3D cartesian uniform mesh
m_3d_cartesian_uniform = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
# 3D cartesian nonuniform mesh
m_3d_cartesian_nonuniform = createMesh3D(x, y, z)
# 3D cylindrical uniform mesh
m_3d_cylindrical_uniform = createMeshCylindrical3D(Nr, Ntheta, Nz, Lr, Ltheta, Lz)
# 3D cylindrical nonuniform mesh
m_3d_cylindrical_nonuniform = createMeshCylindrical3D(r, theta, z)
println("Various mesh types in JFVM!")
The mesh structure can be created using the function createBC. The default boundary conditions are Neumann, i.e., $$\nabla \phi .\mathbf{n} = 0.$$ However, we can easily change it to Dirichlet or Robin (oh, Robin! You are everywhere) by altering the values of $a$, $b$, and $c$ in the following form of boundary condition: $$a \nabla \phi .\mathbf{n} + b \phi= c.$$ You can change these values on every cell face on each boundary separately. I'll demonstrate it in a moment. Right now, we create a boundary condition structure and change the left and right boundaries to Dirichlet (fixed values):
In [6]:
Nx = 10
Lx = 1.0
m = createMesh1D(Nx, Lx)
BC = createBC(m)
BC.left.a[:]=BC.right.a[:]=0.0
BC.left.b[:]=BC.right.b[:]=1.0
BC.left.c[:]=1.0
BC.right.c[:]=0.0
Out[6]:
I use ghost cells to handle boundary conditions. Therefore, the number of cells in each direction has two more values for the ghost cells. The value of ghost cells depends on the boundary conditions. When we assign the initial condition to the cell values, we need to provide the boundary condition structure to calculate the values of the ghost cells. For other cell values, e.g., transfer coefficients, no boundary condition is required but the ghost cells are still initialized to make some computation more convenient (for instance averaging on the cell faces). The cell values are created by:
In [7]:
c_init = 0.0 # initial value of the variable
c_old = createCellVariable(m, 0.0, BC)
c_old.value
Out[7]:
As you can see, the value of the ghost cell near the left boundary is 2.0, which makes sense. The transfer coefficient, here the diffusion coefficient, is also assigned to each cell (no need to say that it can be heterogeneous, i.e., different values on each cell). However, the transfer coefficients must be averaged on the cell faces. I have implemented a few averaging schemes in the code, that you can see here:
In [8]:
D_val = 1.0 # value of the diffusion coefficient
D_cell = createCellVariable(m, D_val) # assigned to cells
# Harmonic average
D_face = harmonicMean(D_cell)
# Geometric average
D_face = geometricMean(D_cell)
# Arithmetic average
D_face = arithmeticMean(D_cell)
# UpwindMean, tvdMean, and linearMean can also be used for other purposes, e.g., linearizing a PDE
Out[8]:
This is not really related to this code, but let's write a couple of lines about it. What can help you in choosing an effective time step is the length of the domain and the diffusion coefficient. Here, we can choose the time step to be $$ \Delta t = \frac{1}{N_{steps}} \sqrt{\frac{L^2}{D}}$$ Do the dimensional analysis for yourself.
In [9]:
N_steps = 20 # number of time steps
dt= sqrt(Lx^2/D_val)/N_steps # time step
Out[9]:
This part is actually quite short, because the main work is done. The PDE is discretized term by term. I learned this method from FiPy. You can do it as follows:
In [10]:
M_diff = diffusionTerm(D_face) # matrix of coefficient for diffusion term
(M_bc, RHS_bc)=boundaryConditionTerm(BC) # matrix of coefficient and RHS for the BC
for i =1:5
(M_t, RHS_t)=transientTerm(c_old, dt, 1.0)
M=M_t-M_diff+M_bc # add all the [sparse] matrices of coefficient
RHS=RHS_bc+RHS_t # add all the RHS's together
c_old = solveLinearPDE(m, M, RHS) # solve the PDE for the first time step, and replace the old value
end
In [11]:
Nx = 10
Lx = 1.0
m = createMesh1D(Nx, Lx)
BC = createBC(m)
BC.left.a[:]=BC.right.a[:]=0.0
BC.left.b[:]=BC.right.b[:]=1.0
BC.left.c[:]=1.0
BC.right.c[:]=0.0
c_init = 0.0 # initial value of the variable
c_old = createCellVariable(m, 0.0, BC)
D_val = 1.0 # value of the diffusion coefficient
D_cell = createCellVariable(m, D_val) # assigned to cells
# Harmonic average
D_face = harmonicMean(D_cell)
N_steps = 20 # number of time steps
dt= sqrt(Lx^2/D_val)/N_steps # time step
M_diff = diffusionTerm(D_face) # matrix of coefficient for diffusion term
(M_bc, RHS_bc)=boundaryConditionTerm(BC) # matrix of coefficient and RHS for the BC
for i =1:5
(M_t, RHS_t)=transientTerm(c_old, dt, 1.0)
M=M_t-M_diff+M_bc # add all the [sparse] matrices of coefficient
RHS=RHS_bc+RHS_t # add all the RHS's together
c_old = solveLinearPDE(m, M, RHS) # solve the PDE
visualizeCells(c_old)
end
In [12]:
function trans_diff_dirichlet(m::MeshStructure)
# a transient diffusion
BC = createBC(m)
BC.left.a[:]=BC.right.a[:]=0.0
BC.left.b[:]=BC.right.b[:]=1.0
BC.left.c[:]=1.0
BC.right.c[:]=0.0
c_init = 0.0 # initial value of the variable
c_old = createCellVariable(m, 0.0, BC)
D_val = 1.0 # value of the diffusion coefficient
D_cell = createCellVariable(m, D_val) # assigned to cells
# Harmonic average
D_face = harmonicMean(D_cell)
N_steps = 20 # number of time steps
dt= sqrt(m.facecenters.x[end]^2/D_val)/N_steps # time step
M_diff = diffusionTerm(D_face) # matrix of coefficient for diffusion term
(M_bc, RHS_bc)=boundaryConditionTerm(BC) # matrix of coefficient and RHS for the BC
for i =1:5
(M_t, RHS_t)=transientTerm(c_old, dt, 1.0)
M=M_t-M_diff+M_bc # add all the [sparse] matrices of coefficient
RHS=RHS_bc+RHS_t # add all the RHS's together
c_old = solveLinearPDE(m, M, RHS) # solve the PDE
end
return c_old
end
Out[12]:
The function trans_dif_dirichlet accept a mesh structure and solves a transient diffusion equation with dirichlet boundary conditions on the left and right boundaries. Let us call it by sending a 2D cartesian mesh as the input argument:
In [13]:
Nx = 10
Ny = 15
Lx = 1.0
Ly = 2.0
m = createMesh2D(Nx, Ny, Lx, Ly)
c = trans_diff_dirichlet(m)
visualizeCells(c)
Out[13]:
See! As simple as that. Let's got to 3D!
In [14]:
using PyPlot
Nx = 10
Ny = 15
Nz = 20
Lx = 1.0
Ly = 2.0
Lz = 3.0
m = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
c = trans_diff_dirichlet(m)
img = visualizeCells(c)
imshow(img)
Out[14]:
There is indeed one way to capture a screenshot of the mayavi 3d plot and include it in the IJulia notebook. The visualizeCell returns an screenshot of the 3D plot, which can be visualized using the imshow function of PyPlot. Don't forget to load PyPlot first.
One other feature that I like to show here is the boundary condition for radial and cylindrical coordinates. The left boundary in the r direction is in fact axial-symmetry and can not be Dirichlet. We have to shift the coordinate a little bit to the right to be able to have a Dirichlet BC on the left side of a radial axis:
In [15]:
Nr=10
Ntheta = 15
Lr = 1.0
Ltheta = 2*pi
m = createMeshRadial2D(Nr, Ntheta, Lr, Ltheta)
m.facecenters.x+=0.1
m.cellcenters.x+=0.1
c = trans_diff_dirichlet(m)
visualizeCells(c)
Out[15]:
And this diffusion equation would not be complete without a 3D cylinder:
In [17]:
Nr=10
Ntheta = 15
Nz = 12
Lr = 1.0
Ltheta = 2*pi
Lz = 3.0
m = createMeshCylindrical3D(Nr, Ntheta, Nz, Lr, Ltheta, Lz)
m.facecenters.x+=0.1
m.cellcenters.x+=0.1
c = trans_diff_dirichlet(m)
img=visualizeCells(c)
imshow(img)
Out[17]: