mg2d_test: reduce the symmetry of the test problem

Should make for a more robust test.

1 files changed, 45 insertions, 41 deletions

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;
             }