# 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)$$

$$\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}$$

$$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);

## 并行程序

\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}

\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()