Line data Source code
1 : #include "cfd/core/derived_fields.h"
2 : #include "cfd/core/memory.h"
3 : #include "cfd/solvers/navier_stokes_solver.h"
4 :
5 : #include <math.h>
6 : #include <stddef.h>
7 :
8 : #ifdef CFD_ENABLE_OPENMP
9 : #include <omp.h>
10 : #endif
11 :
12 : /* OpenMP version-gated pragmas */
13 : #if defined(_OPENMP) && (_OPENMP >= 201307) /* OMP 4.0 */
14 : # define OMP_FOR_SIMD _Pragma("omp parallel for simd")
15 : #else
16 : # define OMP_FOR_SIMD _Pragma("omp parallel for")
17 : #endif
18 :
19 : // Minimum grid size to benefit from parallelization
20 : // Below this threshold, thread overhead exceeds benefit
21 : #define OMP_THRESHOLD 1000
22 :
23 : //=============================================================================
24 : // DERIVED FIELDS LIFECYCLE
25 : //=============================================================================
26 :
27 37 : derived_fields* derived_fields_create(size_t nx, size_t ny, size_t nz) {
28 37 : derived_fields* derived = (derived_fields*)cfd_calloc(1, sizeof(derived_fields));
29 37 : if (derived) {
30 37 : derived->nx = nx;
31 37 : derived->ny = ny;
32 37 : derived->nz = nz;
33 37 : derived->velocity_magnitude = NULL;
34 : }
35 37 : return derived;
36 : }
37 :
38 38 : void derived_fields_destroy(derived_fields* derived) {
39 38 : if (!derived) {
40 : return;
41 : }
42 :
43 37 : derived_fields_clear(derived);
44 37 : cfd_free(derived);
45 : }
46 :
47 39 : void derived_fields_clear(derived_fields* derived) {
48 39 : if (!derived) {
49 : return;
50 : }
51 :
52 38 : if (derived->velocity_magnitude) {
53 28 : cfd_free(derived->velocity_magnitude);
54 28 : derived->velocity_magnitude = NULL;
55 : }
56 :
57 : // Reset statistics
58 38 : derived->stats_computed = 0;
59 : }
60 :
61 : //=============================================================================
62 : // FIELD STATISTICS
63 : //=============================================================================
64 :
65 146 : field_stats calculate_field_statistics(const double* data, size_t count) {
66 146 : field_stats stats = {0};
67 :
68 146 : if (!data || count == 0) {
69 0 : return stats;
70 : }
71 :
72 : #ifdef CFD_ENABLE_OPENMP
73 146 : if (count >= OMP_THRESHOLD) {
74 : // Parallel reduction for min, max, and sum
75 0 : double min_val = data[0];
76 0 : double max_val = data[0];
77 0 : double sum_val = 0.0;
78 :
79 : // Use signed type for OpenMP compatibility with MSVC
80 0 : long long nn = (long long)count;
81 0 : long long i;
82 : #if _OPENMP >= 201107 /* OMP 3.1: min/max reductions */
83 0 : #pragma omp parallel for reduction(min:min_val) reduction(max:max_val) reduction(+:sum_val)
84 : for (i = 0; i < nn; i++) {
85 : double val = data[i];
86 : if (val < min_val) min_val = val;
87 : if (val > max_val) max_val = val;
88 : sum_val += val;
89 : }
90 : #else
91 : /* OMP < 3.1: only + reduction available; compute min/max serially */
92 : #pragma omp parallel for reduction(+:sum_val)
93 : for (i = 0; i < nn; i++) {
94 : sum_val += data[i];
95 : }
96 : for (i = 0; i < nn; i++) {
97 : if (data[i] < min_val) min_val = data[i];
98 : if (data[i] > max_val) max_val = data[i];
99 : }
100 : #endif
101 :
102 0 : stats.min_val = min_val;
103 0 : stats.max_val = max_val;
104 0 : stats.sum_val = sum_val;
105 : } else {
106 : // Sequential for small arrays
107 146 : stats.min_val = data[0];
108 146 : stats.max_val = data[0];
109 146 : stats.sum_val = 0.0;
110 :
111 14338 : for (size_t i = 0; i < count; i++) {
112 14192 : double val = data[i];
113 14192 : if (val < stats.min_val) {
114 344 : stats.min_val = val;
115 : }
116 14192 : if (val > stats.max_val) {
117 552 : stats.max_val = val;
118 : }
119 14192 : stats.sum_val += val;
120 : }
121 : }
122 : #else
123 : // Sequential implementation when OpenMP not available
124 : stats.min_val = data[0];
125 : stats.max_val = data[0];
126 : stats.sum_val = 0.0;
127 :
128 : for (size_t i = 0; i < count; i++) {
129 : double val = data[i];
130 : if (val < stats.min_val)
131 : stats.min_val = val;
132 : if (val > stats.max_val)
133 : stats.max_val = val;
134 : stats.sum_val += val;
135 : }
136 : #endif
137 :
138 146 : stats.avg_val = stats.sum_val / (double)count;
139 146 : return stats;
140 : }
141 :
142 : //=============================================================================
143 : // DERIVED FIELD COMPUTATION
144 : //=============================================================================
145 :
146 33 : void derived_fields_compute_velocity_magnitude(derived_fields* derived, const flow_field* field) {
147 33 : if (!derived || !field || !field->u || !field->v) {
148 : return;
149 : }
150 :
151 29 : size_t n = derived->nx * derived->ny * derived->nz;
152 :
153 : // Allocate if needed
154 29 : if (!derived->velocity_magnitude) {
155 28 : derived->velocity_magnitude = (double*)cfd_malloc(n * sizeof(double));
156 28 : if (!derived->velocity_magnitude) {
157 : return;
158 : }
159 : }
160 :
161 : // Compute velocity magnitude (embarrassingly parallel)
162 29 : const double* u = field->u;
163 29 : const double* v = field->v;
164 29 : const double* w = field->w;
165 29 : double* vel_mag = derived->velocity_magnitude;
166 :
167 : #ifdef CFD_ENABLE_OPENMP
168 29 : if (n >= OMP_THRESHOLD) {
169 : // Use signed type for OpenMP compatibility with MSVC
170 0 : long long nn = (long long)n;
171 0 : long long i;
172 0 : OMP_FOR_SIMD
173 : for (i = 0; i < nn; i++) {
174 : double w_i = w ? w[i] : 0.0;
175 : vel_mag[i] = sqrt((u[i] * u[i]) + (v[i] * v[i]) + (w_i * w_i));
176 : }
177 : } else {
178 2777 : for (size_t i = 0; i < n; i++) {
179 2748 : double w_i = w ? w[i] : 0.0;
180 2748 : vel_mag[i] = sqrt((u[i] * u[i]) + (v[i] * v[i]) + (w_i * w_i));
181 : }
182 : }
183 : #else
184 : for (size_t i = 0; i < n; i++) {
185 : double w_i = w ? w[i] : 0.0;
186 : vel_mag[i] = sqrt(u[i] * u[i] + v[i] * v[i] + w_i * w_i);
187 : }
188 : #endif
189 : }
190 :
191 21 : void derived_fields_compute_statistics(derived_fields* derived, const flow_field* field) {
192 21 : if (!derived || !field) {
193 : return;
194 : }
195 :
196 21 : size_t count = derived->nx * derived->ny * derived->nz;
197 :
198 : // Compute statistics for primary fields
199 21 : if (field->u) {
200 21 : derived->u_stats = calculate_field_statistics(field->u, count);
201 : }
202 21 : if (field->v) {
203 21 : derived->v_stats = calculate_field_statistics(field->v, count);
204 : }
205 21 : if (field->w) {
206 21 : derived->w_stats = calculate_field_statistics(field->w, count);
207 : }
208 21 : if (field->p) {
209 21 : derived->p_stats = calculate_field_statistics(field->p, count);
210 : }
211 21 : if (field->rho) {
212 21 : derived->rho_stats = calculate_field_statistics(field->rho, count);
213 : }
214 21 : if (field->T) {
215 21 : derived->T_stats = calculate_field_statistics(field->T, count);
216 : }
217 :
218 : // Compute velocity magnitude statistics if available
219 21 : if (derived->velocity_magnitude) {
220 20 : derived->vel_mag_stats = calculate_field_statistics(derived->velocity_magnitude, count);
221 : }
222 :
223 21 : derived->stats_computed = 1;
224 : }
|
