mg2d_test: add a test for fixed-derivative boundary condition

Also test more terms in the equation.

1 files changed, 72 insertions, 29 deletions

diff --git a/mg2d_test.c b/mg2d_test.c
index 6bdbb00..d1d3d67 100644
--- a/mg2d_test.c
+++ b/mg2d_test.c
@@ -36,6 +36,33 @@ static double sol_dxy(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)
+{
+    return cos(M_PI * x) * cos(M_PI * y);
+}
+static double sol_dx(double x, double y)
+{
+    return -M_PI * sin(M_PI * x) * cos(M_PI * y);
+}
+static double sol_dy(double x, double y)
+{
+    return -M_PI * cos(M_PI * x) * sin(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 * sin(M_PI * x) * sin(M_PI * y);
+}
+#define BC_TYPE MG2D_BC_TYPE_FIXDIFF
 #endif
 
 int main(int argc, char **argv)
@@ -76,57 +103,65 @@ int main(int argc, char **argv)
     {
         MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0L];
 
-        bnd->type = MG2D_BC_TYPE_FIXVAL;
+        bnd->type = BC_TYPE;
 
         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;
+        if (bnd->type == MG2D_BC_TYPE_FIXVAL) {
+            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]);
+                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;
+        bnd->type = BC_TYPE;
 
         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;
+        if (bnd->type == MG2D_BC_TYPE_FIXVAL) {
+            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]);
+                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;
+        bnd->type = BC_TYPE;
 
         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;
+        if (bnd->type == MG2D_BC_TYPE_FIXVAL) {
+            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]);
+                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;
+        bnd->type = BC_TYPE;
 
         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;
+        if (bnd->type == MG2D_BC_TYPE_FIXVAL) {
+            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 (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]);
+            }
         }
     }
 
@@ -135,24 +170,32 @@ 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[0]);
+        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];
 
             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[ELL_RELAX_DIFF_COEFF_00][ctx->diff_coeffs_stride * y + x] = 2.0 * M_PI;
+            ctx->diff_coeffs[MG2D_DIFF_COEFF_00][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);
+            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);
         }
-        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);