Poisson Equation and Relaxation

We want to solve the Poisson equation:

\[\nabla^2 \phi = f\]

We’ll do it in 1-d, writing our equation as:

\[\frac{d^2 \phi}{dx^2} = f\]

on a domain \([a, b]\).

Discretization

We will represent \(\phi\) and \(f\) on a 1-d uniform grid:

_images/fd_grid_bnd.png

Fig. 14 Our 1-d finite difference grid

Here we represent the discrete value of \(\phi\) with a subscript:

\[\phi_i = \phi(x_i)\]

and similarly for \(f_i\).

Notice that with this grid we have a point on each boundary. We will set these points to the boundary conditions:

\[\phi(a) = A\]
\[\phi(b) = B\]

We now discretize the Poisson equation. Let’s start with the Taylor expansions for \(\phi_{i+1}\) and \(\phi_{i-1}\):

\[\phi_{i\pm 1} = \phi_i \pm \left . \frac{d\phi}{dx} \right |_i \Delta x + \frac{1}{2} \left . \frac{d^2\phi}{dx^2} \right |_i \Delta x^2 \pm \frac{1}{6} \left . \frac{d^3\phi}{dx^3} \right |_i \Delta x^3 + \mathcal{O}(\Delta x^4)\]

adding \(\phi_{i+1}\) and \(\phi_{i-1}\), we get:

\[\phi_{i+1} + \phi_{i-1} = 2 \phi_i + \left . \frac{d^2\phi}{dx^2} \right |_i (\Delta x^2) + \mathcal{O}(\Delta x^4)\]

which we can solve for:

\[\left . \frac{d^2\phi}{dx^2} \right |_i = \frac{\phi_{i+1} - 2 \phi_i + \phi_{i-1}}{\Delta x^2} + \mathcal{O}(\Delta x^2)\]

This shows that our approximation to the second derivative is second-order accurate. We can then write the discrete Poisson equation as:

\[\frac{\phi_{i+1} - 2 \phi_i + \phi_{i-1}}{\Delta x^2} = f_i\]

Since the boundary values are fixed, for a grid of \(N\) points, we need to update the points \(1 \ldots N-2\) using this expression.

Relaxation

We can solve for the update for a single zone:

\[\phi_i = \frac{1}{2} (\phi_{i+1} + \phi_{i-1} - \Delta x^2 f_i)\]

Our solution procedure is to iteratively apply this, updating the values \(\phi_i\) as we go along. This process is called relaxation or smoothing, and the approach we will use, where we use the new information immediately as it becomes available is called Gauss-Seidel relaxation.

Stopping

We need a measure to determine when we should stop relaxing the solution. We define the residual, \(r\), as:

\[r_i = \frac{\phi_{i+1} - 2 \phi_i + \phi_{i-1}}{\Delta x^2} - f_i\]

This is a measure of how we’ll we satisfy the discrete equation.

We will keep iterating until the norm of the residual is small compared to the norm of the source:

\[\| r \| < \epsilon \| f \|\]

Once that is satisfied, then we have solved the problem as well as we can do with the number of points we are using in the discretization.

This requires us to define a vector norm. We’ll use the L2 norm:

\[\| v \| \equiv \left [ \Delta x \sum_{i=0}^{N-1} |v_i|^2 \right ]^{1/2}\]

Implementation

Here’s a class that implements smoothing:

Listing 93 poisson.H
#ifndef POISSON_H
#define POISSON_H

#include <vector>
#include <cassert>
#include <cmath>
#include <limits>
#include <fstream>
#include <iomanip>
#include <functional>

///
/// Solve the 1-d Poisson problem on a node-centered finite-difference
/// grid with Dirichlet BCs using Gauss-Seidel smoothing
///
class Poisson {

private:

    double xmin;
    double xmax;

    int N;

    std::vector<double> x;
    std::vector<double> phi;
    std::vector<double> f;
    
    double dx;

public:

    Poisson(double xmin_in, double xmax_in, int N_in)
        : xmin{xmin_in}, xmax{xmax_in}, N{N_in},
          x(N_in, 0), phi(N_in, 0), f(N_in, 0)
    {

        // initialize the coordinates
        assert (xmax > xmin);

        dx = (xmax - xmin) / static_cast<double>(N-1);

        for (int i = 0; i < N; ++i) {
            x[i] = xmin + static_cast<double>(i) * dx;
        }

    }

    ///
    /// set the left Dirichlet boundary condition
    ///
    void set_left_bc(double val) {phi[0] = val;}

    ///
    /// set the right Dirichlet boundary condition
    ///
    void set_right_bc(double val) {phi[N-1] = val;}

