Convergence stall when solving 2D Poisson PDE with pure Neumann boundaries (finite differences)

Question

I recently started coding a small library of 2D PDE solvers (time dependent and time independent), and my first attempt was a 2D Poisson equation of the form: $$\nabla(\epsilon\nabla\varphi)=\nabla\epsilon\nabla\varphi+\epsilon\nabla^2\varphi=\rho.$$ The general idea was to write a 2D Poisson solver capable of solving for $\varphi$ in arbitrary 2D geometry and with any combination of Dirichlet, Neumann, and periodic boundary conditions.

However, upon implementation and after initial testing, I noticed that when I simulate a rectangular box surrounded by 4 Neumann boundaries at the box edges, the convergence gets stuck at some high value, and it never reaches the desired level. The convergence criterion is that the square root of the sum of the squares of the residual is smaller than some small number (e.g. 1e-4) (see further down in the code example). During the stall, convergence gets stuck at value larger than 1 (the exact value depends on the number of points and the profile of the right-hand side [RHS]).

I have encountered several Stack Exchange posts discussing pure Neumann boundaries and their impact on the overall convergence progress, but none of the proposed solutions worked for me. I have tried removing the mean value of the RHS of the equation from the RHS during the residual calculation, which did not work. I have also tried to remove mean value of the main variable, which also did not resulted in any improvements.

The first version of the implementation was done in C++, based on the finite differences, where the derivatives were implemented in a central difference manner: $$ \frac{\partial\varphi}{\partial x}\to\frac{\varphi_{i+1,j}-\varphi_{i-1,j}}{2\Delta x} $$ and $$\frac{\partial^2\varphi}{\partial x^2}\to\frac{\varphi_{i+1,j}-2\varphi_{i,j}+\varphi_{i-1,j}}{\Delta x^2}.$$ The solver that I implemented is the iterative successive over-relaxation (SOR) method.

For the ease of use, during the discretization I introduced the following notation:

1. Value $\varphi_{i,j}$ is the center, denoted as $\varphi_{c}$
2. Values $\varphi_{i\pm 1,j}$ are right and left, denoted as $\varphi_{r,l}$, respectively
3. Values $\varphi_{i,j\pm 1}$ are top and bottom, denoted as $\varphi_{t,b}$, respectively
4. Coefficients multiplying $\varphi_{c,l,r,b,t}$ are denoted $C, L, R, B, T$, respectively, where in the simulation region where the calculations are to be done:$$C=2\epsilon_c(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}),$$ $$R,L=\frac{1}{\Delta x^2}[\epsilon_c\pm 0.25(\epsilon_r-\epsilon_l)],$$ $$T,B=\frac{1}{\Delta y^2}[\epsilon_c\pm 0.25(\epsilon_t-\epsilon_b)].$$

I have also merged RHS of the equation with possible Neumann boundary conditions (if present): $mRHS_c=\rho_c+Nbc_c.$ In addition, to save some possible unnecessary multiplications, I also stored "normalized" RHS as $nRHS_c=mRHS_c/C.$

Regarding the boundary conditions, on places where Dirichlet and Neumann BCs are expected, solver will skip those points, as for Dirichlet BCs those points are already solved, and for Neumann BCs I assume a sort of "electrode" or "out of simulation space"-region is present. However, both Dirichlet and Neumann BCs do appear through $\varphi_{r,l,t,b}$. For Dirichlet this is straight-forward, but for Neumann, I had to introduce "ghost" points at the boundary of the simulation region. The Neumann BCs are given as: $$\frac{\partial\varphi(x=X_0)}{\partial x}=G_0\to\frac{\varphi_{1,j}-\varphi_{-1,j}}{2\Delta x}=G_0; \frac{\partial\varphi(x=X_1)}{\partial x}=G_1\to\frac{\varphi_{N+1,j}-\varphi_{N-1,j}}{2\Delta x}=G_1$$ $$\frac{\partial\varphi(y=Y_0)}{\partial y}=H_0\to\frac{\varphi_{i,1}-\varphi_{i,-1}}{2\Delta y}=H_0; \frac{\partial\varphi(y=Y_1)}{\partial y}=H_1\to\frac{\varphi_{i,M+1}-\varphi_{i,M-1}}{2\Delta y}=H_1$$ So, if there is Neumann boundary $G_0$ at starting point $x=X_0$, I will express $\varphi_l$ as $\varphi_r-2\Delta xG_0$, at final point $x=X_1$ $\varphi_r$ becomes $\varphi_l+2\Delta xG_1$ (and similar for $y$ direction). The terms $\pm2\Delta (x,y)(G_{0,1},H_{0,1})$ I store in the $NBC$ variable, which is then included in the $mRHS$, as exemplified here in the code snapshot: $$mRHS_c=\rho_c+2\Delta x(LG_{0c}-RG_{1c})+2\Delta y(BH_{0c}-TH_{1c})$$ Matrices $G_{0c}, G_{1c}, H_{0c}, H_{1c}$ are non-zero at locations of Neumann BCs, and zero everywhere else. Generally, the overview of the simulation region is done with another integer matrix lattice which is equal:

