Vibrating string

Solves the one-dimensional wave equation via finite differences, MPI version.

// solve the one-dimensional wave equation

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

const int N = 1001;  // total global number of grid points
const double L = 1, c = 1, dt = 0.001, t_end = 2.4;
enum { left_copy, right_copy };

// update string elongation
void string(const double *u, const double *u_old, double *u_new, int N, double eps) {
  u_new[0] = u[0];
  for (int i = 1; i < N - 1; ++i)
    u_new[i] = eps * (u[i - 1] + u[i + 1]) + 2.0 * (1.0 - eps) * u[i] - u_old[i];
  u_new[N - 1] = u[N - 1];
}

void *secure_malloc(size_t size) {
  void *p = malloc(size);
  if (p == NULL)
    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
  return p;
}

// initial elongation
static inline double u_0(double x) {
  if (x <= 0 || x >= L)
    return 0;
  return exp(-200.0 * (x - 0.5 * L) * (x - 0.5 * L));
}

// initial velocity
static inline double u_0_dt(double x) {
  return 0.0;
}

int main(int argc, char *argv[]) {
  int C_rank, C_size, *N_l, *N0_l;
  double dx = L / (N - 1), eps = dt * dt * c * c / (dx * dx), *u = NULL, *u_l, *u_old_l,
         *u_new_l, *u_temp;
  MPI_Status statuses[4];
  MPI_Request requests[4];
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &C_size);
  MPI_Comm_rank(MPI_COMM_WORLD, &C_rank);
  N_l = secure_malloc(C_size * sizeof(*N_l));
  N0_l = secure_malloc(C_size * sizeof(*N0_l));
  // calculate size and position of local grids
  // the local grids include one mirror grid point at each end or the boundary condition
  for (int i = 0; i < C_size; ++i) {
    // number of local grid points is N_l[C_rank]
    N_l[i] = (i + 1) * (N - 2) / C_size - i * (N - 2) / C_size + 2;
    // position of local
    N0_l[i] = i * (N - 2) / C_size;
  }
  u_old_l = secure_malloc((N_l[C_rank]) * sizeof(*u_old_l));
  u_l = secure_malloc((N_l[C_rank]) * sizeof(*u_l));
  u_new_l = secure_malloc((N_l[C_rank]) * sizeof(*u_new_l));
  // 1st time step uses elongation and velocity
  for (int i = 0; i < N_l[C_rank]; ++i) {
    double x = (i + N0_l[C_rank]) * dx;
    u_old_l[i] = u_0(x);
    u_l[i] = 0.5 * eps * (u_0(x - dx) + u_0(x + dx)) + (1.0 - eps) * u_0(x) + dt * u_0_dt(x);
  }
  // solve wave equation by using elongation at current time and one step before
  for (double t = 2 * dt; t <= t_end; t += dt) {
    string(u_l, u_old_l, u_new_l, N_l[C_rank], eps);
    MPI_Isend(&u_new_l[N_l[C_rank] - 2], 1, MPI_DOUBLE,
              C_rank + 1 < C_size ? C_rank + 1 : MPI_PROC_NULL, right_copy, MPI_COMM_WORLD,
              &requests[0]);
    MPI_Isend(&u_new_l[1], 1, MPI_DOUBLE, C_rank - 1 >= 0 ? C_rank - 1 : MPI_PROC_NULL,
              left_copy, MPI_COMM_WORLD, &requests[1]);
    MPI_Irecv(&u_new_l[0], 1, MPI_DOUBLE, C_rank - 1 >= 0 ? C_rank - 1 : MPI_PROC_NULL,
              right_copy, MPI_COMM_WORLD, &requests[2]);
    MPI_Irecv(&u_new_l[N_l[C_rank] - 1], 1, MPI_DOUBLE,
              C_rank + 1 < C_size ? C_rank + 1 : MPI_PROC_NULL, left_copy, MPI_COMM_WORLD,
              &requests[3]);
    MPI_Waitall(4, requests, statuses);
    u_temp = u_old_l;
    u_old_l = u_l;
    u_l = u_new_l;
    u_new_l = u_temp;
  }
  // exclude overlapping grid points
  for (int i = 0; i < C_size; ++i) {
    N_l[i] -= 2;
    N0_l[i] += 1;
  }
  // gather and print data
  if (C_rank == 0)
    u = secure_malloc(N * sizeof(*u));
  MPI_Gatherv(u_l + 1, N_l[C_rank], MPI_DOUBLE, u, N_l, N0_l, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  if (C_rank == 0) {
    u[0] = u[N - 1] = 0;  // boundary condition
    for (int i = 0; i < N; ++i)
      printf("%g\n", u[i]);
  }
  free(u);
  free(u_old_l);
  free(u_l);
  free(u_new_l);
  MPI_Finalize();
  return EXIT_SUCCESS;
}

