# Choosing the right cell averaging method for linearized PDEs

## Choosing the right averaging method¶

In the previous post, I discussed a linearization procedure that could be used to solve a nonlinear diffusion equation. We saw that the linearization converted the nonlinear diffusion equation to a PDE with linear diffusion and convection terms. Here, I'm going to solve that equation again, but this time on a nonuniform mesh. My objective is to see what happens if we choose an averaging term that is not consistent with the chosen method for the discretization of the convection term.

### Another look at the equations¶

The equation that we solved was $$\frac{\partial\phi}{\partial t}+\nabla . (-D(\phi)\nabla \phi) =0$$ where we can linearize the diffusion term as $$D\left(\phi\right)\nabla \phi=D\left(\phi_{0}\right)\nabla \phi+\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0}}\left(\phi-\phi_{0}\right).$$ Let's put it back in the original diffusion equation: $$\frac{\partial\phi}{\partial t}+\nabla . \left(-D\left(\phi_{0}\right)\nabla \phi-\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0}}\left(\phi-\phi_{0}\right)\right) =0$$ Of course this is an approximation. The idea is that we gues a $\phi_0$, then solve the equation for $\phi$, and replace it in $\phi_0$ until the absolute or relative difference between them is less than a chosen small value. I write the above equation in a form that is actually used in a finite volume scheme with all the averaging on the cell faces. It reads $$\frac{\partial\phi}{\partial t}+\nabla . \left(-D\left(\phi_{0}\right)\nabla \phi\right)+\nabla.\left(-\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0,face}}\phi\right) =\nabla.\left(-\nabla \phi_{0}\left(\frac{\partial D}{\partial \phi}\right)_{\phi_{0,face}}\phi_{0,face}\right)$$ Here's the trick: on the left hand side, we have, from left to right, a transient term, a diffusion term, a convection term, and on the right hand side, a divergence term. We can use any type of discretization for the convective term, as long as the same averaging term to calculate the face value on the right hand side, that would be:

Discretization Averaging
convectionTerm linearMean
convectionUpwindTerm upwindMean
convectionTvdTerm tvdMean

In this example, I'm going to show what will happen if you use an averaging technique that is not consistent with your chosen discretization method for the convective term.

### Implementation¶

I'm repeating the nonlinear diffusion model of the previous post, this time on a nonuniform grid. Two functions are created: the direct substitution, which is assumed to give the right solution, and a linearized equation, with the possibility of using central difference or upwind convection discretization and linear, upwind, and arithmetic averaging methods.

### Definition of functions in JFVM¶

In [1]:
using JFVM, PyPlot

INFO: Loading help data...


### Linearized PDE (Newton method)¶

In [2]:
function diffusion_newton(conv_method, avg_method)
L= 1.0 # domain length
Nx= 100
dx_min=L/Nx
x=[0.0, dx_min]
while x[end]<L
push!(x, x[end]+1.05(x[end]-x[end-1]))
end
x[end]=L
m= createMesh1D(x) # create a nonuniform 1D mesh
D0= 1.0 # diffusion coefficient constant
# Define the diffusion coefficientand its derivative
D(phi)=D0*(1.0+phi.^2)
dD(phi)=2.0*D0*phi
# create boundary condition
BC = createBC(m)
BC.left.a[:]=0.0
BC.left.b[:]=1.0
BC.left.c[:]=5.0
BC.right.a[:]=0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0
(Mbc, RHSbc)= boundaryConditionTerm(BC)
# define the initial condition
phi0= 0.0 # initial value
phi_old= createCellVariable(m, phi0, BC) # create initial cells
phi_val= copyCell(phi_old)
# initialize the solution for upwind case
phi_face= linearMean(phi_val)
# define the time steps
dt= 0.01*L*L/D0 # a proper time step for diffusion process
for i=1:10
err=100
(Mt, RHSt) = transientTerm(phi_old, dt)
while err>1e-10
# calculate the diffusion coefficient
Dface= harmonicMean(cellEval(D,phi_val))
# calculate the face value of phi_0
if avg_method==1
phi_face= linearMean(phi_val)
elseif avg_method==2
phi_face= arithmeticMean(phi_val)
elseif avg_method==3
phi_face= upwindMean(phi_val, u)
end
# calculate the velocity for convection term
# diffusion term
Mdif= diffusionTerm(Dface)
# convection term
if conv_method==1
Mconv= convectionTerm(u)
elseif conv_method==2
Mconv= convectionUpwindTerm(u)
end
# divergence term on the RHS
RHSlin= divergenceTerm(u.*phi_face)
# matrix of coefficients
M= Mt-Mdif+Mconv+Mbc
# RHS vector
RHS= RHSbc+RHSt+RHSlin
# call the linear solver
phi_new= solveLinearPDE(m, M, RHS)
# calculate the error
err= maximum(abs(phi_new.value-phi_val.value))
# assign the new phi to the old phi
phi_val=copyCell(phi_new)
end
phi_old= copyCell(phi_val)
#visualizeCells(phi_val)
end
return phi_val
end

