Implement the multigrid scheme.

1 files changed, 244 insertions, 0 deletions

diff --git a/mg2d.h b/mg2d.h
new file mode 100644
index 0000000..3f9cbe7
--- /dev/null
+++ b/mg2d.h
@@ -0,0 +1,244 @@
+/*
+ * Multigrid solver for a 2nd order 2D linear PDE.
+ * Copyright 2018 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef MG2D_H
+#define MG2D_H
+
+#include <stddef.h>
+#include <stdint.h>
+
+/**
+ * Type of the boundary condition on a given boundary.
+ */
+enum MG2DBCType {
+    /**
+     * The value of the unknown function is fixed on the boundary.
+     */
+    MG2D_BC_TYPE_FIXVAL,
+    /**
+     * The normal derivative of the unkown function is fixed on the boundary.
+     *
+     * TODO: the only supported value of the derivative is currently zero, i.e.
+     * periodic boundary conditions.
+     */
+    MG2D_BC_TYPE_FIXDIFF,
+};
+
+/**
+ * Location of the boundary.
+ */
+enum MG2DBoundaryLoc {
+    /**
+     * Lower coord0, lower coord1.
+     */
+    MG2D_BOUNDARY_00,
+    /**
+     * Lower coord0, upper coord1.
+     */
+    MG2D_BOUNDARY_01,
+    /**
+     * Upper coord0, lower coord1.
+     */
+    MG2D_BOUNDARY_10,
+    /**
+     * Upper coord0, upper coord1.
+     */
+    MG2D_BOUNDARY_11,
+};
+
+/**
+ * Derivative specifier.
+ */
+enum MG2DDiffCoeff {
+    /**
+     * Zeroth derivative, i.e. function value.
+     */
+    MG2D_DIFF_COEFF_00,
+    /**
+     * First derivative wrt coord1.
+     */
+    MG2D_DIFF_COEFF_01,
+    /**
+     * First derivative wrt coord0.
+     */
+    MG2D_DIFF_COEFF_10,
+    /**
+     * Second derivative, once wrt coord0 and once wrt coord1
+     */
+    MG2D_DIFF_COEFF_11,
+    /**
+     * Second derivative wrt coord1
+     */
+    MG2D_DIFF_COEFF_02,
+    /**
+     * Second derivative wrt coord0
+     */
+    MG2D_DIFF_COEFF_20,
+    /**
+     * Total number of allowed derivatives.
+     */
+    MG2D_DIFF_COEFF_NB,
+};
+
+typedef struct MG2DInternal MG2DInternal;
+
+/**
+ * Boundary condition definition.
+ */
+typedef struct MG2DBoundary {
+    /**
+     * Type of the boundary condition.
+     */
+    enum MG2DBCType type;
+    /**
+     * For type = MG2D_BC_TYPE_FIXVAL:
+     *     Values of the unknown function on the boundary.
+     *
+     * Ignored otherwise.
+     */
+    double *val;
+    /**
+     * Number of elements in val.
+     */
+    size_t  val_len;
+} MG2DBoundary;
+
+typedef struct MG2DContext {
+    /**
+     * Solver private data, not to be accessed in any way by the caller.
+     */
+    MG2DInternal *priv;
+
+    /**
+     * A callback that will be used to print diagnostic messages.
+     *
+     * Defaults to fprintf(stderr, ...)
+     */
+    void (*log_callback)(const struct MG2DContext *ctx, int level,
+                         const char *fmt, va_list);
+
+    /**
+     * Arbitrary user data, e.g. to be used by the log callback.
+     */
+    void *opaque;
+
+    /**
+     * Number of grid points in each direction, set by mg2d_solver_alloc().
+     * Read-only for the caller.
+     */
+    size_t domain_size;
+
+    /**
+     * Distance between the neighbouring grid points.
+     * Must be set by the caller.
+     */
+    double step[2];
+
+    /**
+     * Order of the finite difference operators used for approximating derivatives.
+     * Must be set by the caller.
+     */
+    size_t fd_stencil;
+
+    /**
+     * Boundary specification, indexed by MG2DBoundaryLoc.
+     * To be filled by the caller.
+     */
+    MG2DBoundary *boundaries[4];
+
+    /**
+     * Maximum number of solver iterations.
+     */
+    unsigned int maxiter;
+
+    /**
+     * Value of the residual at which to consider the equation solved.
+     */
+    double tol;
+
+    /**
+     * Number of cycles on each subgrid.
+     */
+    unsigned int nb_cycles;
+
+    /**
+     * Number of relaxation steps before and after each coarser-grid solve.
+     */
+    unsigned int nb_relax_pre;
+    unsigned int nb_relax_post;
+
+    /**
+     * Values of the unknown function.
+     *
+     * Allocated and initialized to zero by the solver in mg2d_solver_alloc(),
+     * owned by the solver.
+     * May be filled by the caller before solving to set the initial guess.
+     * Afterwards updated in mg2d_solve().
+     */
+    double   *u;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t u_stride;
+
+    /**
+     * Values of the right-hand side.
+     *
+     * Allocated by the solver in mg2d_solver_alloc(), owned by the solver.
+     * Must be filled by the caller before solving.
+     */
+    double    *rhs;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t  rhs_stride;
+
+    /**
+     * Coefficients C_{*} that define the differential equation.
+     *
+     * Allocated by the solver in mg2d_solver_alloc(), owned by the solver.
+     * Must be filled by the caller before solving.
+     */
+    double   *diff_coeffs[MG2D_DIFF_COEFF_NB];
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t diff_coeffs_stride;
+} MG2DContext;
+
+/**
+ * Allocate the solver for the given domain size.
+ *
+ * @return The solver context on success, NULL on failure.
+ */
+MG2DContext *mg2d_solver_alloc(size_t log2_gridsize);
+/**
+ * Solve the equation, after all the required fields have been filled by the
+ * caller.
+ *
+ * This function may be called more than once.
+ *
+ * @return 0 on success, a negative error code on failure.
+ */
+int mg2d_solve(MG2DContext *ctx);
+/**
+ * Free the solver and write NULL to the provided pointer.
+ */
+void mg2d_solver_free(MG2DContext **ctx);
+
+#endif /* MG2D_H */