aboutsummaryrefslogtreecommitdiff
path: root/mg2d_mpi_test.c
diff options
context:
space:
mode:
Diffstat (limited to 'mg2d_mpi_test.c')
-rw-r--r--mg2d_mpi_test.c283
1 files changed, 283 insertions, 0 deletions
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 <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <string.h>
+#include <unistd.h>
+
+#include <mpi.h>
+
+#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>\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;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-05 03:02:47 +0000