1. 0 (denoted as EXT) at places external to the simulation region
2. 1 (denoted as SIM) at places where simulation occurs
3. 2 (denoted as PBC) at periodic boundaries
4. -1 (denoted as DBC) at Dirichlet boundaries
5. -2 (denoted as NBC) at Neumann boundaries.

For dummy example of 4x4 points (+ ghost points with Neumann BC at edges, Height times Width of the simulation box is 0.4x0.4 [a.u.]), lattice looks like:

-2 -2 -2 -2 -2 -2
-2  1  1  1  1 -2
-2  1  1  1  1 -2
-2  1  1  1  1 -2
-2  1  1  1  1 -2
-2 -2 -2 -2 -2 -2 

Finally, the discretized equation looks like: $$R\varphi_r+L\varphi_l+T\varphi_t+B\varphi_b-C\varphi_c=mRHS_c$$ and in SOR scheme, for iteration $k$: $$\varphi_c^{(k)}=\varphi_c^{(k-1)}+\omega_{SOR}(-\varphi_c^{(k-1)}-nRHS_c+[R\varphi_r^{(k-1)}+L\varphi_l^{(k-1)}+T\varphi_t^{(k-1)}+B\varphi_b^{(k-1)}]/C)=\varphi_c^{(k-1)}+\omega_{SOR}\frac{Residual_c^{(k-1)}}{C}$$ The parameter $\omega_{SOR}$ is the over-relaxation coefficient in range $(1,2)$. The initial guess at iteration $k=0$ is $\varphi_c=0$ (for all i,j in valid simulation region).

Typically, the convergence stall problem arises when I set $G_0=H_0=-G_1=-H_1=1$. Also, for the sake of simplicity, let $\rho_{i,j}=0$ for all i,j. Then matrix $mRHS_{i,j}$ is

0  0  0  0  0 0 
0 30 15 15 30 0 
0 15  0  0 15 0
0 15  0  0 15 0
0 30 15 15 30 0
0  0  0  0  0 0 

Typically, during the realistic test, the number of points and the box dimensions are larger, e.g. Nx=101; Ny=201, Lx=10; Ly=20;, respectively, with dx,y=Lx,y/(Nx,y-1). Here is the example of the code (for the sake of the argument, let matrix $\epsilon_{i,j}=1$ for all i,j):

void mergeNnorm_rhs(double* nRHS, double* mRHS, int Nx, int Ny, double dx, double dy, int* lattice, double* rho, double* epsilon, double* NBCl, double* NBCr, double* NBCb, double* NBCt, int ghost_region = 0)
{
    // Define local variables:
    int i, j, k;
    double h2x = 1.0 / dx / dx,
           h2y = 1.0 / dy / dy,
           hxy = 2.0 * (dx * dx + dy * dy) / dx / dx / dy / dy;
    double L, R, B, T, C;
    double eps_x, eps_y;

#pragma omp parallel for schedule(static) shared(lattice, nRHS, mRHS, rho, epsilon, NBCl, NBCr, NBCb, NBCt) private(i, j, k, L, R, B, T, C, eps_x, eps_y) firstprivate(Nx, Ny, dx, dy, h2x, h2y, hxy, ghost_region) default(none)
    for (i = ghost_region; i < Nx + ghost_region; i++)
        for (j = ghost_region; j < Ny + ghost_region; j++) {
            k = i * (Ny + 2 * ghost_region) + j;
            if (lattice[k] > EXT) {
                // Calculate matrix coefficients:
                eps_x = 0.25 * (epsilon[k + Ny] - epsilon[k - Ny]);
                eps_y = 0.25 * (epsilon[k + 1] - epsilon[k - 1]);
                // Calculate coefficients of the outter 2D stencil components:
                R = h2x * (epsilon[k] + eps_x);
                L = h2x * (epsilon[k] - eps_x);
                T = h2y * (epsilon[k] + eps_y);
                B = h2y * (epsilon[k] - eps_y);
                C = epsilon[k] * hxy;
                // Populate RHS matrices:
                mRHS[k] = rho[k] + 2.0 * dx * (L * NBCl[k] - R * NBCr[k]) + 2.0 * dy * (B * NBCb[k] - T * NBCt[k]);
                nRHS[k] = mRHS[k] / C;
            }
        }
}