Out[2]:
diffusion_newton (generic function with 1 method)

### Direct substitution¶

In [3]:
function diffusion_direct()
L= 1.0 # domain length
Nx= 100
dx_min=L/Nx
x=[0.0, dx_min]
while x[end]<L
push!(x, x[end]+1.05(x[end]-x[end-1]))
end
x[end]=L
m= createMesh1D(x) # create a nonuniform 1D mesh
D0= 1.0 # diffusion coefficient constant
# Define the diffusion coefficientand its derivative
D(phi)=D0*(1.0+phi.^2)
dD(phi)=2.0*D0*phi
# create boundary condition
BC = createBC(m)
BC.left.a[:]=0.0
BC.left.b[:]=1.0
BC.left.c[:]=5.0
BC.right.a[:]=0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0
(Mbc, RHSbc)= boundaryConditionTerm(BC)
# define the initial condition
phi0= 0.0 # initial value
phi_old= createCellVariable(m, phi0, BC) # create initial cells
phi_val= copyCell(phi_old)
# define the time steps
dt= 0.01*L*L/D0 # a proper time step for diffusion process
for i=1:10
err=100
(Mt, RHSt) = transientTerm(phi_old, dt)
while err>1e-10
# calculate the diffusion coefficient
Dface= harmonicMean(cellEval(D,phi_val))
# diffusion term
Mdif= diffusionTerm(Dface)
# matrix of coefficients
M= Mt-Mdif+Mbc
# RHS vector
RHS= RHSbc+RHSt
# call the linear solver
phi_new= solveLinearPDE(m, M, RHS)
# calculate the error
err= maximum(abs(phi_new.value-phi_val.value))
# assign the new phi to the old phi
phi_val=copyCell(phi_new)
end
phi_old= copyCell(phi_val)
#visualizeCells(phi_val)
end
return phi_val
end

Out[3]:
diffusion_direct (generic function with 1 method)

### Case 1: upwind convection with arithmetic average¶

In [8]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(2,2)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r")
legend(["direct subs", "Newton"])
title("Upwind/Arithmetic")
xlabel("x")
ylabel(L"\phi", fontsize=14);


As you can see, the converged solution of the Newton method, which uses inconsistent convection/averaging terms does not match the solution of the direct substitution solver.

### Case 2: central convection with arithmetic average¶

In [9]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(1,2)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r")
legend(["direct subs", "Newton"])
title("Central/arithmetic")
xlabel("x")
ylabel(L"\phi", fontsize=14);


Again, choosing the inconsistent methods results in an erroneous solution (for the Newton method).

### Case 3: central convection with linear average¶

In [10]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(1,1)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r")
legend(["direct subs", "Newton"])
title("Central/linear")
xlabel("x")
ylabel(L"\phi", fontsize=14);


Beautiful! The results match perfectly.

### Case 4: upwind convection with upwind average¶

