使用OpenMPI的二维fdtd

二维FDTD（有限差分时域）方法是求解电磁波传播问题的一种常见方法。使用OpenMPI进行并行计算可以加速FDTD方法的求解过程。

下面是一个使用OpenMPI实现的二维FDTD程序：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>

#define NX 100
#define NY 100
#define TMAX 200

int main(int argc, char **argv)
{
    int size, rank, i, j, t, istart, iend;
    double dx = 0.01, dy = 0.01, dt = 0.01;
    double Ez[NX][NY] = {0}, Hx[NX][NY] = {0}, Hy[NX][NY] = {0};
    double Sc = 1 / sqrt(2);
    double epsz = 8.854e-12, mu = 1.256e-6, c = 1.0 / sqrt(epsz * mu);
    double omega = 2 * M_PI * c / (dx * sqrt(2));

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    istart = rank * (NX - 1) / size + 1;
    iend = (rank + 1) * (NX - 1) / size;

    for (t = 0; t < TMAX; t++)
    {
        for (i = istart; i <= iend; i++)
        {
            for (j = 1; j < NY; j++)
            {
                Hx[i][j] -= dt / (mu * dy) * (Ez[i][j] - Ez[i][j - 1]);
            }
        }

        MPI_Status status;
        if (rank > 0)
        {
            MPI_Send(&Hx[istart][1], NY - 1, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&Hx[istart - 1][1], NY - 1, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD, &status);
        }
        if (rank < size - 1)
        {
            MPI_Send(&Hx[iend - 1][1], NY - 1, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&Hx[iend][1], NY - 1, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD, &status);
        }

        for (i = istart; i <= iend; i++)
        {
            for (j = 0; j < NY - 1; j++)
            {
                Hy[i][j] += dt / (mu * dx) * (Ez[i][j] - Ez[i - 1][j]);
            }
        }

        if (rank > 0)
        {
            MPI_Send(&Hy[istart][1], NY - 1, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&Hy[istart - 1][1], NY - 1, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD, &status);
        }
        if (rank < size - 1)
        {
            MPI_Send(&Hy[iend - 1][1], NY - 1, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&Hy[iend][1], NY - 1, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD, &status);
        }

        for (i = istart; i <= iend; i++)
        {
            for (j = 1; j < NY - 1; j++)
            {
                Ez[i][j] += dt / epsz * Sc * (Hy[i][j] - Hy[i - 1][j] - Hx[i][j] + Hx[i][j - 1]);
            }
        }

        if (rank > 0)
        {
            MPI_Send(&Ez[istart][1], NY - 2, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&Ez[istart - 1][1], NY - 2, MPI_DOUBLE, rank - 1, 0, MPI_COMM_WORLD, &status);
        }
        if (rank < size - 1)
        {
            MPI_Send(&Ez[iend - 1][1], NY - 2, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD);
            MPI_Recv(&Ez[iend][1], NY - 2, MPI_DOUBLE, rank + 1, 0, MPI_COMM_WORLD, &status);
        }
        if (rank == 0)
        {
            Ez[0][1] = sin(omega * t * dt);
        }
    }

    MPI_Finalize();

    return 0;
}

这个程序使用了MPI进行并行计算。每个进程计算一部分网格点，然后使用MPI_Send和MPI_Recv函数进行通信，将自己计算的边界值传递给相邻的进程。

在程序中，istart和iend表示当前进程计算的网格点范围。MPI_Send和MPI_Recv函数用于进程之间的通信，status参数用于存储通信的状态。

在主循环中，程序首先更新Hx和Hy的值，然后进行边界通信。接下来更新Ez的值，并再次进行边界通信。最后，程序在rank = 0的进程上更新Ez的初始值。

这个程序可以使用mpicc编译器进行编译：

mpicc -o fdtd fdtd.c -lm

然后使用mpirun命令运行程序：

mpirun -n 4 ./fdtd

这个命令将使用4个进程运行程序。在程序运行时，每个进程将计算一部分网格点，并与相邻的进程进行通信，以获得边界值。程序将在每个时间步长打印Ez的值。