Line data Source code
1 : /**
2 : * Inlet Boundary Conditions - Common Definitions
3 : *
4 : * Shared helper functions and data structures used by all inlet BC
5 : * backend implementations (scalar, OMP, AVX2, NEON).
6 : *
7 : * This header is internal and should only be included by inlet BC
8 : * implementation files.
9 : */
10 :
11 : #ifndef CFD_BOUNDARY_CONDITIONS_INLET_COMMON_H
12 : #define CFD_BOUNDARY_CONDITIONS_INLET_COMMON_H
13 :
14 : #include "boundary_conditions_internal.h"
15 : #include "cfd/core/indexing.h"
16 : #include <math.h>
17 :
18 : /* ============================================================================
19 : * Edge validation and index conversion
20 : *
21 : * bc_edge_t uses bit flags (0x01..0x20) for potential combining.
22 : * We convert these to sequential indices (0-5) for array lookup.
23 : * ============================================================================ */
24 :
25 : /**
26 : * Check if edge value is valid (exactly one of the six valid edges).
27 : */
28 150 : static inline bool bc_inlet_is_valid_edge(bc_edge_t edge) {
29 74 : return edge == BC_EDGE_LEFT || edge == BC_EDGE_RIGHT ||
30 : edge == BC_EDGE_BOTTOM || edge == BC_EDGE_TOP ||
31 : edge == BC_EDGE_FRONT || edge == BC_EDGE_BACK;
32 : }
33 :
34 : /**
35 : * Convert bc_edge_t bit flag to sequential array index (0-5).
36 : * BC_EDGE_LEFT (0x01) -> 0, BC_EDGE_RIGHT (0x02) -> 1,
37 : * BC_EDGE_BOTTOM (0x04) -> 2, BC_EDGE_TOP (0x08) -> 3,
38 : * BC_EDGE_FRONT (0x10) -> 4, BC_EDGE_BACK (0x20) -> 5
39 : *
40 : * Uses bit manipulation to find the position of the least significant set bit
41 : * by counting trailing zeros in the bit flag (clamped to the range 0-5).
42 : * The behavior is only defined for valid, single-bit edges; callers must
43 : * ensure validity via bc_inlet_is_valid_edge() before calling.
44 : */
45 48 : static inline int bc_inlet_edge_to_index(bc_edge_t edge) {
46 48 : int idx = 0;
47 48 : unsigned int e = (unsigned int)edge;
48 173 : while ((e & 1) == 0 && idx < 5) {
49 65 : e >>= 1;
50 65 : idx++;
51 : }
52 60 : return idx;
53 : }
54 :
55 : /* ============================================================================
56 : * Table-driven edge configuration (indexed 0-5)
57 : * ============================================================================ */
58 :
59 : /**
60 : * Mass flow velocity direction by edge.
61 : */
62 : typedef struct {
63 : int u_sign; /* Sign of u velocity for mass flow (+1 or -1, 0 if not u-direction) */
64 : int v_sign; /* Sign of v velocity for mass flow (+1 or -1, 0 if not v-direction) */
65 : int w_sign; /* Sign of w velocity for mass flow (+1 or -1, 0 if not w-direction) */
66 : } bc_inlet_mass_flow_dir_t;
67 :
68 : /**
69 : * Mass flow velocity direction by edge index (0=left, 1=right, 2=bottom, 3=top, 4=front, 5=back).
70 : */
71 : static const bc_inlet_mass_flow_dir_t bc_inlet_mass_flow_dir[6] = {
72 : { .u_sign = +1, .v_sign = 0, .w_sign = 0 }, /* LEFT: Flow into domain (+x) */
73 : { .u_sign = -1, .v_sign = 0, .w_sign = 0 }, /* RIGHT: Flow into domain (-x) */
74 : { .u_sign = 0, .v_sign = +1, .w_sign = 0 }, /* BOTTOM: Flow into domain (+y) */
75 : { .u_sign = 0, .v_sign = -1, .w_sign = 0 }, /* TOP: Flow into domain (-y) */
76 : { .u_sign = 0, .v_sign = 0, .w_sign = -1 }, /* FRONT: Flow into domain (-z) */
77 : { .u_sign = 0, .v_sign = 0, .w_sign = +1 }, /* BACK: Flow into domain (+z) */
78 : };
79 :
80 : /* ============================================================================
81 : * Index computation functions for each edge (2D: within a single xy-plane)
82 : * ============================================================================ */
83 :
84 : typedef size_t (*bc_inlet_idx_func_t)(size_t i, size_t nx, size_t ny);
85 :
86 562 : static inline size_t bc_inlet_idx_left(size_t j, size_t nx, size_t ny) {
87 562 : (void)ny;
88 562 : return IDX_2D(0, j, nx);
89 : }
90 :
91 28 : static inline size_t bc_inlet_idx_right(size_t j, size_t nx, size_t ny) {
92 28 : (void)ny;
93 28 : return IDX_2D(nx - 1, j, nx);
94 : }
95 :
96 157 : static inline size_t bc_inlet_idx_bottom(size_t i, size_t nx, size_t ny) {
97 157 : (void)nx; (void)ny;
98 157 : return i;
99 : }
100 :
101 29 : static inline size_t bc_inlet_idx_top(size_t i, size_t nx, size_t ny) {
102 29 : return IDX_2D(i, ny - 1, nx);
103 : }
104 :
105 : /**
106 : * Edge loop configuration for table-driven boundary application.
107 : * For 2D x/y edges: use_ny_for_count selects the loop count.
108 : * For 3D z-face edges: is_z_face is set and count is nx*ny (full plane).
109 : */
110 : typedef struct {
111 : bc_inlet_idx_func_t idx_fn; /* Index computation function (NULL for z-faces) */
112 : int use_ny_for_count; /* 1 if loop count is ny, 0 if nx (ignored for z-faces) */
113 : int is_z_face; /* 1 for front/back z-face edges */
114 : } bc_inlet_edge_loop_t;
115 :
116 : /**
117 : * Edge loop configuration indexed 0-5 (left, right, bottom, top, front, back).
118 : */
119 : static const bc_inlet_edge_loop_t bc_inlet_edge_loops[6] = {
120 : { .idx_fn = bc_inlet_idx_left, .use_ny_for_count = 1, .is_z_face = 0 }, /* LEFT */
121 : { .idx_fn = bc_inlet_idx_right, .use_ny_for_count = 1, .is_z_face = 0 }, /* RIGHT */
122 : { .idx_fn = bc_inlet_idx_bottom, .use_ny_for_count = 0, .is_z_face = 0 }, /* BOTTOM */
123 : { .idx_fn = bc_inlet_idx_top, .use_ny_for_count = 0, .is_z_face = 0 }, /* TOP */
124 : { .idx_fn = NULL, .use_ny_for_count = 0, .is_z_face = 1 }, /* FRONT */
125 : { .idx_fn = NULL, .use_ny_for_count = 0, .is_z_face = 1 }, /* BACK */
126 : };
127 :
128 : /* ============================================================================
129 : * Velocity computation helpers
130 : * ============================================================================ */
131 :
132 : /**
133 : * Get base velocity components from inlet configuration.
134 : */
135 761 : static inline void bc_inlet_get_base_velocity(const bc_inlet_config_t* config,
136 : double* u_base, double* v_base) {
137 761 : switch (config->spec_type) {
138 677 : case BC_INLET_SPEC_VELOCITY:
139 677 : *u_base = config->spec.velocity.u;
140 677 : *v_base = config->spec.velocity.v;
141 677 : break;
142 :
143 8 : case BC_INLET_SPEC_MAGNITUDE_DIR:
144 8 : *u_base = config->spec.magnitude_dir.magnitude * cos(config->spec.magnitude_dir.direction);
145 8 : *v_base = config->spec.magnitude_dir.magnitude * sin(config->spec.magnitude_dir.direction);
146 8 : break;
147 :
148 76 : case BC_INLET_SPEC_MASS_FLOW: {
149 : /* Validate edge before computing direction-dependent velocity */
150 152 : if (!bc_inlet_is_valid_edge(config->edge)) {
151 0 : *u_base = 0.0;
152 0 : *v_base = 0.0;
153 0 : return;
154 : }
155 : /* For 2D per unit depth: velocity = mass_flow / (density * inlet_length)
156 : * where mass_flow is kg/(s*m) and density*length gives kg/m^2 */
157 76 : double rho_L = config->spec.mass_flow.density * config->spec.mass_flow.inlet_length;
158 76 : if (rho_L <= 0.0) {
159 : /* Invalid density or inlet_length - return zero velocity */
160 28 : *u_base = 0.0;
161 28 : *v_base = 0.0;
162 28 : return;
163 : }
164 48 : double avg_velocity = config->spec.mass_flow.mass_flow_rate / rho_L;
165 48 : int idx = bc_inlet_edge_to_index(config->edge);
166 48 : *u_base = avg_velocity * bc_inlet_mass_flow_dir[idx].u_sign;
167 48 : *v_base = avg_velocity * bc_inlet_mass_flow_dir[idx].v_sign;
168 48 : break;
169 : }
170 :
171 0 : default:
172 0 : *u_base = 0.0;
173 0 : *v_base = 0.0;
174 0 : break;
175 : }
176 : }
177 :
178 : /**
179 : * Apply profile shape to base velocity.
180 : */
181 761 : static inline void bc_inlet_apply_profile(const bc_inlet_config_t* config,
182 : double u_base, double v_base,
183 : double position,
184 : double* u_out, double* v_out) {
185 761 : switch (config->profile) {
186 398 : case BC_INLET_PROFILE_UNIFORM:
187 398 : *u_out = u_base;
188 398 : *v_out = v_base;
189 398 : break;
190 :
191 331 : case BC_INLET_PROFILE_PARABOLIC: {
192 331 : double profile_factor = 4.0 * position * (1.0 - position);
193 331 : *u_out = u_base * profile_factor;
194 331 : *v_out = v_base * profile_factor;
195 331 : break;
196 : }
197 :
198 32 : case BC_INLET_PROFILE_CUSTOM:
199 32 : if (config->custom_profile != NULL) {
200 16 : config->custom_profile(position, u_out, v_out, config->custom_profile_user_data);
201 : } else {
202 16 : *u_out = u_base;
203 16 : *v_out = v_base;
204 : }
205 : break;
206 :
207 0 : default:
208 0 : *u_out = u_base;
209 0 : *v_out = v_base;
210 0 : break;
211 : }
212 761 : }
213 :
214 : /**
215 : * Compute inlet velocity at given position (combines base velocity + profile).
216 : */
217 761 : static inline void bc_inlet_compute_velocity(const bc_inlet_config_t* config,
218 : double position,
219 : double* u_out, double* v_out) {
220 761 : double u_base, v_base;
221 761 : bc_inlet_get_base_velocity(config, &u_base, &v_base);
222 761 : bc_inlet_apply_profile(config, u_base, v_base, position, u_out, v_out);
223 761 : }
224 :
225 : /**
226 : * Compute w-velocity for z-face inlets from mass flow specification.
227 : * For non-mass-flow specs on z-faces, returns 0.0 (u,v are set by compute_velocity).
228 : */
229 1 : static inline double bc_inlet_compute_w(const bc_inlet_config_t* config) {
230 1 : if (config->spec_type == BC_INLET_SPEC_MASS_FLOW) {
231 0 : double rho_L = config->spec.mass_flow.density * config->spec.mass_flow.inlet_length;
232 0 : if (rho_L <= 0.0) return 0.0;
233 0 : double avg_velocity = config->spec.mass_flow.mass_flow_rate / rho_L;
234 0 : int idx = bc_inlet_edge_to_index(config->edge);
235 0 : return avg_velocity * bc_inlet_mass_flow_dir[idx].w_sign;
236 : }
237 : return 0.0;
238 : }
239 :
240 : #endif /* CFD_BOUNDARY_CONDITIONS_INLET_COMMON_H */
|
