From 5bd1bccffd411384b02ca772822d87fac126e67f Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Thu, 23 May 2019 11:39:59 +0200 Subject: mg2d: add support for MPI-based multi-component solves For the moment, only the finest component is distributed, any coarser levels are gathered to rank 0. That should change in the future. --- mg2d_mpi_test.c | 283 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 283 insertions(+) create mode 100644 mg2d_mpi_test.c (limited to 'mg2d_mpi_test.c') diff --git a/mg2d_mpi_test.c b/mg2d_mpi_test.c new file mode 100644 index 0000000..1bd2784 --- /dev/null +++ b/mg2d_mpi_test.c @@ -0,0 +1,283 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include "mg2d.h" +#include "mg2d_boundary.h" +#include "mg2d_constants.h" + +#include "components.h" + +#define ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x)) +#define MIN(x, y) ((x) > (y) ? (y) : (x)) + +#define MAXITER 64 +#define TOL 5e-15 + +#define DOMAIN_SIZE 1.0 +#define FD_STENCIL 2 + +static const double pde_coeffs[MG2D_DIFF_COEFF_NB] = { + [MG2D_DIFF_COEFF_00] = 1.0, + [MG2D_DIFF_COEFF_10] = 0.9, + [MG2D_DIFF_COEFF_01] = 1.1, + [MG2D_DIFF_COEFF_20] = 1.2, + [MG2D_DIFF_COEFF_02] = 0.8, + [MG2D_DIFF_COEFF_11] = 0.7, +}; + +#if 1 +static double sol_00(double x, double y) +{ + return sin(M_PI * x) * sin(M_PI * y); +} +static double sol_10(double x, double y) +{ + return M_PI * cos(M_PI * x) * sin(M_PI * y); +} +static double sol_01(double x, double y) +{ + return M_PI * sin(M_PI * x) * cos(M_PI * y); +} +static double sol_20(double x, double y) +{ + return -M_PI * M_PI * sol_00(x, y); +} +static double sol_02(double x, double y) +{ + return -M_PI * M_PI * sol_00(x, y); +} +static double sol_11(double x, double y) +{ + return M_PI * M_PI * cos(M_PI * x) * cos(M_PI * y); +} +#define BC_TYPE MG2D_BC_TYPE_FIXVAL +#else +static double sol_00(double x, double y) +{ + return cos(M_PI * x) * cos(M_PI * y); +} +static double sol_10(double x, double y) +{ + return -M_PI * sin(M_PI * x) * cos(M_PI * y); +} +static double sol_01(double x, double y) +{ + return -M_PI * cos(M_PI * x) * sin(M_PI * y); +} +static double sol_20(double x, double y) +{ + return -M_PI * M_PI * sol_00(x, y); +} +static double sol_02(double x, double y) +{ + return -M_PI * M_PI * sol_00(x, y); +} +static double sol_11(double x, double y) +{ + return M_PI * M_PI * sin(M_PI * x) * sin(M_PI * y); +} +#define BC_TYPE MG2D_BC_TYPE_REFLECT +#endif + +static double (*sol[MG2D_DIFF_COEFF_NB])(double x, double y) = { + [MG2D_DIFF_COEFF_00] = sol_00, + [MG2D_DIFF_COEFF_10] = sol_10, + [MG2D_DIFF_COEFF_01] = sol_01, + [MG2D_DIFF_COEFF_20] = sol_20, + [MG2D_DIFF_COEFF_02] = sol_02, + [MG2D_DIFF_COEFF_11] = sol_11, +}; + +int main(int argc, char **argv) +{ + MG2DContext *ctx = NULL; + long int gridsize; + int ret = 0; + + DomainGeometry *dg = NULL; + DomainComponent *dc = NULL; + + size_t patch_start, patch_end, patch_size_y; + + char processor_name[MPI_MAX_PROCESSOR_NAME]; + int nb_processes, rank, processor_name_len; + + if (argc < 2) { + fprintf(stderr, "Usage: %s \n", argv[0]); + return 1; + } + gridsize = strtol(argv[1], NULL, 0); + if (gridsize <= 0) { + fprintf(stderr, "Invalid parameters: %ld\n", gridsize); + return 1; + } + + MPI_Init(NULL, NULL); + + MPI_Comm_size(MPI_COMM_WORLD, &nb_processes); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Get_processor_name(processor_name, &processor_name_len); + + fprintf(stderr, "This is process %d out of %d, running on %s\n", + rank, nb_processes, processor_name); + + dg = mg2di_dg_alloc(nb_processes); + if (!dg) { + fprintf(stderr, "Error allocating domain geometry\n"); + ret = 1; + goto fail; + } + + dg->domain_size[0] = gridsize; + dg->domain_size[1] = gridsize; + + for (unsigned int i = 0; i < dg->nb_components; i++) { + size_t patch_start, patch_end, patch_size_y; + + patch_size_y = (gridsize + nb_processes - 1) / nb_processes; + patch_start = i * patch_size_y; + patch_end = MIN((i + 1) * patch_size_y, gridsize); + patch_size_y = patch_end - patch_start; + if (patch_size_y <= 0) { + fprintf(stderr, "Too many processes for grid size %ld: %d\n", gridsize, nb_processes); + ret = 1; + goto fail; + } + + dg->components[i].interior.start[0] = 0; + dg->components[i].interior.start[1] = patch_start; + dg->components[i].interior.size[0] = gridsize; + dg->components[i].interior.size[1] = patch_size_y; + + dg->components[i].exterior.start[0] = -FD_STENCIL; + dg->components[i].exterior.start[1] = i ? patch_start : -FD_STENCIL; + dg->components[i].exterior.size[0] = gridsize + 2 * FD_STENCIL; + dg->components[i].exterior.size[1] = patch_size_y + ((i == 0) * FD_STENCIL) + ((i == nb_processes - 1) * FD_STENCIL); + + dg->components[i].bnd_is_outer[MG2D_BOUNDARY_0L] = 1; + dg->components[i].bnd_is_outer[MG2D_BOUNDARY_0U] = 1; + dg->components[i].bnd_is_outer[MG2D_BOUNDARY_1L] = i == 0; + dg->components[i].bnd_is_outer[MG2D_BOUNDARY_1U] = i == nb_processes - 1; + } + dc = &dg->components[rank]; + + patch_size_y = (gridsize + nb_processes - 1) / nb_processes; + patch_start = rank * patch_size_y; + patch_end = MIN((rank + 1) * patch_size_y, gridsize); + patch_size_y = patch_end - patch_start; + if (patch_size_y <= 0) { + fprintf(stderr, "Too many processes for grid size %ld: %d\n", gridsize, nb_processes); + ret = 1; + goto fail; + } + + //while (ret == 0) + // sleep(1); + + ctx = mg2d_solver_alloc_mpi(MPI_COMM_WORLD, (size_t [2]){dc->interior.start[0], dc->interior.start[1]}, + dc->interior.size); + if (!ctx) { + fprintf(stderr, "Error allocating the solver context\n"); + return 1; + } + + ctx->step[0] = DOMAIN_SIZE / (gridsize - 1); + ctx->step[1] = DOMAIN_SIZE / (gridsize - 1); + + ctx->fd_stencil = FD_STENCIL; + + ctx->maxiter = MAXITER; + ctx->nb_relax_pre = 2; + ctx->nb_cycles = 1; + ctx->nb_relax_post = 2; + ctx->tol = TOL / (ctx->step[0] * ctx->step[1]); + ctx->nb_threads = 1; + ctx->log_level = MG2D_LOG_INFO; + + for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) { + MG2DBoundary *bnd = ctx->boundaries[bnd_loc]; + const int ci = mg2d_bnd_coord_idx(bnd_loc); + const int bnd_dir = mg2d_bnd_out_dir(bnd_loc); + + double coord[2]; + + if (!dc->bnd_is_outer[bnd_loc]) + continue; + + bnd->type = BC_TYPE; + + memset(bnd->val, 0, dc->interior.size[!ci] * sizeof(*bnd->val)); + + if (bnd->type == MG2D_BC_TYPE_FIXVAL) { + for (int j = 1; j < ctx->fd_stencil; j++) { + double *dst = bnd->val + j * bnd->val_stride; + + coord[ci] = mg2d_bnd_is_upper(bnd_loc) * DOMAIN_SIZE + bnd_dir * j * ctx->step[ci]; + + for (ptrdiff_t k = -j; k < (ptrdiff_t)dc->interior.size[!ci] + j; k++) { + coord[!ci] = (k + dc->interior.start[!ci]) * ctx->step[!ci]; + dst[k] = sol[MG2D_DIFF_COEFF_00](coord[0], coord[1]); + } + } + } + } + + for (size_t y = 0; y < dc->interior.size[1]; y++) { + const double y_coord = (y + dc->interior.start[1]) * ctx->step[1]; + + memset(ctx->u + y * ctx->u_stride, 0, sizeof(*ctx->u) * dc->interior.size[0]); + + for (size_t x = 0; x < dc->interior.size[0]; x++) { + const double x_coord = x * ctx->step[0]; + double rhs = 0.0; + + for (int i = 0; i < MG2D_DIFF_COEFF_NB; i++) { + ctx->diff_coeffs[i][ctx->diff_coeffs_stride * y + x] = pde_coeffs[i]; + rhs += pde_coeffs[i] * sol[i](x_coord, y_coord); + } + + ctx->rhs[y * ctx->rhs_stride + x] = rhs; + } + } + + ret = mg2d_solve(ctx); + if (ret < 0) { + fprintf(stderr, "Error solving the equation\n"); + ret = 1; + goto fail; + } + + mg2d_print_stats(ctx, NULL); + + { + double max_err = 0.0; + + for (size_t y = 0; y < dc->interior.size[1]; y++) { + const double y_coord = (y + dc->interior.start[1]) * ctx->step[1]; + + for (size_t x = 0; x < dc->interior.size[0]; x++) { + const double x_coord = x * ctx->step[0]; + double err = fabs(ctx->u[y * ctx->u_stride + x] - sol[MG2D_DIFF_COEFF_00](x_coord, y_coord)); + if (err > max_err) + max_err = err; + } + } + MPI_Reduce(rank ? &max_err : MPI_IN_PLACE, &max_err, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); + if (rank == 0) { + fprintf(stderr, "max(|solution - exact|): %g\n", max_err); + fprintf(stdout, "%ld %g\n", gridsize, max_err); + } + } + +fail: + mg2d_solver_free(&ctx); + mg2di_dg_free(&dg); + MPI_Finalize(); + return ret; +} -- cgit v1.2.3