Line data Source code
1 : #include "cfd/api/simulation_api.h"
2 : #include "cfd/core/cfd_status.h"
3 : #include "cfd/core/derived_fields.h"
4 : #include "cfd/core/grid.h"
5 : #include "cfd/core/memory.h"
6 : #include "cfd/io/output_registry.h"
7 : #include "cfd/solvers/navier_stokes_solver.h"
8 :
9 :
10 : #include "cfd/core/cfd_init.h"
11 : #include <math.h>
12 : #include <stdio.h>
13 : #include <stdlib.h>
14 : #include <string.h>
15 :
16 : // simulation_data struct is defined in simulation_api.h
17 :
18 : // Default solver type used when none is specified
19 : #define DEFAULT_SOLVER_TYPE NS_SOLVER_TYPE_EXPLICIT_EULER
20 :
21 : // Ensure registry is initialized
22 : // Internal helper to create simulation with a specific solver
23 74 : static simulation_data* create_simulation_with_solver(size_t nx, size_t ny, size_t nz,
24 : double xmin, double xmax,
25 : double ymin, double ymax,
26 : double zmin, double zmax,
27 : const char* solver_type) {
28 : // Lazy initialization of library
29 74 : if (!cfd_is_initialized()) {
30 7 : if (cfd_init() != CFD_SUCCESS) {
31 0 : cfd_set_error(CFD_ERROR, "Failed to initialize CFD library");
32 0 : return NULL;
33 : }
34 : }
35 :
36 74 : if (nx == 0 || ny == 0 || nz == 0) {
37 2 : cfd_set_error(CFD_ERROR_INVALID, "Simulation grid dimensions must be positive");
38 2 : return NULL;
39 : }
40 72 : if (xmax <= xmin || ymax <= ymin || (nz > 1 && zmax <= zmin)) {
41 2 : cfd_set_error(CFD_ERROR_INVALID, "Simulation bounds invalid");
42 2 : return NULL;
43 : }
44 :
45 70 : simulation_data* sim_data = (simulation_data*)cfd_malloc(sizeof(simulation_data));
46 70 : if (!sim_data) {
47 : return NULL;
48 : }
49 :
50 : // Initialize base output directory
51 70 : snprintf(sim_data->output_base_dir, sizeof(sim_data->output_base_dir), "../../artifacts");
52 :
53 : // Create and initialize grid
54 70 : sim_data->grid = grid_create(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax);
55 70 : if (!sim_data->grid) {
56 0 : cfd_free(sim_data);
57 0 : return NULL;
58 : }
59 70 : grid_initialize_uniform(sim_data->grid);
60 :
61 : // Create flow_field
62 70 : sim_data->field = flow_field_create(nx, ny, nz);
63 70 : if (!sim_data->field) {
64 0 : grid_destroy(sim_data->grid);
65 0 : cfd_free(sim_data);
66 0 : return NULL;
67 : }
68 70 : initialize_flow_field(sim_data->field, sim_data->grid);
69 :
70 : // Initialize solver parameters with defaults
71 70 : sim_data->params = ns_solver_params_default();
72 70 : sim_data->params.dt = 0.001;
73 70 : sim_data->params.cfl = 0.2;
74 70 : sim_data->params.mu = 0.01;
75 70 : sim_data->params.max_iter = 1;
76 :
77 : // Initialize stats
78 70 : sim_data->last_stats = ns_solver_stats_default();
79 :
80 : // Create output registry
81 70 : sim_data->outputs = output_registry_create();
82 70 : if (!sim_data->outputs) {
83 0 : flow_field_destroy(sim_data->field);
84 0 : grid_destroy(sim_data->grid);
85 0 : cfd_free(sim_data);
86 0 : return NULL;
87 : }
88 70 : sim_data->run_prefix = NULL;
89 70 : sim_data->current_time = 0.0;
90 :
91 : // Initialize NSSolver Registry
92 70 : sim_data->registry = cfd_registry_create();
93 70 : if (!sim_data->registry) {
94 0 : output_registry_destroy(sim_data->outputs);
95 0 : flow_field_destroy(sim_data->field);
96 0 : grid_destroy(sim_data->grid);
97 0 : cfd_free(sim_data);
98 0 : return NULL;
99 : }
100 70 : cfd_registry_register_defaults(sim_data->registry);
101 :
102 : // Create and initialize the solver
103 70 : sim_data->solver = cfd_solver_create(sim_data->registry, solver_type);
104 70 : if (!sim_data->solver) {
105 : // Failed to create solver - cleanup and return NULL
106 1 : cfd_registry_destroy(sim_data->registry);
107 1 : output_registry_destroy(sim_data->outputs);
108 1 : flow_field_destroy(sim_data->field);
109 1 : grid_destroy(sim_data->grid);
110 1 : cfd_free(sim_data);
111 1 : return NULL;
112 : }
113 :
114 69 : solver_init(sim_data->solver, sim_data->grid, &sim_data->params);
115 :
116 69 : return sim_data;
117 : }
118 :
119 : // Initialize simulation with default solver
120 70 : simulation_data* init_simulation(size_t nx, size_t ny, size_t nz,
121 : double xmin, double xmax,
122 : double ymin, double ymax,
123 : double zmin, double zmax) {
124 70 : return create_simulation_with_solver(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax,
125 : DEFAULT_SOLVER_TYPE);
126 : }
127 :
128 : // Initialize simulation with a specific solver type
129 4 : simulation_data* init_simulation_with_solver(size_t nx, size_t ny, size_t nz,
130 : double xmin, double xmax,
131 : double ymin, double ymax,
132 : double zmin, double zmax,
133 : const char* solver_type) {
134 4 : if (!solver_type) {
135 1 : solver_type = DEFAULT_SOLVER_TYPE;
136 : }
137 4 : return create_simulation_with_solver(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax,
138 : solver_type);
139 : }
140 :
141 : // Set the solver for an existing simulation
142 3 : void simulation_set_solver(simulation_data* sim_data, struct NSSolver* solver) {
143 3 : if (!sim_data || !solver) {
144 2 : cfd_set_error(CFD_ERROR_INVALID, "Invalid arguments for simulation_set_solver");
145 2 : return;
146 : }
147 :
148 : // Destroy existing solver
149 1 : if (sim_data->solver) {
150 1 : solver_destroy(sim_data->solver);
151 : }
152 :
153 1 : sim_data->solver = solver;
154 1 : solver_init(solver, sim_data->grid, &sim_data->params);
155 : }
156 :
157 : // Set the solver by type name
158 6 : int simulation_set_solver_by_name(simulation_data* sim_data, const char* solver_type) {
159 6 : if (!sim_data || !solver_type) {
160 4 : cfd_set_error(CFD_ERROR_INVALID, "Invalid arguments for simulation solver");
161 4 : return -1;
162 : }
163 :
164 2 : struct NSSolver* solver = cfd_solver_create(sim_data->registry, solver_type);
165 2 : if (!solver) {
166 : return -1;
167 : }
168 :
169 1 : simulation_set_solver(sim_data, solver);
170 1 : return 0;
171 : }
172 :
173 : // Get the current solver
174 2 : struct NSSolver* simulation_get_solver(simulation_data* sim_data) {
175 2 : return sim_data ? sim_data->solver : NULL;
176 : }
177 :
178 : // Get statistics from the last solve
179 3 : const ns_solver_stats_t* simulation_get_stats(const simulation_data* sim_data) {
180 3 : return sim_data ? &sim_data->last_stats : NULL;
181 : }
182 :
183 : // Run simulation step
184 21 : cfd_status_t run_simulation_step(simulation_data* sim_data) {
185 21 : if (!sim_data || !sim_data->solver) {
186 : return CFD_ERROR_INVALID;
187 : }
188 :
189 : // Use fixed time step for animation stability
190 20 : sim_data->params.dt = 0.005;
191 :
192 40 : cfd_status_t status = solver_step(sim_data->solver, sim_data->field, sim_data->grid,
193 20 : &sim_data->params, &sim_data->last_stats);
194 20 : if (status != CFD_SUCCESS) {
195 : return status;
196 : }
197 :
198 : // Accumulate simulation time
199 20 : sim_data->current_time += sim_data->params.dt;
200 20 : return CFD_SUCCESS;
201 : }
202 :
203 0 : cfd_status_t run_simulation_solve(simulation_data* sim_data) {
204 0 : if (!sim_data || !sim_data->solver) {
205 : return CFD_ERROR_INVALID;
206 : }
207 :
208 : // Use fixed time step for animation stability
209 0 : sim_data->params.dt = 0.005;
210 :
211 0 : cfd_status_t status = solver_solve(sim_data->solver, sim_data->field, sim_data->grid,
212 0 : &sim_data->params, &sim_data->last_stats);
213 :
214 : // Accumulate simulation time based on iterations performed
215 0 : sim_data->current_time += sim_data->params.dt * sim_data->last_stats.iterations;
216 :
217 0 : return status;
218 : }
219 :
220 : // Free simulation data
221 : // Free simulation data
222 69 : void free_simulation(simulation_data* sim_data) {
223 69 : if (!sim_data) {
224 : return;
225 : }
226 :
227 69 : if (sim_data->solver) {
228 69 : solver_destroy(sim_data->solver);
229 69 : sim_data->solver = NULL;
230 : }
231 :
232 69 : if (sim_data->registry) {
233 69 : cfd_registry_destroy(sim_data->registry);
234 69 : sim_data->registry = NULL;
235 : }
236 :
237 69 : if (sim_data->outputs) {
238 69 : output_registry_destroy(sim_data->outputs);
239 69 : sim_data->outputs = NULL;
240 : }
241 :
242 69 : if (sim_data->run_prefix) {
243 6 : cfd_free(sim_data->run_prefix);
244 6 : sim_data->run_prefix = NULL;
245 : }
246 :
247 69 : flow_field_destroy(sim_data->field);
248 69 : grid_destroy(sim_data->grid);
249 69 : cfd_free(sim_data);
250 : }
251 :
252 : // Static list of available solver names (const strings have static lifetime)
253 : // This avoids use-after-free issues with dynamically created registries
254 : static const char* const s_solver_names[] = {
255 : NS_SOLVER_TYPE_EXPLICIT_EULER,
256 : NS_SOLVER_TYPE_EXPLICIT_EULER_OPTIMIZED,
257 : NS_SOLVER_TYPE_PROJECTION,
258 : NS_SOLVER_TYPE_PROJECTION_OPTIMIZED,
259 : NS_SOLVER_TYPE_EXPLICIT_EULER_GPU,
260 : NS_SOLVER_TYPE_PROJECTION_JACOBI_GPU,
261 : #ifdef CFD_ENABLE_OPENMP
262 : NS_SOLVER_TYPE_EXPLICIT_EULER_OMP,
263 : NS_SOLVER_TYPE_PROJECTION_OMP,
264 : #endif
265 : };
266 :
267 : static const int s_solver_count = (int)(sizeof(s_solver_names) / sizeof(s_solver_names[0]));
268 :
269 : // NSSolver discovery and listing
270 6 : int simulation_list_solvers(const char** names, int max_count) {
271 6 : if (names && max_count > 0) {
272 5 : int fill_count = (s_solver_count < max_count) ? s_solver_count : max_count;
273 39 : for (int i = 0; i < fill_count; i++) {
274 34 : names[i] = s_solver_names[i];
275 : }
276 : }
277 6 : return s_solver_count;
278 : }
279 :
280 3 : int simulation_has_solver(const char* solver_type) {
281 3 : if (!solver_type) {
282 : return 0;
283 : }
284 13 : for (int i = 0; i < s_solver_count; i++) {
285 12 : if (strcmp(s_solver_names[i], solver_type) == 0) {
286 : return 1;
287 : }
288 : }
289 : return 0;
290 : }
291 :
292 :
293 : //=============================================================================
294 : // OUTPUT REGISTRY API
295 : //=============================================================================
296 :
297 : // Register output for automatic generation
298 15 : void simulation_register_output(simulation_data* sim_data, output_field_type field_type,
299 : int interval, const char* prefix) {
300 15 : if (!sim_data || !sim_data->outputs) {
301 2 : cfd_set_error(CFD_ERROR_INVALID, "Invalid simulation data");
302 2 : return;
303 : }
304 13 : output_registry_add(sim_data->outputs, field_type, interval, prefix);
305 : }
306 :
307 : // Clear all registered outputs
308 8 : void simulation_clear_outputs(simulation_data* sim_data) {
309 8 : if (!sim_data || !sim_data->outputs) {
310 : return;
311 : }
312 7 : output_registry_clear(sim_data->outputs);
313 : }
314 :
315 : // Set base output directory
316 2 : void simulation_set_output_dir(simulation_data* sim_data, const char* base_dir) {
317 2 : if (sim_data && base_dir && strlen(base_dir) > 0) {
318 2 : snprintf(sim_data->output_base_dir, sizeof(sim_data->output_base_dir), "%s", base_dir);
319 : }
320 2 : }
321 :
322 : // Set run name prefix
323 10 : void simulation_set_run_prefix(simulation_data* sim_data, const char* prefix) {
324 10 : if (!sim_data) {
325 : return;
326 : }
327 :
328 : // Free existing prefix
329 9 : if (sim_data->run_prefix) {
330 2 : cfd_free(sim_data->run_prefix);
331 2 : sim_data->run_prefix = NULL;
332 : }
333 :
334 : // Set new prefix
335 9 : if (prefix) {
336 8 : size_t len = strlen(prefix) + 1;
337 8 : sim_data->run_prefix = (char*)cfd_malloc(len);
338 8 : if (sim_data->run_prefix) {
339 8 : snprintf(sim_data->run_prefix, len, "%s", prefix);
340 : }
341 : }
342 : }
343 :
344 : //=============================================================================
345 : // AUTOMATIC OUTPUT GENERATION
346 : //=============================================================================
347 :
348 : // Check if any output type that uses velocity magnitude is registered
349 9 : static int needs_velocity_magnitude(const output_registry* outputs) {
350 17 : return output_registry_has_type(outputs, OUTPUT_VELOCITY_MAGNITUDE) ||
351 15 : output_registry_has_type(outputs, OUTPUT_CSV_TIMESERIES) ||
352 23 : output_registry_has_type(outputs, OUTPUT_CSV_CENTERLINE) ||
353 7 : output_registry_has_type(outputs, OUTPUT_CSV_STATISTICS);
354 : }
355 :
356 : // Check if any output type that uses statistics is registered
357 9 : static int needs_statistics(const output_registry* outputs) {
358 17 : return output_registry_has_type(outputs, OUTPUT_CSV_TIMESERIES) ||
359 8 : output_registry_has_type(outputs, OUTPUT_CSV_STATISTICS);
360 : }
361 :
362 : // Automatically write all registered outputs for current step
363 11 : void simulation_write_outputs(simulation_data* sim_data, int step) {
364 11 : if (!sim_data || !sim_data->outputs) {
365 2 : cfd_set_error(CFD_ERROR_INVALID, "Invalid arguments for simulation_write_outputs");
366 2 : return;
367 : }
368 :
369 : // Get run directory (creates it if needed)
370 9 : const char* run_dir =
371 18 : output_registry_get_run_dir(sim_data->outputs, sim_data->output_base_dir,
372 9 : sim_data->run_prefix, sim_data->grid->nx, sim_data->grid->ny);
373 :
374 : // Compute derived fields only when needed
375 9 : derived_fields* derived = NULL;
376 9 : int compute_vel_mag = needs_velocity_magnitude(sim_data->outputs);
377 9 : int compute_stats = needs_statistics(sim_data->outputs);
378 :
379 9 : if (compute_vel_mag || compute_stats) {
380 14 : derived = derived_fields_create(sim_data->grid->nx, sim_data->grid->ny,
381 7 : sim_data->grid->nz);
382 7 : if (derived) {
383 : // Compute velocity magnitude if needed (CSV outputs)
384 7 : if (compute_vel_mag) {
385 7 : derived_fields_compute_velocity_magnitude(derived, sim_data->field);
386 : }
387 : // Compute statistics if needed (CSV timeseries or statistics outputs)
388 7 : if (compute_stats) {
389 6 : derived_fields_compute_statistics(derived, sim_data->field);
390 : }
391 : }
392 : }
393 :
394 : // Write all registered outputs with pre-computed derived fields
395 9 : output_registry_write_outputs(sim_data->outputs, run_dir, step, sim_data->current_time,
396 9 : sim_data->field, derived, sim_data->grid, &sim_data->params,
397 9 : &sim_data->last_stats);
398 :
399 : // Clean up derived fields
400 9 : if (derived) {
401 7 : derived_fields_destroy(derived);
402 : }
403 : }
|
