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