Line data Source code
1 : /**
2 : * @file boundary_copy_utils.h
3 : * @brief Shared boundary copying utilities for projection method solvers
4 : *
5 : * This header provides helper functions used across all projection solver
6 : * backends (CPU, AVX2, OpenMP) to copy boundary velocity values.
7 : *
8 : * Include from backend subdirectories as: #include "../boundary_copy_utils.h"
9 : */
10 :
11 : #ifndef BOUNDARY_COPY_UTILS_H
12 : #define BOUNDARY_COPY_UTILS_H
13 :
14 : #include "cfd/core/indexing.h"
15 : #include <stddef.h>
16 :
17 : /**
18 : * Copy boundary values between velocity arrays
19 : *
20 : * Used to preserve caller-set boundary conditions (e.g., lid velocity for
21 : * cavity flow) during the projection method's predictor and corrector steps.
22 : *
23 : * Copies ALL four boundaries (top, bottom, left, right). For flows with
24 : * outlet boundaries (Neumann BCs), use copy_dirichlet_boundaries instead.
25 : *
26 : * @param dst_u Destination u-velocity array
27 : * @param dst_v Destination v-velocity array
28 : * @param src_u Source u-velocity array
29 : * @param src_v Source v-velocity array
30 : * @param nx Grid width
31 : * @param ny Grid height
32 : */
33 : static inline void copy_boundary_velocities(double* dst_u, double* dst_v,
34 : const double* src_u, const double* src_v,
35 : size_t nx, size_t ny) {
36 : // Bottom and top boundaries (j = 0 and j = ny-1)
37 : for (size_t i = 0; i < nx; i++) {
38 : dst_u[i] = src_u[i];
39 : dst_v[i] = src_v[i];
40 : dst_u[IDX_2D(i, ny - 1, nx)] = src_u[IDX_2D(i, ny - 1, nx)];
41 : dst_v[IDX_2D(i, ny - 1, nx)] = src_v[IDX_2D(i, ny - 1, nx)];
42 : }
43 : // Left and right boundaries (i = 0 and i = nx-1)
44 : for (size_t j = 1; j < ny - 1; j++) {
45 : dst_u[IDX_2D(0, j, nx)] = src_u[IDX_2D(0, j, nx)];
46 : dst_v[IDX_2D(0, j, nx)] = src_v[IDX_2D(0, j, nx)];
47 : dst_u[IDX_2D(nx - 1, j, nx)] = src_u[IDX_2D(nx - 1, j, nx)];
48 : dst_v[IDX_2D(nx - 1, j, nx)] = src_v[IDX_2D(nx - 1, j, nx)];
49 : }
50 : }
51 :
52 : /**
53 : * Copy only Dirichlet boundary values (walls + inlet), skip outlet
54 : *
55 : * Used for flows with outlet boundaries (e.g., Poiseuille channel flow).
56 : * Copies top, bottom, and left boundaries, but NOT the right boundary
57 : * to allow the pressure correction to update the outlet velocity.
58 : *
59 : * @param dst_u Destination u-velocity array
60 : * @param dst_v Destination v-velocity array
61 : * @param src_u Source u-velocity array
62 : * @param src_v Source v-velocity array
63 : * @param nx Grid width
64 : * @param ny Grid height
65 : */
66 : static inline void copy_dirichlet_boundaries(double* dst_u, double* dst_v,
67 : const double* src_u, const double* src_v,
68 : size_t nx, size_t ny) {
69 : // Bottom and top boundaries (j = 0 and j = ny-1)
70 : for (size_t i = 0; i < nx; i++) {
71 : dst_u[i] = src_u[i];
72 : dst_v[i] = src_v[i];
73 : dst_u[IDX_2D(i, ny - 1, nx)] = src_u[IDX_2D(i, ny - 1, nx)];
74 : dst_v[IDX_2D(i, ny - 1, nx)] = src_v[IDX_2D(i, ny - 1, nx)];
75 : }
76 : // Left boundary only (i = 0) - skip right boundary to preserve outlet
77 : for (size_t j = 1; j < ny - 1; j++) {
78 : dst_u[IDX_2D(0, j, nx)] = src_u[IDX_2D(0, j, nx)];
79 : dst_v[IDX_2D(0, j, nx)] = src_v[IDX_2D(0, j, nx)];
80 : }
81 : }
82 :
83 : /**
84 : * Copy boundary values for 3D velocity arrays (u, v, w on all 6 faces)
85 : *
86 : * When nz <= 1, copies only the 4 xy-boundary edges (same as 2D) and
87 : * ignores w. When nz > 1, also copies the front/back z-faces.
88 : *
89 : * @param dst_u, dst_v, dst_w Destination velocity arrays
90 : * @param src_u, src_v, src_w Source velocity arrays
91 : * @param nx, ny, nz Grid dimensions (nz=1 for 2D)
92 : */
93 23296 : static inline void copy_boundary_velocities_3d(
94 : double* dst_u, double* dst_v, double* dst_w,
95 : const double* src_u, const double* src_v, const double* src_w,
96 : size_t nx, size_t ny, size_t nz)
97 : {
98 23296 : size_t plane = nx * ny;
99 :
100 63118 : for (size_t k = 0; k < nz; k++) {
101 39822 : size_t base = k * plane;
102 :
103 : /* Bottom and top y-boundaries (j = 0 and j = ny-1) */
104 1467836 : for (size_t i = 0; i < nx; i++) {
105 1428014 : size_t bot = base + i;
106 1428014 : size_t top = base + IDX_2D(i, ny - 1, nx);
107 1428014 : dst_u[bot] = src_u[bot];
108 1428014 : dst_v[bot] = src_v[bot];
109 1428014 : dst_u[top] = src_u[top];
110 1428014 : dst_v[top] = src_v[top];
111 1428014 : if (nz > 1) {
112 269952 : dst_w[bot] = src_w[bot];
113 269952 : dst_w[top] = src_w[top];
114 : }
115 : }
116 : /* Left and right x-boundaries (i = 0 and i = nx-1) */
117 1385004 : for (size_t j = 1; j < ny - 1; j++) {
118 1345182 : size_t left = base + IDX_2D(0, j, nx);
119 1345182 : size_t right = base + IDX_2D(nx - 1, j, nx);
120 1345182 : dst_u[left] = src_u[left];
121 1345182 : dst_v[left] = src_v[left];
122 1345182 : dst_u[right] = src_u[right];
123 1345182 : dst_v[right] = src_v[right];
124 1345182 : if (nz > 1) {
125 234464 : dst_w[left] = src_w[left];
126 234464 : dst_w[right] = src_w[right];
127 : }
128 : }
129 : }
130 :
131 : /* Front and back z-faces (k = 0 and k = nz-1) */
132 23296 : if (nz > 1) {
133 1218 : size_t back_base = (nz - 1) * plane;
134 16526 : for (size_t j = 1; j < ny - 1; j++) {
135 219156 : for (size_t i = 1; i < nx - 1; i++) {
136 203848 : size_t off = IDX_2D(i, j, nx);
137 : /* k = 0 face */
138 203848 : dst_u[off] = src_u[off];
139 203848 : dst_v[off] = src_v[off];
140 203848 : dst_w[off] = src_w[off];
141 : /* k = nz-1 face */
142 203848 : dst_u[back_base + off] = src_u[back_base + off];
143 203848 : dst_v[back_base + off] = src_v[back_base + off];
144 203848 : dst_w[back_base + off] = src_w[back_base + off];
145 : }
146 : }
147 : }
148 23296 : }
149 :
150 : #endif /* BOUNDARY_COPY_UTILS_H */
|