Solves the one-dimensional wave equation via finite differences, MPL version.

// solve the time-dependent one-dimensional wave equation
// via a finite difference discretization and explicit time stepping

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <vector>
#include <mpl/mpl.hpp>

const int n{1001};        // total number of grid points
const double l{1};        // lengths of domain
const double c{1};        // speed of sound
const double dt{0.001};   // temporal step width
const double t_end{2.4};  // simulation time

enum class copy : int { left, right };

// update grid points
void string(const std::vector<double> &u, const std::vector<double> &u_old,
            std::vector<double> &u_new, double eps) {
  using size_type = std::vector<double>::size_type;
  const size_type N{u.size()};
  u_new[0] = u[0];
  for (size_type i{1}; i < N - 1; ++i)
    u_new[i] = eps * (u[i - 1] + u[i + 1]) + 2.0 * (1.0 - eps) * u[i] - u_old[i];
  u_new[N - 1] = u[N - 1];
}

// initial elongation of string
inline double u_0(double x) {
  if (x <= 0 or x >= l)
    return 0;
  return std::exp(-200.0 * (x - 0.5 * l) * (x - 0.5 * l));
}

// initial velocity of string
inline double u_0_dt([[maybe_unused]] double x) {
  return 0.0;
}

int main() {
  const double dx{l / (n - 1)};  // grid spacing
  const double eps{dt * dt * c * c / (dx * dx)};
  const mpl::communicator &comm_world{mpl::environment::comm_world()};
  const int c_size{comm_world.size()};
  const int c_rank{comm_world.rank()};
  std::vector<int> n_l, n_0_l;
  for (int i{0}; i < c_size; ++i) {
    // number of local grid points of process i
    n_l.push_back((i + 1) * (n - 2) / c_size - i * (n - 2) / c_size + 2);
    // position of local grid of process i within the global grid
    n_0_l.push_back(i * (n - 2) / c_size);
  }
  // grid data for times (t-dt), t and t+dt
  std::vector<double> u_old_l(n_l[c_rank]);
  std::vector<double> u_l(n_l[c_rank]);
  std::vector<double> u_new_l(n_l[c_rank]);
  // 1st propagation step uses current elongation and velocity
  // calculate all grid points including overlapping border data
  for (int i{0}; i < n_l[c_rank]; ++i) {
    double x = (i + n_0_l[c_rank]) * dx;
    u_old_l[i] = u_0(x);
    u_l[i] = 0.5 * eps * (u_0(x - dx) + u_0(x + dx)) + (1.0 - eps) * u_0(x) + dt * u_0_dt(x);
  }
  // propagate
  double t{2 * dt};
  while (t <= t_end) {
    // make one time step to get elongation
    string(u_l, u_old_l, u_new_l, eps);
    // update border data
    mpl::irequest_pool r;
    r.push(comm_world.isend(u_new_l[n_l[c_rank] - 2],
                            c_rank + 1 < c_size ? c_rank + 1 : mpl::proc_null, copy::right));
    r.push(comm_world.isend(u_new_l[1], c_rank - 1 >= 0 ? c_rank - 1 : mpl::proc_null,
                            copy::left));
    r.push(comm_world.irecv(u_new_l[0], c_rank - 1 >= 0 ? c_rank - 1 : mpl::proc_null,
                            copy::right));
    r.push(comm_world.irecv(u_new_l[n_l[c_rank] - 1],
                            c_rank + 1 < c_size ? c_rank + 1 : mpl::proc_null, copy::left));
    r.waitall();
    std::swap(u_l, u_old_l);
    std::swap(u_new_l, u_l);
    t += dt;
  }
  // finally, gather all the data at rank 0 and print result
  mpl::layouts<double> layouts;
  for (int i{0}; i < c_size; ++i)
    layouts.push_back(mpl::indexed_layout<double>({{n_l[i] - 2, n_0_l[i] + 1}}));
  mpl::contiguous_layout<double> layout(n_l[c_rank] - 2);
  if (c_rank == 0) {
    std::vector<double> u(n, 0);
    comm_world.gatherv(0, u_l.data() + 1, layout, u.data(), layouts);
    for (int i{0}; i < n; ++i)
      std::cout << dx * i << '\t' << u[i] << '\n';
  } else
    comm_world.gatherv(0, u_l.data() + 1, layout);
  return EXIT_SUCCESS;
}