In [11]:
phi_ds=diffusion_direct()
phi_n=diffusion_newton(2,3)
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r")
legend(["direct subs", "Newton"])
title("Upwind/upwind")
xlabel("x")
ylabel(L"\phi", fontsize=14);


Again, we get the correct result by choosing a right pair of discretization/averaging methods. Depending on the Peclet number (convective to diffusive rate ratio), using the upwind method gives a numerically more stable solution.

# Solving nonlinear PDEs with FVM

There are many cases where the transfer coefficient in the convection-diffusion equations is a function of the independent variable. This makes the PDE nonlinear. While working on my thesis, I realized that it is not that difficult to use the Taylor expansion to linearize the nonlinear terms in a PDE. Here, I'm going to give a simple example. In another post, I will explain how the same idea can be used to numerically solve the Buckley-Leverett equation.

### The PDE¶

I'm going to try a simple one here. I say try because I'm actually making it up right now. The equation reads $$\frac{\partial\phi}{\partial t}+\nabla . (-D\nabla \phi) =0$$ where $$D=D_0(1+\phi^2).$$ It means that the diffusion coefficient is increased when the concentration of the solute increases.
To linearize this equation, I assume that $\phi$ and $\nabla \phi$ are independent variables. Using Taylor expansion around a point $\phi_0$, we can obtain the following general form: $$f\left(\phi\right)\nabla \phi=f\left(\phi_{0}\right)\nabla \phi+\nabla \phi_{0}\left(\frac{\partial f}{\partial \phi}\right)_{\phi_{0}}\left(\phi-\phi_{0}\right).$$ In our case, this equation simplifies to $$D_0(1+\phi^2)\nabla \phi=D_0(1+\phi_0^2)\nabla \phi+2D_0\phi_0\nabla \phi_0(\phi-\phi_0).$$ Putting back everythin together, we obtain this linearized form for our PDE $$\frac{\partial\phi}{\partial t}+\nabla.(-D_{0}(1+\phi_{0}^{2})\nabla\phi)+\nabla.(-2D_{0}\phi_{0}\nabla\phi_{0}\phi)=\nabla.(-2D_{0}\phi_{0}\nabla\phi_{0}\phi_{0}).$$ Don't freak out! This equation is actually quite simple. By linearizing, we have added a linear convection term to our nonlinear diffusion equation. This equation is still an approximation of the real PDE. We have to solve the linear equation for $\phi$ by initializing $\phi_0$. Then, we assign the new value of $\phi$ to $\phi_0$ until it converges to a solution. The question that still remains is that "is this really a good idea to make a diffusion equation more complicated by adding a convection term to it?".

### The solution procedure¶

Solving this equation is done in two loops. The outer loop is for the time steps, and the inner loop is for the nonlinear equation. We will assemple an equation which looks like this: $$\text{Transient}(\Delta t,\phi)-\text{Diffusion}(D(\phi_0),\phi)+\text{Convection}(\mathbf{u}(\phi_0),\phi)=\text{divergence}(D(\phi_0), \nabla\phi_0)$$

### Implementation in JFVM¶

Here, I'll implement the procedure in JFVM. The same procedure can be done in FVtool.

In [1]:
using JFVM, PyPlot

INFO: Loading help data...

