I wrote a code which solves the 2D Poisson equation with homogeneous Dirichlet BC everywhere and a source term of -1. I am using the classical Jacobi iteration method. The grid is $N_x \, \mathrm{x} \, N_y$, and all the arrays are 0 based.

The code yields good results with a 5 points stencil (2$^{\mathrm{nd}}$ order finite difference method). With this stencil, the code loops all interior points (i.e. from $1$ to $N_x-2$ in the $x$ direction), while the boundary points are never touched.

However, the code diverges when using a 9 points stentil, 4$^\mathrm{th}$ order central difference method. Again, only the interior points are iterated on. Since the discretization uses two points on each side, I am extrapolating past the boundary points (ghost points) with a first order approximation, i.e. $U(-1) = U(0), U(Nx) = U(Nx-1)$.

Do I have to use special discretization methods near the boundary? I could use non-symmetric finite difference schemes, but I would expect my method (with ghost points) to work. Below is a sample of my code. For this case, the grid is equidistant ($\Delta_x = \Delta_y$, or in the code $h(1) = h(2)$).

! Jacobi scheme                                                                                          
if(method_order == second_order) then                                                                    
    ! 2nd order method: 5 points stencil                                                                 
    do jj=min_y,max_y                                                                                    
        do ii=min_x,max_x                                                                                
            Ujacob(ii,jj) = 0.25_dp                                                     &               
                &   * (Uold(ii-1,jj) + Uold(ii+1,jj) + Uold(ii,jj-1) + Uold(ii,jj+1)    &               
                &   - h(1)**2 * Urhs(ii,jj))                                                             
        end do                                                                                           
    end do                                                                                               
else if(method_order == fourth_order) then                                                               

    ! Fill buffer points                                                                                  
    Ujacob(     -1, 0:grid_ny-1) = -Ujacob(        1, 0:grid_ny-1)                                        
    Ujacob(grid_nx, 0:grid_ny-1) = -Ujacob(grid_nx-2, 0:grid_ny-1)                                        

    Ujacob(0:grid_nx-1,      -1) = -Ujacob(0:grid_nx-1,         1)                                        
    Ujacob(0:grid_nx-1, grid_ny) = -Ujacob(0:grid_nx-1, grid_ny-2)  

    ! 4th order method: 9 points stencil                            
    do jj=min_y,max_y                                                                                    
        do ii=min_x,max_x                                                                                
            Ujacob(ii,jj) = 1./60._dp                                                   &               
        &   *(-1._dp * (Uold(ii-2,jj) + Uold(ii+2,jj) + Uold(ii,jj-2) + Uold(ii,jj+2))  &               
        &   + 16._dp * (Uold(ii-1,jj) + Uold(ii+1,jj) + Uold(ii,jj-1) + Uold(ii,jj+1))  &               
        &   - 12._dp * h(1)**2 * Urhs(ii,jj))                                                            
        end do                                                                                           
    end do                                                                                               
end if


Following Wolfgang's advice, I fill my ghost points by mirroring what's happening next to the boundary. However, that does not help and the code still diverges.

  • $\begingroup$ When you say that $U(-1)=U(0)$, you assume that the function is constant zero outside the domain, and consequently that $U'(0)=0$. In other words, you are imposing a second (and wrong) boundary condition. I bet you'd get better results if you used $U(-1)=-U(1)$. $\endgroup$ Sep 5, 2016 at 21:26
  • $\begingroup$ @WolfgangBangerth thanks for the advice. See edit (I have changed the code section as well) $\endgroup$
    – solalito
    Sep 6, 2016 at 13:57
  • $\begingroup$ Can u provide full code ?, since we need extra equation for interior cells next to boundary, u can consider following form, here i'm using 1D $u_{i,j} =(u_{i+1,j} + u_{i-1,j})/2-h^2*(Urhs)$, since we know $u_{i,j}$ at boundary, we can subs. that value and use these equations. I'm not sure whether it will work or not . I believe it is better than applying gradient BC $\endgroup$
    – AGN
    Sep 6, 2016 at 15:36

1 Answer 1


Common Ghost cell treatment in CFD

Ghost cell treatment is an art in CFD, we can't arbitrarily choose ghost cell for all the problems. Let's assume $0$ is boundary point and $-1$ is ghost point and $u$ is variable then commonly used ghost cell procedure is, u(-1) =u(0), this is some sense equivalent $\frac{du}{dx}=0$ another procedure is u(-1) = -u(0), this kind of procedure is followed if we want normal velocity is zero at wall ($v_n=0$) commonly used in Euler equation.

Here I'm restricting myself to 1D heat equation. From the physics behind it, we can extend that to other elliptic equations. All the above B.C described may fail for heat equation.

  • If we impose u(-1) =u(0) that is equivalent to insulated wall B.C, but we not sure whether the wall is insulated or not. Maintaining same temperature at two points without B.C is non-physical
  • If we use u(-1) = -u(0), this also violate physics, lets say u(0)=1000 K, u(-1)=-1000 K. Is it physical? -1000 K doesn't exist. Even if exist we never allow this much jump in numerical scheme unless there is shock. Since it is elliptic solution, expected to be smooth.

Steady state 1-D Heat equation is $$\frac{d^2u}{dx^2} = 0$$ $2^{nd}$ order FD eqivalent is $$\frac{d^2u}{dx^2} = (u_{i-1}-2u_i+u_{i+1})/h^2 =0$$ then $$u_{i-1}=2*u_i-u_{i+1}$$ for i=0 and nx this become here $-1$ and $nx+1$ are ghost cell points $$u_{-1}=2*u_0-u_{1}$$ $$u_{nx+1}=2*u_{nx}-u_{nx-1}$$

these equations can be used for ghost cell points, please note that these equations are second order accurate and may have impact on solution.Because interior points are fouth-order scheme. Here I have attached a sample matlab code, you can play with it.

clear all
close all
u=0*x; %I.C
u(2)=10; %bc because x(1) is ghost cell
u(nx-1)=10; %bc because x(nx) is ghost cell
for j1=1:1:itmax
for i1=1:nx
    if i1==1;
        u(i1)=u(2)*2-u(i1+2);  % Ghost cell values are framed 
                                    %from GDE in second order

    if i1==nx;
    if i1==2
    if i1==nx-1
    if i1>2 && i1<nx-1
for i1=3:1:nx-2
    u_res(i1)=(-u(i1-2)+16*u(i1-1)+16*u(i1+1)-u(i1+2))/30-u(i1); %residue calc.
if max(abs(u_res))<0.0001 %convergence creteria


I will be happy if you do following things

  1. You can check the effect of second order ghost cell treatment for this problem on convergence
  2. You can set one $4^{th}$ order biased ghost cell interpolation scehme inside the code and check the effect of it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.