Line data Source code
1 : #include "cfd/core/cfd_status.h"
2 : #include "cfd/core/cpu_features.h"
3 : #include "cfd/core/filesystem.h"
4 : #include "cfd/core/gpu_device.h"
5 : #include "cfd/core/grid.h"
6 : #include "cfd/core/logging.h"
7 : #include "cfd/core/memory.h"
8 : #include "cfd/solvers/navier_stokes_solver.h"
9 : #include "cfd/solvers/poisson_solver.h"
10 :
11 :
12 : #ifdef _WIN32
13 : #define WIN32_LEAN_AND_MEAN
14 : #include <windows.h>
15 : #endif
16 :
17 : #include <math.h>
18 : #include <stddef.h>
19 : #include <stdio.h>
20 : #include <stdlib.h>
21 : #include <string.h>
22 :
23 : #ifdef CFD_ENABLE_OPENMP
24 : #include <omp.h>
25 : #endif
26 :
27 : /* Compute max T for stats reporting */
28 8121 : static double compute_max_temperature(const flow_field* field) {
29 8121 : double max_t = 0.0;
30 8121 : if (field->T) {
31 8121 : size_t n = field->nx * field->ny * field->nz;
32 8121 : max_t = field->T[0];
33 11490115 : for (size_t i = 1; i < n; i++) {
34 11481994 : if (field->T[i] > max_t) max_t = field->T[i];
35 : }
36 : }
37 8121 : return max_t;
38 : }
39 :
40 : #ifdef CFD_HAS_CUDA
41 : /* Reject energy-equation params on the GPU backend, which does not yet
42 : * implement the energy equation. The CPU, OMP, and SIMD backends now support
43 : * it, so only the CUDA wrappers (below) still need this guard. */
44 : static cfd_status_t check_energy_unsupported(const ns_solver_params_t* params) {
45 : if (params->alpha > 0.0 || params->beta != 0.0) {
46 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
47 : "Energy equation (alpha/beta) not supported by the GPU backend; "
48 : "use a CPU, OMP, or SIMD solver");
49 : return CFD_ERROR_UNSUPPORTED;
50 : }
51 : return CFD_SUCCESS;
52 : }
53 : #endif /* CFD_HAS_CUDA */
54 :
55 : // Forward declarations for internal solver implementations
56 : // These are not part of the public API
57 : cfd_status_t explicit_euler_impl(flow_field* field, const grid* grid, const ns_solver_params_t* params);
58 : cfd_status_t rk2_impl(flow_field* field, const grid* grid, const ns_solver_params_t* params);
59 : void explicit_euler_optimized_impl(flow_field* field, const grid* grid,
60 : const ns_solver_params_t* params);
61 : #ifdef CFD_ENABLE_OPENMP
62 : cfd_status_t explicit_euler_omp_impl(flow_field* field, const grid* grid,
63 : const ns_solver_params_t* params);
64 : cfd_status_t rk2_omp_impl(flow_field* field, const grid* grid,
65 : const ns_solver_params_t* params);
66 : cfd_status_t solve_projection_method_omp(flow_field* field, const grid* grid,
67 : const ns_solver_params_t* params);
68 : #endif
69 :
70 : // GPU solver functions
71 : cfd_status_t solve_projection_method_gpu(flow_field* field, const grid* grid,
72 : const ns_solver_params_t* params,
73 : const gpu_config_t* config);
74 : int gpu_is_available(void);
75 :
76 : // SIMD solver functions
77 :
78 : cfd_status_t explicit_euler_simd_init(struct NSSolver* solver, const grid* grid,
79 : const ns_solver_params_t* params);
80 : void explicit_euler_simd_destroy(struct NSSolver* solver);
81 :
82 : cfd_status_t explicit_euler_simd_step(struct NSSolver* solver, flow_field* field, const grid* grid,
83 : const ns_solver_params_t* params, ns_solver_stats_t* stats);
84 :
85 :
86 : cfd_status_t projection_simd_init(ns_solver_t* solver, const grid* grid, const ns_solver_params_t* params);
87 : void projection_simd_destroy(ns_solver_t* solver);
88 :
89 : cfd_status_t projection_simd_step(ns_solver_t* solver, flow_field* field, const grid* grid,
90 : const ns_solver_params_t* params, ns_solver_stats_t* stats);
91 :
92 : cfd_status_t rk2_avx2_init(ns_solver_t* solver, const grid* grid,
93 : const ns_solver_params_t* params);
94 : void rk2_avx2_destroy(ns_solver_t* solver);
95 : cfd_status_t rk2_avx2_step(ns_solver_t* solver, flow_field* field, const grid* grid,
96 : const ns_solver_params_t* params, ns_solver_stats_t* stats);
97 : cfd_status_t rk2_avx2_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
98 : const ns_solver_params_t* params, ns_solver_stats_t* stats);
99 :
100 : #ifdef _WIN32
101 : #else
102 : #include <sys/time.h>
103 : #endif
104 :
105 : // Maximum number of registered solver types
106 : #define MAX_REGISTERED_SOLVERS 32
107 :
108 : // Registry entry
109 : typedef struct {
110 : char name[64];
111 : ns_solver_factory_func factory;
112 : const char* description;
113 : ns_solver_backend_t backend; // Backend type for efficient filtering
114 : } solver_registry_entry;
115 :
116 : // ns_solver_registry_t structure
117 : struct NSSolverRegistry {
118 : solver_registry_entry entries[MAX_REGISTERED_SOLVERS];
119 : int count;
120 : };
121 :
122 : // Forward declarations for built-in solver factories
123 : static ns_solver_t* create_explicit_euler_solver(void);
124 : static ns_solver_t* create_rk2_solver(void);
125 : static ns_solver_t* create_rk2_optimized_solver(void);
126 : static ns_solver_t* create_explicit_euler_optimized_solver(void);
127 : static ns_solver_t* create_projection_solver(void);
128 : static ns_solver_t* create_projection_optimized_solver(void);
129 : #ifdef CFD_HAS_CUDA
130 : static ns_solver_t* create_explicit_euler_gpu_solver(void);
131 : static ns_solver_t* create_projection_gpu_solver(void);
132 : #endif
133 : #ifdef CFD_ENABLE_OPENMP
134 : static ns_solver_t* create_explicit_euler_omp_solver(void);
135 : static ns_solver_t* create_projection_omp_solver(void);
136 : static ns_solver_t* create_rk2_omp_solver(void);
137 : #endif
138 :
139 : // External projection method solver functions
140 : extern cfd_status_t solve_projection_method(flow_field* field, const grid* grid,
141 : const ns_solver_params_t* params);
142 : extern void solve_projection_method_optimized(flow_field* field, const grid* grid,
143 : const ns_solver_params_t* params);
144 :
145 : // External GPU solver functions (from solver_gpu.cu or solver_gpu_stub.c)
146 : #include "cfd/core/cfd_status.h"
147 : #include "cfd/core/gpu_device.h"
148 :
149 :
150 : // Helper to get current time in milliseconds
151 16854 : static double get_time_ms(void) {
152 : #ifdef _WIN32
153 : LARGE_INTEGER freq, counter;
154 : QueryPerformanceFrequency(&freq);
155 : QueryPerformanceCounter(&counter);
156 : return (double)counter.QuadPart * 1000.0 / (double)freq.QuadPart;
157 : #else
158 16854 : struct timeval tv;
159 16854 : gettimeofday(&tv, NULL);
160 16854 : return tv.tv_sec * 1000.0 + tv.tv_usec / 1000.0;
161 : #endif
162 : }
163 :
164 : /**
165 : * solver Registry Implementation
166 : */
167 :
168 274 : ns_solver_registry_t* cfd_registry_create(void) {
169 274 : ns_solver_registry_t* registry = (ns_solver_registry_t*)cfd_calloc(1, sizeof(ns_solver_registry_t));
170 274 : return registry;
171 : }
172 :
173 274 : void cfd_registry_destroy(ns_solver_registry_t* registry) {
174 274 : if (registry) {
175 274 : cfd_free(registry);
176 : }
177 274 : }
178 :
179 271 : void cfd_registry_register_defaults(ns_solver_registry_t* registry) {
180 271 : if (!registry) {
181 : return;
182 : }
183 :
184 : // Register built-in solvers
185 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_EXPLICIT_EULER, create_explicit_euler_solver);
186 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_RK2, create_rk2_solver);
187 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_RK2_OPTIMIZED, create_rk2_optimized_solver);
188 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_EXPLICIT_EULER_OPTIMIZED,
189 : create_explicit_euler_optimized_solver);
190 :
191 : // Register projection method solvers
192 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_PROJECTION, create_projection_solver);
193 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_PROJECTION_OPTIMIZED,
194 : create_projection_optimized_solver);
195 :
196 : // Register GPU solvers (requires CUDA)
197 : #ifdef CFD_HAS_CUDA
198 : cfd_registry_register(registry, NS_SOLVER_TYPE_EXPLICIT_EULER_GPU,
199 : create_explicit_euler_gpu_solver);
200 : cfd_registry_register(registry, NS_SOLVER_TYPE_PROJECTION_JACOBI_GPU,
201 : create_projection_gpu_solver);
202 : #endif
203 :
204 : // Register OpenMP solvers
205 : #ifdef CFD_ENABLE_OPENMP
206 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_EXPLICIT_EULER_OMP,
207 : create_explicit_euler_omp_solver);
208 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_PROJECTION_OMP, create_projection_omp_solver);
209 271 : cfd_registry_register(registry, NS_SOLVER_TYPE_RK2_OMP, create_rk2_omp_solver);
210 : #endif
211 : }
212 :
213 : /**
214 : * Infer backend from solver type name based on naming convention.
215 : * Returns the likely backend, or NS_SOLVER_BACKEND_SCALAR as default.
216 : * Note: Defined here (before first use) and also used by cfd_solver_create_checked.
217 : */
218 2481 : static ns_solver_backend_t infer_backend_from_type(const char* type_name) {
219 2481 : if (!type_name) {
220 : return NS_SOLVER_BACKEND_SCALAR;
221 : }
222 :
223 : /* Check for GPU suffix first (most specific) */
224 2481 : if (strstr(type_name, "_gpu") != NULL) {
225 : return NS_SOLVER_BACKEND_CUDA;
226 : }
227 :
228 : /* Check for OMP suffix */
229 2479 : if (strstr(type_name, "_omp") != NULL) {
230 : return NS_SOLVER_BACKEND_OMP;
231 : }
232 :
233 : /* Check for optimized suffix (SIMD) */
234 1664 : if (strstr(type_name, "_optimized") != NULL) {
235 816 : return NS_SOLVER_BACKEND_SIMD;
236 : }
237 :
238 : /* Default to scalar */
239 : return NS_SOLVER_BACKEND_SCALAR;
240 : }
241 :
242 2475 : int cfd_registry_register(ns_solver_registry_t* registry, const char* type_name,
243 : ns_solver_factory_func factory) {
244 2475 : if (!registry || !type_name || !factory) {
245 2 : cfd_set_error(CFD_ERROR_INVALID, "Invalid arguments for solver registration");
246 2 : return -1;
247 : }
248 2473 : if (strlen(type_name) == 0) {
249 1 : cfd_set_error(CFD_ERROR_INVALID, "solver type name cannot be empty");
250 1 : return -1;
251 : }
252 2472 : if (registry->count >= MAX_REGISTERED_SOLVERS) {
253 1 : cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Max registered solvers limit reached");
254 1 : return -1;
255 : }
256 :
257 : /* Infer backend from type name for efficient filtering later */
258 2471 : ns_solver_backend_t backend = infer_backend_from_type(type_name);
259 :
260 : // Check if already registered
261 12723 : for (int i = 0; i < registry->count; i++) {
262 10252 : if (strcmp(registry->entries[i].name, type_name) == 0) {
263 : // Update existing entry
264 0 : registry->entries[i].factory = factory;
265 0 : registry->entries[i].backend = backend;
266 0 : return 0;
267 : }
268 : }
269 :
270 : // Add new entry
271 2471 : snprintf(registry->entries[registry->count].name,
272 : sizeof(registry->entries[registry->count].name), "%s", type_name);
273 2471 : registry->entries[registry->count].factory = factory;
274 2471 : registry->entries[registry->count].backend = backend;
275 2471 : registry->count++;
276 :
277 2471 : return 0;
278 : }
279 :
280 0 : int cfd_registry_unregister(ns_solver_registry_t* registry, const char* type_name) {
281 0 : if (!registry || !type_name) {
282 : return -1;
283 : }
284 :
285 0 : for (int i = 0; i < registry->count; i++) {
286 0 : if (strcmp(registry->entries[i].name, type_name) == 0) {
287 : // Shift remaining entries
288 0 : for (int j = i; j < registry->count - 1; j++) {
289 0 : registry->entries[j] = registry->entries[j + 1];
290 : }
291 0 : registry->count--;
292 0 : return 0;
293 : }
294 : }
295 : return -1;
296 : }
297 :
298 1 : int cfd_registry_list(ns_solver_registry_t* registry, const char** names, int max_count) {
299 1 : if (!registry) {
300 : return 0;
301 : }
302 :
303 1 : int count = (registry->count < max_count) ? registry->count : max_count;
304 1 : if (names) {
305 10 : for (int i = 0; i < count; i++) {
306 9 : names[i] = registry->entries[i].name;
307 : }
308 : }
309 : return registry->count;
310 : }
311 :
312 0 : int cfd_registry_has(ns_solver_registry_t* registry, const char* type_name) {
313 0 : if (!registry || !type_name) {
314 : return 0;
315 : }
316 :
317 0 : for (int i = 0; i < registry->count; i++) {
318 0 : if (strcmp(registry->entries[i].name, type_name) == 0) {
319 : return 1;
320 : }
321 : }
322 : return 0;
323 : }
324 :
325 0 : const char* cfd_registry_get_description(ns_solver_registry_t* registry, const char* type_name) {
326 0 : if (!registry || !type_name) {
327 : return NULL;
328 : }
329 :
330 : // Create a temporary solver to get its description
331 0 : ns_solver_t* solver = cfd_solver_create(registry, type_name);
332 0 : if (solver) {
333 0 : const char* desc = solver->description;
334 0 : solver_destroy(solver);
335 0 : return desc;
336 : }
337 : return NULL;
338 : }
339 :
340 : /**
341 : * solver Creation and Management
342 : */
343 :
344 291 : ns_solver_t* cfd_solver_create(ns_solver_registry_t* registry, const char* type_name) {
345 291 : if (!registry || !type_name) {
346 1 : cfd_set_error(CFD_ERROR_INVALID, "Invalid arguments for solver creation");
347 1 : return NULL;
348 : }
349 :
350 1073 : for (int i = 0; i < registry->count; i++) {
351 1062 : if (strcmp(registry->entries[i].name, type_name) == 0) {
352 279 : ns_solver_t* solver = registry->entries[i].factory();
353 279 : if (!solver) {
354 : /* Factory returned NULL - check if error was already set by the factory.
355 : * Factories should set specific errors:
356 : * - CFD_ERROR_UNSUPPORTED: backend not available (GPU, SIMD, etc.)
357 : * - CFD_ERROR_NOMEM: allocation failed
358 : * If no error was set, assume memory allocation failure. */
359 0 : if (cfd_get_last_status() == CFD_SUCCESS) {
360 0 : cfd_set_error(CFD_ERROR_NOMEM, "Failed to allocate solver");
361 : }
362 : }
363 279 : return solver;
364 : }
365 : }
366 :
367 : /* Solver type not found in registry - use CFD_ERROR_NOT_FOUND
368 : * to distinguish from CFD_ERROR_UNSUPPORTED (backend unavailable) */
369 11 : char error_msg[128];
370 11 : snprintf(error_msg, sizeof(error_msg), "Solver type '%s' not registered", type_name);
371 11 : cfd_set_error(CFD_ERROR_NOT_FOUND, error_msg);
372 11 : return NULL;
373 : }
374 :
375 280 : void solver_destroy(ns_solver_t* solver) {
376 280 : if (!solver) {
377 : return;
378 : }
379 :
380 279 : if (solver->destroy) {
381 279 : solver->destroy(solver);
382 : }
383 279 : cfd_free(solver);
384 : }
385 :
386 :
387 225 : cfd_status_t solver_init(ns_solver_t* solver, const grid* grid, const ns_solver_params_t* params) {
388 225 : if (!solver) {
389 : return CFD_ERROR_INVALID;
390 : }
391 224 : if (!solver->init) {
392 : return CFD_SUCCESS; // Optional
393 : }
394 :
395 224 : return solver->init(solver, grid, params);
396 : }
397 :
398 :
399 8421 : cfd_status_t solver_step(ns_solver_t* solver, flow_field* field, const grid* grid,
400 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
401 8421 : if (!solver || !field || !grid || !params) {
402 : return CFD_ERROR_INVALID;
403 : }
404 8417 : if (!solver->step) {
405 : return CFD_ERROR;
406 : }
407 :
408 8417 : double start_time = get_time_ms();
409 :
410 8417 : cfd_status_t status = solver->step(solver, field, grid, params, stats);
411 8417 : double end_time = get_time_ms();
412 :
413 8417 : if (stats) {
414 8417 : stats->elapsed_time_ms = end_time - start_time;
415 8417 : stats->status = status;
416 : }
417 :
418 : return status;
419 : }
420 :
421 :
422 11 : cfd_status_t solver_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
423 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
424 11 : if (!solver || !field || !grid || !params) {
425 : return CFD_ERROR_INVALID;
426 : }
427 10 : if (!solver->solve) {
428 : return CFD_ERROR;
429 : }
430 :
431 10 : double start_time = get_time_ms();
432 :
433 10 : cfd_status_t status = solver->solve(solver, field, grid, params, stats);
434 10 : double end_time = get_time_ms();
435 :
436 10 : if (stats) {
437 10 : stats->elapsed_time_ms = end_time - start_time;
438 10 : stats->status = status;
439 : }
440 :
441 : return status;
442 : }
443 :
444 1 : void solver_apply_boundary(ns_solver_t* solver, flow_field* field, const grid* grid) {
445 1 : if (!solver || !field || !grid) {
446 : return;
447 : }
448 :
449 0 : if (solver->apply_boundary) {
450 0 : solver->apply_boundary(solver, field, grid);
451 : } else {
452 : // Fall back to default boundary conditions
453 0 : apply_boundary_conditions(field, grid);
454 : }
455 : }
456 :
457 1 : double solver_compute_dt(ns_solver_t* solver, const flow_field* field, const grid* grid,
458 : const ns_solver_params_t* params) {
459 1 : if (!solver || !field || !grid || !params) {
460 : return 0.0;
461 : }
462 :
463 0 : if (solver->compute_dt) {
464 0 : return solver->compute_dt(solver, field, grid, params);
465 : }
466 :
467 : // Default implementation
468 0 : double max_vel = 0.0;
469 0 : double min_dx = grid->dx[0];
470 0 : double min_dy = grid->dy[0];
471 :
472 0 : for (size_t i = 0; i < grid->nx - 1; i++) {
473 0 : if (grid->dx[i] < min_dx) {
474 0 : min_dx = grid->dx[i];
475 : }
476 : }
477 0 : for (size_t j = 0; j < grid->ny - 1; j++) {
478 0 : if (grid->dy[j] < min_dy) {
479 0 : min_dy = grid->dy[j];
480 : }
481 : }
482 :
483 0 : for (size_t i = 0; i < field->nx * field->ny; i++) {
484 0 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
485 0 : if (vel > max_vel) {
486 0 : max_vel = vel;
487 : }
488 : }
489 :
490 0 : if (max_vel < 1e-10) {
491 0 : max_vel = 1.0;
492 : }
493 :
494 0 : double dt = params->cfl * fmin(min_dx, min_dy) / max_vel;
495 0 : return fmin(fmax(dt, 1e-6), 0.01);
496 : }
497 :
498 : /**
499 : * Built-in solver: Explicit Euler
500 : * Wraps the existing solve_navier_stokes function
501 : */
502 :
503 : typedef struct {
504 : int initialized;
505 : } explicit_euler_context;
506 :
507 139 : static cfd_status_t explicit_euler_init(ns_solver_t* solver, const grid* grid,
508 : const ns_solver_params_t* params) {
509 139 : (void)grid;
510 139 : (void)params;
511 :
512 139 : explicit_euler_context* ctx =
513 139 : (explicit_euler_context*)cfd_malloc(sizeof(explicit_euler_context));
514 139 : if (!ctx) {
515 : return CFD_ERROR;
516 : }
517 :
518 139 : ctx->initialized = 1;
519 139 : solver->context = ctx;
520 139 : return CFD_SUCCESS;
521 : }
522 :
523 175 : static void explicit_euler_destroy(ns_solver_t* solver) {
524 175 : if (solver->context) {
525 139 : cfd_free(solver->context);
526 139 : solver->context = NULL;
527 : }
528 175 : }
529 :
530 6206 : static cfd_status_t explicit_euler_step(ns_solver_t* solver, flow_field* field, const grid* grid,
531 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
532 6206 : (void)solver;
533 :
534 6206 : if (field->nx < 3 || field->ny < 3) {
535 : return CFD_ERROR_INVALID;
536 : }
537 :
538 : // Create params with single iteration
539 6206 : ns_solver_params_t step_params = *params;
540 6206 : step_params.max_iter = 1;
541 :
542 6206 : explicit_euler_impl(field, grid, &step_params);
543 :
544 6206 : if (stats) {
545 848 : stats->iterations = 1;
546 :
547 : // Compute max velocity
548 848 : double max_vel = 0.0;
549 848 : double max_p = 0.0;
550 1047079 : for (size_t i = 0; i < field->nx * field->ny; i++) {
551 1046231 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
552 1046231 : if (vel > max_vel) {
553 8160 : max_vel = vel;
554 : }
555 1046231 : if (fabs(field->p[i]) > max_p) {
556 15552 : max_p = fabs(field->p[i]);
557 : }
558 : }
559 848 : stats->max_velocity = max_vel;
560 848 : stats->max_pressure = max_p;
561 1696 : stats->max_temperature = compute_max_temperature(field);
562 : }
563 :
564 : return CFD_SUCCESS;
565 : }
566 :
567 2 : static cfd_status_t explicit_euler_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
568 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
569 2 : (void)solver;
570 :
571 2 : if (field->nx < 3 || field->ny < 3) {
572 : return CFD_ERROR_INVALID;
573 : }
574 :
575 2 : explicit_euler_impl(field, grid, params);
576 :
577 2 : if (stats) {
578 2 : stats->iterations = params->max_iter;
579 :
580 2 : double max_vel = 0.0;
581 2 : double max_p = 0.0;
582 164 : for (size_t i = 0; i < field->nx * field->ny; i++) {
583 162 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
584 162 : if (vel > max_vel) {
585 24 : max_vel = vel;
586 : }
587 162 : if (fabs(field->p[i]) > max_p) {
588 14 : max_p = fabs(field->p[i]);
589 : }
590 : }
591 2 : stats->max_velocity = max_vel;
592 2 : stats->max_pressure = max_p;
593 4 : stats->max_temperature = compute_max_temperature(field);
594 : }
595 :
596 : return CFD_SUCCESS;
597 : }
598 :
599 116 : static ns_solver_t* create_explicit_euler_solver(void) {
600 116 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
601 116 : if (!s) {
602 : return NULL;
603 : }
604 :
605 116 : s->name = NS_SOLVER_TYPE_EXPLICIT_EULER;
606 116 : s->description = "Basic explicit Euler finite difference solver for 2D Navier-Stokes";
607 116 : s->version = "1.0.0";
608 116 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT;
609 116 : s->backend = NS_SOLVER_BACKEND_SCALAR;
610 :
611 116 : s->init = explicit_euler_init;
612 116 : s->destroy = explicit_euler_destroy;
613 116 : s->step = explicit_euler_step;
614 116 : s->solve = explicit_euler_solve;
615 116 : s->apply_boundary = NULL; // Use default
616 116 : s->compute_dt = NULL; // Use default
617 :
618 116 : return s;
619 : }
620 :
621 : /**
622 : * Built-in solver: RK2 (Heun's Method)
623 : * Second-order Runge-Kutta time integration
624 : */
625 :
626 3992 : static cfd_status_t rk2_step(ns_solver_t* solver, flow_field* field, const grid* grid,
627 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
628 3992 : (void)solver;
629 :
630 3992 : if (field->nx < 3 || field->ny < 3) {
631 : return CFD_ERROR_INVALID;
632 : }
633 :
634 3992 : ns_solver_params_t step_params = *params;
635 3992 : step_params.max_iter = 1;
636 :
637 3992 : cfd_status_t status = rk2_impl(field, grid, &step_params);
638 :
639 3992 : if (stats) {
640 634 : stats->iterations = 1;
641 634 : double max_vel = 0.0;
642 634 : double max_p = 0.0;
643 866362 : for (size_t i = 0; i < field->nx * field->ny; i++) {
644 865728 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
645 865728 : if (vel > max_vel) max_vel = vel;
646 865728 : if (fabs(field->p[i]) > max_p) max_p = fabs(field->p[i]);
647 : }
648 634 : stats->max_velocity = max_vel;
649 634 : stats->max_pressure = max_p;
650 1268 : stats->max_temperature = compute_max_temperature(field);
651 : }
652 :
653 : return status;
654 : }
655 :
656 2 : static cfd_status_t rk2_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
657 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
658 2 : (void)solver;
659 :
660 2 : if (field->nx < 3 || field->ny < 3) {
661 : return CFD_ERROR_INVALID;
662 : }
663 :
664 2 : cfd_status_t status = rk2_impl(field, grid, params);
665 :
666 2 : if (stats) {
667 2 : stats->iterations = params->max_iter;
668 2 : double max_vel = 0.0;
669 2 : double max_p = 0.0;
670 164 : for (size_t i = 0; i < field->nx * field->ny; i++) {
671 162 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
672 162 : if (vel > max_vel) max_vel = vel;
673 162 : if (fabs(field->p[i]) > max_p) max_p = fabs(field->p[i]);
674 : }
675 2 : stats->max_velocity = max_vel;
676 2 : stats->max_pressure = max_p;
677 4 : stats->max_temperature = compute_max_temperature(field);
678 : }
679 :
680 : return status;
681 : }
682 :
683 30 : static ns_solver_t* create_rk2_solver(void) {
684 30 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
685 30 : if (!s) {
686 : return NULL;
687 : }
688 :
689 30 : s->name = NS_SOLVER_TYPE_RK2;
690 30 : s->description = "RK2 (Heun's method) - 2nd order explicit time integration";
691 30 : s->version = "1.0.0";
692 30 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT;
693 30 : s->backend = NS_SOLVER_BACKEND_SCALAR;
694 :
695 30 : s->init = explicit_euler_init;
696 30 : s->destroy = explicit_euler_destroy;
697 30 : s->step = rk2_step;
698 30 : s->solve = rk2_solve;
699 30 : s->apply_boundary = NULL;
700 30 : s->compute_dt = NULL;
701 :
702 30 : return s;
703 : }
704 :
705 7 : static ns_solver_t* create_rk2_optimized_solver(void) {
706 7 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
707 7 : if (!s) {
708 : return NULL;
709 : }
710 :
711 7 : s->name = NS_SOLVER_TYPE_RK2_OPTIMIZED;
712 7 : s->description = "AVX2/SIMD + OpenMP RK2 (Heun's method)";
713 7 : s->version = "1.0.0";
714 7 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT |
715 : NS_SOLVER_CAP_SIMD;
716 7 : s->backend = NS_SOLVER_BACKEND_SIMD;
717 :
718 7 : s->init = rk2_avx2_init;
719 7 : s->destroy = rk2_avx2_destroy;
720 7 : s->step = rk2_avx2_step;
721 7 : s->solve = rk2_avx2_solve;
722 7 : s->apply_boundary = NULL;
723 7 : s->compute_dt = NULL;
724 :
725 7 : return s;
726 : }
727 :
728 305 : static cfd_status_t explicit_euler_simd_step_guarded(ns_solver_t* solver, flow_field* field,
729 : const grid* grid, const ns_solver_params_t* params,
730 : ns_solver_stats_t* stats) {
731 305 : return explicit_euler_simd_step(solver, field, grid, params, stats);
732 : }
733 :
734 1 : static cfd_status_t explicit_euler_simd_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
735 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
736 1 : if (!solver || !field || !grid || !params) {
737 : return CFD_ERROR_INVALID;
738 : }
739 21 : for (int i = 0; i < params->max_iter; i++) {
740 20 : cfd_status_t status = explicit_euler_simd_step(solver, field, grid, params, NULL);
741 20 : if (status != CFD_SUCCESS) {
742 0 : return status;
743 : }
744 : }
745 :
746 1 : if (stats) {
747 1 : stats->iterations = params->max_iter;
748 1 : double max_vel = 0.0;
749 1 : double max_p = 0.0;
750 82 : for (size_t i = 0; i < field->nx * field->ny; i++) {
751 81 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
752 81 : if (vel > max_vel) {
753 12 : max_vel = vel;
754 : }
755 81 : if (fabs(field->p[i]) > max_p) {
756 7 : max_p = fabs(field->p[i]);
757 : }
758 : }
759 1 : stats->max_velocity = max_vel;
760 1 : stats->max_pressure = max_p;
761 2 : stats->max_temperature = compute_max_temperature(field);
762 : }
763 : return CFD_SUCCESS;
764 : }
765 :
766 15 : static ns_solver_t* create_explicit_euler_optimized_solver(void) {
767 15 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
768 15 : if (!s) {
769 : return NULL;
770 : }
771 :
772 15 : s->name = NS_SOLVER_TYPE_EXPLICIT_EULER_OPTIMIZED;
773 15 : s->description = "SIMD-optimized explicit Euler solver (AVX2)";
774 15 : s->version = "1.0.0";
775 15 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_SIMD;
776 15 : s->backend = NS_SOLVER_BACKEND_SIMD;
777 :
778 15 : s->init = explicit_euler_simd_init;
779 15 : s->destroy = explicit_euler_simd_destroy;
780 15 : s->step = explicit_euler_simd_step_guarded;
781 15 : s->solve = explicit_euler_simd_solve;
782 15 : s->apply_boundary = NULL;
783 15 : s->compute_dt = NULL;
784 :
785 15 : return s;
786 : }
787 :
788 : /**
789 : * Built-in solver: Projection Method (Chorin's Method)
790 : */
791 :
792 : typedef struct {
793 : int initialized;
794 : } projection_context;
795 :
796 43 : static cfd_status_t projection_init(ns_solver_t* solver, const grid* grid, const ns_solver_params_t* params) {
797 43 : (void)grid;
798 43 : (void)params;
799 43 : projection_context* ctx = (projection_context*)cfd_malloc(sizeof(projection_context));
800 43 : if (!ctx) {
801 : return CFD_ERROR;
802 : }
803 43 : ctx->initialized = 1;
804 43 : solver->context = ctx;
805 43 : return CFD_SUCCESS;
806 : }
807 :
808 63 : static void projection_destroy(ns_solver_t* solver) {
809 63 : if (solver->context) {
810 52 : cfd_free(solver->context);
811 52 : solver->context = NULL;
812 : }
813 63 : }
814 :
815 6083 : static cfd_status_t projection_step(ns_solver_t* solver, flow_field* field, const grid* grid,
816 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
817 6083 : (void)solver;
818 6083 : if (field->nx < 3 || field->ny < 3) {
819 : return CFD_ERROR_INVALID;
820 : }
821 :
822 6083 : ns_solver_params_t step_params = *params;
823 6083 : step_params.max_iter = 1;
824 :
825 6083 : cfd_status_t status = solve_projection_method(field, grid, &step_params);
826 6083 : if (status != CFD_SUCCESS) {
827 : return status;
828 : }
829 :
830 6083 : if (stats) {
831 6083 : stats->iterations = 1;
832 : // Compute max velocity/pressure
833 6083 : double max_vel = 0.0;
834 6083 : double max_p = 0.0;
835 5123077 : for (size_t i = 0; i < field->nx * field->ny; i++) {
836 5116994 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
837 5116994 : if (vel > max_vel) {
838 391455 : max_vel = vel;
839 : }
840 5116994 : if (fabs(field->p[i]) > max_p) {
841 93869 : max_p = fabs(field->p[i]);
842 : }
843 : }
844 6083 : stats->max_velocity = max_vel;
845 6083 : stats->max_pressure = max_p;
846 12166 : stats->max_temperature = compute_max_temperature(field);
847 : }
848 : return CFD_SUCCESS;
849 : }
850 :
851 2 : static cfd_status_t projection_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
852 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
853 2 : (void)solver;
854 2 : if (field->nx < 3 || field->ny < 3) {
855 : return CFD_ERROR_INVALID;
856 : }
857 :
858 2 : cfd_status_t status = solve_projection_method(field, grid, params);
859 2 : if (status != CFD_SUCCESS) {
860 : return status;
861 : }
862 :
863 2 : if (stats) {
864 2 : stats->iterations = params->max_iter;
865 : // Compute max velocity/pressure
866 2 : double max_vel = 0.0;
867 2 : double max_p = 0.0;
868 164 : for (size_t i = 0; i < field->nx * field->ny; i++) {
869 162 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
870 162 : if (vel > max_vel) {
871 14 : max_vel = vel;
872 : }
873 162 : if (fabs(field->p[i]) > max_p) {
874 2 : max_p = fabs(field->p[i]);
875 : }
876 : }
877 2 : stats->max_velocity = max_vel;
878 2 : stats->max_pressure = max_p;
879 4 : stats->max_temperature = compute_max_temperature(field);
880 : }
881 : return CFD_SUCCESS;
882 : }
883 :
884 47 : static ns_solver_t* create_projection_solver(void) {
885 47 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
886 47 : if (!s) {
887 : return NULL;
888 : }
889 :
890 47 : s->name = NS_SOLVER_TYPE_PROJECTION;
891 47 : s->description = "Projection method (Chorin's method)";
892 47 : s->version = "1.0.0";
893 47 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT;
894 47 : s->backend = NS_SOLVER_BACKEND_SCALAR;
895 :
896 47 : s->init = projection_init;
897 47 : s->destroy = projection_destroy;
898 47 : s->step = projection_step;
899 47 : s->solve = projection_solve;
900 47 : s->apply_boundary = NULL;
901 47 : s->compute_dt = NULL;
902 :
903 47 : return s;
904 : }
905 :
906 1 : static cfd_status_t projection_simd_step_guarded(ns_solver_t* solver, flow_field* field,
907 : const grid* grid, const ns_solver_params_t* params,
908 : ns_solver_stats_t* stats) {
909 1 : return projection_simd_step(solver, field, grid, params, stats);
910 : }
911 :
912 0 : static cfd_status_t projection_simd_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
913 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
914 0 : if (!solver || !field || !grid || !params) {
915 : return CFD_ERROR_INVALID;
916 : }
917 : // Use the step function which utilizes the persistent context
918 0 : for (int i = 0; i < params->max_iter; i++) {
919 0 : cfd_status_t status = projection_simd_step(solver, field, grid, params,
920 : NULL); // Pass NULL stats for individual steps
921 0 : if (status != CFD_SUCCESS) {
922 0 : return status;
923 : }
924 : }
925 :
926 0 : if (stats) {
927 0 : stats->iterations = params->max_iter;
928 : // Compute max velocity/pressure
929 0 : double max_vel = 0.0;
930 0 : double max_p = 0.0;
931 0 : for (size_t i = 0; i < field->nx * field->ny; i++) {
932 0 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
933 0 : if (vel > max_vel) {
934 0 : max_vel = vel;
935 : }
936 0 : if (fabs(field->p[i]) > max_p) {
937 0 : max_p = fabs(field->p[i]);
938 : }
939 : }
940 0 : stats->max_velocity = max_vel;
941 0 : stats->max_pressure = max_p;
942 0 : stats->max_temperature = compute_max_temperature(field);
943 : }
944 : return CFD_SUCCESS;
945 : }
946 :
947 19 : static ns_solver_t* create_projection_optimized_solver(void) {
948 19 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
949 19 : if (!s) {
950 : return NULL;
951 : }
952 :
953 19 : s->name = NS_SOLVER_TYPE_PROJECTION_OPTIMIZED;
954 19 : s->description = "SIMD-optimized Projection solver (AVX2)";
955 19 : s->version = "1.0.0";
956 19 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_SIMD;
957 19 : s->backend = NS_SOLVER_BACKEND_SIMD;
958 :
959 19 : s->init = projection_simd_init;
960 19 : s->destroy = projection_simd_destroy;
961 19 : s->step = projection_simd_step_guarded;
962 19 : s->solve = projection_simd_solve;
963 19 : s->apply_boundary = NULL;
964 19 : s->compute_dt = NULL;
965 :
966 19 : return s;
967 : }
968 :
969 : #ifdef CFD_HAS_CUDA
970 : /**
971 : * Built-in solver: GPU-Accelerated Explicit Euler
972 : * Uses CUDA for GPU acceleration (requires CUDA support)
973 : */
974 :
975 : typedef struct {
976 : gpu_solver_context_t* gpu_ctx;
977 : gpu_config_t gpu_config_t;
978 : int use_gpu;
979 : } gpu_solver_wrapper_context;
980 :
981 : static cfd_status_t gpu_euler_init(ns_solver_t* solver, const grid* grid, const ns_solver_params_t* params) {
982 : gpu_solver_wrapper_context* ctx =
983 : (gpu_solver_wrapper_context*)cfd_malloc(sizeof(gpu_solver_wrapper_context));
984 : if (!ctx) {
985 : return CFD_ERROR;
986 : }
987 :
988 : ctx->gpu_config_t = gpu_config_default();
989 : ctx->use_gpu = gpu_should_use(&ctx->gpu_config_t, grid->nx, grid->ny, grid->nz, params->max_iter);
990 : ctx->gpu_ctx = NULL;
991 :
992 : if (ctx->use_gpu) {
993 : ctx->gpu_ctx = gpu_solver_create(grid->nx, grid->ny, grid->nz, &ctx->gpu_config_t);
994 : if (!ctx->gpu_ctx) {
995 : cfd_free(ctx);
996 : CFD_LOG_ERROR("gpu", "GPU Euler init: Failed to create GPU context");
997 : return CFD_ERROR_UNSUPPORTED;
998 : }
999 : }
1000 :
1001 : solver->context = ctx;
1002 : return CFD_SUCCESS;
1003 : }
1004 :
1005 : static void gpu_euler_destroy(ns_solver_t* solver) {
1006 : if (solver->context) {
1007 : gpu_solver_wrapper_context* ctx = (gpu_solver_wrapper_context*)solver->context;
1008 : if (ctx->gpu_ctx) {
1009 : gpu_solver_destroy(ctx->gpu_ctx);
1010 : }
1011 : cfd_free(ctx);
1012 : solver->context = NULL;
1013 : }
1014 : }
1015 :
1016 : static cfd_status_t gpu_euler_step(ns_solver_t* solver, flow_field* field, const grid* grid,
1017 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1018 : gpu_solver_wrapper_context* ctx = (gpu_solver_wrapper_context*)solver->context;
1019 :
1020 : if (field->nx < 3 || field->ny < 3) {
1021 : return CFD_ERROR_INVALID;
1022 : }
1023 : cfd_status_t rc = check_energy_unsupported(params);
1024 : if (rc != CFD_SUCCESS) return rc;
1025 :
1026 : ns_solver_params_t step_params = *params;
1027 : step_params.max_iter = 1;
1028 :
1029 : if (ctx && ctx->use_gpu && ctx->gpu_ctx) {
1030 : // Upload, step, download
1031 : if (gpu_solver_upload(ctx->gpu_ctx, field) == 0) {
1032 : gpu_solver_stats_t gpu_stats;
1033 : if (gpu_solver_step(ctx->gpu_ctx, grid, &step_params, &gpu_stats) == 0) {
1034 : gpu_solver_download(ctx->gpu_ctx, field);
1035 :
1036 : if (stats) {
1037 : stats->iterations = 1;
1038 : stats->elapsed_time_ms = gpu_stats.kernel_time_ms;
1039 : // Compute max velocity/pressure
1040 : double max_vel = 0.0, max_p = 0.0;
1041 : for (size_t i = 0; i < field->nx * field->ny; i++) {
1042 : double vel =
1043 : sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1044 : if (vel > max_vel) {
1045 : max_vel = vel;
1046 : }
1047 : if (fabs(field->p[i]) > max_p) {
1048 : max_p = fabs(field->p[i]);
1049 : }
1050 : }
1051 : stats->max_velocity = max_vel;
1052 : stats->max_pressure = max_p;
1053 : stats->max_temperature = compute_max_temperature(field);
1054 : }
1055 : return CFD_SUCCESS;
1056 : }
1057 : }
1058 : // GPU operation failed
1059 : CFD_LOG_ERROR("gpu", "GPU Euler step: GPU operation failed");
1060 : return CFD_ERROR_INVALID;
1061 : }
1062 :
1063 : // GPU not available - could be: CUDA not compiled in, gpu_should_use() returned false,
1064 : // or GPU initialization failed
1065 : CFD_LOG_ERROR("gpu", "GPU Euler step: GPU solver not initialized");
1066 : return CFD_ERROR_UNSUPPORTED;
1067 : }
1068 :
1069 : static cfd_status_t gpu_euler_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
1070 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1071 : gpu_solver_wrapper_context* ctx = (gpu_solver_wrapper_context*)solver->context;
1072 :
1073 : if (field->nx < 3 || field->ny < 3) {
1074 : return CFD_ERROR_INVALID;
1075 : }
1076 : cfd_status_t rc = check_energy_unsupported(params);
1077 : if (rc != CFD_SUCCESS) return rc;
1078 :
1079 : if (ctx && ctx->use_gpu && ctx->gpu_ctx) {
1080 : // Use full GPU solver
1081 : solve_navier_stokes_gpu(field, grid, params, &ctx->gpu_config_t);
1082 :
1083 : if (stats) {
1084 : stats->iterations = params->max_iter;
1085 : gpu_solver_stats_t gpu_stats = gpu_solver_get_stats(ctx->gpu_ctx);
1086 : stats->elapsed_time_ms = gpu_stats.kernel_time_ms;
1087 :
1088 : double max_vel = 0.0, max_p = 0.0;
1089 : for (size_t i = 0; i < field->nx * field->ny; i++) {
1090 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1091 : if (vel > max_vel) {
1092 : max_vel = vel;
1093 : }
1094 : if (fabs(field->p[i]) > max_p) {
1095 : max_p = fabs(field->p[i]);
1096 : }
1097 : }
1098 : stats->max_velocity = max_vel;
1099 : stats->max_pressure = max_p;
1100 : stats->max_temperature = compute_max_temperature(field);
1101 : }
1102 : return CFD_SUCCESS;
1103 : }
1104 :
1105 : // GPU not available - could be: CUDA not compiled in, gpu_should_use() returned false,
1106 : // or GPU initialization failed
1107 : CFD_LOG_ERROR("gpu", "GPU Euler solve: GPU solver not initialized");
1108 : return CFD_ERROR_UNSUPPORTED;
1109 : }
1110 :
1111 : static ns_solver_t* create_explicit_euler_gpu_solver(void) {
1112 : /* Check if GPU is available at runtime before creating the solver */
1113 : if (!gpu_is_available()) {
1114 : cfd_set_error(CFD_ERROR_UNSUPPORTED, "CUDA GPU not available at runtime");
1115 : return NULL;
1116 : }
1117 :
1118 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
1119 : if (!s) {
1120 : return NULL;
1121 : }
1122 :
1123 : s->name = NS_SOLVER_TYPE_EXPLICIT_EULER_GPU;
1124 : s->description = "GPU-accelerated explicit Euler solver (CUDA)";
1125 : s->version = "1.0.0";
1126 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_GPU;
1127 : s->backend = NS_SOLVER_BACKEND_CUDA;
1128 :
1129 : s->init = gpu_euler_init;
1130 : s->destroy = gpu_euler_destroy;
1131 : s->step = gpu_euler_step;
1132 : s->solve = gpu_euler_solve;
1133 : s->apply_boundary = NULL;
1134 : s->compute_dt = NULL;
1135 :
1136 : return s;
1137 : }
1138 :
1139 : /**
1140 : * Built-in solver: GPU-Accelerated Projection Method
1141 : */
1142 :
1143 : static cfd_status_t gpu_projection_step(ns_solver_t* solver, flow_field* field, const grid* grid,
1144 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1145 : (void)solver;
1146 : /* The GPU projection backend implements the energy equation (advection-
1147 : * diffusion + Boussinesq buoyancy + per-face thermal BCs). Heat-source
1148 : * callbacks and unsupported thermal BC types are rejected inside
1149 : * solve_projection_method_gpu. */
1150 : ns_solver_params_t step_params = *params;
1151 : step_params.max_iter = 1;
1152 : /* Override thresholds to allow single-step GPU execution on small grids */
1153 : gpu_config_t cfg = gpu_config_default();
1154 : cfg.min_grid_size = 1;
1155 : cfg.min_steps = 1;
1156 : cfd_status_t rc = solve_projection_method_gpu(field, grid, &step_params, &cfg);
1157 : if (rc == CFD_SUCCESS && stats) {
1158 : stats->max_temperature = compute_max_temperature(field);
1159 : }
1160 : return rc;
1161 : }
1162 :
1163 : static cfd_status_t gpu_projection_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
1164 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1165 : (void)solver;
1166 : /* Override thresholds to allow GPU execution on small grids */
1167 : gpu_config_t cfg = gpu_config_default();
1168 : cfg.min_grid_size = 1;
1169 : cfg.min_steps = 1;
1170 : cfd_status_t rc = solve_projection_method_gpu(field, grid, params, &cfg);
1171 : if (rc == CFD_SUCCESS && stats) {
1172 : stats->max_temperature = compute_max_temperature(field);
1173 : }
1174 : return rc;
1175 : }
1176 :
1177 : static ns_solver_t* create_projection_gpu_solver(void) {
1178 : /* Check if GPU is available at runtime before creating the solver */
1179 : if (!gpu_is_available()) {
1180 : cfd_set_error(CFD_ERROR_UNSUPPORTED, "CUDA GPU not available at runtime");
1181 : return NULL;
1182 : }
1183 :
1184 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
1185 : if (!s) {
1186 : return NULL;
1187 : }
1188 :
1189 : s->name = NS_SOLVER_TYPE_PROJECTION_JACOBI_GPU;
1190 : s->description = "GPU-accelerated projection method with Jacobi iteration (CUDA)";
1191 : s->version = "1.0.0";
1192 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_GPU;
1193 : s->backend = NS_SOLVER_BACKEND_CUDA;
1194 :
1195 : s->init = NULL;
1196 : s->destroy = NULL;
1197 : s->step = gpu_projection_step;
1198 : s->solve = gpu_projection_solve;
1199 : s->apply_boundary = NULL;
1200 : s->compute_dt = NULL;
1201 :
1202 : return s;
1203 : }
1204 : #endif /* CFD_HAS_CUDA */
1205 :
1206 : #ifdef CFD_ENABLE_OPENMP
1207 : /**
1208 : * Built-in solver: Explicit Euler OpenMP
1209 : */
1210 :
1211 212 : static cfd_status_t explicit_euler_omp_step(ns_solver_t* solver, flow_field* field, const grid* grid,
1212 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1213 212 : (void)solver;
1214 212 : if (field->nx < 3 || field->ny < 3) {
1215 : return CFD_ERROR_INVALID;
1216 : }
1217 212 : ns_solver_params_t step_params = *params;
1218 212 : step_params.max_iter = 1;
1219 :
1220 212 : explicit_euler_omp_impl(field, grid, &step_params);
1221 :
1222 212 : if (stats) {
1223 212 : stats->iterations = 1;
1224 212 : double max_vel = 0.0, max_p = 0.0;
1225 1751380 : for (size_t i = 0; i < field->nx * field->ny; i++) {
1226 1751168 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1227 1751168 : if (vel > max_vel) {
1228 4080 : max_vel = vel;
1229 : }
1230 1751168 : if (fabs(field->p[i]) > max_p) {
1231 6785 : max_p = fabs(field->p[i]);
1232 : }
1233 : }
1234 212 : stats->max_velocity = max_vel;
1235 212 : stats->max_pressure = max_p;
1236 424 : stats->max_temperature = compute_max_temperature(field);
1237 : }
1238 : return CFD_SUCCESS;
1239 : }
1240 :
1241 1 : static cfd_status_t explicit_euler_omp_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
1242 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1243 1 : (void)solver;
1244 1 : if (field->nx < 3 || field->ny < 3) {
1245 : return CFD_ERROR_INVALID;
1246 : }
1247 1 : explicit_euler_omp_impl(field, grid, params);
1248 :
1249 1 : if (stats) {
1250 1 : stats->iterations = params->max_iter;
1251 1 : double max_vel = 0.0, max_p = 0.0;
1252 82 : for (size_t i = 0; i < field->nx * field->ny; i++) {
1253 81 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1254 81 : if (vel > max_vel) {
1255 12 : max_vel = vel;
1256 : }
1257 81 : if (fabs(field->p[i]) > max_p) {
1258 7 : max_p = fabs(field->p[i]);
1259 : }
1260 : }
1261 1 : stats->max_velocity = max_vel;
1262 1 : stats->max_pressure = max_p;
1263 2 : stats->max_temperature = compute_max_temperature(field);
1264 : }
1265 : return CFD_SUCCESS;
1266 : }
1267 :
1268 16 : static ns_solver_t* create_explicit_euler_omp_solver(void) {
1269 16 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
1270 16 : if (!s) {
1271 : return NULL;
1272 : }
1273 :
1274 16 : s->name = NS_SOLVER_TYPE_EXPLICIT_EULER_OMP;
1275 16 : s->description = "OpenMP-parallelized explicit Euler solver";
1276 16 : s->version = "1.0.0";
1277 16 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_PARALLEL;
1278 16 : s->backend = NS_SOLVER_BACKEND_OMP;
1279 :
1280 16 : s->init = explicit_euler_init; // Can reuse existing init
1281 16 : s->destroy = explicit_euler_destroy; // Can reuse existing destroy
1282 16 : s->step = explicit_euler_omp_step;
1283 16 : s->solve = explicit_euler_omp_solve;
1284 16 : s->apply_boundary = NULL;
1285 16 : s->compute_dt = NULL;
1286 :
1287 16 : return s;
1288 : }
1289 :
1290 : /**
1291 : * Built-in NSSolver: RK2 OpenMP
1292 : */
1293 :
1294 212 : static cfd_status_t rk2_omp_step(ns_solver_t* solver, flow_field* field, const grid* grid,
1295 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1296 212 : (void)solver;
1297 212 : if (field->nx < 3 || field->ny < 3) {
1298 : return CFD_ERROR_INVALID;
1299 : }
1300 212 : ns_solver_params_t step_params = *params;
1301 212 : step_params.max_iter = 1;
1302 :
1303 212 : cfd_status_t status = rk2_omp_impl(field, grid, &step_params);
1304 :
1305 212 : if (stats) {
1306 212 : stats->iterations = 1;
1307 212 : double max_vel = 0.0, max_p = 0.0;
1308 212 : ptrdiff_t i;
1309 212 : ptrdiff_t n = (ptrdiff_t)(field->nx * field->ny);
1310 : #if _OPENMP >= 201107
1311 212 : #pragma omp parallel for reduction(max: max_vel, max_p) schedule(static)
1312 : #endif
1313 : for (i = 0; i < n; i++) {
1314 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1315 : if (vel > max_vel) {
1316 : max_vel = vel;
1317 : }
1318 : double ap = fabs(field->p[i]);
1319 : if (ap > max_p) {
1320 : max_p = ap;
1321 : }
1322 : }
1323 212 : stats->max_velocity = max_vel;
1324 212 : stats->max_pressure = max_p;
1325 424 : stats->max_temperature = compute_max_temperature(field);
1326 : }
1327 : return status;
1328 : }
1329 :
1330 1 : static cfd_status_t rk2_omp_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
1331 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1332 1 : (void)solver;
1333 1 : if (field->nx < 3 || field->ny < 3) {
1334 : return CFD_ERROR_INVALID;
1335 : }
1336 1 : cfd_status_t status = rk2_omp_impl(field, grid, params);
1337 :
1338 1 : if (stats) {
1339 1 : stats->iterations = params->max_iter;
1340 1 : double max_vel = 0.0, max_p = 0.0;
1341 1 : ptrdiff_t i;
1342 1 : ptrdiff_t n = (ptrdiff_t)(field->nx * field->ny);
1343 : #if _OPENMP >= 201107
1344 1 : #pragma omp parallel for reduction(max: max_vel, max_p) schedule(static)
1345 : #endif
1346 : for (i = 0; i < n; i++) {
1347 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1348 : if (vel > max_vel) {
1349 : max_vel = vel;
1350 : }
1351 : double ap = fabs(field->p[i]);
1352 : if (ap > max_p) {
1353 : max_p = ap;
1354 : }
1355 : }
1356 1 : stats->max_velocity = max_vel;
1357 1 : stats->max_pressure = max_p;
1358 2 : stats->max_temperature = compute_max_temperature(field);
1359 : }
1360 : return status;
1361 : }
1362 :
1363 13 : static ns_solver_t* create_rk2_omp_solver(void) {
1364 13 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
1365 13 : if (!s) {
1366 : return NULL;
1367 : }
1368 :
1369 13 : s->name = NS_SOLVER_TYPE_RK2_OMP;
1370 13 : s->description = "OpenMP-parallelized RK2 (Heun's method)";
1371 13 : s->version = "1.0.0";
1372 13 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_PARALLEL;
1373 13 : s->backend = NS_SOLVER_BACKEND_OMP;
1374 :
1375 13 : s->init = explicit_euler_init;
1376 13 : s->destroy = explicit_euler_destroy;
1377 13 : s->step = rk2_omp_step;
1378 13 : s->solve = rk2_omp_solve;
1379 13 : s->apply_boundary = NULL;
1380 13 : s->compute_dt = NULL;
1381 :
1382 13 : return s;
1383 : }
1384 :
1385 : /**
1386 : * Built-in NSSolver: Projection OpenMP
1387 : */
1388 :
1389 9 : static cfd_status_t projection_omp_init(ns_solver_t* solver, const grid* grid,
1390 : const ns_solver_params_t* params) {
1391 9 : (void)params;
1392 9 : if (!solver || !grid) {
1393 : return CFD_ERROR_INVALID;
1394 : }
1395 :
1396 : /* OMP projection requires OMP CG Poisson solver.
1397 : * Never fall back to scalar CG — it would serialize the Poisson solve. */
1398 9 : poisson_solver_t* test_solver = poisson_solver_create(
1399 : POISSON_METHOD_CG, POISSON_BACKEND_OMP);
1400 9 : if (!test_solver) {
1401 0 : CFD_LOG_WARNING("projection", "OMP CG Poisson solver not available");
1402 0 : return CFD_ERROR_UNSUPPORTED;
1403 : }
1404 9 : poisson_solver_destroy(test_solver);
1405 :
1406 9 : projection_context* ctx = (projection_context*)cfd_malloc(sizeof(projection_context));
1407 9 : if (!ctx) {
1408 : return CFD_ERROR_NOMEM;
1409 : }
1410 9 : ctx->initialized = 1;
1411 9 : solver->context = ctx;
1412 9 : return CFD_SUCCESS;
1413 : }
1414 :
1415 122 : static cfd_status_t projection_omp_step(ns_solver_t* solver, flow_field* field, const grid* grid,
1416 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1417 122 : (void)solver;
1418 122 : if (field->nx < 3 || field->ny < 3) {
1419 : return CFD_ERROR_INVALID;
1420 : }
1421 122 : ns_solver_params_t step_params = *params;
1422 122 : step_params.max_iter = 1;
1423 :
1424 122 : cfd_status_t status = solve_projection_method_omp(field, grid, &step_params);
1425 122 : if (status != CFD_SUCCESS) {
1426 : return status;
1427 : }
1428 :
1429 122 : if (stats) {
1430 122 : stats->iterations = 1;
1431 122 : double max_vel = 0.0, max_p = 0.0;
1432 122 : ptrdiff_t i;
1433 122 : ptrdiff_t n = (ptrdiff_t)(field->nx * field->ny);
1434 : #if _OPENMP >= 201107
1435 122 : #pragma omp parallel for reduction(max: max_vel, max_p) schedule(static)
1436 : #endif
1437 : for (i = 0; i < n; i++) {
1438 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1439 : if (vel > max_vel) {
1440 : max_vel = vel;
1441 : }
1442 : double ap = fabs(field->p[i]);
1443 : if (ap > max_p) {
1444 : max_p = ap;
1445 : }
1446 : }
1447 122 : stats->max_velocity = max_vel;
1448 122 : stats->max_pressure = max_p;
1449 244 : stats->max_temperature = compute_max_temperature(field);
1450 : }
1451 : return CFD_SUCCESS;
1452 : }
1453 :
1454 1 : static cfd_status_t projection_omp_solve(ns_solver_t* solver, flow_field* field, const grid* grid,
1455 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
1456 1 : (void)solver;
1457 1 : if (field->nx < 3 || field->ny < 3) {
1458 : return CFD_ERROR_INVALID;
1459 : }
1460 1 : cfd_status_t status = solve_projection_method_omp(field, grid, params);
1461 1 : if (status != CFD_SUCCESS) {
1462 : return status;
1463 : }
1464 :
1465 1 : if (stats) {
1466 1 : stats->iterations = params->max_iter;
1467 1 : double max_vel = 0.0, max_p = 0.0;
1468 82 : for (size_t i = 0; i < field->nx * field->ny; i++) {
1469 81 : double vel = sqrt((field->u[i] * field->u[i]) + (field->v[i] * field->v[i]));
1470 81 : if (vel > max_vel) {
1471 7 : max_vel = vel;
1472 : }
1473 81 : if (fabs(field->p[i]) > max_p) {
1474 1 : max_p = fabs(field->p[i]);
1475 : }
1476 : }
1477 1 : stats->max_velocity = max_vel;
1478 1 : stats->max_pressure = max_p;
1479 2 : stats->max_temperature = compute_max_temperature(field);
1480 : }
1481 : return CFD_SUCCESS;
1482 : }
1483 :
1484 16 : static ns_solver_t* create_projection_omp_solver(void) {
1485 16 : ns_solver_t* s = (ns_solver_t*)cfd_calloc(1, sizeof(*s));
1486 16 : if (!s) {
1487 : return NULL;
1488 : }
1489 :
1490 16 : s->name = NS_SOLVER_TYPE_PROJECTION_OMP;
1491 16 : s->description = "OpenMP-parallelized Projection solver";
1492 16 : s->version = "1.0.0";
1493 16 : s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_PARALLEL;
1494 :
1495 16 : s->init = projection_omp_init; // Checks OMP CG availability
1496 16 : s->destroy = projection_destroy;
1497 16 : s->step = projection_omp_step;
1498 16 : s->solve = projection_omp_solve;
1499 16 : s->apply_boundary = NULL;
1500 16 : s->compute_dt = NULL;
1501 16 : s->backend = NS_SOLVER_BACKEND_OMP;
1502 :
1503 16 : return s;
1504 : }
1505 : #endif
1506 :
1507 : //=============================================================================
1508 : // Backend Availability API
1509 : //=============================================================================
1510 :
1511 30 : int cfd_backend_is_available(ns_solver_backend_t backend) {
1512 30 : switch (backend) {
1513 : case NS_SOLVER_BACKEND_SCALAR:
1514 : return 1; // Always available
1515 :
1516 12 : case NS_SOLVER_BACKEND_SIMD:
1517 12 : return cfd_has_simd();
1518 :
1519 : case NS_SOLVER_BACKEND_OMP:
1520 : #ifdef CFD_ENABLE_OPENMP
1521 : return 1;
1522 : #else
1523 : return 0;
1524 : #endif
1525 :
1526 4 : case NS_SOLVER_BACKEND_CUDA:
1527 4 : return gpu_is_available();
1528 :
1529 1 : default:
1530 1 : return 0;
1531 : }
1532 : }
1533 :
1534 8 : const char* cfd_backend_get_name(ns_solver_backend_t backend) {
1535 8 : switch (backend) {
1536 : case NS_SOLVER_BACKEND_SCALAR:
1537 : return "scalar";
1538 2 : case NS_SOLVER_BACKEND_SIMD:
1539 2 : return "simd";
1540 1 : case NS_SOLVER_BACKEND_OMP:
1541 1 : return "openmp";
1542 3 : case NS_SOLVER_BACKEND_CUDA:
1543 3 : return "cuda";
1544 1 : default:
1545 1 : return "unknown";
1546 : }
1547 : }
1548 :
1549 13 : int cfd_registry_list_by_backend(ns_solver_registry_t* registry, ns_solver_backend_t backend,
1550 : const char** names, int max_count) {
1551 13 : if (!registry) {
1552 : return 0;
1553 : }
1554 :
1555 : int count = 0;
1556 :
1557 : /* If names is NULL, just count total matches (discovery mode).
1558 : * Otherwise, fill the array up to max_count and return number filled. */
1559 111 : for (int i = 0; i < registry->count; i++) {
1560 : /* Use the stored backend instead of creating temporary solvers.
1561 : * This is much more efficient and avoids side effects from factory calls. */
1562 100 : if (registry->entries[i].backend == backend) {
1563 28 : if (names == NULL) {
1564 : /* Discovery mode: count all matches */
1565 3 : count++;
1566 25 : } else if (count < max_count) {
1567 : /* Fill mode: add to array if space available */
1568 24 : names[count] = registry->entries[i].name;
1569 24 : count++;
1570 : } else {
1571 : /* Array full, stop filling but don't count more */
1572 : break;
1573 : }
1574 : }
1575 : }
1576 :
1577 : return count;
1578 : }
1579 :
1580 12 : ns_solver_t* cfd_solver_create_checked(ns_solver_registry_t* registry, const char* type_name) {
1581 12 : if (!registry || !type_name) {
1582 2 : cfd_set_error(CFD_ERROR_INVALID, "Invalid arguments for solver creation");
1583 2 : return NULL;
1584 : }
1585 :
1586 : /* Check backend availability BEFORE creating the solver */
1587 10 : ns_solver_backend_t expected_backend = infer_backend_from_type(type_name);
1588 10 : if (!cfd_backend_is_available(expected_backend)) {
1589 2 : const char* backend_name = cfd_backend_get_name(expected_backend);
1590 2 : char error_msg[128];
1591 2 : snprintf(error_msg, sizeof(error_msg),
1592 : "Backend '%s' is not available on this system", backend_name);
1593 2 : cfd_set_error(CFD_ERROR_UNSUPPORTED, error_msg);
1594 2 : return NULL;
1595 : }
1596 :
1597 : /* Now create the solver - backend is available */
1598 8 : ns_solver_t* solver = cfd_solver_create(registry, type_name);
1599 8 : if (!solver) {
1600 : /* Error already set by cfd_solver_create */
1601 : return NULL;
1602 : }
1603 :
1604 : return solver;
1605 : }
|