In [16]:
function diffusion_newton()
L= 1.0 # domain length
Nx= 100 # number of cells
m= createMesh1D(Nx, L) # create a 1D mesh
D0= 1.0 # diffusion coefficient constant
# Define the diffusion coefficientand its derivative
D(phi)=D0*(1.0+phi.^2)
dD(phi)=2.0*D0*phi
# create boundary condition
BC = createBC(m)
BC.left.a[:]=0.0
BC.left.b[:]=1.0
BC.left.c[:]=5.0
BC.right.a[:]=0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0
(Mbc, RHSbc)= boundaryConditionTerm(BC)
# define the initial condition
phi0= 0.0 # initial value
phi_old= createCellVariable(m, phi0, BC) # create initial cells
phi_val= copyCell(phi_old)
# define the time steps
dt= 0.001*L*L/D0 # a proper time step for diffusion process
for i=1:10
err=100
(Mt, RHSt) = transientTerm(phi_old, dt)
while err>1e-10
# calculate the diffusion coefficient
Dface= harmonicMean(cellEval(D,phi_val))
# calculate the face value of phi_0
phi_face= linearMean(phi_val)
# calculate the velocity for convection term
# diffusion term
Mdif= diffusionTerm(Dface)
# convection term
Mconv= convectionTerm(u)
# divergence term on the RHS
RHSlin= divergenceTerm(u.*phi_face)
# matrix of coefficients
M= Mt-Mdif-Mconv+Mbc
# RHS vector
RHS= RHSbc+RHSt-RHSlin
# call the linear solver
phi_new= solveLinearPDE(m, M, RHS)
# calculate the error
err= maximum(abs(phi_new.value-phi_val.value))
# assign the new phi to the old phi
phi_val=copyCell(phi_new)
end
phi_old= copyCell(phi_val)
visualizeCells(phi_val)
end
return phi_val
end

Out[16]:
diffusion_newton (generic function with 1 method)
In [33]:
@time phi_n=diffusion_newton()

elapsed time: 0.065438686 seconds (15208324 bytes allocated)

Out[33]:
CellValue{Float64}(MeshStructure(1,[100],CellSize{Float64}([0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01  …  0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01],[0.0],[0.0]),CellLocation{Float64}([0.005,0.015,0.025,0.035,0.045,0.055,0.065,0.075,0.085,0.095  …  0.905,0.915,0.925,0.935,0.945,0.955,0.965,0.975,0.985,0.995],[0.0],[0.0]),FaceLocation{Float64}([0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09  …  0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1.0],[0.0],[0.0]),[1],[1]),[5.0196,4.9804,4.94059,4.90018,4.85914,4.81745,4.77511,4.7321,4.68839,4.64398  …  0.000323843,0.000250452,0.000192389,0.000146146,0.000108911,7.84098e-5,5.27679e-5,3.04028e-5,9.92833e-6,-9.92833e-6])

The solution looks quite good. The diffusion coefficient is higher at higher concentrations, and it can be clearly observed in the above concentration profile (if you compare it with a diffusion problem with constant diffusion coefficient, in your head).

### Some questions¶

The convection term that appears after linearization of the problem has a velocity term, which reads $$\mathbf{u}(\phi_0)=−2.0D_0\phi_0\nabla\phi_0$$ I have written a function, called gradientTerm which gives you the value of the gradient on the cell faces. This gradient is multiplied by a term that needs to be averaged on each cell face. Here, I have calculated the arithmetic average of $\phi_0$ using a function called suprise suprise linearMean. This is of course one way of doing it. The other ways are geometricMean, harmonicMean, upwindMean, tvdMean, and aritheticMean. Here, we will discuss them.

### Averaging techniques for the face values in FVM¶

The following figure shows two neigbor cells with different length and values.
The averaging terms that are used in JFVM are as follows:

#### Linear mean¶

It is equivalent to a linear interpolation or the lever rule. It reads $$\phi_{face}=\frac{\Delta x_1\phi_2+\Delta x_2\phi_1}{\Delta x_1+\Delta x_2}$$

#### Arithmetic mean¶

The value of each cell is weighted by its length: $$\phi_{face}=\frac{\Delta x_1\phi_1+\Delta x_2\phi_2}{\Delta x_1+\Delta x_2}$$

#### Geometric mean¶

It reads $$\phi_{face}=\exp\left(\frac{\Delta x_1\log(\phi_1)+\Delta x_2\log(\phi_2)}{\Delta x_1+\Delta x_2}\right)$$

#### Harmonic mean¶

It reads $$\phi_{face}=\frac{\Delta x_1+\Delta x_2}{\Delta x_1/\phi_1+\Delta x_2/\phi_2}$$ Note that if the value of one of the cells is zero, the geometric and harmonic means are also equal to zero.
These methods have different applications for different ptoblems. We will discuss them in future posts.