/**/
double residual_2D_PDE_cfd_direct(int Nx, int Ny, double dx, double dy, int* lattice, double* phi, double* mRHS, double* epsilon)
{
    // Define local variables:
    int i, j, k;
    double diff;
    double phi_l, phi_r, phi_b, phi_t;
    double L, R, B, T, C;
    double eps_x, eps_y;
    double res = 0.0;
    double h2x = 1.0 / dx / dx,
           h2y = 1.0 / dy / dy,
           hxy = 2.0 * (dx * dx + dy * dy) / dx / dx / dy / dy;

#pragma omp parallel for schedule(static) shared(lattice, phi, mRHS, epsilon) private(i, j, k, diff, phi_l, phi_r, phi_b, phi_t, L, R, B, T, C, eps_x, eps_y) firstprivate(Nx, Ny, h2x, h2y, hxy) reduction(+:res) default(none)
    for (i = 1; i < Nx - 1; i++)
        for (j = 1; j < Ny - 1; j++) {
            k = i * Ny + j;
            if (lattice[k] > EXT) {
                // Calculate outter 2D stencil components:
                phi_r = lattice[k + Ny] == PBC ? phi[Ny + j]            : (lattice[k + Ny] == NBC ? phi[k - Ny] : phi[k + Ny]);
                phi_l = lattice[k - Ny] == PBC ? phi[Ny * (Nx - 2) + j] : (lattice[k - Ny] == NBC ? phi[k + Ny] : phi[k - Ny]);
                phi_t = lattice[k + 1]  == PBC ? phi[Ny * i + 1]        : (lattice[k + 1]  == NBC ? phi[k - 1]  : phi[k + 1]);
                phi_b = lattice[k - 1]  == PBC ? phi[Ny * i + Ny - 2]   : (lattice[k - 1]  == NBC ? phi[k + 1]  : phi[k - 1]);
                // Calculate matrix coefficients:
                eps_x = 0.25 * (epsilon[k + Ny] - epsilon[k - Ny]);
                eps_y = 0.25 * (epsilon[k + 1] - epsilon[k - 1]);
                // Calculate coefficients of the outter 2D stencil components:
                R = h2x * (epsilon[k] + eps_x);
                L = h2x * (epsilon[k] - eps_x);
                T = h2y * (epsilon[k] + eps_y);
                B = h2y * (epsilon[k] - eps_y);
                C = epsilon[k] * hxy;
                // Calculate current deviation between left- and right-hand sides and update residual calculations:
                diff = R * phi_r + L * phi_l + T * phi_t + B * phi_b - C * phi[k] - mRHS[k];
                res += diff * diff;
            }
        }

    // Represent residual as square root of the sum of the squares of the LHS and RHS difference:
    return sqrt(res);
}
/**/
void OddEven_2D_PDE_cfd_direct1D(void (*solver)(int, int, double, double, int, double, int*, double*, double*, double*), int Nx, int Ny, double dx, double dy, int* lattice, double* phi, double* nRHS, double* mRHS, double* epsilon, double Rconv, int SORtype = MT0, int verbose = 0)
{
    // Define local variables:
    int iter = 0;
    int Nmax = 100 * Nx * Ny;
    int Ncheck = 1;
    int Nstage = 2;
    double eps = 1e-4;
    double residual = 0.0;
    double w = 1.0;

    do {
        // Calculate residual every Ncheck iterations:
        residual = residual_2D_PDE_cfd_direct(Nx, Ny, dx, dy, lattice, phi, mRHS, epsilon);

        // Do Ncheck iterations of solving the 2D elliptic PDE with designated solver:
        for (int loc_it = 0; loc_it < Ncheck; loc_it++) {
            iter++;
            for (int stage = 0; stage < Nstage; stage++) {
                switch (SORtype) {
                case MT0:
                    w = Rconv;
                    break;
                case NR1:
                    w = (iter == 1 && stage == 0) ? 1.0 : ((iter == 1) ? 1.0 / (1.0 - 0.5 * Rconv * Rconv) : 1.0 / (1.0 - 0.25 * Rconv * Rconv * w));
                    break;
                case NR2:
                    w = (iter == 1 && stage == 0) ? 1.0 : ((iter == 1) ? 1.0 / (1.0 - 0.5 * Rconv * Rconv) : 1.0 / (1.0 - 0.25 * Rconv * Rconv * w));
                    break;
                }
                (*solver)(Nx, Ny, dx, dy, stage, w, lattice, phi, nRHS, epsilon);
            }
        }
        //cout << scientific << setw(16) << setprecision(12) << setfill(' ') << residual << "\n";
    } while (residual > eps && iter <= Nmax);

    // Calculate residual of the final iteration:
    residual = residual_2D_PDE_cfd_direct(Nx, Ny, dx, dy, lattice, phi, mRHS, epsilon);

    // Report to user the status of the solution:
    if (verbose) {
        if (iter >= Nmax) printf("Warning: Maximum number of iterations reached! (%d)\nFinal convergence achieved:\t%.10e\n", iter, residual);
        else {
            printf("Solution converged!\n\tTotal number of iterations:\t%d\n\tConvergence achieved:\t%.10e\n", iter, residual);
        }
    }
}

