Diffusion equation 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.$

Analytical solution

The anaytical solution of the diffusion problem in a semi-infinite domain is given by $$c=1-\textrm{erf}(x/2\sqrt{Dt}).$$ We can choose a long domain to mimic a semi-infinite domain.

Start the solver



In [16]:

    
using JFVM
using PyPlot

Domain size and physical constants



In [23]:

    
D_val = 1.0 # [m^2/s]
Lx = 50.0 # [m]
c_init = 0.0
c_left = 1.0
c_right = 0.0









    Out[23]:





0.0

A diffusion function

First, we define a function that accepts a mesh, boundary condition, diffusion coefficient, and final simulation time, and returns a cell value, which is the final concentration profile. I also define a boundary condition function to assign values of one and zero to the left and right boundaries.



In [8]:

    
function trans_diff(m::MeshStructure, BC::BoundaryCondition, D_val::Real, t_final::Real)
    # a transient diffusion 
    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= t_final/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:N_steps
        (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

function dirich_bc(m::MeshStructure)
    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
    return BC
end









    Out[8]:





dirich_bc (generic function with 1 method)

Solution

I define two domains, one with uniform grids and another with a nonuniform grids with smaller cells near the left boundary. The total number of grids are the same for both case. The figure below shows how you can use grid refinement near the more important areas of the domain to get the best out of your numerical model.



In [34]:

    
t_final = 5
Nx = 11 # number of cells
m_uni = createMesh1D(Nx, Lx)
BC_uni = dirich_bc(m_uni)
c_uni = trans_diff(m_uni, BC_uni, D_val, t_final)
x = linspace(0,50,200)
m_nonuni = createMesh1D([linspace(0,10,5), linspace(14,50,5)])
BC_nonuni = dirich_bc(m_nonuni)
c_nonuni = trans_diff(m_nonuni, BC_nonuni, D_val, t_final)
c_analytic = 1-erf(x./(2.*sqrt(D_val*t_final)))
visualizeCells(c_uni)
visualizeCells(c_nonuni)
plot(x, c_analytic)
legend(["uniform", "nonuniform", "analytic"])
axis([0, 10, 0, 1])









    












    Out[34]:





4-element Array{Int64,1}:
  0
 10
  0
  1