### Direct substitution¶

The diffusion problem that we solved above can also be solved using a direct substitution method. Here, we don't use the derivatives and no linearization is required. The code is written here:

In [18]:
function diffusion_direct()
L= 1.0 # domain length
Nx= 100 # number of cells
m= createMesh1D(Nx, L) # create a 1D mesh
D0= 1.0 # diffusion coefficient constant
# Define the diffusion coefficientand its derivative
D(phi)=D0*(1.0+phi.^2)
dD(phi)=2.0*D0*phi
# create boundary condition
BC = createBC(m)
BC.left.a[:]=0.0
BC.left.b[:]=1.0
BC.left.c[:]=5.0
BC.right.a[:]=0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0
(Mbc, RHSbc)= boundaryConditionTerm(BC)
# define the initial condition
phi0= 0.0 # initial value
phi_old= createCellVariable(m, phi0, BC) # create initial cells
phi_val= copyCell(phi_old)
# define the time steps
dt= 0.001*L*L/D0 # a proper time step for diffusion process
for i=1:10
err=100
(Mt, RHSt) = transientTerm(phi_old, dt)
while err>1e-10
# calculate the diffusion coefficient
Dface= harmonicMean(cellEval(D,phi_val))
# diffusion term
Mdif= diffusionTerm(Dface)
# matrix of coefficients
M= Mt-Mdif+Mbc
# RHS vector
RHS= RHSbc+RHSt
# call the linear solver
phi_new= solveLinearPDE(m, M, RHS)
# calculate the error
err= maximum(abs(phi_new.value-phi_val.value))
# assign the new phi to the old phi
phi_val=copyCell(phi_new)
end
phi_old= copyCell(phi_val)
visualizeCells(phi_val)
end
return phi_val
end

Out[18]:
diffusion_direct (generic function with 1 method)
In [34]:
@time phi_ds=diffusion_direct()

elapsed time: 0.115909663 seconds (22458788 bytes allocated, 26.14% gc time)

Out[34]:
CellValue{Float64}(MeshStructure(1,[100],CellSize{Float64}([0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01  …  0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01],[0.0],[0.0]),CellLocation{Float64}([0.005,0.015,0.025,0.035,0.045,0.055,0.065,0.075,0.085,0.095  …  0.905,0.915,0.925,0.935,0.945,0.955,0.965,0.975,0.985,0.995],[0.0],[0.0]),FaceLocation{Float64}([0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09  …  0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1.0],[0.0],[0.0]),[1],[1]),[5.0196,4.9804,4.94059,4.90018,4.85914,4.81745,4.77511,4.7321,4.68839,4.64398  …  0.000323843,0.000250452,0.000192389,0.000146146,0.000108911,7.84098e-5,5.27679e-5,3.04028e-5,9.92833e-6,-9.92833e-6])

### Show and compare the results¶

Here, I'm plotting the results of the Newton and the direct substitution methods for comparison.

In [40]:
plot(phi_ds.domain.cellcenters.x, phi_ds.value[2:end-1], "y",
phi_n.domain.cellcenters.x, phi_n.value[2:end-1], "--r")
legend(["direct subs", "Newton"])

Out[40]:
PyObject <matplotlib.legend.Legend object at 0x7ffa453f94d0>

### Conclusion¶

The final results shows that we can linearize a nonlinear PDE, relatively easily and solve the equations using a linear PDE solver. No difference is observed between the final results of the direct substitution method (with no linearization) and the Newton method (with a linearized PDE). We can also observe that, as expected, the Newton method is relatively faster.

## Update¶

The same problem implemented in Matlab:

# Switched to Nikola

Today, I switched this blog from Octopress to Nikola. I did it mostly because I'm more comfortable with Python, and of course the new exciting version of ipython, which is now called Jupyter. I think I'm going to stick to this format for a while.