ell_relax -> ell_grid_solve

Generalize the API to allow for multiple solver types. This is done in
preparation for the exact linear system inversion solver.

1 files changed, 218 insertions, 0 deletions

diff --git a/ell_grid_solve.h b/ell_grid_solve.h
new file mode 100644
index 0000000..70782c8
--- /dev/null
+++ b/ell_grid_solve.h
@@ -0,0 +1,218 @@
+/*
+ * Solver for a 2nd order 2D linear PDE using finite differencing on a grid.
+ * 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_ELL_GRID_SOLVE_H
+#define MG2D_ELL_GRID_SOLVE_H
+
+/**
+ * The problem being solved is a linear partial differential
+ * equation
+ *
+ *   ∑ C_{ab} ∂_a　∂_b u  + ∑ C_{a} ∂_a u + C u = rhs
+ *  a,b                     a
+ *
+ *  where
+ *      * a and b identify spatial directions and run from 0 to 1
+ *      * u = u(x_a) is the unknown function
+ *      * C_{ab}, C_{a} and C are the coefficients in front of
+ *        the corresponding derivative of unknown function
+ *      * rhs is the right-hand side of the equation
+ * C_{*} and rhs are all (known) functions of space and define the
+ * equation to be solved.
+ */
+
+#include <stddef.h>
+#include <stdint.h>
+#include <threadpool.h>
+
+#include "log.h"
+#include "mg2d_boundary.h"
+#include "mg2d_constants.h"
+
+enum EGSType {
+    /**
+     * Solve the equation using relaxation.
+     *
+     * solver_data is EGSRelaxContext
+     * mg2di_egs_solve() does a single relaxation step
+     */
+    EGS_SOLVER_RELAXATION,
+};
+
+typedef struct EGSInternal EGSInternal;
+
+typedef struct EGSRelaxContext {
+    /**
+     * The time stepping factor in relaxation.
+     */
+    double relax_factor;
+
+    int64_t count_correct;
+    int64_t time_correct;
+} EGSRelaxContext;
+
+typedef struct EGSContext {
+    enum EGSType solver_type;
+
+    /**
+     * Solver type-specific data.
+     *
+     * Should be cast into the struct specified in documentation for this type.
+     */
+    void *solver_data;
+    /**
+     * Solver private data, not to be accessed in any way by the caller.
+     */
+    EGSInternal *priv;
+
+    /**
+     * The logging context, to be filled by the caller before mg2di_egs_init().
+     */
+    MG2DLogger logger;
+
+    /**
+     * The thread pool used for execution. May be set by the caller before
+     * mg2di_egs_init().
+     */
+    TPContext *tp;
+
+    /**
+     * Flags indicating supported CPU features.
+     */
+    int cpuflags;
+
+    /**
+     * Size of the solver grid, set by mg2di_egs_alloc().
+     * Read-only for the caller.
+     */
+    size_t domain_size[2];
+
+    /**
+     * Distance between the neighbouring grid points.
+     * Must be set by the caller before mg2di_egs_init().
+     */
+    double step[2];
+
+    /**
+     * Order of the finite difference operators used for approximating derivatives.
+     * Must be set by the caller before mg2di_egs_init().
+     */
+    size_t fd_stencil;
+
+    /**
+     * Boundary specification, indexed by MG2DBoundaryLoc.
+     * To be filled by the caller before mg2di_egs_init().
+     */
+    MG2DBoundary *boundaries[4];
+
+    /**
+     * Values of the unknown function.
+     *
+     * Allocated by the solver in mg2di_egs_alloc(), owned by the solver.
+     * Must be filled by the caller before mg2di_egs_init() to set the
+     * initial guess.
+     * Afterwards updated in mg2di_egs_step().
+     */
+    double   *u;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t u_stride;
+
+    /**
+     * Values of the right-hand side.
+     *
+     * Allocated by the solver in mg2di_egs_alloc(), owned by the solver.
+     * Must be filled by the caller before mg2di_egs_init().
+     */
+    double    *rhs;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t  rhs_stride;
+
+    /**
+     * Values of the right-hand side.
+     *
+     * Allocated by the solver in mg2di_egs_alloc(), owned by the solver.
+     * Read-only for the caller. Initialized after mg2di_egs_init(),
+     * afterwards updated in mg2di_egs_step().
+     */
+    double   *residual;
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t residual_stride;
+
+    /**
+     * Maximum of the absolute value of residual.
+     */
+    double residual_max;
+
+    /**
+     * Coefficients C_{*} that define the differential equation.
+     *
+     * Allocated by the solver in mg2di_egs_alloc(), owned by the solver.
+     * Must be filled by the caller before mg2di_egs_init().
+     */
+    double   *diff_coeffs[MG2D_DIFF_COEFF_NB];
+    /**
+     * Distance between neighbouring gridpoints along coord1.
+     */
+    ptrdiff_t diff_coeffs_stride;
+
+    /* timings */
+    int64_t time_boundaries;
+    int64_t count_boundaries;
+    int64_t time_res_calc;
+    int64_t count_res;
+} EGSContext;
+
+/**
+ * Allocate the solver for the given domain size.
+ *
+ * @param domain_size number of grid points in each direction.
+ *
+ * @return The solver context on success, NULL on failure.
+ */
+EGSContext *mg2di_egs_alloc(enum EGSType solver_type, size_t domain_size[2]);
+/**
+ * Initialize the solver for use, after all the required fields are filled by
+ * the caller.
+ *
+ * This function may be called multiple times to re-initialize the solver after
+ * certain parameters (e.g. the right-hand side or the equation coefficients)
+ * change.
+ *
+ * @return 0 on success, a negative error code on failure.
+ */
+int              mg2di_egs_init(EGSContext *ctx);
+/**
+ * Free the solver and write NULL to the provided pointer.
+ */
+void             mg2di_egs_free(EGSContext **ctx);
+
+/**
+ * Solve the equation defined by the data filled in the context. Precise
+ * semantics depends on the solver type.
+ *
+ * @return 0 on success, a negative error code on failure.
+ */
+int mg2di_egs_solve(EGSContext *ctx);
+
+#endif /* MG2D_ELL_GRID_SOLVE_H */