Successive over-relaxation not converging (when not done in-place)

Question

I'm trying to find the potential given some boundary conditions using the successive over-relaxation method.

I have 2 solutions:

-One iterates over all elements and applies the formula field[y,x] = (1-alpha)*field[y,x] + (field[max(y-1,0),x] + field[min(y+1,field.shape[0]-1),x] + field[y,max(x-1,0)] + field[y,min(x+1,field.shape[1]-1)]) * alpha/4 in place. This is slow because it doesn't access memory in a nice way.

-The other one, I create 4 matrices shifted in the 4 directions by 1. I Apply the same formula by then adding the matrices up. This however doesn't take into account modifications done during the current iteration. This is significantly faster then the previous one.

With alpha = 1.9 the first algorithm converges while the second one doesn't. For alpha = 1.0 both converge but very slowly.

Can anyone tell me what I'm doing wrong? And how can I fix the fast solution.

Full code:

#! python3

import numpy
import math
import time

def solve_laplace(boundary, mask, file = None, alpha = 1.0, threshold = 0.0001):
    """
    We are using the successive over-relaxation method. We iterate until our solution changes less than some threshold value.

    Vm+1(x,y,...) = alpha*( ((Vm(x-1,y,...) + Vm(x+1,y,...) + Vm(x,y-1,...) + Vm(x,y+1,...) + ...)/(2*nr dimensions) ) + (1-alpha)*Vm(x,y,...)
    """

    dim = boundary.ndim

    threshold = 0.0001
    field = numpy.zeros_like(boundary)
    numpy.copyto(field, boundary, casting = "safe", where = mask)
    last_diff = float("infinity")

    for iter_nr in range(10000):#max number of iterations
        prev = field.copy() #make a copy of the field at the start of the iteration (python always stores pointers unless you explicitly copy something)

        for d in range(dim): #can be scaled to arbitrary dimensions, using 2D for testing

            #these 2 blocks are hard to follow but they work, read the comments
            front = prev[tuple(0 if i==d else slice(None) for i in range(dim))] #select front face of cube/whatever
            front = front[tuple(numpy.newaxis if i==d else slice(None) for i in range(dim))] #prepare it for next step
            front = numpy.concatenate((front,prev),d) #add it the previous iteration's result
            front = front[tuple(slice(-1) if i==d else slice(None) for i in range(dim))] #remove the back side of the previous iteration's result
            #we now have the volume shifted right by 1 pixel, x now corresponds to the x-1 term

            back = prev[tuple(-1 if i==d else slice(None) for i in range(dim))] #select back face of cube/whatever
            back = back[tuple(numpy.newaxis if i==d else slice(None) for i in range(dim))] #prepare it for next step
            back = numpy.concatenate((prev,back),d) #add it the previous iteration's result
            back = back[tuple(slice(1,None) if i==d else slice(None) for i in range(dim))] #remove the front side of the previous iteration's result
            #we now have the volume shifted left by 1 pixel, x now corresponds to the x+1 term

            field += (front + back) * alpha/(2*dim) #this part of the formula: alpha*( ((Vm(x-1,y,...) + Vm(x+1,y,...) + Vm(x,y-1,...) + Vm(x,y+1,...))/(2*nr dimensions)
            #numpy.copyto(field, boundary, casting = "safe", where = mask)

        field -= alpha*prev #this part of the formula: (1-alpha)*Vm(x,y,...)
        #reset values at boundaries
        numpy.copyto(field, boundary, casting = "safe", where = mask) 

        #check if the difference is less than threshold
        average = math.sqrt(numpy.average(field**2)) #sqrt of average of squares, just so i get a positive number
        diff = math.sqrt(numpy.average((field-prev)**2)) #standard deviation

        if last_diff < diff/average:
            print("Solution is diverging.")
            break

        if diff/average < threshold:
            print("Found solution after", iter_nr,"iteratiorn.")
            break

        last_diff = diff/average

    if file is not None:
        numpy.save(file,field)
    return field



