Skip to main content

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)
    u= -faceEval(dD, phi_face).*gradientTerm(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
            u= -faceEval(dD, phi_face).*gradientTerm(phi_val)
            # 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
        u= faceEval(dD, phi_face).*gradientTerm(phi_val)
        # 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.
neighbor cells FVM 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: