Line data Source code
1 : #include "cfd/core/grid.h"
2 : #include "cfd/core/cfd_status.h"
3 : #include "cfd/core/memory.h"
4 :
5 :
6 : #include <math.h>
7 : #include <stddef.h>
8 :
9 333 : grid* grid_create(size_t nx, size_t ny, size_t nz,
10 : double xmin, double xmax,
11 : double ymin, double ymax,
12 : double zmin, double zmax) {
13 333 : if (nx == 0 || ny == 0 || nz == 0) {
14 5 : cfd_set_error(CFD_ERROR_INVALID, "grid dimensions must be positive");
15 5 : return NULL;
16 : }
17 328 : if (xmax <= xmin || ymax <= ymin) {
18 7 : cfd_set_error(CFD_ERROR_INVALID, "grid bounds invalid (max must be > min)");
19 7 : return NULL;
20 : }
21 321 : if (nz > 1 && zmax <= zmin) {
22 2 : cfd_set_error(CFD_ERROR_INVALID, "grid z-bounds invalid (zmax must be > zmin when nz > 1)");
23 2 : return NULL;
24 : }
25 :
26 319 : grid* new_grid = (grid*)cfd_malloc(sizeof(grid));
27 319 : if (new_grid == NULL) {
28 : return NULL;
29 : }
30 :
31 319 : new_grid->nx = nx;
32 319 : new_grid->ny = ny;
33 319 : new_grid->nz = nz;
34 319 : new_grid->xmin = xmin;
35 319 : new_grid->xmax = xmax;
36 319 : new_grid->ymin = ymin;
37 319 : new_grid->ymax = ymax;
38 :
39 : // Allocate memory for x/y grid arrays
40 319 : new_grid->x = (double*)cfd_calloc(nx, sizeof(double));
41 319 : new_grid->y = (double*)cfd_calloc(ny, sizeof(double));
42 319 : new_grid->dx = (double*)cfd_calloc(nx - 1, sizeof(double));
43 319 : new_grid->dy = (double*)cfd_calloc(ny - 1, sizeof(double));
44 :
45 319 : if (!new_grid->x || !new_grid->y || !new_grid->dx || !new_grid->dy) {
46 0 : grid_destroy(new_grid);
47 0 : return NULL;
48 : }
49 :
50 : // Set z-dimension fields and precomputed constants
51 319 : if (nz > 1) {
52 32 : new_grid->zmin = zmin;
53 32 : new_grid->zmax = zmax;
54 32 : new_grid->z = (double*)cfd_calloc(nz, sizeof(double));
55 32 : new_grid->dz = (double*)cfd_calloc(nz - 1, sizeof(double));
56 32 : if (!new_grid->z || !new_grid->dz) {
57 0 : grid_destroy(new_grid);
58 0 : return NULL;
59 : }
60 32 : new_grid->stride_z = nx * ny;
61 : // inv_dz2 is set during grid_initialize_uniform/stretched
62 32 : new_grid->inv_dz2 = 0.0;
63 32 : new_grid->k_start = 1;
64 32 : new_grid->k_end = nz - 1;
65 : } else {
66 : // 2D mode: no z allocation, zero precomputed constants
67 287 : new_grid->zmin = 0.0;
68 287 : new_grid->zmax = 0.0;
69 287 : new_grid->z = NULL;
70 287 : new_grid->dz = NULL;
71 287 : new_grid->stride_z = 0;
72 287 : new_grid->inv_dz2 = 0.0;
73 287 : new_grid->k_start = 0;
74 287 : new_grid->k_end = 1;
75 : }
76 :
77 : return new_grid;
78 : }
79 :
80 320 : void grid_destroy(grid* grid) {
81 320 : if (grid != NULL) {
82 319 : cfd_free(grid->x);
83 319 : cfd_free(grid->y);
84 319 : cfd_free(grid->dx);
85 319 : cfd_free(grid->dy);
86 319 : cfd_free(grid->z);
87 319 : cfd_free(grid->dz);
88 319 : cfd_free(grid);
89 : }
90 320 : }
91 :
92 297 : void grid_initialize_uniform(grid* grid) {
93 297 : double dx = (grid->xmax - grid->xmin) / (grid->nx - 1);
94 297 : double dy = (grid->ymax - grid->ymin) / (grid->ny - 1);
95 :
96 : // Set x-coordinates
97 7824 : for (size_t i = 0; i < grid->nx; i++) {
98 7527 : grid->x[i] = grid->xmin + (i * dx);
99 : }
100 :
101 : // Set y-coordinates
102 7731 : for (size_t j = 0; j < grid->ny; j++) {
103 7434 : grid->y[j] = grid->ymin + (j * dy);
104 : }
105 :
106 : // Set cell sizes
107 7527 : for (size_t i = 0; i < grid->nx - 1; i++) {
108 7230 : grid->dx[i] = dx;
109 : }
110 7434 : for (size_t j = 0; j < grid->ny - 1; j++) {
111 7137 : grid->dy[j] = dy;
112 : }
113 :
114 : // Initialize z-direction if 3D
115 297 : if (grid->nz > 1 && grid->z && grid->dz) {
116 29 : double dz_val = (grid->zmax - grid->zmin) / (grid->nz - 1);
117 :
118 330 : for (size_t k = 0; k < grid->nz; k++) {
119 301 : grid->z[k] = grid->zmin + (k * dz_val);
120 : }
121 301 : for (size_t k = 0; k < grid->nz - 1; k++) {
122 272 : grid->dz[k] = dz_val;
123 : }
124 :
125 29 : grid->inv_dz2 = 1.0 / (dz_val * dz_val);
126 : }
127 297 : }
128 :
129 15 : void grid_initialize_stretched(grid* grid, double beta) {
130 : // Tanh stretching: clusters points near both boundaries
131 : // Higher beta = more clustering near boundaries (useful for boundary layers)
132 : //
133 : // Formula: x[i] = xmin + (xmax - xmin) * (1 + tanh(beta * (2*xi - 1)) / tanh(beta)) / 2
134 : // where xi = i / (n-1) goes from 0 to 1
135 : //
136 : // At xi=0: tanh(-beta)/tanh(beta) = -1, so x = xmin
137 : // At xi=1: tanh(+beta)/tanh(beta) = +1, so x = xmax
138 : // At xi=0.5: tanh(0)/tanh(beta) = 0, so x = (xmin + xmax) / 2
139 :
140 : // Handle beta = 0 or very small beta: fall back to uniform grid
141 : // (avoids division by zero since tanh(0) = 0)
142 15 : if (fabs(beta) < 1e-10) {
143 1 : grid_initialize_uniform(grid);
144 1 : return;
145 : }
146 :
147 14 : double tanh_beta = tanh(beta);
148 :
149 : // Initialize x-direction with stretching
150 274 : for (size_t i = 0; i < grid->nx; i++) {
151 260 : double xi = (double)i / (grid->nx - 1);
152 260 : grid->x[i] = grid->xmin + (grid->xmax - grid->xmin) *
153 260 : (1.0 + tanh(beta * (2.0 * xi - 1.0)) / tanh_beta) / 2.0;
154 : }
155 :
156 : // Initialize y-direction with stretching
157 274 : for (size_t j = 0; j < grid->ny; j++) {
158 260 : double eta = (double)j / (grid->ny - 1);
159 260 : grid->y[j] = grid->ymin + (grid->ymax - grid->ymin) *
160 260 : (1.0 + tanh(beta * (2.0 * eta - 1.0)) / tanh_beta) / 2.0;
161 : }
162 :
163 : // Calculate cell sizes
164 260 : for (size_t i = 0; i < grid->nx - 1; i++) {
165 246 : grid->dx[i] = grid->x[i + 1] - grid->x[i];
166 : }
167 260 : for (size_t j = 0; j < grid->ny - 1; j++) {
168 246 : grid->dy[j] = grid->y[j + 1] - grid->y[j];
169 : }
170 :
171 : // Initialize z-direction with stretching if 3D
172 14 : if (grid->nz > 1 && grid->z && grid->dz) {
173 22 : for (size_t k = 0; k < grid->nz; k++) {
174 21 : double zeta = (double)k / (grid->nz - 1);
175 21 : grid->z[k] = grid->zmin + (grid->zmax - grid->zmin) *
176 21 : (1.0 + (tanh(beta * (2.0 * zeta - 1.0)) / tanh_beta)) / 2.0;
177 : }
178 21 : for (size_t k = 0; k < grid->nz - 1; k++) {
179 20 : grid->dz[k] = grid->z[k + 1] - grid->z[k];
180 : }
181 :
182 : // Use minimum dz for inv_dz2 (conservative for CFL)
183 1 : double dz_min = grid->dz[0];
184 20 : for (size_t k = 1; k < grid->nz - 1; k++) {
185 19 : if (grid->dz[k] < dz_min) {
186 1 : dz_min = grid->dz[k];
187 : }
188 : }
189 1 : grid->inv_dz2 = 1.0 / (dz_min * dz_min);
190 : }
191 : }
|
