# Are there necessary conditions for SOR algorithm to converge when used for 2D discrete Poisson problem, with PBC, besides solvability condition?

I am trying to solve a 2D poisson problem that is supposed to represent diffusion of chemicals on a grid: $$\nabla^2 R_{ij}=f_{ij}$$. I am discretizing the problem with the standard central discretization, and since I have PBC this leads to singularity in the matrix of the linear problem, which means that a solvability condition is induced; to account for this, I am making sure that what I solve is actually $$\nabla^2 R_{ij}=\tilde{f_{ij}}$$ with $$\tilde{f_{ij}}=f_{ij}-f_{avg}$$, so that the source term sums to zero over all the domain. This is the code I'm using:

# define SOR(param,mat,acc) function to implement the SOR algorithm for integration
# on a grid; assumes fixed boundary conditions

def SOR_point(param,iter):
"""
R: current resource concentration (on current grid) matrice n x n x n_R
N: state vector for individuals (complete matrix on the current grid)
param: dictionary containing parameters
mat: dictionary containing matrices
acc: desired accuracy, determins when to stop algorithm

returns the equilibrium matrix n x n x n_r with concentrations and prints execution time

"""

n   = param['n']               # grid size
n_r = 1.                       # numer of nutrients
dx  = param['L0']/n            # step size

# initialization and boundary conditions setting
R_best = np.zeros((n+2,n+2))
avg_r = np.mean(R_best)
R_best= R_best-avg_r

# keep track of size of update
delta = R_best.copy()
delta_max = np.max(R_best)
delta_max_list = [2,1]

# make source
f = f_point(n+2,10)
avg = np.mean(f)
f = f-avg

# SOR algorithm and keep track of time
t1 = time()

it=0
while (it<iter):
it+=1
#impose PBC through ghost points
R_best[0, 1:n+1]  = R_best[-2, 1:n+1]
R_best[-1, 1:n+1] = R_best[1, 1:n+1]
R_best[1:n+1, 0]  = R_best[1:n+1, -2]
R_best[1:n+1, -1] = R_best[1:n+1, 1]
# implement red/blue updates to make starting point ininfluential
# loop on red
for i in np.arange(1,n+1):
start = 2 if i%2==0 else 1
for j in np.arange(start,n+1,2):
# next index for pbc (for the previous one pbc are automatic)
delta[i,j] = 0.25*(R_best[i,j+1]+R_best[i,j-1]+R_best[i+1,j]+R_best[i-1,j]-f[i,j]*dx**2)-R_best[i,j]
R_best[i,j] += delta[i,j]*param['sor']
# loop in blue
for i in np.arange(1,n+1):
start = 1 if i%2==0 else 2
for j in np.arange(start,n+1,2):
# next index for pbc (for the previous one pbc are automatic)
delta[i,j] = 0.25*(R_best[i,j+1]+R_best[i,j-1]+R_best[i+1,j]+R_best[i-1,j]-f[i,j]*dx**2)-R_best[i,j]
R_best[i,j] += delta[i,j]*param['sor']

delta_max = np.max(np.abs(delta[1:n+1,1:n+1]))
delta_max_list.append(delta_max)

print("N_iter %d delta_max %e\r" % (len(delta_max_list), delta_max), end='')

t2 = time()
print("\nTotal running time: %.2f min" % ((t2-t1)/60))
print("Code speed: %.1f iterations per second" %(len(delta_max_list)/(t2-t1)))

return R_best[1:n+1,1:n+1]


where

def f_point(n,c):
i=np.random.randint(n)
j=np.random.randint(n)
mat = np.zeros((n,n))
mat[i,j]=c
return mat


as you can see, I am using red-blue grid iteration. My problem is that the $$\delta$$, i.e. the update of the algorithm stops decreasing at some point, hence I start moving away from convergence, and my solution starts increasing in absolute value. I have tried playing around with the relaxation parameter, but this doesn't help. What could be the issue here? I checked my implementation of periodic boundary condition and I see no mistake, am I missing something that needs to be ensured to guarantee convergence of the algorithm? If so, how can I fix this to numerically solve the problem on this grid? Note that this is just a preliminary try with a point-source, but the specific problem I need to solve has a very inhomogeneous source term that is non-zero at each point of the grid.

• I'm going to analyze this question. @Federica Sibilla, are you still interested or did you already solve it by yourself? May 11 at 8:34
• @Rigel I am now using a different algorithm but I’d be very interested in understanding the problem anyways, if you could help!! May 12 at 16:45

You are trying to solve a Poisson problem with periodic boundary conditions. As you know, it only has solutions when the integral of the source term is zero. However, due to finite precision arithmetic, your source term induces some tiny contaminant inflow into the domain, and that adds up over time. I edited your code (see below) in such a way that you get in a more ideal situation where there is practically no additional contaminant entering/exiting the domain. However, for production code, I suggest to loook for methods for controlling the tiny inflow that you have, in order to be able to solve more generic problems. It should not be difficult to find one and implement it (I had no time tor that right now).

