From 99cd43b5c824309961a0dae2782704b73b02e881 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 28 Dec 2018 10:15:39 +0100 Subject: Add test programs for relaxation and full multigrid. --- relax_test.c | 209 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 209 insertions(+) create mode 100644 relax_test.c (limited to 'relax_test.c') diff --git a/relax_test.c b/relax_test.c new file mode 100644 index 0000000..b750aba --- /dev/null +++ b/relax_test.c @@ -0,0 +1,209 @@ +#include +#include +#include +#include +#include + +#include "ell_relax.h" +#include "log.h" + +#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 + +static double findmax(double *arr, size_t len) +{ + double ret = 0.0; + for (size_t i = 0; i < len; i++) { + double val = fabs(*arr++); + if (val > ret) + ret = val; + } + return ret; +} + +int main(int argc, char **argv) +{ + EllRelaxContext *ctx; + long int log2N, log2maxiter; + long int N, maxiter; + double res_old, res_new; + int ret = 0; + + if (argc < 3) { + fprintf(stderr, "Usage: %s \n", argv[0]); + return 1; + } + log2N = strtol(argv[1], NULL, 0); + log2maxiter = strtol(argv[2], NULL, 0); + if (log2N <= 0 || log2maxiter <= 0 || + log2N >= sizeof(N) * 8 || log2maxiter >= sizeof(maxiter) * 8) { + fprintf(stderr, "Invalid log2N/log2maxiter values: %ld/%ld\n", + log2N, log2maxiter); + return 1; + } + + N = (1L << log2N) + 1; + maxiter = 1L << log2maxiter; + + ctx = mg2di_ell_relax_alloc((size_t [2]){N, N}); + if (!ctx) { + fprintf(stderr, "Error allocating the solver context\n"); + return 1; + } + + ctx->logger.log = mg2di_log_default_callback; + + ctx->step[0] = 1.0 / (N - 1.0); + ctx->step[1] = 1.0 / (N - 1.0); + + ctx->fd_stencil = 2; + + { + EllRelaxBoundary *bnd = &ctx->boundaries[ELL_RELAX_BOUNDARY_0L]; + + bnd->type = ELL_RELAX_BC_TYPE_FIXVAL; + + memset(bnd->val, 0, N * 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[1] + j; k++) + dst[k] = sol(-j * ctx->step[0], k * ctx->step[1]); + } + } + { + EllRelaxBoundary *bnd = &ctx->boundaries[ELL_RELAX_BOUNDARY_0U]; + + bnd->type = ELL_RELAX_BC_TYPE_FIXVAL; + + memset(bnd->val, 0, N * 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[1] + j; k++) + dst[k] = sol((N - 1 + j) * ctx->step[0], k * ctx->step[1]); + } + } + { + EllRelaxBoundary *bnd = &ctx->boundaries[ELL_RELAX_BOUNDARY_1L]; + + bnd->type = ELL_RELAX_BC_TYPE_FIXVAL; + + memset(bnd->val, 0, N * 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[0] + j; k++) + dst[k] = sol(k * ctx->step[0], -j * ctx->step[1]); + } + } + { + EllRelaxBoundary *bnd = &ctx->boundaries[ELL_RELAX_BOUNDARY_1U]; + + bnd->type = ELL_RELAX_BC_TYPE_FIXVAL; + + memset(bnd->val, 0, N * 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[0] + j; k++) + dst[k] = sol(k * ctx->step[0], (N - 1 + j) * ctx->step[1]); + } + } + + for (size_t y = 0; y < ctx->domain_size[1]; y++) { + const double y_coord = y * ctx->step[1]; + + memset(ctx->u + y * ctx->u_stride, 0, sizeof(*ctx->u) * ctx->domain_size[0]); + + for (size_t x = 0; x < ctx->domain_size[0]; x++) { + const double x_coord = x * ctx->step[0]; + + ctx->diff_coeffs[ELL_RELAX_DIFF_COEFF_02][ctx->diff_coeffs_stride * y + x] = 1.0; + ctx->diff_coeffs[ELL_RELAX_DIFF_COEFF_20][ctx->diff_coeffs_stride * y + x] = 1.0; + ctx->diff_coeffs[ELL_RELAX_DIFF_COEFF_11][ctx->diff_coeffs_stride * y + x] = 1.0; + + 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[ELL_RELAX_DIFF_COEFF_00] + y * ctx->diff_coeffs_stride, 0, + sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size[0]); + memset(ctx->diff_coeffs[ELL_RELAX_DIFF_COEFF_01] + y * ctx->diff_coeffs_stride, 0, + sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size[0]); + memset(ctx->diff_coeffs[ELL_RELAX_DIFF_COEFF_10] + y * ctx->diff_coeffs_stride, 0, + sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size[0]); + } + + ret = mg2di_ell_relax_init(ctx); + if (ret < 0) { + fprintf(stderr, "Error initializing the solver context: %d\n", ret); + ret = 1; + goto fail; + } + + res_old = findmax(ctx->residual, ctx->residual_stride * ctx->domain_size[1]); + + for (int i = 0; i < maxiter; i++) { + ret = mg2di_ell_relax_step(ctx); + if (ret < 0) { + fprintf(stderr, "Error during relaxation\n"); + ret = 1; + goto fail; + } + res_new = findmax(ctx->residual, ctx->residual_stride * ctx->domain_size[1]); + if (res_new > 1e3) { + fprintf(stderr, "Diverged at step %d: %g -> %g\n", i, res_old, res_new); + goto fail; + } + fprintf(stderr, "Step %d: %g -> %g (%g)\n", i, res_old, res_new, res_old / res_new); + res_old = res_new; + } + + { + double max_err = 0.0; + + for (size_t y = 0; y < ctx->domain_size[1]; y++) { + const double y_coord = y * ctx->step[1]; + + for (size_t x = 0; x < ctx->domain_size[0]; 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); + } + +fail: + mg2di_ell_relax_free(&ctx); + return ret; +} -- cgit v1.2.3