Add test programs for relaxation and full multigrid.

1 files changed, 184 insertions, 0 deletions

diff --git a/mg2d_test.c b/mg2d_test.c
new file mode 100644
index 0000000..9b55963
--- /dev/null
+++ b/mg2d_test.c
@@ -0,0 +1,184 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <string.h>
+
+#include "mg2d.h"
+
+#define MAXITER 64
+#define TOL 1e-9
+
+#if 1
+static double sol(double x, double y)
+{
+    return sin(M_PI * x) * sin(M_PI * y);
+}
+static double sol_dx(double x, double y)
+{
+    return M_PI * cos(M_PI * x) * sin(M_PI * y);
+}
+static double sol_dy(double x, double y)
+{
+    return M_PI * sin(M_PI * x) * cos(M_PI * y);
+}
+static double sol_dxx(double x, double y)
+{
+    return -M_PI * M_PI * sol(x, y);
+}
+static double sol_dyy(double x, double y)
+{
+    return -M_PI * M_PI * sol(x, y);
+}
+static double sol_dxy(double x, double y)
+{
+    return M_PI * M_PI * cos(M_PI * x) * cos(M_PI * y);
+}
+#endif
+
+int main(int argc, char **argv)
+{
+    MG2DContext *ctx;
+    long int gridsize;
+    int ret = 0;
+
+    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;
+    }
+
+    ctx = mg2d_solver_alloc(gridsize);
+    if (!ctx) {
+        fprintf(stderr, "Error allocating the solver context\n");
+        return 1;
+    }
+
+    ctx->step[0] = 1.0 / (gridsize - 1);
+    ctx->step[1] = 1.0 / (gridsize - 1);
+
+    ctx->fd_stencil = 2;
+
+    ctx->maxiter = MAXITER;
+    ctx->nb_relax_pre  = 2;
+    ctx->nb_cycles     = 1;
+    ctx->nb_relax_post = 2;
+    ctx->tol = TOL;
+    ctx->nb_threads = 1;
+
+    {
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0L];
+
+        bnd->type = MG2D_BC_TYPE_FIXVAL;
+
+        memset(bnd->val, 0, gridsize * sizeof(*bnd->val));
+
+        for (int j = 1; j < ctx->fd_stencil; j++) {
+            double *dst = bnd->val + j * bnd->val_stride;
+
+            for (ptrdiff_t k = -j; k < (ptrdiff_t)ctx->domain_size + j; k++)
+                dst[k] = sol(-j * ctx->step[0], k * ctx->step[1]);
+        }
+    }
+    {
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0U];
+
+        bnd->type = MG2D_BC_TYPE_FIXVAL;
+
+        memset(bnd->val, 0, gridsize * sizeof(*bnd->val));
+
+        for (int j = 1; j < ctx->fd_stencil; j++) {
+            double *dst = bnd->val + j * bnd->val_stride;
+
+            for (ptrdiff_t k = -j; k < (ptrdiff_t)ctx->domain_size + j; k++)
+                dst[k] = sol((gridsize - 1 + j) * ctx->step[0], k * ctx->step[1]);
+        }
+    }
+    {
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1L];
+
+        bnd->type = MG2D_BC_TYPE_FIXVAL;
+
+        memset(bnd->val, 0, gridsize * sizeof(*bnd->val));
+
+        for (int j = 1; j < ctx->fd_stencil; j++) {
+            double *dst = bnd->val + j * bnd->val_stride;
+
+            for (ptrdiff_t k = -j; k < (ptrdiff_t)ctx->domain_size + j; k++)
+                dst[k] = sol(k * ctx->step[0], -j * ctx->step[1]);
+        }
+    }
+    {
+        MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1U];
+
+        bnd->type = MG2D_BC_TYPE_FIXVAL;
+
+        memset(bnd->val, 0, gridsize * sizeof(*bnd->val));
+
+        for (int j = 1; j < ctx->fd_stencil; j++) {
+            double *dst = bnd->val + j * bnd->val_stride;
+
+            for (ptrdiff_t k = -j; k < (ptrdiff_t)ctx->domain_size + j; k++)
+                dst[k] = sol(k * ctx->step[0], (gridsize - 1 + j) * ctx->step[1]);
+        }
+    }
+
+
+    for (size_t y = 0; y < ctx->domain_size; y++) {
+        const double y_coord = y * ctx->step[1];
+
+        memset(ctx->u + y * ctx->u_stride, 0, sizeof(*ctx->u) * ctx->domain_size);
+        //memset(ctx->rhs + y * ctx->rhs_stride, 0, sizeof(*ctx->rhs) * ctx->domain_size[0]);
+
+        for (size_t x = 0; x < ctx->domain_size; x++) {
+            const double x_coord = x * ctx->step[0];
+
+            ctx->diff_coeffs[MG2D_DIFF_COEFF_02][ctx->diff_coeffs_stride * y + x] = 1.0;
+            ctx->diff_coeffs[MG2D_DIFF_COEFF_20][ctx->diff_coeffs_stride * y + x] = 1.0;
+            ctx->diff_coeffs[MG2D_DIFF_COEFF_11][ctx->diff_coeffs_stride * y + x] = 1.0;
+            //ctx->diff_coeffs[ELL_RELAX_DIFF_COEFF_00][ctx->diff_coeffs_stride * y + x] = 2.0 * M_PI;
+
+            ctx->rhs[y * ctx->rhs_stride + x] = sol_dxx(x_coord, y_coord) + sol_dyy(x_coord, y_coord) + sol_dxy(x_coord, y_coord);
+        }
+        memset(ctx->diff_coeffs[MG2D_DIFF_COEFF_00] + y * ctx->diff_coeffs_stride, 0,
+               sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size);
+        memset(ctx->diff_coeffs[MG2D_DIFF_COEFF_01] + y * ctx->diff_coeffs_stride, 0,
+               sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size);
+        memset(ctx->diff_coeffs[MG2D_DIFF_COEFF_10] + y * ctx->diff_coeffs_stride, 0,
+               sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size);
+    }
+
+    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 < ctx->domain_size; y++) {
+            const double y_coord = y * ctx->step[1];
+
+            for (size_t x = 0; x < ctx->domain_size; x++) {
+                const double x_coord = x * ctx->step[0];
+                double err = fabs(ctx->u[y * ctx->u_stride + x] - sol(x_coord, y_coord));
+                if (err > max_err)
+                    max_err = err;
+            }
+        }
+        fprintf(stderr, "max(|solution - exact|): %g\n", max_err);
+        fprintf(stdout, "%ld %g\n", gridsize, max_err);
+    }
+
+fail:
+    mg2d_solver_free(&ctx);
+    return ret;
+}