In addition to that, I suggest repeating the setting of the BCs also right before the blue update since the red update also edits the cells that are affected by the BCs.

In the code here below I also added some printing and debugging stuff that might be useful for you in the future (you'll need the matplotlib package to run it).

Check the output picture "delta.png", you'll see that now the delta drops until machine accuracy (1e-16), and then it stays approximately constant at that level, which is what we expect to see.

This is the main file you should run:

import matplotlib.pyplot as plt
import os, glob

#%%
iter = 200
L0 = 1
n = 14
sor_omega = 1.5

R_best, R_hist, delta_max_list = sor.SOR_point({'L0':L0, 'n':n, 'sor':sor_omega}, iter)

#%%
for R_hist_old_file in glob.glob('R_hist*.png'):
os.remove(R_hist_old_file)

for ii in range(4):
fig, ax = plt.subplots()
ax.imshow(R_hist[ii], vmin=-0.04, vmax=0.04)
ax.set_title(f'{ii}')
fig.savefig(f'R_hist[{ii}].png')

for ii in range(5, iter+1, 20):
fig, ax = plt.subplots()
ax.imshow(R_hist[ii], vmin=-0.04, vmax=0.04)
ax.set_title(f'{ii}')
fig.savefig(f'R_hist[{ii}].png')

fig, ax = plt.subplots()
ax.semilogy(delta_max_list, '-o')
fig.savefig('delta.png')


and it calls this edited version of the code you shared:

import numpy as np

# define SOR(param,mat,acc) function to implement the SOR algorithm for integration
# on a grid; assumes fixed boundary conditions

def SOR_point(param,iter):
"""
R: current resource concentration (on current grid) matrice n x n x n_R
N: state vector for individuals (complete matrix on the current grid)
param: dictionary containing parameters
mat: dictionary containing matrices
acc: desired accuracy, determins when to stop algorithm

returns the equilibrium matrix n x n x n_r with concentrations and prints execution time

"""

n   = param['n']               # grid size
n_r = 1.                       # numer of nutrients
dx  = param['L0']/n            # step size

# initialization and boundary conditions setting
R_best = np.zeros((n+2,n+2))
avg_r = np.mean(R_best)
R_best= R_best-avg_r

# keep track of size of update
delta = R_best.copy()
delta_max = np.max(R_best)
delta_max_list = [2,1]

# make source
f = f_point(n+2,10)
avg = np.mean(f)
f = f-avg

# SOR algorithm and keep track of time
# t1 = time()
# Make a list for recording all the intermediate steps [Rigel]
R_hist = []
it=0
while (it<iter):
it+=1
#impose PBC through ghost points
R_best[0, 1:n+1]  = R_best[-2, 1:n+1]
R_best[-1, 1:n+1] = R_best[1, 1:n+1]
R_best[1:n+1, 0]  = R_best[1:n+1, -2]
R_best[1:n+1, -1] = R_best[1:n+1, 1]

# implement red/blue updates to make starting point ininfluential
# loop on red
for i in np.arange(1,n+1):
start = 2 if i%2==0 else 1
for j in np.arange(start,n+1,2):
# next index for pbc (for the previous one pbc are automatic)
delta[i,j] = 0.25*(R_best[i,j+1]+R_best[i,j-1]+R_best[i+1,j]+R_best[i-1,j]-f[i,j]*dx**2)-R_best[i,j]
R_best[i,j] += delta[i,j]*param['sor']
# loop in blue
# Repeat BC setting here [Rigel]
R_best[0, 1:n+1]  = R_best[-2, 1:n+1]
R_best[-1, 1:n+1] = R_best[1, 1:n+1]
R_best[1:n+1, 0]  = R_best[1:n+1, -2]
R_best[1:n+1, -1] = R_best[1:n+1, 1]
for i in np.arange(1,n+1):
start = 1 if i%2==0 else 2
for j in np.arange(start,n+1,2):
# next index for pbc (for the previous one pbc are automatic)
delta[i,j] = 0.25*(R_best[i,j+1]+R_best[i,j-1]+R_best[i+1,j]+R_best[i-1,j]-f[i,j]*dx**2)-R_best[i,j]
R_best[i,j] += delta[i,j]*param['sor']

delta_max = np.max(np.abs(delta[1:n+1,1:n+1]))
delta_max_list.append(delta_max)
# Record the current update, including bcs [Rigel]
R_hist.append(R_best.copy())

print("N_iter %d delta_max %e\r" % (len(delta_max_list), delta_max), end='')

# t2 = time()
# print("\nTotal running time: %.2f min" % ((t2-t1)/60))
# print("Code speed: %.1f iterations per second" %(len(delta_max_list)/(t2-t1)))

# Returning also the previous values for debugging [Rigel]
return R_best[1:n+1,1:n+1], R_hist, delta_max_list

def f_point(n,c):
# Fix the point so we reduce variablity and better understand the problem [Rigel]
i = 7 # i=np.random.randint(n)
j  = 10 # j=np.random.randint(n)
mat = np.zeros((n,n))
mat[i,j] = c
mat[j,i] = -c
return mat