Line data Source code
1 : /**
2 : * Boundary Conditions - Parameterized Core Implementation
3 : *
4 : * Shared implementation for Neumann, Periodic, and Dirichlet BCs.
5 : * Before including, define:
6 : * BC_CORE_FUNC_PREFIX - function name prefix (scalar or omp)
7 : * BC_CORE_USE_OMP - 1 for OpenMP pragmas, 0 for plain loops
8 : *
9 : * Generates:
10 : * bc_apply_neumann_<prefix>_impl()
11 : * bc_apply_periodic_<prefix>_impl()
12 : * bc_apply_dirichlet_<prefix>_impl()
13 : *
14 : * 3D support: all functions accept (nz, stride_z). When nz <= 1,
15 : * z-face loops are skipped and x/y-face loops run for a single plane.
16 : */
17 :
18 : #include "boundary_conditions_internal.h"
19 : #include "cfd/core/indexing.h"
20 : #include <limits.h>
21 :
22 : /* Token pasting helpers */
23 : #define BC_CORE_PASTE2(a, b) a##b
24 : #define BC_CORE_PASTE(a, b) BC_CORE_PASTE2(a, b)
25 : #define BC_CORE_FN(name) BC_CORE_PASTE(BC_CORE_PASTE(bc_apply_, name), BC_CORE_PASTE(_, BC_CORE_PASTE(BC_CORE_FUNC_PREFIX, _impl)))
26 :
27 : /* Loop variable type and cast for OMP compatibility */
28 : #if BC_CORE_USE_OMP
29 : #define BC_LOOP_VAR int
30 : #define BC_LOOP_LIMIT(n) ((n) > (size_t)INT_MAX ? INT_MAX : (int)(n))
31 : #define BC_OMP_FOR _Pragma("omp parallel for schedule(static)")
32 : #else
33 : #define BC_LOOP_VAR size_t
34 : #define BC_LOOP_LIMIT(n) (n)
35 : #define BC_OMP_FOR
36 : #endif
37 :
38 : /* --------------------------------------------------------------------
39 : * Neumann (zero-gradient)
40 : * -------------------------------------------------------------------- */
41 191824 : void BC_CORE_FN(neumann)(double* field, size_t nx, size_t ny,
42 : size_t nz, size_t stride_z) {
43 191824 : BC_LOOP_VAR j, i;
44 191824 : size_t k;
45 :
46 : /* x-faces (left/right) for each z-plane */
47 383652 : for (k = 0; k < nz; k++) {
48 191828 : size_t base = k * stride_z;
49 4378 : BC_OMP_FOR
50 6762229 : for (j = 0; j < BC_LOOP_LIMIT(ny); j++) {
51 6574779 : field[base + IDX_2D(0, j, nx)] = field[base + IDX_2D(1, j, nx)];
52 6574779 : field[base + IDX_2D(nx - 1, j, nx)] = field[base + IDX_2D(nx - 2, j, nx)];
53 : }
54 : }
55 :
56 : /* y-faces (bottom/top) for each z-plane */
57 383652 : for (k = 0; k < nz; k++) {
58 191828 : size_t base = k * stride_z;
59 191828 : double* bottom_dst = field + base;
60 191828 : double* bottom_src = field + base + nx;
61 191828 : double* top_dst = field + base + ((ny - 1) * nx);
62 191828 : double* top_src = field + base + ((ny - 2) * nx);
63 :
64 4378 : BC_OMP_FOR
65 6757805 : for (i = 0; i < BC_LOOP_LIMIT(nx); i++) {
66 6570355 : bottom_dst[i] = bottom_src[i];
67 6570355 : top_dst[i] = top_src[i];
68 : }
69 : }
70 :
71 : /* z-faces (back/front) - only when nz > 1 */
72 191824 : if (nz > 1) {
73 1 : double* back_dst = field; /* k=0 plane */
74 1 : double* back_src = field + stride_z; /* k=1 plane */
75 1 : double* front_dst = field + ((nz - 1) * stride_z); /* k=nz-1 plane */
76 1 : double* front_src = field + ((nz - 2) * stride_z); /* k=nz-2 plane */
77 1 : size_t plane_size = nx * ny;
78 :
79 1 : BC_OMP_FOR
80 0 : for (i = 0; i < BC_LOOP_LIMIT(plane_size); i++) {
81 0 : back_dst[i] = back_src[i];
82 0 : front_dst[i] = front_src[i];
83 : }
84 : }
85 191824 : }
86 :
87 : /* --------------------------------------------------------------------
88 : * Periodic
89 : * -------------------------------------------------------------------- */
90 26149 : void BC_CORE_FN(periodic)(double* field, size_t nx, size_t ny,
91 : size_t nz, size_t stride_z) {
92 26149 : BC_LOOP_VAR j, i;
93 26149 : size_t k;
94 :
95 : /* x-faces (left/right) for each z-plane */
96 52302 : for (k = 0; k < nz; k++) {
97 26153 : size_t base = k * stride_z;
98 26153 : BC_OMP_FOR
99 0 : for (j = 0; j < BC_LOOP_LIMIT(ny); j++) {
100 0 : field[base + IDX_2D(0, j, nx)] = field[base + IDX_2D(nx - 2, j, nx)];
101 0 : field[base + IDX_2D(nx - 1, j, nx)] = field[base + IDX_2D(1, j, nx)];
102 : }
103 : }
104 :
105 : /* y-faces (bottom/top) for each z-plane */
106 52302 : for (k = 0; k < nz; k++) {
107 26153 : size_t base = k * stride_z;
108 26153 : double* bottom_dst = field + base;
109 26153 : double* bottom_src = field + base + ((ny - 2) * nx);
110 26153 : double* top_dst = field + base + ((ny - 1) * nx);
111 26153 : double* top_src = field + base + nx;
112 :
113 26153 : BC_OMP_FOR
114 0 : for (i = 0; i < BC_LOOP_LIMIT(nx); i++) {
115 0 : bottom_dst[i] = bottom_src[i];
116 0 : top_dst[i] = top_src[i];
117 : }
118 : }
119 :
120 : /* z-faces (back/front) - only when nz > 1 */
121 26149 : if (nz > 1) {
122 1 : double* back_dst = field; /* k=0 plane */
123 1 : double* back_src = field + ((nz - 2) * stride_z); /* k=nz-2 plane */
124 1 : double* front_dst = field + ((nz - 1) * stride_z); /* k=nz-1 plane */
125 1 : double* front_src = field + stride_z; /* k=1 plane */
126 1 : size_t plane_size = nx * ny;
127 :
128 1 : BC_OMP_FOR
129 0 : for (i = 0; i < BC_LOOP_LIMIT(plane_size); i++) {
130 0 : back_dst[i] = back_src[i];
131 0 : front_dst[i] = front_src[i];
132 : }
133 : }
134 26149 : }
135 :
136 : /* --------------------------------------------------------------------
137 : * Dirichlet (fixed value)
138 : * -------------------------------------------------------------------- */
139 7856 : void BC_CORE_FN(dirichlet)(double* field, size_t nx, size_t ny,
140 : size_t nz, size_t stride_z,
141 : const bc_dirichlet_values_t* values) {
142 7856 : BC_LOOP_VAR j, i;
143 7856 : size_t k;
144 7856 : double val_left = values->left;
145 7856 : double val_right = values->right;
146 7856 : double val_bottom = values->bottom;
147 7856 : double val_top = values->top;
148 :
149 : /* x-faces (left/right) for each z-plane */
150 15728 : for (k = 0; k < nz; k++) {
151 7872 : size_t base = k * stride_z;
152 7859 : BC_OMP_FOR
153 425 : for (j = 0; j < BC_LOOP_LIMIT(ny); j++) {
154 412 : field[base + IDX_2D(0, j, nx)] = val_left;
155 412 : field[base + IDX_2D(nx - 1, j, nx)] = val_right;
156 : }
157 : }
158 :
159 : /* y-faces (bottom/top) for each z-plane */
160 15728 : for (k = 0; k < nz; k++) {
161 7872 : size_t base = k * stride_z;
162 7872 : double* bottom_row = field + base;
163 7872 : double* top_row = field + base + ((ny - 1) * nx);
164 :
165 7859 : BC_OMP_FOR
166 425 : for (i = 0; i < BC_LOOP_LIMIT(nx); i++) {
167 412 : bottom_row[i] = val_bottom;
168 412 : top_row[i] = val_top;
169 : }
170 : }
171 :
172 : /* z-faces (back/front) - only when nz > 1 */
173 7856 : if (nz > 1) {
174 4 : double val_back = values->back;
175 4 : double val_front = values->front;
176 4 : double* back_plane = field; /* k=0 plane */
177 4 : double* front_plane = field + ((nz - 1) * stride_z); /* k=nz-1 plane */
178 4 : size_t plane_size = nx * ny;
179 :
180 4 : BC_OMP_FOR
181 0 : for (i = 0; i < BC_LOOP_LIMIT(plane_size); i++) {
182 0 : back_plane[i] = val_back;
183 0 : front_plane[i] = val_front;
184 : }
185 : }
186 7856 : }
187 :
188 : /* Clean up macros to allow re-inclusion with different parameters */
189 : #undef BC_LOOP_VAR
190 : #undef BC_LOOP_LIMIT
191 : #undef BC_OMP_FOR
192 : #undef BC_CORE_FN
193 : #undef BC_CORE_PASTE
194 : #undef BC_CORE_PASTE2
195 : #undef BC_CORE_FUNC_PREFIX
196 : #undef BC_CORE_USE_OMP
|
