From 3d873841dfca19b632525b1db7a33e0f9750dc06 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Wed, 1 Aug 2018 17:12:44 +0200 Subject: Implement the multigrid scheme. --- mg2d.h | 244 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 mg2d.h (limited to 'mg2d.h') 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 + * + * 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 . + */ + +#ifndef MG2D_H +#define MG2D_H + +#include +#include + +/** + * 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 */ -- cgit v1.2.3