/**/
void SORoe_solver_cfd_direct1D(int Nx, int Ny, double dx, double dy, int stage, double w, int* lattice, double* phi, double* nRHS, double* epsilon)
{
    // Define local variables:
    int i, j, k;
    double gx = dx * dx * 0.5 / (dx * dx + dy * dy),
           gy = dy * dy * 0.5 / (dx * dx + dy * dy);
    double r, l, t, b;
    double phi_l, phi_r, phi_b, phi_t;
    double eps_x, eps_y;

#pragma omp parallel for schedule(static) shared(lattice, phi, nRHS, epsilon) private(i, j, k, phi_l, phi_r, phi_b, phi_t, l, r, b, t, eps_x, eps_y) firstprivate(Nx, Ny, gx, gy, stage, w) default(none)
    for (i = 1; i < Nx - 1; i++)
        for (j = 1 + ((i + stage) % 2); j < Ny - 1; j += 2) {
            k = i * Ny + j;
            if (lattice[k] > EXT) {
                // Calculate outter 2D stencil components:
                phi_r = lattice[k + Ny] == PBC ? phi[Ny + j]            : (lattice[k + Ny] == NBC ? phi[k - Ny] : phi[k + Ny]);
                phi_l = lattice[k - Ny] == PBC ? phi[Ny * (Nx - 2) + j] : (lattice[k - Ny] == NBC ? phi[k + Ny] : phi[k - Ny]);
                phi_t = lattice[k + 1]  == PBC ? phi[Ny * i + 1]        : (lattice[k + 1]  == NBC ? phi[k - 1]  : phi[k + 1]);
                phi_b = lattice[k - 1]  == PBC ? phi[Ny * i + Ny - 2]   : (lattice[k - 1]  == NBC ? phi[k + 1]  : phi[k - 1]);
                // Calculate epsilon corrections
                eps_x = 0.25 * (epsilon[k + Ny] - epsilon[k - Ny]) / epsilon[k];
                eps_y = 0.25 * (epsilon[k + 1]  - epsilon[k - 1])  / epsilon[k];
                // Calculate normalized coefficients of the outter 2D stencil components:
                r = gy * (1.0 + eps_x);
                l = gy * (1.0 - eps_x);
                t = gx * (1.0 + eps_y);
                b = gx * (1.0 - eps_y);
                // Update SOR iteration:
                phi[k] += w * (r * phi_r + l * phi_l + t * phi_t + b * phi_b - phi[k] - nRHS[k]);
            }
        }
}

If you can please provide me with any kind of insight or advice on what should I do to resolve the possible convergence issue, please let me know. My guess is that SOR should be able to tackle this kind of pure Neumann problems, but I am stuck for some time with this, and I do not know how to resolve this problem. Many thanks in advance~

Answer 1

For this equation, with the Neumann boundary conditions the problem becomes mathematically ill-posed; if there is a solution $\phi(x,y)$ then $\phi(x,y)+const$ is also a solution. This ill-posedness manifests itself in your solver's inability to converge.