    ///
    /// set the source term, f, in L phi = f
    ///
    void set_source(std::function<double(double)> func) {
    
        for (int i = 0; i < N; ++i) {
            f[i] = func(x[i]);
        }
    }

    ///
    /// return the number of points
    ///
    int npts() {return N;}

    ///
    /// return the coordinate vector
    ///
    const std::vector<double>& xc() {return x;}

    ///
    /// return the source vector
    ///
    std::vector<double>& get_source() {
        return f;
    }

    ///
    /// return the solution vector
    ///
    std::vector<double>& get_phi() {
        return phi;
    }

    ///
    /// do Gauss-Seidel smoothing for n_smooth iterations
    ///
    void smooth(int n_smooth) {

        // perform Gauss-Seidel smoothing

        // we only operate on the interior nodes

        for (int i = 0; i < n_smooth; ++i) {

            for (int j = 1; j < N-1; ++j) {
                phi[j] = 0.5 * (phi[j-1] + phi[j+1] - dx * dx * f[j]);
            }
        }
    }

    ///
    /// solve the Poisson problem via relaxation until the residual
    /// norm is tol compared to the source norm
    ///
    void solve(double tol) {

        double err = std::numeric_limits<double>::max();

        double source_norm = norm(f);

        while (err > tol) {

            smooth(10);

            double r_norm = norm(residual());

            if (source_norm != 0.0) {
                err = r_norm / source_norm;
            } else {
                err = r_norm;
            }
        }

    }

    ///
    /// compute the residual
    ///
    std::vector<double> residual() {
        std::vector<double> r(N, 0);

        for (int i = 1; i < N-1; ++i) {
            r[i] = f[i] - (phi[i+1] - 2.0 * phi[i] + phi[i-1]) / (dx * dx);
        }

        return r;

    }

    ///
    /// given a vector e on our grid, return the L2 norm
    ///
    double norm(const std::vector<double>& e) {
        double l2{0.0};

        for (int i = 0; i < N; ++i) {
            l2 += e[i] * e[i];
        }

        l2 = std::sqrt(dx * l2);

        return l2;
    }

    ///
    /// output the solution to file fname
    ///
    void write_solution(const std::string& fname) {
        auto of = std::ofstream(fname);

        for (int i = 0; i < N; ++i) {
            of << std::setw(20) << x[i] << " " << std::setw(20) << phi[i] << std::setw(20) << f[i] << std::endl;
        }

    }


};
#endif

Let’s try it on the problem:

\[\phi^{\prime\prime} = \sin(x)\]

on \([0, 1]\) with the boundary conditions:

\[\phi(0) = 0\]
\[\phi(1) = 1\]

This has the analytic solution:

\[\phi(x) = -\sin(x) + (1 + \sin(1)) x\]

Here’s the driver that implements this problem:

Listing 94 poisson.cpp
#include "poisson.H"

const double TOL = 1.e-10;

int main() {

    auto p = Poisson(0.0, 1.0, 64);

    p.set_source([] (double x) {return std::sin(x);});

    p.set_left_bc(0.0);
    p.set_right_bc(1.0);

    p.solve(TOL);

    p.write_solution("poisson.txt");
}

Convergence

This method should converge as second-order – i.e., if we double the number of points, the error should go down by a factor of 4.

try it…

Let’s add a method to our class that accepts a function that provides the analytic solution and returns the norm of the error (solution compared to analytic solution).

Then run with 64, 128, 256 points and see how the error changes.

2D implementation

Here’s a 2D implementation:

We can solve the problem:

\[\phi_{xx} + \phi_{yy} = -2 \left [(1-6x^2)y^2(1-y^2) + (1-6y^2)x^2(1-x^2) \right ]\]

on \([0, 1] \times [0, 1]\) with homogeneous Dirichlet BCs. The analytic solution in this case is:

\[\phi(x, y) = (x^2 - x^4) (y^4 - y^2)\]

Here’s the driver:

Listing 96 poisson.cpp
#include "poisson2d.H"

const double TOL = 1.e-10;

int main() {

    auto p = Poisson(0.0, 1.0, 0.0, 1.0, 64, 64);

    p.set_source([] (double x, double y)
                 {return -2.0 * ((1.0 - 6.0 * x * x) * y * y * (1.0 - y * y) +
                                 (1.0 - 6.0 * y * y) * x * x * (1.0 - x * x));});

    p.set_xlo_bc(0.0);
    p.set_xhi_bc(0.0);
    p.set_ylo_bc(0.0);
    p.set_yhi_bc(0.0);

    p.solve(TOL);

    p.write_solution("poisson.txt");
}

This outputs the data in a way that can be visualized with gnuplot using:

set pm3d map
set size square
splot 'poisson.txt'
_images/poisson2d.png

We’ll use this example next week when we do parallel programming.