From 1e3d13b0f885479fca1e02bf4fdbdba219dc0578 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Wed, 6 Feb 2019 14:09:34 +0100 Subject: ell_grid_solve: use the new boundary API to simplify code --- ell_grid_solve.c | 301 +++++++++++++++++++------------------------------------ 1 file changed, 105 insertions(+), 196 deletions(-) diff --git a/ell_grid_solve.c b/ell_grid_solve.c index 75960e7..0cefdfd 100644 --- a/ell_grid_solve.c +++ b/ell_grid_solve.c @@ -77,28 +77,6 @@ struct EGSInternal { TPContext *tp_internal; }; -static const struct { - unsigned int stride_idx; - unsigned int is_upper; -} boundary_def[] = { - [MG2D_BOUNDARY_0L] = { - .stride_idx = 1, - .is_upper = 0, - }, - [MG2D_BOUNDARY_0U] = { - .stride_idx = 1, - .is_upper = 1, - }, - [MG2D_BOUNDARY_1L] = { - .stride_idx = 0, - .is_upper = 0 - }, - [MG2D_BOUNDARY_1U] = { - .stride_idx = 0, - .is_upper = 1, - }, -}; - static const double fd_denoms[][MG2D_DIFF_COEFF_NB] = { { [MG2D_DIFF_COEFF_00] = 1.0, @@ -170,15 +148,14 @@ static void boundaries_apply(EGSContext *ctx) start = gettime(); for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) { - const int si = boundary_def[i].stride_idx; - const ptrdiff_t stride_boundary = strides[si]; - const ptrdiff_t stride_offset = strides[!si]; - const size_t size_boundary = ctx->domain_size[si]; - const size_t size_offset = ctx->domain_size[!si]; + 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_offset = ctx->domain_size[ci]; - double *dst = ctx->u + boundary_def[i].is_upper * ((size_offset - 1) * stride_offset); - const ptrdiff_t dst_strides[] = { stride_boundary, - (boundary_def[i].is_upper ? 1 : -1) * stride_offset }; + double *dst = ctx->u + mg2d_bnd_is_upper(i) * ((size_offset - 1) * stride_offset); + const ptrdiff_t dst_strides[] = { stride_boundary, mg2d_bnd_out_dir(i) * stride_offset }; switch (ctx->boundaries[i]->type) { case MG2D_BC_TYPE_FIXVAL: @@ -193,73 +170,35 @@ static void boundaries_apply(EGSContext *ctx) } /* fill in the corner ghosts */ - if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF || - ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF) { - double *dst = ctx->u; - int fact_x = -1, fact_z = -1; - - if (ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF) - fact_z *= -1; - else - fact_x *= -1; - - for (int i = 1; i <= FD_STENCIL_MAX; i++) - for (int j = 1; j <= FD_STENCIL_MAX; j++) { - const ptrdiff_t idx_dst = -j * strides[1] - i; - const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i; - dst[idx_dst] = dst[idx_src]; - } - } - if (ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF || - ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF) { - double *dst = ctx->u + ctx->domain_size[0] - 1; - int fact_x = 1, fact_z = -1; - - if (ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF) - fact_z *= -1; - else - fact_x *= -1; - - for (int i = 1; i <= FD_STENCIL_MAX; i++) - for (int j = 1; j <= FD_STENCIL_MAX; j++) { - const ptrdiff_t idx_dst = -j * strides[1] + i; - const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i; - dst[idx_dst] = dst[idx_src]; - } - } - if (ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF || - ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) { - double *dst = ctx->u + strides[1] * (ctx->domain_size[1] - 1); - int fact_x = -1, fact_z = 1; - - if (ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) - fact_z *= -1; - else - fact_x *= -1; - - for (int i = 1; i <= FD_STENCIL_MAX; i++) - for (int j = 1; j <= FD_STENCIL_MAX; j++) { - const ptrdiff_t idx_dst = j * strides[1] - i; - const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i; - dst[idx_dst] = dst[idx_src]; - } - } - if (ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF || - ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) { - double *dst = ctx->u + strides[1] * (ctx->domain_size[1] - 1) + ctx->domain_size[0] - 1; - int fact_x = 1, fact_z = 1; - - if (ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) - fact_z *= -1; - else - fact_x *= -1; - - for (int i = 1; i <= FD_STENCIL_MAX; i++) - for (int j = 1; j <= FD_STENCIL_MAX; j++) { - const ptrdiff_t idx_dst = j * strides[1] + i; - const ptrdiff_t idx_src = fact_z * j * strides[1] + fact_x * i; - dst[idx_dst] = dst[idx_src]; - } + 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); + + 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 = ctx->u + + mg2d_bnd_is_upper(loc_y) * ((ctx->domain_size[1] - 1) * priv->stride) + + mg2d_bnd_is_upper(loc_x) * (ctx->domain_size[0] - 1); + + if (bnd_y->type == MG2D_BC_TYPE_FIXDIFF) + fact_y *= -1; + else if (bnd_x->type == MG2D_BC_TYPE_FIXDIFF) + fact_x *= -1; + else + continue; + + 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]; + } + } } ctx->time_boundaries += gettime() - start; ctx->count_boundaries++; @@ -297,19 +236,6 @@ static int solve_relax_step(EGSContext *ctx) return 0; } -static int bnd_id(int x_lower, int x_upper, int y_lower, int y_upper) -{ - if (x_lower) - return MG2D_BOUNDARY_0L; - if (x_upper) - return MG2D_BOUNDARY_0U; - if (y_lower) - return MG2D_BOUNDARY_1L; - if (y_upper) - return MG2D_BOUNDARY_1U; - return -1; -} - static void fill_mat_s1(double *mat_row, double **diff_coeffs, double *fd_factors, ptrdiff_t idx_src, ptrdiff_t row_stride, ptrdiff_t row_offset) { @@ -397,15 +323,12 @@ static int solve_exact(EGSContext *ctx) memset(e->mat, 0, SQR(e->N) * sizeof(*e->mat)); for (ptrdiff_t idx1 = 0; idx1 < ctx->domain_size[1]; idx1++) { - const int boundary_y_lower = idx1 < ctx->fd_stencil; - const int boundary_y_upper = idx1 >= ctx->domain_size[1] - ctx->fd_stencil; - - for (ptrdiff_t idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) { - const int boundary_x_lower = idx0 < ctx->fd_stencil; - const int boundary_x_upper = idx0 >= ctx->domain_size[0] - ctx->fd_stencil; + int is_bnd[4], boundary; - const int boundary = boundary_y_lower || boundary_y_upper || boundary_x_lower || boundary_x_upper; + is_bnd[MG2D_BOUNDARY_1L] = idx1 < ctx->fd_stencil; + is_bnd[MG2D_BOUNDARY_1U] = idx1 >= ctx->domain_size[1] - ctx->fd_stencil; + for (ptrdiff_t idx0 = 0; idx0 < ctx->domain_size[0]; idx0++) { const ptrdiff_t idx_src = idx1 * priv->stride + idx0; const ptrdiff_t mat_row_idx = idx1 * ctx->domain_size[0] + idx0; @@ -414,6 +337,11 @@ static int solve_exact(EGSContext *ctx) 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; + + boundary = is_bnd[0] || is_bnd[1] || is_bnd[2] || is_bnd[3]; + if (boundary) { memset(e->scratch_line, 0, e->N_ghosts * sizeof(*e->scratch_line)); row_stride = ctx->domain_size[0] + 2 * ctx->fd_stencil; @@ -430,92 +358,71 @@ static int solve_exact(EGSContext *ctx) e->rhs[mat_row_idx] = ctx->rhs[idx_src]; if (boundary) { - const int bnd_fixval = (idx0 == 0 && ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL) || - (idx1 == 0 && ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXVAL) || - (idx0 == ctx->domain_size[1] - 1 && ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXVAL) || - (idx1 == ctx->domain_size[0] - 1 && ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXVAL); - double *mat_dst = e->mat + e->N * mat_row_idx; - ptrdiff_t bnd_idx = (boundary_x_lower || boundary_x_upper) ? idx0 : idx1; - if (bnd_fixval) { - const MG2DBoundary *bnd = ctx->boundaries[bnd_id(boundary_x_lower, boundary_x_upper, - boundary_y_lower, boundary_y_upper)]; + const MG2DBoundary *bnd_fixval = NULL; - mat_dst[mat_row_idx] = 1.0; - e->rhs[mat_row_idx] = bnd->val[bnd_idx]; + 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] = 1.0; + e->rhs[mat_row_idx] = bnd_fixval->val[idx[!ci]]; memset(e->scratch_line, 0, e->N_ghosts * sizeof(*e->scratch_line)); - } else { + break; + } + + if (!bnd_fixval) { /* apply the boundary conditions to eliminate the ghostpoint values */ // fixdiff - if (boundary_y_lower && ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXDIFF) { - for (ptrdiff_t idx1_bnd = 0; idx1_bnd < ctx->fd_stencil; idx1_bnd++) - for (ptrdiff_t idx0_col = -(ptrdiff_t)ctx->fd_stencil; idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil); idx0_col++) { - const ptrdiff_t idx_ghost = -(idx1_bnd + 1) * row_stride + idx0_col; - mat_row[(idx1_bnd + 1) * row_stride + idx0_col] += mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; - } - } - if (boundary_y_upper && ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXDIFF) { - for (ptrdiff_t idx1_bnd = 0; idx1_bnd < ctx->fd_stencil; idx1_bnd++) - for (ptrdiff_t idx0_col = -(ptrdiff_t)ctx->fd_stencil; idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil); idx0_col++) { - const ptrdiff_t idx_ghost = (ctx->domain_size[1] - 1 + idx1_bnd + 1) * row_stride + idx0_col; - mat_row[(ctx->domain_size[1] - 1 - (idx1_bnd + 1)) * row_stride + idx0_col] += mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; - } - } - if (boundary_x_lower && ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXDIFF) { - for (ptrdiff_t idx1_col = -(ptrdiff_t)ctx->fd_stencil; idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil); idx1_col++) - for (ptrdiff_t idx0_bnd = 0; idx0_bnd < ctx->fd_stencil; idx0_bnd++) { - const ptrdiff_t idx_ghost = idx1_col * row_stride - (idx0_bnd + 1); - mat_row[idx1_col * row_stride + idx0_bnd + 1] += mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; - } - } - if (boundary_x_upper && ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXDIFF) { - for (ptrdiff_t idx1_col = -(ptrdiff_t)ctx->fd_stencil; idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil); idx1_col++) - for (ptrdiff_t idx0_bnd = 0; idx0_bnd < ctx->fd_stencil; idx0_bnd++) { - const ptrdiff_t idx_ghost = idx1_col * row_stride + ctx->domain_size[0] - 1 + idx0_bnd + 1; - mat_row[idx1_col * row_stride + ctx->domain_size[0] - 1 - (idx0_bnd + 1)] += mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; + 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_FIXDIFF) + 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 - if (boundary_y_lower && ctx->boundaries[MG2D_BOUNDARY_1L]->type == MG2D_BC_TYPE_FIXVAL) { - const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1L]; - for (ptrdiff_t idx1_bnd = 1; idx1_bnd < ctx->fd_stencil; idx1_bnd++) - for (ptrdiff_t idx0_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil - 1); idx0_col++) { - const ptrdiff_t idx_ghost = -idx1_bnd * row_stride + idx0_col; - e->rhs[mat_row_idx] -= bnd->val[idx1_bnd * bnd->val_stride + idx0_col] * mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; - } - } - if (boundary_y_upper && ctx->boundaries[MG2D_BOUNDARY_1U]->type == MG2D_BC_TYPE_FIXVAL) { - const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_1U]; - for (ptrdiff_t idx1_bnd = 1; idx1_bnd < ctx->fd_stencil; idx1_bnd++) - for (ptrdiff_t idx0_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx0_col < (ptrdiff_t)(ctx->domain_size[0] + ctx->fd_stencil - 1); idx0_col++) { - const ptrdiff_t idx_ghost = (ctx->domain_size[1] - 1 + idx1_bnd) * row_stride + idx0_col; - e->rhs[mat_row_idx] -= bnd->val[idx1_bnd * bnd->val_stride + idx0_col] * mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; - } - } - if (boundary_x_lower && ctx->boundaries[MG2D_BOUNDARY_0L]->type == MG2D_BC_TYPE_FIXVAL) { - const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0L]; - for (ptrdiff_t idx1_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil - 1); idx1_col++) - for (ptrdiff_t idx0_bnd = 1; idx0_bnd < ctx->fd_stencil; idx0_bnd++) { - const ptrdiff_t idx_ghost = idx1_col * row_stride - idx0_bnd; - e->rhs[mat_row_idx] -= bnd->val[idx0_bnd * bnd->val_stride + idx1_col] * mat_row[idx_ghost]; - mat_row[idx1_col * row_stride - idx0_bnd] = 0.0; - } - } - if (boundary_x_upper && ctx->boundaries[MG2D_BOUNDARY_0U]->type == MG2D_BC_TYPE_FIXVAL) { - const MG2DBoundary *bnd = ctx->boundaries[MG2D_BOUNDARY_0U]; - for (ptrdiff_t idx1_col = -((ptrdiff_t)ctx->fd_stencil - 1); idx1_col < (ptrdiff_t)(ctx->domain_size[1] + ctx->fd_stencil - 1); idx1_col++) - for (ptrdiff_t idx0_bnd = 1; idx0_bnd < ctx->fd_stencil; idx0_bnd++) { - const ptrdiff_t idx_ghost = idx1_col * row_stride + ctx->domain_size[0] - 1 + idx0_bnd; - e->rhs[mat_row_idx] -= bnd->val[idx0_bnd * bnd->val_stride + idx1_col] * mat_row[idx_ghost]; - mat_row[idx_ghost] = 0.0; + 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; bnd_idx < (ptrdiff_t)(bnd_len[0] + ctx->fd_stencil); bnd_idx++) { + e->rhs[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; } } @@ -637,14 +544,15 @@ int mg2di_egs_init(EGSContext *ctx) for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) { MG2DBoundary *bnd = ctx->boundaries[i]; + const int ci = mg2d_bnd_coord_idx(i); if (bnd->type == MG2D_BC_TYPE_FIXDIFF) { - for (int k = 0; k < ctx->domain_size[boundary_def[i].stride_idx]; k++) + for (int k = 0; k < ctx->domain_size[!ci]; k++) if (bnd->val[k] != 0.0) { mg2di_log(&ctx->logger, 0, "Only zero boundary derivative supported\n"); return -ENOSYS; } } else if (bnd->type == MG2D_BC_TYPE_FIXVAL) { - priv->residual_calc_size[!boundary_def[i].stride_idx]--; + priv->residual_calc_size[ci]--; } } @@ -709,7 +617,8 @@ static int arrays_alloc(EGSContext *ctx, const size_t domain_size[2]) priv->stride = stride; for (int i = 0; i < ARRAY_ELEMS(ctx->boundaries); i++) { - ctx->boundaries[i] = mg2di_bc_alloc(domain_size[boundary_def[i].stride_idx]); + const int ci = mg2d_bnd_coord_idx(i); + ctx->boundaries[i] = mg2di_bc_alloc(domain_size[!ci]); if (!ctx->boundaries[i]) return -ENOMEM; } -- cgit v1.2.3