From d519d82a3e4b32944b77b1ae26cfefa45ec29d71 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 26 Jan 2019 10:19:30 +0100 Subject: 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. --- ell_grid_solve.h | 218 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 ell_grid_solve.h (limited to 'ell_grid_solve.h') 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 + * + * 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_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 +#include +#include + +#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 */ -- cgit v1.2.3