Line data Source code
1 : /**
2 : * @file linear_solver.c
3 : * @brief Core linear solver implementation
4 : *
5 : * Implements:
6 : * - Default parameter functions
7 : * - Backend selection
8 : * - Solver lifecycle (create, init, destroy)
9 : * - Common solve loop
10 : * - Legacy poisson_solve() wrapper
11 : */
12 :
13 : #include "cfd/solvers/poisson_solver.h"
14 : #include "linear_solver_internal.h"
15 :
16 : #include "cfd/boundary/boundary_conditions.h"
17 : #include "cfd/core/indexing.h"
18 : #include "cfd/core/logging.h"
19 : #include "cfd/core/memory.h"
20 :
21 : #include <math.h>
22 : #include <stdio.h>
23 : #include <stdlib.h>
24 : #include <string.h>
25 :
26 : #ifdef _WIN32
27 : #define WIN32_LEAN_AND_MEAN
28 : #include <windows.h>
29 : #else
30 : #include <sys/time.h>
31 : #endif
32 :
33 : /* ============================================================================
34 : * DEFAULT PARAMETERS
35 : * ============================================================================ */
36 :
37 465 : poisson_solver_params_t poisson_solver_params_default(void) {
38 465 : poisson_solver_params_t params;
39 465 : params.tolerance = 1e-6;
40 465 : params.absolute_tolerance = 1e-10;
41 465 : params.max_iterations = 5000; /* Increased from 1000 for CG on fine grids */
42 465 : params.omega = 0.0; /* Auto-compute optimal omega for grid dimensions */
43 465 : params.check_interval = 1;
44 465 : params.verbose = false;
45 465 : params.preconditioner = POISSON_PRECOND_NONE;
46 465 : return params;
47 : }
48 :
49 6478 : poisson_solver_stats_t poisson_solver_stats_default(void) {
50 6478 : poisson_solver_stats_t stats;
51 6478 : stats.status = POISSON_ERROR;
52 6478 : stats.iterations = 0;
53 6478 : stats.initial_residual = 0.0;
54 6478 : stats.final_residual = 0.0;
55 6478 : stats.elapsed_time_ms = 0.0;
56 6478 : return stats;
57 : }
58 :
59 : /* ============================================================================
60 : * TIMING
61 : * ============================================================================ */
62 :
63 25818 : double poisson_solver_get_time_ms(void) {
64 : #ifdef _WIN32
65 : LARGE_INTEGER freq, counter;
66 : QueryPerformanceFrequency(&freq);
67 : QueryPerformanceCounter(&counter);
68 : return (double)counter.QuadPart * 1000.0 / (double)freq.QuadPart;
69 : #else
70 25818 : struct timeval tv;
71 25818 : gettimeofday(&tv, NULL);
72 25818 : return tv.tv_sec * 1000.0 + tv.tv_usec / 1000.0;
73 : #endif
74 : }
75 :
76 : /* ============================================================================
77 : * BACKEND SELECTION
78 : * ============================================================================ */
79 :
80 : static poisson_solver_backend_t g_default_backend = POISSON_BACKEND_AUTO;
81 :
82 2 : poisson_solver_backend_t poisson_solver_get_backend(void) {
83 2 : return g_default_backend;
84 : }
85 :
86 2 : const char* poisson_solver_get_backend_name(void) {
87 2 : switch (g_default_backend) {
88 : case POISSON_BACKEND_SCALAR: return "scalar";
89 0 : case POISSON_BACKEND_OMP: return "omp";
90 0 : case POISSON_BACKEND_SIMD: return "simd";
91 0 : case POISSON_BACKEND_GPU: return "gpu";
92 1 : case POISSON_BACKEND_AUTO:
93 1 : default: return "auto";
94 : }
95 : }
96 :
97 2 : bool poisson_solver_set_backend(poisson_solver_backend_t backend) {
98 2 : if (!poisson_solver_backend_available(backend)) {
99 : return false;
100 : }
101 2 : g_default_backend = backend;
102 2 : return true;
103 : }
104 :
105 43 : bool poisson_solver_backend_available(poisson_solver_backend_t backend) {
106 43 : switch (backend) {
107 : case POISSON_BACKEND_AUTO:
108 : case POISSON_BACKEND_SCALAR:
109 : return true;
110 :
111 : case POISSON_BACKEND_OMP:
112 : #ifdef CFD_ENABLE_OPENMP
113 : return true;
114 : #else
115 : return false;
116 : #endif
117 :
118 16 : case POISSON_BACKEND_SIMD:
119 16 : return poisson_solver_simd_backend_available();
120 :
121 : case POISSON_BACKEND_GPU:
122 : #ifdef CFD_HAS_CUDA
123 : return true;
124 : #else
125 : return false;
126 : #endif
127 :
128 : default:
129 : return false;
130 : }
131 : }
132 :
133 : /* ============================================================================
134 : * SOLVER LIFECYCLE
135 : * ============================================================================ */
136 :
137 : /**
138 : * Auto-select best available backend
139 : *
140 : * Priority: SIMD (runtime detection) > Scalar
141 : */
142 2 : static poisson_solver_backend_t select_best_backend(void) {
143 : /* Prefer SIMD with runtime detection */
144 2 : if (poisson_solver_simd_backend_available()) {
145 0 : return POISSON_BACKEND_SIMD;
146 : }
147 : return POISSON_BACKEND_SCALAR;
148 : }
149 :
150 319 : poisson_solver_t* poisson_solver_create(
151 : poisson_solver_method_t method,
152 : poisson_solver_backend_t backend)
153 : {
154 : /* Auto-select backend if requested */
155 319 : if (backend == POISSON_BACKEND_AUTO) {
156 2 : backend = select_best_backend();
157 : }
158 :
159 : /* Create appropriate solver - no silent fallbacks */
160 319 : switch (method) {
161 28 : case POISSON_METHOD_JACOBI:
162 28 : switch (backend) {
163 1 : case POISSON_BACKEND_SIMD:
164 1 : return create_jacobi_simd_solver();
165 27 : case POISSON_BACKEND_SCALAR:
166 27 : return create_jacobi_scalar_solver();
167 : default:
168 : return NULL; /* Requested backend not available for Jacobi */
169 : }
170 :
171 25 : case POISSON_METHOD_SOR:
172 : case POISSON_METHOD_GAUSS_SEIDEL:
173 25 : switch (backend) {
174 1 : case POISSON_BACKEND_SIMD:
175 1 : return create_sor_simd_solver();
176 24 : case POISSON_BACKEND_SCALAR:
177 24 : return create_sor_scalar_solver();
178 : default:
179 : return NULL; /* SOR not available for OMP/GPU */
180 : }
181 :
182 28 : case POISSON_METHOD_REDBLACK_SOR:
183 28 : switch (backend) {
184 1 : case POISSON_BACKEND_SIMD:
185 1 : return create_redblack_simd_solver();
186 : #ifdef CFD_ENABLE_OPENMP
187 5 : case POISSON_BACKEND_OMP:
188 5 : return create_redblack_omp_solver();
189 : #endif
190 22 : case POISSON_BACKEND_SCALAR:
191 22 : return create_redblack_scalar_solver();
192 : default:
193 : return NULL; /* Requested backend not available for Red-Black */
194 : }
195 :
196 222 : case POISSON_METHOD_CG:
197 222 : switch (backend) {
198 17 : case POISSON_BACKEND_SIMD:
199 17 : return create_cg_simd_solver();
200 : #ifdef CFD_ENABLE_OPENMP
201 19 : case POISSON_BACKEND_OMP:
202 19 : return create_cg_omp_solver();
203 : #endif
204 186 : case POISSON_BACKEND_SCALAR:
205 186 : return create_cg_scalar_solver();
206 : default:
207 : return NULL; /* Requested backend not available for CG */
208 : }
209 :
210 15 : case POISSON_METHOD_BICGSTAB:
211 15 : switch (backend) {
212 1 : case POISSON_BACKEND_SIMD:
213 1 : return create_bicgstab_simd_solver();
214 12 : case POISSON_BACKEND_SCALAR:
215 12 : return create_bicgstab_scalar_solver();
216 : default:
217 : return NULL; /* Requested backend not available for BiCGSTAB */
218 : }
219 :
220 : case POISSON_METHOD_MULTIGRID:
221 : /* Not yet implemented */
222 : return NULL;
223 :
224 : default:
225 : return NULL;
226 : }
227 : }
228 :
229 279 : cfd_status_t poisson_solver_init(
230 : poisson_solver_t* solver,
231 : size_t nx, size_t ny, size_t nz,
232 : double dx, double dy, double dz,
233 : const poisson_solver_params_t* params)
234 : {
235 279 : if (!solver) {
236 : return CFD_ERROR_INVALID;
237 : }
238 :
239 : /* Require at least one interior cell in each active dimension.
240 : * For 2D (nz <= 1): nx, ny >= 3. For 3D (nz > 1): nx, ny, nz >= 3.
241 : * Rejects degenerate grids (e.g. nz==2) where k_start==k_end. */
242 278 : if (nx < 3 || ny < 3 || (nz > 1 && nz < 3)) {
243 : return CFD_ERROR_INVALID;
244 : }
245 :
246 275 : solver->nx = nx;
247 275 : solver->ny = ny;
248 275 : solver->nz = nz;
249 275 : solver->dx = dx;
250 275 : solver->dy = dy;
251 275 : solver->dz = dz;
252 :
253 275 : if (params) {
254 240 : solver->params = *params;
255 : } else {
256 35 : solver->params = poisson_solver_params_default();
257 : }
258 :
259 : /* Adjust max_iterations for Jacobi (needs more iterations) */
260 275 : if (solver->method == POISSON_METHOD_JACOBI && params == NULL) {
261 2 : solver->params.max_iterations = 2000;
262 : }
263 :
264 : /* Call solver-specific init if provided */
265 275 : if (solver->init) {
266 275 : return solver->init(solver, nx, ny, nz, dx, dy, dz, &solver->params);
267 : }
268 :
269 : return CFD_SUCCESS;
270 : }
271 :
272 297 : void poisson_solver_destroy(poisson_solver_t* solver) {
273 297 : if (!solver) {
274 : return;
275 : }
276 :
277 295 : if (solver->destroy) {
278 295 : solver->destroy(solver);
279 : }
280 :
281 295 : cfd_free(solver);
282 : }
283 :
284 : /* ============================================================================
285 : * SOLVER OPERATIONS
286 : * ============================================================================ */
287 :
288 185532 : double poisson_solver_compute_residual(
289 : poisson_solver_t* solver,
290 : const double* x,
291 : const double* rhs)
292 : {
293 185532 : if (!solver || !x || !rhs) {
294 : return -1.0;
295 : }
296 :
297 185529 : size_t nx = solver->nx;
298 185529 : size_t ny = solver->ny;
299 185529 : double dx2 = solver->dx * solver->dx;
300 185529 : double dy2 = solver->dy * solver->dy;
301 185529 : double inv_dz2 = poisson_solver_compute_inv_dz2(solver->dz);
302 :
303 185529 : size_t stride_z, k_start, k_end;
304 185529 : poisson_solver_compute_3d_bounds(solver->nz, nx, ny,
305 : &stride_z, &k_start, &k_end);
306 :
307 185529 : double max_residual = 0.0;
308 :
309 408522 : for (size_t k = k_start; k < k_end; k++) {
310 6968316 : for (size_t j = 1; j < ny - 1; j++) {
311 259652064 : for (size_t i = 1; i < nx - 1; i++) {
312 252906741 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
313 :
314 : /* Compute Laplacian: d^2x/dx^2 + d^2x/dy^2 + d^2x/dz^2 */
315 252906741 : double laplacian =
316 252906741 : (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) / dx2
317 252906741 : + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) / dy2
318 252906741 : + (x[idx + stride_z] + x[idx - stride_z]
319 252906741 : - 2.0 * x[idx]) * inv_dz2;
320 :
321 252906741 : double residual = fabs(laplacian - rhs[idx]);
322 252906741 : if (residual > max_residual) {
323 4235945 : max_residual = residual;
324 : }
325 : }
326 : }
327 : }
328 :
329 : return max_residual;
330 : }
331 :
332 198228 : void poisson_solver_apply_bc(
333 : poisson_solver_t* solver,
334 : double* x)
335 : {
336 198228 : if (!solver || !x) {
337 : return;
338 : }
339 :
340 198227 : if (solver->apply_bc) {
341 10329 : solver->apply_bc(solver, x);
342 10329 : return;
343 : }
344 :
345 : /* Default: Neumann BCs (zero gradient) on all faces */
346 187898 : size_t nx = solver->nx;
347 187898 : size_t ny = solver->ny;
348 187898 : size_t nz = solver->nz;
349 187898 : size_t plane_size = nx * ny;
350 :
351 : /* Z-face Neumann: copy adjacent interior plane to boundary planes */
352 187898 : if (nz > 1) {
353 2 : memcpy(x, x + plane_size, plane_size * sizeof(double));
354 2 : memcpy(x + (nz - 1) * plane_size,
355 2 : x + (nz - 2) * plane_size,
356 : plane_size * sizeof(double));
357 : }
358 :
359 : /* Apply 2D Neumann BCs on each z-plane using solver's own backend.
360 : * This avoids BC_BACKEND_AUTO selecting a different backend (e.g. OMP)
361 : * than the solver itself, which could spawn unexpected threads. */
362 375810 : for (size_t k = 0; k < nz; k++) {
363 187912 : double* plane = x + k * plane_size;
364 187912 : switch (solver->backend) {
365 462 : case POISSON_BACKEND_OMP:
366 462 : bc_apply_scalar_omp(plane, nx, ny, BC_TYPE_NEUMANN);
367 462 : break;
368 0 : case POISSON_BACKEND_SIMD:
369 0 : bc_apply_scalar_simd(plane, nx, ny, BC_TYPE_NEUMANN);
370 0 : break;
371 187450 : default:
372 187450 : bc_apply_scalar_cpu(plane, nx, ny, BC_TYPE_NEUMANN);
373 187450 : break;
374 : }
375 : }
376 : }
377 :
378 : /**
379 : * Common solve loop used by all solvers
380 : */
381 45 : cfd_status_t poisson_solver_solve_common(
382 : poisson_solver_t* solver,
383 : double* x,
384 : double* x_temp,
385 : const double* rhs,
386 : poisson_solver_stats_t* stats)
387 : {
388 45 : if (!solver || !x || !rhs) {
389 : return CFD_ERROR_INVALID;
390 : }
391 :
392 45 : if (!solver->iterate) {
393 : return CFD_ERROR_UNSUPPORTED;
394 : }
395 :
396 45 : poisson_solver_params_t* params = &solver->params;
397 45 : double start_time = poisson_solver_get_time_ms();
398 :
399 : /* Compute initial residual */
400 45 : double initial_res = poisson_solver_compute_residual(solver, x, rhs);
401 45 : double tolerance = params->tolerance * initial_res;
402 :
403 : /* Ensure minimum absolute tolerance */
404 45 : if (tolerance < params->absolute_tolerance) {
405 9 : tolerance = params->absolute_tolerance;
406 : }
407 :
408 45 : if (stats) {
409 45 : stats->initial_residual = initial_res;
410 : }
411 :
412 : /* Already converged? */
413 45 : if (initial_res < params->absolute_tolerance) {
414 9 : if (stats) {
415 9 : stats->status = POISSON_CONVERGED;
416 9 : stats->iterations = 0;
417 9 : stats->final_residual = initial_res;
418 9 : stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time;
419 : }
420 9 : return CFD_SUCCESS;
421 : }
422 :
423 20076 : int converged = 0;
424 : int iter;
425 : double res = initial_res;
426 :
427 20076 : for (iter = 0; iter < params->max_iterations; iter++) {
428 20074 : double new_res = 0.0;
429 :
430 : /* Perform one iteration */
431 20074 : cfd_status_t status = solver->iterate(solver, x, x_temp, rhs,
432 20074 : (iter % params->check_interval == 0) ? &new_res : NULL);
433 :
434 20074 : if (status != CFD_SUCCESS) {
435 0 : if (stats) {
436 0 : stats->status = POISSON_ERROR;
437 0 : stats->iterations = iter + 1;
438 0 : stats->final_residual = res;
439 0 : stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time;
440 : }
441 0 : return status;
442 : }
443 :
444 : /* Check convergence at intervals */
445 20074 : if (iter % params->check_interval == 0) {
446 20074 : res = new_res;
447 :
448 20074 : if (params->verbose) {
449 0 : CFD_LOG_DEBUG("poisson", "Iter %d: residual = %.6e", iter, res);
450 : }
451 :
452 20074 : if (res < tolerance || res < params->absolute_tolerance) {
453 34 : converged = 1;
454 34 : break;
455 : }
456 : }
457 : }
458 :
459 36 : double end_time = poisson_solver_get_time_ms();
460 :
461 36 : if (stats) {
462 36 : stats->iterations = iter + 1;
463 36 : stats->final_residual = res;
464 36 : stats->elapsed_time_ms = end_time - start_time;
465 36 : stats->status = converged ? POISSON_CONVERGED : POISSON_MAX_ITER;
466 : }
467 :
468 36 : return converged ? CFD_SUCCESS : CFD_ERROR_MAX_ITER;
469 : }
470 :
471 6477 : cfd_status_t poisson_solver_solve(
472 : poisson_solver_t* solver,
473 : double* x,
474 : double* x_temp,
475 : const double* rhs,
476 : poisson_solver_stats_t* stats)
477 : {
478 6477 : if (!solver) {
479 : return CFD_ERROR_INVALID;
480 : }
481 :
482 : /* Use solver-specific solve if provided, otherwise common loop */
483 6477 : if (solver->solve) {
484 6432 : double start_time = poisson_solver_get_time_ms();
485 6432 : cfd_status_t status = solver->solve(solver, x, x_temp, rhs, stats);
486 6432 : if (stats) {
487 6432 : stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time;
488 : }
489 6432 : return status;
490 : }
491 :
492 45 : return poisson_solver_solve_common(solver, x, x_temp, rhs, stats);
493 : }
494 :
495 165400 : cfd_status_t poisson_solver_iterate(
496 : poisson_solver_t* solver,
497 : double* x,
498 : double* x_temp,
499 : const double* rhs,
500 : double* residual)
501 : {
502 165400 : if (!solver || !x || !rhs) {
503 : return CFD_ERROR_INVALID;
504 : }
505 :
506 165400 : if (!solver->iterate) {
507 : return CFD_ERROR_UNSUPPORTED;
508 : }
509 :
510 165400 : return solver->iterate(solver, x, x_temp, rhs, residual);
511 : }
512 :
513 : /* ============================================================================
514 : * CACHED SOLVER INSTANCES
515 : * ============================================================================ */
516 :
517 : /*
518 : * Cached solver instances for poisson_solve() convenience API.
519 : * Avoids creating/destroying solvers on each call.
520 : */
521 : static poisson_solver_t* g_cached_jacobi_simd = NULL;
522 : static poisson_solver_t* g_cached_sor = NULL;
523 : static poisson_solver_t* g_cached_sor_simd = NULL;
524 : static poisson_solver_t* g_cached_redblack_simd = NULL;
525 : static poisson_solver_t* g_cached_redblack_omp = NULL;
526 : static poisson_solver_t* g_cached_redblack_scalar = NULL;
527 : static poisson_solver_t* g_cached_cg_scalar = NULL;
528 : static poisson_solver_t* g_cached_cg_omp = NULL;
529 : static poisson_solver_t* g_cached_cg_simd = NULL;
530 :
531 : /**
532 : * Cleanup cached solvers (called at program exit)
533 : */
534 10 : static void cleanup_cached_solvers(void) {
535 10 : if (g_cached_jacobi_simd) {
536 0 : poisson_solver_destroy(g_cached_jacobi_simd);
537 0 : g_cached_jacobi_simd = NULL;
538 : }
539 10 : if (g_cached_sor) {
540 1 : poisson_solver_destroy(g_cached_sor);
541 1 : g_cached_sor = NULL;
542 : }
543 10 : if (g_cached_sor_simd) {
544 0 : poisson_solver_destroy(g_cached_sor_simd);
545 0 : g_cached_sor_simd = NULL;
546 : }
547 10 : if (g_cached_redblack_simd) {
548 0 : poisson_solver_destroy(g_cached_redblack_simd);
549 0 : g_cached_redblack_simd = NULL;
550 : }
551 10 : if (g_cached_redblack_omp) {
552 0 : poisson_solver_destroy(g_cached_redblack_omp);
553 0 : g_cached_redblack_omp = NULL;
554 : }
555 10 : if (g_cached_redblack_scalar) {
556 0 : poisson_solver_destroy(g_cached_redblack_scalar);
557 0 : g_cached_redblack_scalar = NULL;
558 : }
559 10 : if (g_cached_cg_scalar) {
560 9 : poisson_solver_destroy(g_cached_cg_scalar);
561 9 : g_cached_cg_scalar = NULL;
562 : }
563 10 : if (g_cached_cg_omp) {
564 5 : poisson_solver_destroy(g_cached_cg_omp);
565 5 : g_cached_cg_omp = NULL;
566 : }
567 10 : if (g_cached_cg_simd) {
568 0 : poisson_solver_destroy(g_cached_cg_simd);
569 0 : g_cached_cg_simd = NULL;
570 : }
571 10 : }
572 :
573 6268 : int poisson_solve_3d(
574 : double* p, double* p_temp, const double* rhs,
575 : size_t nx, size_t ny, size_t nz,
576 : double dx, double dy, double dz,
577 : poisson_solver_type solver_type)
578 : {
579 6268 : poisson_solver_t** solver_ptr;
580 6268 : poisson_solver_method_t method;
581 6268 : poisson_solver_backend_t backend;
582 :
583 6268 : switch (solver_type) {
584 : case POISSON_SOLVER_JACOBI_SIMD:
585 : solver_ptr = &g_cached_jacobi_simd;
586 : method = POISSON_METHOD_JACOBI;
587 : backend = POISSON_BACKEND_SIMD;
588 : break;
589 :
590 1 : case POISSON_SOLVER_REDBLACK_SIMD:
591 1 : solver_ptr = &g_cached_redblack_simd;
592 1 : method = POISSON_METHOD_REDBLACK_SOR;
593 1 : backend = POISSON_BACKEND_SIMD;
594 1 : break;
595 :
596 0 : case POISSON_SOLVER_REDBLACK_OMP:
597 0 : solver_ptr = &g_cached_redblack_omp;
598 0 : method = POISSON_METHOD_REDBLACK_SOR;
599 0 : backend = POISSON_BACKEND_OMP;
600 0 : break;
601 :
602 1 : case POISSON_SOLVER_SOR_SCALAR:
603 1 : solver_ptr = &g_cached_sor;
604 1 : method = POISSON_METHOD_SOR;
605 1 : backend = POISSON_BACKEND_SCALAR;
606 1 : break;
607 :
608 0 : case POISSON_SOLVER_REDBLACK_SCALAR:
609 0 : solver_ptr = &g_cached_redblack_scalar;
610 0 : method = POISSON_METHOD_REDBLACK_SOR;
611 0 : backend = POISSON_BACKEND_SCALAR;
612 0 : break;
613 :
614 6123 : case POISSON_SOLVER_CG_SCALAR:
615 6123 : solver_ptr = &g_cached_cg_scalar;
616 6123 : method = POISSON_METHOD_CG;
617 6123 : backend = POISSON_BACKEND_SCALAR;
618 6123 : break;
619 :
620 0 : case POISSON_SOLVER_CG_SIMD:
621 0 : solver_ptr = &g_cached_cg_simd;
622 0 : method = POISSON_METHOD_CG;
623 0 : backend = POISSON_BACKEND_SIMD;
624 0 : break;
625 :
626 142 : case POISSON_SOLVER_CG_OMP:
627 142 : solver_ptr = &g_cached_cg_omp;
628 142 : method = POISSON_METHOD_CG;
629 142 : backend = POISSON_BACKEND_OMP;
630 142 : break;
631 :
632 0 : case POISSON_SOLVER_SOR_SIMD:
633 0 : solver_ptr = &g_cached_sor_simd;
634 0 : method = POISSON_METHOD_SOR;
635 0 : backend = POISSON_BACKEND_SIMD;
636 0 : break;
637 :
638 0 : default:
639 0 : CFD_LOG_ERROR("poisson", "poisson_solve_3d: Unknown solver type %d", solver_type);
640 0 : return -1;
641 : }
642 :
643 : /* Recreate solver if grid dimensions or spacing changed.
644 : * Compare against the solver's own stored dimensions, not shared globals,
645 : * because each solver type is cached independently. */
646 6268 : if (*solver_ptr == NULL
647 6251 : || (*solver_ptr)->nx != nx || (*solver_ptr)->ny != ny
648 6235 : || (*solver_ptr)->nz != nz
649 6235 : || (*solver_ptr)->dx != dx || (*solver_ptr)->dy != dy
650 6235 : || (*solver_ptr)->dz != dz) {
651 : /* Register cleanup on first use */
652 33 : static int cleanup_registered = 0;
653 33 : if (!cleanup_registered) {
654 10 : atexit(cleanup_cached_solvers);
655 10 : cleanup_registered = 1;
656 : }
657 :
658 : /* Destroy old solver if exists */
659 33 : if (*solver_ptr) {
660 16 : poisson_solver_destroy(*solver_ptr);
661 16 : *solver_ptr = NULL;
662 : }
663 :
664 : /* Create new solver */
665 33 : *solver_ptr = poisson_solver_create(method, backend);
666 :
667 33 : if (*solver_ptr) {
668 31 : poisson_solver_init(*solver_ptr, nx, ny, nz, dx, dy, dz, NULL);
669 : }
670 : }
671 :
672 6268 : if (!*solver_ptr) {
673 : return -1;
674 : }
675 :
676 6266 : poisson_solver_stats_t stats = poisson_solver_stats_default();
677 6266 : cfd_status_t status = poisson_solver_solve(*solver_ptr, p, p_temp, rhs, &stats);
678 :
679 6266 : if (status == CFD_SUCCESS && stats.status == POISSON_CONVERGED) {
680 6266 : return stats.iterations;
681 : }
682 :
683 : /* Return iterations even if not converged (legacy behavior) */
684 : if (stats.iterations > 0) {
685 : return -1; /* Signal non-convergence */
686 : }
687 :
688 : return -1;
689 : }
690 :
691 3 : int poisson_solve(
692 : double* p, double* p_temp, const double* rhs,
693 : size_t nx, size_t ny, double dx, double dy,
694 : poisson_solver_type solver_type)
695 : {
696 3 : return poisson_solve_3d(p, p_temp, rhs, nx, ny, 1, dx, dy, 0.0, solver_type);
697 : }
698 :
699 : /* Direct solver functions - delegate to unified interface */
700 0 : int poisson_solve_sor_scalar(
701 : double* p, const double* rhs,
702 : size_t nx, size_t ny, double dx, double dy)
703 : {
704 : /* SOR doesn't need temp buffer, pass NULL */
705 0 : return poisson_solve(p, NULL, rhs, nx, ny, dx, dy, POISSON_SOLVER_SOR_SCALAR);
706 : }
707 :
708 : /* SIMD functions with runtime CPU detection */
709 0 : int poisson_solve_jacobi_simd(
710 : double* p, double* p_temp, const double* rhs,
711 : size_t nx, size_t ny, double dx, double dy)
712 : {
713 0 : return poisson_solve(p, p_temp, rhs, nx, ny, dx, dy, POISSON_SOLVER_JACOBI_SIMD);
714 : }
715 :
716 0 : int poisson_solve_redblack_simd(
717 : double* p, double* p_temp, const double* rhs,
718 : size_t nx, size_t ny, double dx, double dy)
719 : {
720 0 : return poisson_solve(p, p_temp, rhs, nx, ny, dx, dy, POISSON_SOLVER_REDBLACK_SIMD);
721 : }
|
