From 8633d53539cdcab72607b755de840a20e674d8c2 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 1 Feb 2019 19:39:27 +0100 Subject: mg2d_test: reduce the symmetry of the test problem Should make for a more robust test. --- mg2d_test.c | 86 ++++++++++++++++++++++++++++++++----------------------------- 1 file changed, 45 insertions(+), 41 deletions(-) (limited to 'mg2d_test.c') diff --git a/mg2d_test.c b/mg2d_test.c index d1d3d67..9c3ac17 100644 --- a/mg2d_test.c +++ b/mg2d_test.c @@ -11,60 +11,78 @@ #define MAXITER 64 #define TOL 1e-9 +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(double x, double y) +static double sol_00(double x, double y) { return sin(M_PI * x) * sin(M_PI * y); } -static double sol_dx(double x, double y) +static double sol_10(double x, double y) { return M_PI * cos(M_PI * x) * sin(M_PI * y); } -static double sol_dy(double x, double y) +static double sol_01(double x, double y) { return M_PI * sin(M_PI * x) * cos(M_PI * y); } -static double sol_dxx(double x, double y) +static double sol_20(double x, double y) { - return -M_PI * M_PI * sol(x, y); + return -M_PI * M_PI * sol_00(x, y); } -static double sol_dyy(double x, double y) +static double sol_02(double x, double y) { - return -M_PI * M_PI * sol(x, y); + return -M_PI * M_PI * sol_00(x, y); } -static double sol_dxy(double x, double 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(double x, double y) +static double sol_00(double x, double y) { return cos(M_PI * x) * cos(M_PI * y); } -static double sol_dx(double x, double y) +static double sol_10(double x, double y) { return -M_PI * sin(M_PI * x) * cos(M_PI * y); } -static double sol_dy(double x, double y) +static double sol_01(double x, double y) { return -M_PI * cos(M_PI * x) * sin(M_PI * y); } -static double sol_dxx(double x, double y) +static double sol_20(double x, double y) { - return -M_PI * M_PI * sol(x, y); + return -M_PI * M_PI * sol_00(x, y); } -static double sol_dyy(double x, double y) +static double sol_02(double x, double y) { - return -M_PI * M_PI * sol(x, y); + return -M_PI * M_PI * sol_00(x, y); } -static double sol_dxy(double x, double 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_FIXDIFF #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; @@ -112,7 +130,7 @@ int main(int argc, char **argv) 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]); + dst[k] = sol[MG2D_DIFF_COEFF_00](-j * ctx->step[0], k * ctx->step[1]); } } } @@ -128,7 +146,7 @@ int main(int argc, char **argv) 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]); + dst[k] = sol[MG2D_DIFF_COEFF_00]((gridsize - 1 + j) * ctx->step[0], k * ctx->step[1]); } } } @@ -144,7 +162,7 @@ int main(int argc, char **argv) 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]); + dst[k] = sol[MG2D_DIFF_COEFF_00](k * ctx->step[0], -j * ctx->step[1]); } } } @@ -160,7 +178,7 @@ int main(int argc, char **argv) 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]); + dst[k] = sol[MG2D_DIFF_COEFF_00](k * ctx->step[0], (gridsize - 1 + j) * ctx->step[1]); } } } @@ -170,31 +188,17 @@ int main(int argc, char **argv) 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); - 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); - memset(ctx->diff_coeffs[MG2D_DIFF_COEFF_20] + y * ctx->diff_coeffs_stride, 0, - sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size); - memset(ctx->diff_coeffs[MG2D_DIFF_COEFF_02] + y * ctx->diff_coeffs_stride, 0, - sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size); - memset(ctx->diff_coeffs[MG2D_DIFF_COEFF_11] + y * ctx->diff_coeffs_stride, 0, - sizeof(*ctx->diff_coeffs[0]) * ctx->domain_size); for (size_t x = 0; x < ctx->domain_size; x++) { const double x_coord = x * ctx->step[0]; + double rhs = 0.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_01][ctx->diff_coeffs_stride * y + x] = 1.0; - ctx->diff_coeffs[MG2D_DIFF_COEFF_10][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[MG2D_DIFF_COEFF_00][ctx->diff_coeffs_stride * y + x] = 1.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] = sol_dxx(x_coord, y_coord) + sol_dyy(x_coord, y_coord) + sol(x_coord, y_coord) + sol_dxy(x_coord, y_coord) + sol_dx(x_coord, y_coord) + sol_dy(x_coord, y_coord); + ctx->rhs[y * ctx->rhs_stride + x] = rhs; } } @@ -215,7 +219,7 @@ int main(int argc, char **argv) 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)); + 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; } -- cgit v1.2.3