# Boundary treatment with higher order methods

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


EDIT:

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.

• 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)$. Sep 5, 2016 at 21:26
• @WolfgangBangerth thanks for the advice. See edit (I have changed the code section as well) Sep 6, 2016 at 13:57
• 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
– AGN
Sep 6, 2016 at 15:36

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.

clc
clear all
close all
nx=101;
x=linspace(0,1,nx);
dx=1/(nx-1);
u=0*x; %I.C
u_res=u;
u(2)=10; %bc because x(1) is ghost cell
u(nx-1)=10; %bc because x(nx) is ghost cell
itmax=10000;
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
end

if i1==nx;
u(i1)=u(nx-1)*2-u(nx-2);
end
if i1==2
u(i1)=10;
end
if i1==nx-1
u(i1)=100;
end
if i1>2 && i1<nx-1
u(i1)=(-u(i1-2)+16*u(i1-1)+16*u(i1+1)-u(i1+2))/30;
end
end
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.
end
if max(abs(u_res))<0.0001 %convergence creteria
break
end

end
plot(x(2:nx-1),u(2:nx-1))


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.