parallel jacobi for Five-point difference format

格式分析

Jacobi迭代的基本形式
$$
x_i^{k+1}=\frac{1}{a_{ii}}\left(b_i-\sum_{j\neq i}a_{ij}x_j^k\right)
$$
对于二维Laplace方程，$x,y$方向采用相同的步长，在$x,y$方向采用中心差分离散可得
$$
\frac{\partial^2u(x,y)}{\partial x^2}\big|_{(x_i,y_j)}\approx\frac{2u(x_i,y_j)-u(x_{i-1},y_j)-u(x_{i+1},y_j)}{h^2}
$$

$$
\frac{\partial^2u(x,y)}{\partial y^2}\big|_{(x_i,y_j)}\approx\frac{2u(x_i,y_j)-u(x_{i},y_{j-1})-u(x_{i},y_{j+1})}{h^2}
$$

带入二维Laplace方程即可得Laplace方程在$(x_i,y_i)$点的近似离散方程
$$
4U_{i,j}-U_{i-1,j}-U_{i+1,j}-U_{i,j-1}-U_{i,j+1}=h^2f_{i,j}
$$
伪代码格式

while (not converged) {
    for (i,j)
        xnew[i][j] = (x[i+1][j] + x[i-1][j] + x[i][j+1] + x[i][j-1]+h*h*f)/4;
    for (i,j)
        x[i][j] = xnew[i][j];
}

误差范数估计采用如下

diffnorm = 0;
for (i,j)
    diffnorm += (xnew[i][j] - x[i][j]) * (xnew[i][j] - x[i][j]);
diffnorm = sqrt(diffnorm);

并行程序

为了程序便于实现，本程序未处理边界，可以globMat中处理好边界然后传递给localMat，这样程序会好很多。但是写好了不想改了，就选择了边界为0的一个例子
$$
\begin{aligned}\Delta u &=f \quad \Omega\in[0,1]\times[0,1] \\ u(x,y)&=0 \quad [x,y]\in \partial \Omega\end{aligned}
$$
这里$u=x^2(x-1)^2y^2(y-1)^2$​，
$$
\begin{aligned}f&=2(x^4(6y^2-6y+1)-2x^3(6y^2-6y+1)+\\&x^2(6y^4-12y^3+12y^2-6y+1)-6x(y-1)^2y^2+(y-1)^2y^2)\end{aligned}
$$
。

/*------------------------------------------------

JACOBI PARALLEL CODE
SOLVE \nabla u=f

------------------------------------------------*/

#include <math.h>
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#define N 8
#define eps  1e-9
const double pi = acos(-1.0);
const double T = 1.0;
const long long int maxIter = 1e10;

double uFun(double x,double y){
    return x*x*(x-1)*(x-1)*y*y*(y-1)*(y-1);
}

int main(int argc, char* argv[]) {
    double h = T / (N - 1);

    int rank, size;
    MPI_Status status;
    MPI_Init(&argc, &argv);

    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    if (N % size != 0) {
        MPI_Abort(MPI_COMM_WORLD, 1);
    }
    MPI_Barrier(MPI_COMM_WORLD);
    int rows = N / size;
    double(*localMat)[N] = new double[rows + 2][N];
    double(*globMat)[N] = new double[rows][N];

    for (int i = 0; i < rows + 2; i++) {
        for (int j = 0; j < N; j++) {
            *(*(localMat + i) + j) = 0;
        }
    }
    int iterFrist = rank == 0 ? 2 : 1;
    int iterLast = rank == size - 1 ? rows - 1 : rows;

    double diffLoc = 1, diffGlob = 1;
    int cnt = 0;
    double startTimes = MPI_Wtime();
    while (diffGlob > eps && cnt < maxIter) {
        if (rank < size - 1) {
            MPI_Send(*(localMat + rows), N, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD);
        }
        if (rank > 0) {
            MPI_Recv(*localMat, N, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD, &status);
        }
        if (rank > 0) {
            MPI_Send(*(localMat + 1), N, MPI_DOUBLE, rank - 1, 1, MPI_COMM_WORLD);
        }
        if (rank < size - 1) {
            MPI_Recv(*(localMat + rows + 1), N, MPI_DOUBLE, rank + 1, 1, MPI_COMM_WORLD, &status);
        }
        diffLoc = 0;
        for (int i = iterFrist; i < iterLast + 1; i++) {
            for (int j = 1; j < N - 1; j++) {
                //f(x,y) ---------
                double xx = h * (rank * rows + i - 1);
                double yy = h * j;
                double f = 2 * (pow(xx, 4) * (6 * yy * yy - 6 * yy + 1) -
                                2 * pow(xx, 3) * (6 * yy * yy - 6 * yy + 1) +
                                xx * xx *
                                    (6 * pow(yy, 4) - 12 * pow(yy, 3) +
                                     12 * pow(yy, 2) - 6 * yy + 1) -
                                6 * xx * (yy - 1) * (yy - 1) * yy * yy +
                                (yy - 1) * (yy - 1) * yy * yy); 
                //----------------
                *(*(globMat + i - 1) + j) =
                    (localMat[i][j + 1] + localMat[i][j - 1] + localMat[i - 1][j] +
                     localMat[i + 1][j] - f * h * h) /
                    4.0;
                diffLoc += (globMat[i - 1][j] - localMat[i][j]) * (globMat[i - 1][j] - localMat[i][j]);
            }
        }
        for (int i = iterFrist; i < iterLast + 1; i++) {
            for (int j = 1; j < N - 1; j++) {
                *(*(localMat + i) + j) = *(*(globMat + i - 1) + j);
            }
        }
        MPI_Allreduce(&diffLoc, &diffGlob, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        diffGlob = sqrt(diffGlob);
        cnt++;
    }

    double endTimes = MPI_Wtime();
    double Final[N][N];

    MPI_Gather(*(localMat+1), N * rows, MPI_DOUBLE, Final, N * rows, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    delete [] localMat;
    delete [] globMat;
    double infMax=-1.0;
    if (rank == 0) {
        FILE* f = fopen("data1.txt", "w");
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                fprintf(f, "%.12lf,", Final[i][j]);
                infMax = fmax(infMax,fabs(Final[i][j]-uFun(i*h,j*h)));
            }
            fprintf(f, "\n");
        }
        fclose(f);
    }
    MPI_Barrier(MPI_COMM_WORLD);
    if (rank == 0) {
        fprintf(stdout, "Total time is: %lf\n %d \n %.12lf\n", endTimes - startTimes, cnt,infMax);
    }
    MPI_Finalize();
    return 0;
}

64点差分三维图

import numpy as np 
import matplotlib.pyplot as plt

filename='./data1.txt'
Z = np.genfromtxt(filename,delimiter=',',dtype=np.float)
Z=np.delete(Z,-1,axis=1)

x = np.arange(0,1,1/64)
y = np.arange(0,1,1/64)
X,Y = np.meshgrid(x,y)

fig = plt.figure('himmelblau')
ax = fig.gca(projection='3d')
ax.plot_surface(X,Y,Z,cmap='rainbow')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()

误差分析图