/* * 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 . */ #include "config.h" #include #include #include #include #include #include #include #include #include #include #include "bicgstab.h" #include "common.h" #include "components.h" #include "cpu.h" #include "ell_grid_solve.h" #include "log.h" #include "mg2d_boundary.h" #include "mg2d_boundary_internal.h" #include "mg2d_constants.h" #include "residual_calc.h" #include "timer.h" static const double relax_factors[FD_STENCIL_MAX] = { [0] = 1.0 / 5, [1] = 1.0 / 7, }; typedef struct EGSRelaxInternal { void (*line_add)(ptrdiff_t, double *, const double *, double); double relax_factor; } EGSRelaxInternal; typedef struct EGSExactInternal { void (*fill_mat)(double *mat_row, ptrdiff_t mat_stride, ptrdiff_t row_stride, NDArray **diff_coeffs, ptrdiff_t idx_src); size_t N; size_t N_ghosts; int arrays_alloced; double *mat; double *mat_f; double *rhs; double *rhs_mat; double *rhs_factor; double *x; int *ipiv; double *scratch_lines; ptrdiff_t scratch_lines_allocated; int mat_filled; BiCGStabContext *bicgstab; } EGSExactInternal; struct EGSInternal { NDArray *rhs_base; NDArray *residual_base; NDArray *u_base; NDArray *u_next_base; NDArray *u_next_exterior; NDArray *u_next; int u_next_valid; /* all the diff coeffs concatenated */ NDArray *diff_coeffs_public_base; NDArray *diff_coeffs_base; NDArray *diff_coeffs[MG2D_DIFF_COEFF_NB]; ptrdiff_t (*residual_calc_offset)[2]; size_t (*residual_calc_size)[2]; ptrdiff_t bnd_calc_offset[2]; size_t bnd_calc_size[2]; int bnd_zero[4]; int reflect_skip; ResidualCalcContext *rescalc; EGSRelaxInternal r; EGSExactInternal e; TPContext *tp_internal; MPI_Comm comm; DomainGeometry *dg; unsigned int local_component; size_t fd_stencil; unsigned int steps_per_sync; unsigned int steps_since_sync; int *sync_sendcounts; int *sync_senddispl; MPI_Datatype *sync_sendtypes; int *sync_recvcounts; int *sync_recvdispl; MPI_Datatype *sync_recvtypes; }; static const double fd_denoms[][MG2D_DIFF_COEFF_NB] = { { [MG2D_DIFF_COEFF_00] = 1.0, [MG2D_DIFF_COEFF_10] = 2.0, [MG2D_DIFF_COEFF_01] = 2.0, [MG2D_DIFF_COEFF_20] = 1.0, [MG2D_DIFF_COEFF_02] = 1.0, [MG2D_DIFF_COEFF_11] = 4.0, }, { [MG2D_DIFF_COEFF_00] = 1.0, [MG2D_DIFF_COEFF_10] = 12.0, [MG2D_DIFF_COEFF_01] = 12.0, [MG2D_DIFF_COEFF_20] = 12.0, [MG2D_DIFF_COEFF_02] = 12.0, [MG2D_DIFF_COEFF_11] = 144.0, }, }; static void boundaries_sync(EGSContext *ctx, NDArray *a_dst) { EGSInternal *priv = ctx->priv; if (priv->dg->nb_components > 1) { mg2di_timer_start(&ctx->timer_mpi_sync); MPI_Alltoallw(a_dst->data, priv->sync_sendcounts, priv->sync_senddispl, priv->sync_sendtypes, a_dst->data, priv->sync_recvcounts, priv->sync_recvdispl, priv->sync_recvtypes, priv->comm); mg2di_timer_stop(&ctx->timer_mpi_sync); } } static void residual_calc(EGSContext *ctx, int export_res) { EGSInternal *priv = ctx->priv; ptrdiff_t *offset = priv->residual_calc_offset[priv->steps_since_sync]; size_t *size = priv->residual_calc_size[priv->steps_since_sync]; const double *diff_coeffs0 = NDA_PTR2D(priv->diff_coeffs[0], offset[0], offset[1]); const double *diff_coeffs1 = NDA_PTR2D(priv->diff_coeffs[1], offset[0], offset[1]); NDArray *dst = export_res ? ctx->residual : priv->u_next; int reflect_flags = 0; mg2di_timer_start(&ctx->timer_res_calc); for (int bnd_idx = 0; bnd_idx < 4; bnd_idx++) if (ctx->boundaries[bnd_idx]->type == MG2D_BC_TYPE_REFLECT && priv->dg->components[priv->local_component].bnd_is_outer[bnd_idx]) reflect_flags |= 1 << bnd_idx; mg2di_residual_calc(priv->rescalc, size, &ctx->residual_max, NDA_PTR2D(dst, offset[0], offset[1]), dst->stride[0], NDA_PTR2D(ctx->u, offset[0], offset[1]), ctx->u->stride[0], NDA_PTR2D(ctx->rhs, offset[0], offset[1]), ctx->rhs->stride[0], diff_coeffs0, priv->diff_coeffs[0]->stride[0], diff_coeffs1 - diff_coeffs0, export_res ? 0.0 : 1.0, export_res ? 1.0 : priv->r.relax_factor, reflect_flags, FD_STENCIL_MAX); mg2di_timer_stop(&ctx->timer_res_calc); if (priv->dg->nb_components > 1) { mg2di_timer_start(&ctx->timer_mpi_sync); MPI_Allreduce(MPI_IN_PLACE, &ctx->residual_max, 1, MPI_DOUBLE, MPI_MAX, priv->comm); mg2di_timer_stop(&ctx->timer_mpi_sync); } priv->reflect_skip = ~0; } static void boundaries_apply_fixval(double *dst, const ptrdiff_t dst_stride[2], const double *src, ptrdiff_t src_stride, size_t boundary_size) { for (int j = 0; j < FD_STENCIL_MAX; j++) { for (ptrdiff_t i = -j; i < (ptrdiff_t)boundary_size + j; i++) dst[i * dst_stride[0]] = src[i]; dst += dst_stride[1]; src += src_stride; } } static void boundaries_apply_reflect(double *dst, const ptrdiff_t dst_stride[2], size_t boundary_size) { if (dst_stride[0] == 1) { for (int j = 1; j <= FD_STENCIL_MAX; j++) memcpy(dst + j * dst_stride[1], dst - j * dst_stride[1], sizeof(*dst) * boundary_size); } else { for (size_t i = 0; i < boundary_size; i++) { for (int j = 1; j <= FD_STENCIL_MAX; j++) dst[dst_stride[1] * j] = dst[-dst_stride[1] * j]; dst += dst_stride[0]; } } } /* FIXME we currently use 1st-order forward differences in all cases, * since centered and/or higher-order FDs seem to be unstable * this may be suboptimal and should be investigated closer */ static const double falloff_bnd_fd_factors[][FD_STENCIL_MAX * 2 + 1] = { { 1.0 / 2.0, -1.0 / 2.0, }, { 1.0 / 2.0, -1.0 / 2.0, }, //{ -1.0, 8.0, 0.0, -8.0, 1.0 }, }; static void boundaries_apply_falloff(double *dst, const ptrdiff_t dst_stride[2], size_t boundary_size, EGSContext *ctx, int ci) { const ptrdiff_t offset = ctx->priv->dg->components[ctx->priv->local_component].interior.start[!ci]; for (ptrdiff_t bnd_layer = 1; bnd_layer <= ctx->fd_stencil; bnd_layer++) for (ptrdiff_t bnd_idx = -(ptrdiff_t)ctx->fd_stencil; bnd_idx < (ptrdiff_t)(boundary_size + ctx->fd_stencil); bnd_idx++) { const double x = ctx->step[0] * (ctx->priv->dg->domain_size[ci] - 1 + bnd_layer - ctx->fd_stencil); const double y = ctx->step[0] * (bnd_idx + offset); const double r = sqrt(SQR(x) + SQR(y)); double *row = dst + bnd_layer * dst_stride[1] + bnd_idx * dst_stride[0]; double val = row[-(ptrdiff_t)ctx->fd_stencil * dst_stride[1]] / r * (x / r); for (int i = 1; i < ctx->fd_stencil * 2 + 1; i++) val += dst[(bnd_layer - i) * dst_stride[1] + bnd_idx * dst_stride[0]] * falloff_bnd_fd_factors[ctx->fd_stencil - 1][i] / ctx->step[0]; *row = -val * (ctx->step[0] / falloff_bnd_fd_factors[ctx->fd_stencil - 1][0]); } } static void boundaries_apply(EGSContext *ctx, NDArray *a_dst, int init) { static const enum MG2DBCType bnd_type_order[] = { MG2D_BC_TYPE_FIXVAL, MG2D_BC_TYPE_REFLECT, MG2D_BC_TYPE_FALLOFF }; EGSInternal *priv = ctx->priv; const ptrdiff_t strides[2] = { 1, a_dst->stride[0] }; mg2di_timer_start(&ctx->timer_bnd); for (int order_idx = 0; order_idx < ARRAY_ELEMS(bnd_type_order); order_idx++) { for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) { MG2DBoundary *bnd = ctx->boundaries[i]; const int ci = mg2d_bnd_coord_idx(i); const ptrdiff_t stride_boundary = strides[!ci]; const ptrdiff_t stride_offset = strides[ci]; const size_t size_boundary = ctx->domain_size[!ci]; //const size_t size_boundary = priv->bnd_calc_size[!ci]; const size_t size_offset = ctx->domain_size[ci]; double *dst = a_dst->data + mg2d_bnd_is_upper(i) * ((size_offset - 1) * stride_offset);// + priv->bnd_calc_offset[!ci] * stride_boundary; const ptrdiff_t dst_strides[] = { stride_boundary, mg2d_bnd_out_dir(i) * stride_offset }; if (bnd->type != bnd_type_order[order_idx] || !priv->dg->components[priv->local_component].bnd_is_outer[i]) continue; switch (bnd->type) { case MG2D_BC_TYPE_FIXVAL: if (init && !priv->bnd_zero[i]) { mg2di_timer_start(&ctx->timer_bnd_fixval); boundaries_apply_fixval(dst, dst_strides, bnd->val, bnd->val_stride, size_boundary); mg2di_timer_stop(&ctx->timer_bnd_fixval); } break; case MG2D_BC_TYPE_REFLECT: if (!(priv->reflect_skip & (1 << i))) { mg2di_timer_start(&ctx->timer_bnd_reflect); boundaries_apply_reflect(dst, dst_strides, size_boundary); mg2di_timer_stop(&ctx->timer_bnd_reflect); } break; case MG2D_BC_TYPE_FALLOFF: mg2di_timer_start(&ctx->timer_bnd_falloff); boundaries_apply_falloff(dst, dst_strides, size_boundary, ctx, ci); mg2di_timer_stop(&ctx->timer_bnd_falloff); break; } } } /* fill in the corner ghosts */ mg2di_timer_start(&ctx->timer_bnd_corners); for (int pos_y = 0; pos_y < 2; pos_y++) { enum MG2DBoundaryLoc loc_y = mg2d_bnd_id(1, pos_y); MG2DBoundary *bnd_y = ctx->boundaries[loc_y]; const int dir_y = mg2d_bnd_out_dir(loc_y); if (!priv->dg->components[priv->local_component].bnd_is_outer[loc_y]) continue; for (int pos_x = 0; pos_x < 2; pos_x++) { enum MG2DBoundaryLoc loc_x = mg2d_bnd_id(0, pos_x); MG2DBoundary *bnd_x = ctx->boundaries[loc_x]; const int dir_x = mg2d_bnd_out_dir(loc_x); int fact_x = dir_x, fact_y = dir_y; double *dst = a_dst->data + mg2d_bnd_is_upper(loc_y) * ((ctx->domain_size[1] - 1) * a_dst->stride[0]) + mg2d_bnd_is_upper(loc_x) * (ctx->domain_size[0] - 1); if (!priv->dg->components[priv->local_component].bnd_is_outer[loc_x]) continue; if (bnd_y->type == MG2D_BC_TYPE_REFLECT) fact_y *= -1; else if (bnd_x->type == MG2D_BC_TYPE_REFLECT) fact_x *= -1; else continue; dst = a_dst->data + mg2d_bnd_is_upper(loc_y) * ((ctx->domain_size[1] - 1) * a_dst->stride[0]) + mg2d_bnd_is_upper(loc_x) * (ctx->domain_size[0] - 1); for (int j = 1; j <= FD_STENCIL_MAX; j++) for (int i = 1; i <= FD_STENCIL_MAX; i++) { const ptrdiff_t idx_dst = dir_y * j * strides[1] + dir_x * i; const ptrdiff_t idx_src = fact_y * j * strides[1] + fact_x * i; dst[idx_dst] = dst[idx_src]; } } } mg2di_timer_stop(&ctx->timer_bnd_corners); mg2di_timer_stop(&ctx->timer_bnd); } #if HAVE_NASM void mg2di_line_madd_fma3(ptrdiff_t linesize, double *dst, const double *src, double c); #endif static void line_madd_c(ptrdiff_t linesize, double *dst, const double *src, double c) { for (ptrdiff_t i = 0; i < linesize; i++) dst[i] += c * src[i]; } static int residual_add_task(void *arg, unsigned int job_idx, unsigned int thread_idx) { #if 1 EGSContext *ctx = arg; EGSInternal *priv = ctx->priv; priv->r.line_add(ctx->domain_size[0], ctx->u->data + ctx->u->stride[0] * job_idx, ctx->residual->data + ctx->residual->stride[0] * job_idx, priv->r.relax_factor); #else for (int idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) { ptrdiff_t idx = job_idx * ctx->u->stride[0] + idx0; ctx->u->data[idx] += priv->r.relax_factor * ctx->residual->data[idx]; } #endif return 0; } static int solve_relax_step(EGSContext *ctx, int export_res) { EGSInternal *priv = ctx->priv; EGSRelaxContext *r = ctx->relax; int u_next_valid = priv->u_next_valid; if (u_next_valid) { NDArray *tmp = ctx->u; NDArray *tmp_ext = ctx->u_exterior; NDArray *tmp_base = priv->u_base; ctx->u = priv->u_next; ctx->u_exterior = priv->u_next_exterior; priv->u_base = priv->u_next_base; priv->u_next = tmp; priv->u_next_exterior = tmp_ext; priv->u_next_base = tmp_base; priv->u_next_valid = 0; } if (export_res) { if (!u_next_valid) { mg2di_timer_start(&r->timer_correct); tp_execute(ctx->tp, ctx->domain_size[1], residual_add_task, ctx); mg2di_timer_stop(&r->timer_correct); priv->reflect_skip = 0; boundaries_apply(ctx, ctx->u, 0); priv->steps_since_sync++; if (priv->steps_since_sync >= priv->steps_per_sync) { boundaries_sync(ctx, ctx->u); priv->steps_since_sync = 0; } } residual_calc(ctx, 1); boundaries_sync(ctx, ctx->residual); } else { mg2di_assert(u_next_valid); residual_calc(ctx, 0); boundaries_apply(ctx, priv->u_next, 0); priv->steps_since_sync++; if (priv->steps_since_sync >= priv->steps_per_sync) { boundaries_sync(ctx, priv->u_next); priv->steps_since_sync = 0; } priv->u_next_valid = 1; } return 0; } static void exact_arrays_free(EGSContext *ctx) { EGSInternal *priv = ctx->priv; EGSExactInternal *e = &priv->e; free(ctx->priv->e.scratch_lines); free(ctx->priv->e.mat); free(ctx->priv->e.mat_f); free(ctx->priv->e.rhs); free(ctx->priv->e.rhs_mat); free(ctx->priv->e.rhs_factor); free(ctx->priv->e.x); free(ctx->priv->e.ipiv); mg2di_bicgstab_context_free(&e->bicgstab); e->arrays_alloced = 0; } static int exact_arrays_alloc(EGSContext *ctx) { EGSInternal *priv = ctx->priv; EGSExactInternal *e = &priv->e; int ret; if (e->arrays_alloced) return 0; e->N = ctx->domain_size[0] * ctx->domain_size[1]; e->N_ghosts = (ctx->domain_size[0] + 2 * FD_STENCIL_MAX) * (ctx->domain_size[1] + 2 * FD_STENCIL_MAX); e->mat = calloc(SQR(e->N), sizeof(*e->mat)); e->mat_f = calloc(SQR(e->N), sizeof(*e->mat_f)); e->rhs = calloc(e->N, sizeof(*e->rhs)); e->rhs_mat = calloc(e->N, sizeof(*e->rhs_mat)); e->rhs_factor = calloc(e->N, sizeof(*e->rhs_factor)); e->x = calloc(e->N, sizeof(*e->x)); e->ipiv = calloc(e->N, sizeof(*e->ipiv)); if (!e->mat || !e->mat_f || !e->rhs || !e->rhs_mat || !e->rhs_factor || !e->x || !e->ipiv) { ret = -ENOMEM; goto fail; } ret = mg2di_bicgstab_context_alloc(&e->bicgstab, e->N, 4); if (ret < 0) goto fail; for (int i = 0; i < e->N; i++) e->mat[i * e->N + i] = 1.0; ret = mg2di_bicgstab_init(e->bicgstab, e->mat, e->x); if (ret < 0) goto fail; e->arrays_alloced = 1; return 0; fail: exact_arrays_free(ctx); return ret; } static void fill_mat_s1(double *mat_row, ptrdiff_t mat_stride, ptrdiff_t row_stride, NDArray **diff_coeffs, ptrdiff_t idx_src) { mat_row[0] += diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_src]; mat_row[ 1 * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_src]; mat_row[-1 * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_src]; mat_row[ 1 * row_stride * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_src]; mat_row[-1 * row_stride * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_src]; mat_row[ 1 * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[ 0 * mat_stride] += -2.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[-1 * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[ 1 * row_stride * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[ 0 * row_stride * mat_stride] += -2.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[-1 * row_stride * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[( 1 + 1 * row_stride) * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 1 - 1 * row_stride) * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-1 + 1 * row_stride) * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-1 - 1 * row_stride) * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; } static void fill_mat_s2(double *mat_row, ptrdiff_t mat_stride, ptrdiff_t row_stride, NDArray **diff_coeffs, ptrdiff_t idx_src) { mat_row[0] += diff_coeffs[MG2D_DIFF_COEFF_00]->data[idx_src]; mat_row[ 2 * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_src]; mat_row[ 1 * mat_stride] += 8.0 * diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_src]; mat_row[-1 * mat_stride] += -8.0 * diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_src]; mat_row[-2 * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_10]->data[idx_src]; mat_row[ 2 * row_stride * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_src]; mat_row[ 1 * row_stride * mat_stride] += 8.0 * diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_src]; mat_row[-1 * row_stride * mat_stride] += -8.0 * diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_src]; mat_row[-2 * row_stride * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_01]->data[idx_src]; mat_row[ 2 * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[ 1 * mat_stride] += 16.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[ 0 * mat_stride] += -30.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[-1 * mat_stride] += 16.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[-2 * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_20]->data[idx_src]; mat_row[ 2 * row_stride * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[ 1 * row_stride * mat_stride] += 16.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[ 0 * row_stride * mat_stride] += -30.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[-1 * row_stride * mat_stride] += 16.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[-2 * row_stride * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_02]->data[idx_src]; mat_row[( 2 + 2 * row_stride) * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 2 + 1 * row_stride) * mat_stride] += -8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 2 - 1 * row_stride) * mat_stride] += 8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 2 - 2 * row_stride) * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 1 + 2 * row_stride) * mat_stride] += -8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 1 + 1 * row_stride) * mat_stride] += 64.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 1 - 1 * row_stride) * mat_stride] += -64.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[( 1 - 2 * row_stride) * mat_stride] += 8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-1 + 2 * row_stride) * mat_stride] += 8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-1 + 1 * row_stride) * mat_stride] += -64.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-1 - 1 * row_stride) * mat_stride] += 64.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-1 - 2 * row_stride) * mat_stride] += -8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-2 + 2 * row_stride) * mat_stride] += -1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-2 + 1 * row_stride) * mat_stride] += 8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-2 - 1 * row_stride) * mat_stride] += -8.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; mat_row[(-2 - 2 * row_stride) * mat_stride] += 1.0 * diff_coeffs[MG2D_DIFF_COEFF_11]->data[idx_src]; } static void mat_fill_row(EGSContext *ctx, double *scratch_line, ptrdiff_t idx1, ptrdiff_t idx0) { EGSInternal *priv = ctx->priv; EGSExactInternal *e = &priv->e; const ptrdiff_t idx_src_dc = NDA_IDX2D(priv->diff_coeffs[0], idx0, idx1); const ptrdiff_t mat_row_idx = idx1 * ctx->domain_size[0] + idx0; int is_bnd[4], boundary; ptrdiff_t row_stride, mat_stride; ptrdiff_t row_offset; double *mat_row; is_bnd[MG2D_BOUNDARY_0L] = idx0 < ctx->fd_stencil; is_bnd[MG2D_BOUNDARY_0U] = idx0 >= ctx->domain_size[0] - ctx->fd_stencil; is_bnd[MG2D_BOUNDARY_1L] = idx1 < ctx->fd_stencil; is_bnd[MG2D_BOUNDARY_1U] = idx1 >= ctx->domain_size[1] - ctx->fd_stencil; boundary = is_bnd[0] || is_bnd[1] || is_bnd[2] || is_bnd[3]; if (boundary) { memset(scratch_line, 0, e->N_ghosts * sizeof(*scratch_line)); row_stride = ctx->domain_size[0] + 2 * ctx->fd_stencil; mat_row = scratch_line + row_stride * ctx->fd_stencil + ctx->fd_stencil; mat_stride = 1; } else { mat_row = e->mat + mat_row_idx; row_stride = ctx->domain_size[0]; mat_stride = e->N; } row_offset = idx1 * row_stride + idx0; e->fill_mat(mat_row + row_offset * mat_stride, mat_stride, row_stride, priv->diff_coeffs, idx_src_dc); e->rhs_factor[mat_row_idx] = 1.0; if (boundary) { double *mat_dst = e->mat + mat_row_idx; ptrdiff_t mat_stride_dst = e->N; const MG2DBoundary *bnd_fixval = NULL; for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) { const int ci = mg2d_bnd_coord_idx(bnd_loc); const ptrdiff_t idx[] = { idx0, idx1 }; if (!(idx[ci] == (ctx->domain_size[ci] - 1) * mg2d_bnd_is_upper(bnd_loc) && ctx->boundaries[bnd_loc]->type == MG2D_BC_TYPE_FIXVAL)) continue; bnd_fixval = ctx->boundaries[bnd_loc]; mat_dst[mat_row_idx * mat_stride_dst] = 1.0; e->rhs_mat[mat_row_idx] = bnd_fixval->val[idx[!ci]]; e->rhs_factor[mat_row_idx] = 0.0; memset(scratch_line, 0, e->N_ghosts * sizeof(*scratch_line)); break; } if (!bnd_fixval) { /* apply the boundary conditions to eliminate the ghostpoint values */ // reflect for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) { MG2DBoundary *bnd = ctx->boundaries[bnd_loc]; const int ci = mg2d_bnd_coord_idx(bnd_loc); const int dir = mg2d_bnd_out_dir(bnd_loc); const ptrdiff_t dst_stride[2] = { 1, row_stride }; const ptrdiff_t bnd_stride[2] = { dst_stride[!ci], dst_stride[ci] }; const ptrdiff_t bnd_len[2] = { ctx->domain_size[!ci], ctx->domain_size[ci] }; double *dst = mat_row + mg2d_bnd_is_upper(bnd_loc) * ((bnd_len[1] - 1) * bnd_stride[1]); if (!is_bnd[bnd_loc] || bnd->type != MG2D_BC_TYPE_REFLECT) continue; for (ptrdiff_t bnd_layer = 1; bnd_layer <= ctx->fd_stencil; bnd_layer++) for (ptrdiff_t bnd_idx = -(ptrdiff_t)ctx->fd_stencil; bnd_idx < (ptrdiff_t)(bnd_len[0] + ctx->fd_stencil); bnd_idx++) { dst[-dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] += dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]]; dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] = 0.0; } } // fixval for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) { MG2DBoundary *bnd = ctx->boundaries[bnd_loc]; const int ci = mg2d_bnd_coord_idx(bnd_loc); const int dir = mg2d_bnd_out_dir(bnd_loc); const ptrdiff_t dst_stride[2] = { 1, row_stride }; const ptrdiff_t bnd_stride[2] = { dst_stride[!ci], dst_stride[ci] }; const ptrdiff_t bnd_len[2] = { ctx->domain_size[!ci], ctx->domain_size[ci] }; double *dst = mat_row + mg2d_bnd_is_upper(bnd_loc) * ((bnd_len[1] - 1) * bnd_stride[1]); if (!is_bnd[bnd_loc] || bnd->type != MG2D_BC_TYPE_FIXVAL) continue; for (ptrdiff_t bnd_layer = 1; bnd_layer < ctx->fd_stencil; bnd_layer++) for (ptrdiff_t bnd_idx = -((ptrdiff_t)ctx->fd_stencil - 1); bnd_idx < (ptrdiff_t)(bnd_len[0] + ctx->fd_stencil - 1); bnd_idx++) { e->rhs_mat[mat_row_idx] -= bnd->val[bnd_layer * bnd->val_stride + bnd_idx] * dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]]; dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] = 0.0; } } // falloff for (int bnd_loc = 0; bnd_loc < ARRAY_ELEMS(ctx->boundaries); bnd_loc++) { MG2DBoundary *bnd = ctx->boundaries[bnd_loc]; const int ci = mg2d_bnd_coord_idx(bnd_loc); const int dir = mg2d_bnd_out_dir(bnd_loc); const ptrdiff_t dst_stride[2] = { 1, row_stride }; const ptrdiff_t bnd_stride[2] = { dst_stride[!ci], dst_stride[ci] }; const ptrdiff_t bnd_len[2] = { ctx->domain_size[!ci], ctx->domain_size[ci] }; double *dst = mat_row + mg2d_bnd_is_upper(bnd_loc) * ((bnd_len[1] - 1) * bnd_stride[1]); if (!is_bnd[bnd_loc] || bnd->type != MG2D_BC_TYPE_FALLOFF) continue; for (ptrdiff_t bnd_layer = ctx->fd_stencil; bnd_layer > 0; bnd_layer--) for (ptrdiff_t bnd_idx = -(ptrdiff_t)ctx->fd_stencil; bnd_idx < (ptrdiff_t)(bnd_len[0] + ctx->fd_stencil); bnd_idx++) { const double x = ctx->step[0] * (ctx->domain_size[0] - 1 + bnd_layer - ctx->fd_stencil); const double y = ctx->step[0] * bnd_idx; const double r = sqrt(SQR(x) + SQR(y)); double val = dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]]; for (int i = 1; i < ctx->fd_stencil * 2 + 1; i++) { dst[dir * (bnd_layer - i) * bnd_stride[1] + bnd_idx * bnd_stride[0]] -= val * falloff_bnd_fd_factors[ctx->fd_stencil - 1][i] / falloff_bnd_fd_factors[ctx->fd_stencil - 1][0]; } dst[dir * (bnd_layer - (ptrdiff_t)ctx->fd_stencil) * bnd_stride[1] + bnd_idx * bnd_stride[0]] -= val * (x / SQR(r)) * (ctx->step[0] / falloff_bnd_fd_factors[ctx->fd_stencil - 1][0]); dst[dir * bnd_layer * bnd_stride[1] + bnd_idx * bnd_stride[0]] = 0.0; } } /* copy the interior values */ for (ptrdiff_t idx1_col = 0; idx1_col < ctx->domain_size[1]; idx1_col++) for (ptrdiff_t idx0_col = 0; idx0_col < ctx->domain_size[0]; idx0_col++) { mat_dst[(idx1_col * ctx->domain_size[0] + idx0_col) * mat_stride_dst] = mat_row[idx1_col * row_stride + idx0_col]; mat_row[idx1_col * row_stride + idx0_col] = 0.0; } } /* make sure all the values from the scratch line have been accounted for */ for (ptrdiff_t idx_scratch = 0; idx_scratch < e->N_ghosts; idx_scratch++) if (scratch_line[idx_scratch] != 0.0) abort(); } } static int mat_fill_row_task(void *arg, unsigned int job_idx, unsigned int thread_idx) { EGSContext *ctx = arg; EGSExactInternal *e = &ctx->priv->e; const ptrdiff_t idx1 = job_idx / ctx->domain_size[0]; const ptrdiff_t idx0 = job_idx % ctx->domain_size[0]; mat_fill_row(ctx, e->scratch_lines + e->N_ghosts * thread_idx, idx1, idx0); return 0; } static int solve_exact(EGSContext *ctx) { EGSInternal *priv = ctx->priv; EGSExactContext *ec = ctx->exact; EGSExactInternal *e = &priv->e; unsigned int nb_threads = tp_get_nb_threads(ctx->tp); int ret; if (!e->arrays_alloced) { ret = exact_arrays_alloc(ctx); if (ret < 0) return ret; } if (e->scratch_lines_allocated < nb_threads) { free(e->scratch_lines); e->scratch_lines = NULL; e->scratch_lines_allocated = 0; e->scratch_lines = calloc(e->N_ghosts * nb_threads, sizeof(*e->scratch_lines)); if (!e->scratch_lines) return -ENOMEM; e->scratch_lines_allocated = nb_threads; } mg2di_timer_start(&ec->timer_mat_construct); if (!e->mat_filled) { memset(e->mat, 0, SQR(e->N) * sizeof(*e->mat)); memset(e->rhs_mat, 0, e->N * sizeof(*e->rhs_mat)); tp_execute(ctx->tp, ctx->domain_size[1] * ctx->domain_size[0], mat_fill_row_task, ctx); e->mat_filled = 1; } for (ptrdiff_t idx1 = 0; idx1 < ctx->domain_size[1]; idx1++) for (ptrdiff_t idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) { const ptrdiff_t mat_row_idx = idx1 * ctx->domain_size[0] + idx0; e->rhs[mat_row_idx] = e->rhs_factor[mat_row_idx] * (*NDA_PTR2D(ctx->rhs, idx0, idx1)) + e->rhs_mat[mat_row_idx]; } mg2di_timer_stop(&ec->timer_mat_construct); mg2di_timer_start(&ec->timer_bicgstab); ret = mg2di_bicgstab_solve(e->bicgstab, e->mat, e->rhs, e->x); if (ret >= 0) ec->bicgstab_iterations += ret; mg2di_timer_stop(&ec->timer_bicgstab); if (ret < 0) { char equed = 'N'; double cond, ferr, berr, rpivot; mg2di_timer_start(&ec->timer_lu_solve); ret = LAPACKE_dgesvx(LAPACK_COL_MAJOR, 'N', 'N', e->N, 1, e->mat, e->N, e->mat_f, e->N, e->ipiv, &equed, NULL, NULL, e->rhs, e->N, e->x, e->N, &cond, &ferr, &berr, &rpivot); if (ret == 0) ret = LAPACKE_dgetri(LAPACK_COL_MAJOR, e->N, e->mat_f, e->N, e->ipiv); mg2di_timer_stop(&ec->timer_lu_solve); if (ret != 0) { mg2di_log(&ctx->logger, MG2D_LOG_ERROR, "Error solving the linear system: %d\n", ret); return -EDOM; } mg2di_log(&ctx->logger, MG2D_LOG_DEBUG, "LU factorization solution to a %zdx%zd matrix: " "condition number %16.16g; forward error %16.16g backward error %16.16g\n", e->N, e->N, cond, ferr, berr); ret = mg2di_bicgstab_init(e->bicgstab, e->mat_f, e->x); if (ret < 0) return ret; } mg2di_timer_start(&ec->timer_export); for (size_t idx1 = 0; idx1 < ctx->domain_size[1]; idx1++) memcpy(ctx->u->data + idx1 * ctx->u->stride[0], e->x + idx1 * ctx->domain_size[0], ctx->domain_size[0] * sizeof(*e->x)); mg2di_timer_stop(&ec->timer_export); priv->reflect_skip = 0; boundaries_apply(ctx, ctx->u, 0); residual_calc(ctx, 1); return 0; } int mg2di_egs_solve(EGSContext *ctx, enum EGSType solve_type, int export_res) { int ret; mg2di_timer_start(&ctx->timer_solve); switch (solve_type) { case EGS_SOLVE_RELAXATION: ret = solve_relax_step(ctx, export_res); break; case EGS_SOLVE_EXACT: ret = solve_exact(ctx); break; default: ret = -EINVAL; } mg2di_timer_stop(&ctx->timer_solve); return ret; } static int init_diff_coeffs_task(void *arg, unsigned int job_idx, unsigned int thread_idx) { struct { EGSContext *ctx; enum MG2DDiffCoeff dc; double fact; } *a = arg; EGSContext *ctx = a->ctx; EGSInternal *priv = ctx->priv; NDArray *dst = priv->diff_coeffs[a->dc]; const NDArray *src = ctx->diff_coeffs[a->dc]; const double factor = a->fact; for (ptrdiff_t idx1 = 0; idx1 < dst->shape[1]; idx1++) *NDA_PTR2D(dst, idx1, job_idx) = *NDA_PTR2D(src, idx1, job_idx) * factor; return 0; } static int init_mpi_sync(EGSContext *ctx) { EGSInternal *priv = ctx->priv; const DomainGeometry *dg = priv->dg; Rect *overlaps_recv = NULL, *overlaps_send = NULL; ptrdiff_t *lo; int ret = 0; overlaps_recv = calloc(dg->nb_components, sizeof(*overlaps_recv)); overlaps_send = calloc(dg->nb_components, sizeof(*overlaps_send)); if (!overlaps_recv || !overlaps_send) { ret = -ENOMEM; goto fail; } ret = mg2di_dg_edge_overlaps(overlaps_recv, overlaps_send, dg, priv->local_component, ctx->fd_stencil * priv->steps_per_sync); if (ret < 0) goto fail; priv->sync_sendtypes = calloc(dg->nb_components, sizeof(*priv->sync_sendtypes)); priv->sync_senddispl = calloc(dg->nb_components, sizeof(*priv->sync_senddispl)); priv->sync_sendcounts = calloc(dg->nb_components, sizeof(*priv->sync_sendcounts)); priv->sync_recvtypes = calloc(dg->nb_components, sizeof(*priv->sync_recvtypes)); priv->sync_recvdispl = calloc(dg->nb_components, sizeof(*priv->sync_recvdispl)); priv->sync_recvcounts = calloc(dg->nb_components, sizeof(*priv->sync_recvcounts)); if (!priv->sync_sendtypes || !priv->sync_senddispl || !priv->sync_sendcounts || !priv->sync_recvtypes || !priv->sync_recvdispl || !priv->sync_recvcounts) goto fail; lo = dg->components[priv->local_component].interior.start; /* construct the send/receive parameters */ for (unsigned int i = 0; i < dg->nb_components; i++) { if (i == priv->local_component) { MPI_Type_dup(MPI_INT, &priv->sync_sendtypes[i]); MPI_Type_dup(MPI_INT, &priv->sync_recvtypes[i]); priv->sync_sendcounts[i] = 0; priv->sync_recvcounts[i] = 0; priv->sync_senddispl[i] = 0; priv->sync_recvdispl[i] = 0; continue; } /* receive */ MPI_Type_vector(overlaps_recv[i].size[1], overlaps_recv[i].size[0], ctx->u->stride[0], MPI_DOUBLE, &priv->sync_recvtypes[i]); MPI_Type_commit(&priv->sync_recvtypes[i]); priv->sync_recvcounts[i] = 1; priv->sync_recvdispl[i] = ((overlaps_recv[i].start[1] - lo[1]) * ctx->u->stride[0] + (overlaps_recv[i].start[0] - lo[0])) * sizeof(*ctx->u->data); /* send */ MPI_Type_vector(overlaps_send[i].size[1], overlaps_send[i].size[0], ctx->u->stride[0], MPI_DOUBLE, &priv->sync_sendtypes[i]); MPI_Type_commit(&priv->sync_sendtypes[i]); priv->sync_sendcounts[i] = 1; priv->sync_senddispl[i] = ((overlaps_send[i].start[1] - lo[1]) * ctx->u->stride[0] + (overlaps_send[i].start[0] - lo[0])) * sizeof(*ctx->u->data); } fail: free(overlaps_recv); free(overlaps_send); return ret; } int mg2di_egs_init(EGSContext *ctx, int flags) { EGSInternal *priv = ctx->priv; EGSRelaxContext *r = ctx->relax; DomainComponent *dc = &priv->dg->components[priv->local_component]; int ret = 0; mg2di_timer_start(&ctx->timer_solve); mg2di_timer_start(&ctx->timer_init); if (ctx->step[0] <= DBL_EPSILON || ctx->step[1] <= DBL_EPSILON) { mg2di_log(&ctx->logger, 0, "Spatial step size too small\n"); ret = -EINVAL; goto finish; } /* initialize the FD stencil-dependent state */ if (!priv->fd_stencil) { if (priv->dg->nb_components > 1) { ret = init_mpi_sync(ctx); if (ret < 0) goto finish; } priv->fd_stencil = ctx->fd_stencil; } if (priv->fd_stencil != ctx->fd_stencil) { mg2di_log(&ctx->logger, MG2D_LOG_ERROR, "FD stencil changed\n"); ret = -EINVAL; goto finish; } if (r->relax_factor == 0.0) priv->r.relax_factor = relax_factors[ctx->fd_stencil - 1]; else priv->r.relax_factor = r->relax_factor; priv->r.relax_factor *= ctx->step[0] * ctx->step[0]; if (r->relax_multiplier > 0.0) priv->r.relax_factor *= r->relax_multiplier; priv->r.line_add = line_madd_c; #if HAVE_NASM if (ctx->cpuflags & MG2DI_CPU_FLAG_FMA3) priv->r.line_add = mg2di_line_madd_fma3; #endif if (!ctx->tp) { ret = tp_init(&priv->tp_internal, 1); if (ret < 0) goto finish; ctx->tp = priv->tp_internal; } if (!(flags & EGS_INIT_FLAG_SAME_DIFF_COEFFS)) { struct { EGSContext *ctx; enum MG2DDiffCoeff dc; double fact; } arg = { .ctx = ctx }; arg.dc = MG2D_DIFF_COEFF_00; arg.fact = 1.0 / fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_00]; tp_execute(ctx->tp, ctx->domain_size[1], init_diff_coeffs_task, &arg); arg.dc = MG2D_DIFF_COEFF_10; arg.fact = 1.0 / (fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_10] * ctx->step[0]); tp_execute(ctx->tp, ctx->domain_size[1], init_diff_coeffs_task, &arg); arg.dc = MG2D_DIFF_COEFF_01; arg.fact = 1.0 / (fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_01] * ctx->step[1]); tp_execute(ctx->tp, ctx->domain_size[1], init_diff_coeffs_task, &arg); arg.dc = MG2D_DIFF_COEFF_20; arg.fact = 1.0 / (fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_20] * SQR(ctx->step[0])); tp_execute(ctx->tp, ctx->domain_size[1], init_diff_coeffs_task, &arg); arg.dc = MG2D_DIFF_COEFF_02; arg.fact = 1.0 / (fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_02] * SQR(ctx->step[1])); tp_execute(ctx->tp, ctx->domain_size[1], init_diff_coeffs_task, &arg); arg.dc = MG2D_DIFF_COEFF_11; arg.fact = 1.0 / (fd_denoms[ctx->fd_stencil - 1][MG2D_DIFF_COEFF_11] * ctx->step[0] * ctx->step[1]); tp_execute(ctx->tp, ctx->domain_size[1], init_diff_coeffs_task, &arg); } if (!(flags & EGS_INIT_FLAG_SAME_DIFF_COEFFS)) priv->e.mat_filled = 0; switch (ctx->fd_stencil) { case 1: priv->e.fill_mat = fill_mat_s1; break; case 2: priv->e.fill_mat = fill_mat_s2; break; default: mg2di_log(&ctx->logger, 0, "Invalid finite difference stencil: %zd\n", ctx->fd_stencil); ret = -EINVAL; goto finish; } priv->residual_calc_size[priv->steps_per_sync - 1][0] = ctx->domain_size[0]; priv->residual_calc_size[priv->steps_per_sync - 1][1] = ctx->domain_size[1]; priv->residual_calc_offset[priv->steps_per_sync - 1][0] = ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL && dc->bnd_is_outer[MG2D_BOUNDARY_0L]; priv->residual_calc_offset[priv->steps_per_sync - 1][1] = ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXVAL && dc->bnd_is_outer[MG2D_BOUNDARY_1L]; for (int bnd_idx = 0; bnd_idx < ARRAY_ELEMS(ctx->boundaries); bnd_idx++) { MG2DBoundary *bnd = ctx->boundaries[bnd_idx]; const int ci = mg2d_bnd_coord_idx(bnd_idx); if (!dc->bnd_is_outer[bnd_idx]) continue; if (bnd->type == MG2D_BC_TYPE_FIXVAL) { double maxval = 0.0; for (int j = 0; j < ctx->fd_stencil; j++) { for (ptrdiff_t i = -j; i < (ptrdiff_t)ctx->domain_size[!ci] + j; i++) { double val = fabs(bnd->val[j * bnd->val_stride + i]); maxval = MAX(maxval, val); } } priv->bnd_zero[bnd_idx] = maxval == 0.0; priv->residual_calc_size[priv->steps_per_sync - 1][ci]--; } } for (int step = priv->steps_per_sync - 2; step >= 0; step--) { const int buffer_size = ctx->fd_stencil; int ic_0l = !dc->bnd_is_outer[MG2D_BOUNDARY_0L]; int ic_0u = !dc->bnd_is_outer[MG2D_BOUNDARY_0U]; int ic_1l = !dc->bnd_is_outer[MG2D_BOUNDARY_1L]; int ic_1u = !dc->bnd_is_outer[MG2D_BOUNDARY_1U]; priv->residual_calc_offset[step][0] = priv->residual_calc_offset[step + 1][0] - buffer_size * ic_0l; priv->residual_calc_offset[step][1] = priv->residual_calc_offset[step + 1][1] - buffer_size * ic_1l; priv->residual_calc_size[step][0] = priv->residual_calc_size[step + 1][0] + buffer_size * (ic_0l + ic_0u); priv->residual_calc_size[step][1] = priv->residual_calc_size[step + 1][1] + buffer_size * (ic_1l + ic_1u); if (step == 0) { priv->bnd_calc_offset[0] = -buffer_size * ic_0l; priv->bnd_calc_offset[1] = -buffer_size * ic_1l; priv->bnd_calc_size[0] = ctx->domain_size[0] + buffer_size * (ic_0l + ic_0u); priv->bnd_calc_size[1] = ctx->domain_size[1] + buffer_size * (ic_1l + ic_1u); } } priv->rescalc->tp = ctx->tp; priv->rescalc->fd_stencil = ctx->fd_stencil; priv->rescalc->cpuflags = ctx->cpuflags; ret = mg2di_residual_calc_init(priv->rescalc); if (ret < 0) goto finish; finish: mg2di_timer_stop(&ctx->timer_init); if (ret >= 0) { priv->reflect_skip = 0; if (!(flags & EGS_INIT_FLAG_SAME_DIFF_COEFFS)) { for (int dc_idx = 0; dc_idx < ARRAY_ELEMS(priv->diff_coeffs); dc_idx++) boundaries_sync(ctx, priv->diff_coeffs[dc_idx]); } boundaries_sync(ctx, ctx->rhs); boundaries_apply(ctx, ctx->u, 1); boundaries_sync(ctx, ctx->u); priv->steps_since_sync = 0; boundaries_apply(ctx, priv->u_next, 1); boundaries_sync(ctx, priv->u_next); residual_calc(ctx, 0); boundaries_apply(ctx, priv->u_next, 0); priv->steps_since_sync++; if (priv->steps_since_sync >= priv->steps_per_sync) { boundaries_sync(ctx, priv->u_next); priv->steps_since_sync = 0; } priv->u_next_valid = 1; } mg2di_timer_stop(&ctx->timer_solve); return ret; } static int arrays_alloc(EGSContext *ctx, const size_t domain_size[2]) { EGSInternal *priv = ctx->priv; const DomainComponent *dc = &priv->dg->components[priv->local_component]; NDASlice slice_interior[2]; NDASlice slice_exterior[2]; size_t size_alloc[2]; size_t size_dc[2]; size_t extra_points[4]; int ret; for (int bnd_idx = 0; bnd_idx < 4; bnd_idx++) extra_points[bnd_idx] = dc->bnd_is_outer[bnd_idx] ? 0 : FD_STENCIL_MAX * priv->steps_per_sync; size_alloc[0] = dc->exterior.size[1] + extra_points[MG2D_BOUNDARY_1L] + extra_points[MG2D_BOUNDARY_1U]; size_alloc[1] = dc->exterior.size[0] + extra_points[MG2D_BOUNDARY_0L] + extra_points[MG2D_BOUNDARY_0U]; slice_exterior[0].start = extra_points[MG2D_BOUNDARY_1L]; slice_exterior[0].end = slice_exterior[0].start + dc->exterior.size[1]; slice_exterior[0].step = 1; slice_exterior[1].start = extra_points[MG2D_BOUNDARY_0L]; slice_exterior[1].end = slice_exterior[1].start + dc->exterior.size[0]; slice_exterior[1].step = 1; slice_interior[0].start = slice_exterior[0].start + (dc->interior.start[1] - dc->exterior.start[1]); slice_interior[0].end = slice_interior[0].start + dc->interior.size[1]; slice_interior[0].step = 1; slice_interior[1].start = slice_exterior[1].start + (dc->interior.start[0] - dc->exterior.start[0]); slice_interior[1].end = slice_interior[1].start + dc->interior.size[0]; slice_interior[1].step = 1; ret = ndarray_alloc(&priv->u_base, 2, size_alloc, NDARRAY_ALLOC_ZERO); if (ret < 0) return ret; ret = ndarray_alloc(&priv->u_next_base, 2, size_alloc, NDARRAY_ALLOC_ZERO); if (ret < 0) return ret; ret = ndarray_slice(&ctx->u, priv->u_base, slice_interior); if (ret < 0) return ret; ret = ndarray_slice(&priv->u_next, priv->u_next_base, slice_interior); if (ret < 0) return ret; ret = ndarray_slice(&ctx->u_exterior, priv->u_base, slice_exterior); if (ret < 0) return ret; ret = ndarray_slice(&priv->u_next_exterior, priv->u_next_base, slice_exterior); if (ret < 0) return ret; ret = ndarray_alloc(&priv->rhs_base, 2, size_alloc, 0); if (ret < 0) return ret; ret = ndarray_slice(&ctx->rhs, priv->rhs_base, slice_interior); if (ret < 0) return ret; ret = ndarray_alloc(&priv->residual_base, 2, size_alloc, NDARRAY_ALLOC_ZERO); if (ret < 0) return ret; ret = ndarray_slice(&ctx->residual, priv->residual_base, slice_interior); if (ret < 0) return ret; size_dc[0] = size_alloc[0] * MG2D_DIFF_COEFF_NB; size_dc[1] = size_alloc[1]; ret = ndarray_alloc(&priv->diff_coeffs_public_base, 2, size_dc, NDARRAY_ALLOC_ZERO); if (ret < 0) return ret; ret = ndarray_alloc(&priv->diff_coeffs_base, 2, size_dc, NDARRAY_ALLOC_ZERO); if (ret < 0) return ret; for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) { const NDASlice slice_dc[2] = { NDASLICE(i * size_alloc[0], (i + 1) * size_alloc[0], 1), NDASLICE_NULL }; NDArray *tmp; ret = ndarray_slice(&tmp, priv->diff_coeffs_public_base, slice_dc); if (ret < 0) return ret; ret = ndarray_slice(&ctx->diff_coeffs[i], tmp, slice_interior); ndarray_free(&tmp); if (ret < 0) return ret; ret = ndarray_slice(&tmp, priv->diff_coeffs_base, slice_dc); if (ret < 0) return ret; ret = ndarray_slice(&priv->diff_coeffs[i], tmp, slice_interior); ndarray_free(&tmp); if (ret < 0) return ret; } for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) { const int ci = mg2d_bnd_coord_idx(i); ctx->boundaries[i] = mg2di_bc_alloc(domain_size[!ci]); if (!ctx->boundaries[i]) return -ENOMEM; } return 0; } static EGSContext *egs_alloc(const DomainGeometry *dg, unsigned int local_component) { EGSContext *ctx; int ret; ctx = calloc(1, sizeof(*ctx)); if (!ctx) return NULL; ctx->priv = calloc(1, sizeof(*ctx->priv)); if (!ctx->priv) goto fail; ctx->relax = calloc(1, sizeof(*ctx->relax)); if (!ctx->relax) goto fail; mg2di_timer_init(&ctx->relax->timer_correct); ctx->exact = calloc(1, sizeof(*ctx->exact)); if (!ctx->exact) goto fail; ret = mg2di_dg_copy(&ctx->priv->dg, dg); if (ret < 0) goto fail; ctx->priv->local_component = local_component; ctx->priv->comm = MPI_COMM_NULL; if (dg->nb_components > 1) { char *steps = getenv("MG2D_STEPS_PER_SYNC"); if (steps) ctx->priv->steps_per_sync = strtol(steps, NULL, 0); if (ctx->priv->steps_per_sync <= 0) ctx->priv->steps_per_sync = 2; for (int i = 0; i < dg->nb_components; i++) { if (dg->components[i].interior.size[0] < ctx->priv->steps_per_sync * FD_STENCIL_MAX || dg->components[i].interior.size[1] < ctx->priv->steps_per_sync * FD_STENCIL_MAX) { ctx->priv->steps_per_sync = 1; break; } } } else ctx->priv->steps_per_sync = 1; ctx->priv->residual_calc_offset = calloc(ctx->priv->steps_per_sync, sizeof(*ctx->priv->residual_calc_offset)); ctx->priv->residual_calc_size = calloc(ctx->priv->steps_per_sync, sizeof(*ctx->priv->residual_calc_size)); mg2di_timer_init(&ctx->exact->timer_mat_construct); mg2di_timer_init(&ctx->exact->timer_bicgstab); mg2di_timer_init(&ctx->exact->timer_lu_solve); mg2di_timer_init(&ctx->exact->timer_export); if (!dg->domain_size[0] || !dg->domain_size[1] || dg->domain_size[0] > SIZE_MAX / dg->domain_size[1]) goto fail; ret = arrays_alloc(ctx, dg->components[local_component].interior.size); if (ret < 0) goto fail; ctx->domain_size[0] = dg->components[local_component].interior.size[0]; ctx->domain_size[1] = dg->components[local_component].interior.size[1]; ctx->priv->rescalc = mg2di_residual_calc_alloc(); if (!ctx->priv->rescalc) goto fail; mg2di_timer_init(&ctx->timer_bnd); mg2di_timer_init(&ctx->timer_bnd_fixval); mg2di_timer_init(&ctx->timer_bnd_falloff); mg2di_timer_init(&ctx->timer_bnd_reflect); mg2di_timer_init(&ctx->timer_bnd_corners); mg2di_timer_init(&ctx->timer_res_calc); mg2di_timer_init(&ctx->timer_init); mg2di_timer_init(&ctx->timer_solve); mg2di_timer_init(&ctx->timer_mpi_sync); return ctx; fail: mg2di_egs_free(&ctx); return NULL; } EGSContext *mg2di_egs_alloc(const size_t domain_size[2]) { EGSContext *ctx; DomainGeometry *dg; dg = mg2di_dg_alloc(1); if (!dg) return NULL; dg->domain_size[0] = domain_size[0]; dg->domain_size[1] = domain_size[1]; dg->components[0].interior.start[0] = 0; dg->components[0].interior.start[1] = 0; dg->components[0].interior.size[0] = domain_size[0]; dg->components[0].interior.size[1] = domain_size[1]; dg->components[0].exterior.start[0] = -FD_STENCIL_MAX; dg->components[0].exterior.start[1] = -FD_STENCIL_MAX; dg->components[0].exterior.size[0] = domain_size[0] + 2 * FD_STENCIL_MAX; dg->components[0].exterior.size[1] = domain_size[1] + 2 * FD_STENCIL_MAX; for (int i = 0; i < ARRAY_ELEMS(dg->components[0].bnd_is_outer); i++) dg->components[0].bnd_is_outer[i] = 1; ctx = egs_alloc(dg, 0); mg2di_dg_free(&dg); if (!ctx) return NULL; return ctx; } EGSContext *mg2di_egs_alloc_mpi(MPI_Comm comm, const DomainGeometry *dg) { EGSContext *ctx = NULL; int local_component; MPI_Comm_rank(comm, &local_component); ctx = egs_alloc(dg, local_component); if (!ctx) return NULL; ctx->priv->comm = comm; return ctx; } void mg2di_egs_free(EGSContext **pctx) { EGSContext *ctx = *pctx; if (!ctx) return; mg2di_residual_calc_free(&ctx->priv->rescalc); free(ctx->relax); free(ctx->exact); exact_arrays_free(ctx); ndarray_free(&ctx->u); ndarray_free(&ctx->u_exterior); ndarray_free(&ctx->priv->u_base); ndarray_free(&ctx->priv->u_next); ndarray_free(&ctx->priv->u_next_exterior); ndarray_free(&ctx->priv->u_next_base); ndarray_free(&ctx->rhs); ndarray_free(&ctx->priv->rhs_base); ndarray_free(&ctx->residual); ndarray_free(&ctx->priv->residual_base); for (int i = 0; i < ARRAY_ELEMS(ctx->diff_coeffs); i++) { ndarray_free(&ctx->diff_coeffs[i]); } ndarray_free(&ctx->priv->diff_coeffs_public_base); for (int i = 0; i < ARRAY_ELEMS(ctx->priv->diff_coeffs); i++) { ndarray_free(&ctx->priv->diff_coeffs[i]); } ndarray_free(&ctx->priv->diff_coeffs_base); for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) mg2di_bc_free(&ctx->boundaries[i]); if (ctx->priv->sync_sendtypes) { for (int i = 0; i < ctx->priv->dg->nb_components; i++) { if (ctx->priv->sync_sendtypes[i]) MPI_Type_free(&ctx->priv->sync_sendtypes[i]); } } if (ctx->priv->sync_recvtypes) { for (int i = 0; i < ctx->priv->dg->nb_components; i++) { if (ctx->priv->sync_recvtypes[i]) MPI_Type_free(&ctx->priv->sync_recvtypes[i]); } } free(ctx->priv->sync_sendtypes); free(ctx->priv->sync_recvtypes); free(ctx->priv->sync_sendcounts); free(ctx->priv->sync_senddispl); free(ctx->priv->sync_recvcounts); free(ctx->priv->sync_recvdispl); free(ctx->priv->residual_calc_offset); free(ctx->priv->residual_calc_size); mg2di_dg_free(&ctx->priv->dg); tp_free(&ctx->priv->tp_internal); free(ctx->priv); free(ctx); *pctx = NULL; }