def solve_laplace_slow_2D(boundary, mask, file = None, alpha = 1.9,threshold = 0.0001):
    """
    We are using the successive over-relaxation method. We iterate until our solution changes less than some threshold value.

    Vm+1(x,y,...) = alpha*( ((Vm(x-1,y,...) + Vm(x+1,y,...) + Vm(x,y-1,...) + Vm(x,y+1,...) + ...)/(2*nr dimensions) ) + (1-alpha)*Vm(x,y,...)
    """

    assert boundary.ndim == 2

    field = numpy.zeros_like(boundary)
    numpy.copyto(field, boundary, casting = "safe", where = mask) 
    last_diff = float("infinity")
    start_time = time.time()

    for iter_nr in range(10000):#max number of iterations
        prev = field.copy()
        for y in range(field.shape[0]):
            for x in range(field.shape[1]):
                if not mask[y,x]:
                    field[y,x] = (1-alpha)*field[y,x] + (field[max(y-1,0),x] + field[min(y+1,field.shape[0]-1),x] + field[y,max(x-1,0)] + field[y,min(x+1,field.shape[1]-1)]) * alpha/4

        #check if the difference is less than threshold
        average = math.sqrt(numpy.average(field**2)) #sqrt of average of squares, just so i get a positive number
        diff = math.sqrt(numpy.average((field-prev)**2)) #standard deviation

        if last_diff < diff/average:
            print("Solution is diverging.")
            break

        if diff/average < threshold:
            print("Found solution after the", iter_nr,"iteratiorn.")
            break

        if time.time() - start_time > 3600:
            print("Completed in an hour time at iteration:", iter_nr)
            break

        last_diff = diff/average

        #print(time.time() - start_time, iter_nr, last_diff)

    if file is not None:
        numpy.save(file,field)
    return field

def test():
    boundary = numpy.zeros((51,51))
    boundary[25,25] = 1
    for i in range(51):
        boundary[0,i] = -1
        boundary[50,i] = -1
        boundary[i,0] = -1
        boundary[i,50] = -1
    mask = (boundary != 0)

    print("Trying fast method:")
    solve_laplace(boundary,mask,alpha = 1.5) #diverges
    print("Trying slow method:")
    solve_laplace_slow_2D(boundary,mask,alpha = 1.5) #converges but is very slow

Answer 1

The "fast" version is (damped/over-relaxed) Jacobi instead of SOR, which is a multiplicative method. Even with damping factor of 1.0, Jacobi is not guaranteed to converge, as you would see if you applied it to a 1D problem, or added extreme anisotropy to your multi-dimensional problem.

Your SOR ("slow") implementation is slow because it is written in pure Python rather than because of the memory access. Indeed, the vectorized code accesses memory more times, and thus is actually worse in terms of memory bandwidth. But the interpreter is so much slower that it hides the bandwidth effects.

It is a fact of life that writing numerical kernels in a dynamic language without a JIT will give slow performance, the variability of which may not correlate with machine characteristics. You should be able to speed up the looped code using Cython or Numba, or just write it in a compiled language.

Now it is true that a different sort of local vectorization can be applied to Jacobi, but not to SOR. That can be partly accommodated, especially for SPD problems, by using polynomial acceleration such as Chebyshev (for known spectral bounds) or with Krylov methods (which are spectrally adaptive). Alternatively, you can use multicolor relaxation (e.g., "red-black Gauss-Seidel") to have a method with relatively fast guaranteed convergence. Incorporating into a multigrid algorithm will also greatly speed up your convergence (to a small constant number of iterations independent of problem size).

Answer 2

The second algorithm does look like a Jacobi iteration. In that case, you should try using a relaxation coefficient of less than one (2/3 is the typically prescribed value according to wikipedia's entry on the subject).

If you want to go with Jacobi iteration, scheduled relaxation may help a little. (Yang, Xiang; Mittal, Rajat (June 27, 2014). "Acceleration of the Jacobi iterative method by factors exceeding 100 using scheduled relaxation". Journal of Computational Physics. doi:10.1016/j.jcp.2014.06.010) Essentially, it involves using multiple dampening coffecients with corresponding prescribed numbers of iterations to form a larger cycle. The paper gives several out of the box relaxation schedules, so it should be easy to get running quickly with a minimal amount of overhead... personally, I would go with Conjugate Gradient. It seems to be faster for me.

I've tested out the SRJ method (using a clunky csc_matrix based implementation) to solve a finite difference discretization of poisson's equation with periodic boundaries on systems as large as 512x512 (2D) and 64x64x64 (3D). SRJ seems to converge with about 1/5 as many total iterations as the typically suggested singular over-relaxation coefficient of 2/3 requires...

But, conjugate gradient got far better convergence properties, even unconditioned.

Geometric Multigrid and Incomplete Cholesky preconditioned CG are also definitely good choices, particularly for problems in higher dimensions, and seem to give the best performance of the iterative linear solver's I've seen so far (although I'm still working on comparing with lean algebraic and combinatorial multigrid).

If you are simply solving laplace's equation with spatially constant coefficients, fourier transform based methods (via fftw) are blazingly fast.