PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
vtk_io.h File Reference
#include "variables.h"
#include "logging.h"
#include "io.h"
Include dependency graph for vtk_io.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

PetscInt CreateVTKFileFromMetadata (const char *filename, const VTKMetaData *meta, MPI_Comm comm)
 Creates and writes a VTK file (either .vts or .vtp) from a populated metadata struct.
 
PetscErrorCode PrepareOutputCoordinates (UserCtx *user, PetscScalar **out_coords, PetscInt *out_nx, PetscInt *out_ny, PetscInt *out_nz, PetscInt *out_npoints)
 Creates a C array of coordinates corresponding to a subsampled (legacy-style) grid.
 
PetscErrorCode PrepareOutputEulerianFieldData (UserCtx *user, Vec field_vec, PetscInt num_components, PetscScalar **out_data)
 Creates a C array of field data corresponding to a subsampled (legacy-style) grid.
 
PetscErrorCode PrepareOutputParticleData (UserCtx *user, PostProcessParams *pps, VTKMetaData *meta, PetscInt *p_n_total)
 Gathers, subsamples, and prepares all particle data for VTK output.
 

Function Documentation

◆ CreateVTKFileFromMetadata()

PetscInt CreateVTKFileFromMetadata ( const char *  filename,
const VTKMetaData meta,
MPI_Comm  comm 
)

Creates and writes a VTK file (either .vts or .vtp) from a populated metadata struct.

Parameters
[in]filenameThe output file name.
[in]metaPointer to a VTKMetaData structure containing all necessary fields.
[in]commThe MPI communicator.
Returns
0 on success, non-zero on failure.

Creates and writes a VTK file (either .vts or .vtp) from a populated metadata struct.

Creates a VTK file from prepared metadata and field payloads.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also
CreateVTKFileFromMetadata()

Definition at line 153 of file vtk_io.c.

154{
155 PetscMPIInt rank;
156 MPI_Comm_rank(comm, &rank);
157 PetscErrorCode ierr = 0;
158
159 if (!rank) {
160 char tmp_filename[PETSC_MAX_PATH_LEN];
161 FILE *fp = NULL;
162
163 ierr = PetscSNPrintf(tmp_filename, sizeof(tmp_filename), "%s.tmp", filename);CHKERRQ(ierr);
164 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 writing combined VTK file '%s'.\n", filename);
165 fp = fopen(tmp_filename, "wb");
166 if (!fp) {
167 LOG_ALLOW(GLOBAL, LOG_ERROR, "fopen failed for %s.\n", tmp_filename);
168 return PETSC_ERR_FILE_OPEN;
169 }
170
171 PetscInt boffset = 0;
172
173 ierr = WriteVTKFileHeader(fp, meta, &boffset);
174 if (ierr) {
175 fclose(fp);
176 remove(tmp_filename);
177 return ierr;
178 }
179
180 if (meta->coords) {
181 ierr = WriteVTKAppendedBlock(fp, meta->coords, 3 * meta->npoints, sizeof(PetscScalar));
182 if (ierr) {
183 fclose(fp);
184 remove(tmp_filename);
185 return ierr;
186 }
187 }
188
189 for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
190 const VTKFieldInfo *field = &meta->point_data_fields[i];
191 if (field->data) {
192 ierr = WriteVTKAppendedBlock(fp, field->data, field->num_components * meta->npoints, sizeof(PetscScalar));
193 if (ierr) {
194 fclose(fp);
195 remove(tmp_filename);
196 return ierr;
197 }
198 }
199 }
200 if (meta->fileType == VTK_POLYDATA) {
201 if (meta->connectivity) {
202 ierr = WriteVTKAppendedBlock(fp, meta->connectivity, meta->npoints, sizeof(PetscInt));
203 if (ierr) {
204 fclose(fp);
205 remove(tmp_filename);
206 return ierr;
207 }
208 }
209 if (meta->offsets) {
210 ierr = WriteVTKAppendedBlock(fp, meta->offsets, meta->npoints, sizeof(PetscInt));
211 if (ierr) {
212 fclose(fp);
213 remove(tmp_filename);
214 return ierr;
215 }
216 }
217 }
218 ierr = WriteVTKFileFooter(fp, meta);
219 if (ierr) {
220 fclose(fp);
221 remove(tmp_filename);
222 return ierr;
223 }
224 if (fclose(fp) != 0) {
225 remove(tmp_filename);
226 return PETSC_ERR_FILE_WRITE;
227 }
228 if (rename(tmp_filename, filename) != 0) {
229 remove(tmp_filename);
230 return PETSC_ERR_FILE_WRITE;
231 }
232 LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 finished writing VTK file '%s'.\n", filename);
233 }
234 return 0;
235}
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
@ LOG_ERROR
Critical errors that may halt the program.
Definition logging.h:28
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:30
PetscInt npoints
Definition variables.h:604
PetscInt num_components
Definition variables.h:591
PetscInt num_point_data_fields
Definition variables.h:607
PetscInt * connectivity
Definition variables.h:608
PetscInt * offsets
Definition variables.h:609
VTKFileType fileType
Definition variables.h:602
PetscScalar * data
Definition variables.h:592
PetscScalar * coords
Definition variables.h:605
VTKFieldInfo point_data_fields[20]
Definition variables.h:606
@ VTK_POLYDATA
Definition variables.h:598
Stores all necessary information for a single data array in a VTK file.
Definition variables.h:589
static PetscErrorCode WriteVTKFileFooter(FILE *fp, const VTKMetaData *meta)
Internal helper implementation: WriteVTKFileFooter().
Definition vtk_io.c:136
static PetscErrorCode WriteVTKAppendedBlock(FILE *fp, const void *data, PetscInt num_elements, size_t element_size)
Internal helper implementation: WriteVTKAppendedBlock().
Definition vtk_io.c:21
static PetscErrorCode WriteVTKFileHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset)
Internal helper implementation: WriteVTKFileHeader().
Definition vtk_io.c:122

◆ PrepareOutputCoordinates()

PetscErrorCode PrepareOutputCoordinates ( UserCtx user,
PetscScalar **  out_coords,
PetscInt *  out_nx,
PetscInt *  out_ny,
PetscInt *  out_nz,
PetscInt *  out_npoints 
)

Creates a C array of coordinates corresponding to a subsampled (legacy-style) grid.

This function gathers the full, distributed grid coordinates onto rank 0. On rank 0, it then allocates a new, smaller C array and copies only the coordinates for the nodes within the range [0..IM-2, 0..JM-2, 0..KM-2]. This produces a contiguous array of points for a grid of size (IM-1)x(JM-1)x(KM-1), matching the legacy output. The output arrays are only allocated and valid on rank 0.

Parameters
[in]userThe UserCtx containing the grid information (DM, IM/JM/KM).
[out]out_coordsOn rank 0, a pointer to the newly allocated C array for coordinate data. NULL on other ranks.
[out]out_nxThe number of points in the x-dimension for the new grid (IM-1).
[out]out_nyThe number of points in the y-dimension for the new grid (JM-1).
[out]out_nzThe number of points in the z-dimension for the new grid (KM-1).
[out]out_npointsThe total number of points in the new grid.
Returns
PetscErrorCode

Creates a C array of coordinates corresponding to a subsampled (legacy-style) grid.

Local to this translation unit.

Definition at line 245 of file vtk_io.c.

246{
247 PetscErrorCode ierr;
248 Vec coords_global;
249 PetscInt IM, JM, KM;
250 PetscInt N_full_coords;
251 PetscScalar *full_coords_arr = NULL;
252
253 PetscFunctionBeginUser;
254 ierr = DMDAGetInfo(user->da, NULL, &IM, &JM, &KM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
255
256 // Set the dimensions of the output grid
257 *out_nx = IM - 1;
258 *out_ny = JM - 1;
259 *out_nz = KM - 1;
260 *out_npoints = (*out_nx) * (*out_ny) * (*out_nz);
261
262 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing subsampled coordinates for a %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT " grid.\n", *out_nx, *out_ny, *out_nz);
263
264 // --- Step 1: Gather the full coordinate vector to rank 0 ---
265 ierr = DMGetCoordinates(user->da, &coords_global); CHKERRQ(ierr);
266 // Reuse of your existing, proven utility function's logic
267 ierr = VecToArrayOnRank0(coords_global, &N_full_coords, &full_coords_arr); CHKERRQ(ierr);
268
269 // --- Step 2: On rank 0, subsample the gathered array ---
270 if (user->simCtx->rank == 0) {
271 // We expect N_full_coords to be IM * JM * KM * 3
272 if (N_full_coords != IM * JM * KM * 3) {
273 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Gathered coordinate array has wrong size. Expected %" PetscInt_FMT ", got %" PetscInt_FMT, IM * JM * KM * 3, N_full_coords);
274 }
275
276 // Allocate the smaller output C array
277 ierr = PetscMalloc1(3 * (*out_npoints), out_coords); CHKERRQ(ierr);
278
279 // Loop over the smaller grid dimensions and copy data
280 PetscInt p_out = 0; // Index for the small output array
281 for (PetscInt k = 0; k < *out_nz; k++) {
282 for (PetscInt j = 0; j < *out_ny; j++) {
283 for (PetscInt i = 0; i < *out_nx; i++) {
284 // Calculate the index in the full, 1D source array
285 PetscInt p_in = 3 * (k * (JM * IM) + j * IM + i);
286
287 (*out_coords)[p_out++] = full_coords_arr[p_in + 0]; // x
288 (*out_coords)[p_out++] = full_coords_arr[p_in + 1]; // y
289 (*out_coords)[p_out++] = full_coords_arr[p_in + 2]; // z
290 }
291 }
292 }
293 // Free the temporary full array that was allocated by VecToArrayOnRank0
294 ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
295 } else {
296 // On other ranks, the output pointer is NULL
297 *out_coords = NULL;
298 }
299
300 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Subsampled coordinates prepared on rank 0 with total %" PetscInt_FMT " points.\n", *out_npoints);
301
302 PetscFunctionReturn(0);
303}
PetscErrorCode VecToArrayOnRank0(Vec inVec, PetscInt *N, double **arrayOut)
Gathers the contents of a distributed PETSc Vec into a single array on rank 0.
Definition io.c:1868
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
PetscMPIInt rank
Definition variables.h:646
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PrepareOutputEulerianFieldData()

PetscErrorCode PrepareOutputEulerianFieldData ( UserCtx user,
Vec  field_vec,
PetscInt  num_components,
PetscScalar **  out_data 
)

Creates a C array of field data corresponding to a subsampled (legacy-style) grid.

This function gathers a full, distributed PETSc vector to rank 0. On rank 0, it then allocates a new, smaller C array and copies only the data components for nodes within the range [0..IM-2, 0..JM-2, 0..KM-2]. This produces a contiguous data array that perfectly matches the point ordering of the subsampled coordinates. The output array is only allocated and valid on rank 0.

Parameters
[in]userThe UserCtx for grid information.
[in]field_vecThe full-sized PETSc vector containing the field data (e.g., user->P_nodal).
[in]num_componentsThe number of components for this field (1 for scalar, 3 for vector).
[out]out_dataOn rank 0, a pointer to the newly allocated C array for the field data. NULL on other ranks.
Returns
PetscErrorCode

Creates a C array of field data corresponding to a subsampled (legacy-style) grid.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also
PrepareOutputEulerianFieldData()

Definition at line 313 of file vtk_io.c.

314{
315 PetscErrorCode ierr;
316 PetscInt IM, JM, KM, nx, ny, nz, npoints;
317 PetscInt N_full_field = 0;
318 PetscScalar *full_field_arr = NULL;
319
320 PetscFunctionBeginUser;
321 ierr = DMDAGetInfo(user->da, NULL, &IM, &JM, &KM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
322 nx = IM - 1; ny = JM - 1; nz = KM - 1;
323 npoints = nx * ny * nz;
324
326 "Preparing subsampled field data (dof=%" PetscInt_FMT ") for a %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT " grid.\n",
327 num_components, nx, ny, nz);
328
329 // --- Step 1: Gather the full field vector to rank 0 ---
330 ierr = VecToArrayOnRank0(field_vec, &N_full_field, &full_field_arr); CHKERRQ(ierr);
331
332 // --- Step 2: On rank 0, subsample the gathered array ---
333 if (user->simCtx->rank == 0) {
334 // Sanity check the size of the gathered array
335 if (N_full_field != IM * JM * KM * num_components) {
336 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
337 "Gathered field array has wrong size. Expected %" PetscInt_FMT ", got %" PetscInt_FMT,
338 IM * JM * KM * num_components, N_full_field);
339 }
340
341 // Allocate the smaller output C array (AoS layout expected by your writer)
342 ierr = PetscMalloc1(num_components * npoints, out_data); CHKERRQ(ierr);
343
344 /*
345 // ======== Layout diagnostics (rank 0 only, read-only) ========
346 if (num_components > 1 && full_field_arr) {
347 const PetscInt n_full = IM * JM * KM;
348
349 // Choose a safe interior probe point and its +i neighbor
350 const PetscInt ic = (IM > 4) ? IM/2 : 1;
351 const PetscInt jc = (JM > 4) ? JM/2 : 1;
352 const PetscInt kc = (KM > 4) ? KM/2 : 1;
353 const PetscInt ic1 = (ic + 1 < IM-1) ? (ic + 1) : ic; // stay interior
354
355 const PetscInt t_full_c = (kc * JM + jc) * IM + ic;
356 const PetscInt t_full_c1 = (kc * JM + jc) * IM + ic1;
357
358 // If source were AoS: full[num_components * t + comp]
359 PetscScalar aos_c0 = full_field_arr[num_components * t_full_c + 0];
360 PetscScalar aos_c1 = full_field_arr[num_components * t_full_c1 + 0];
361 PetscScalar aos_c0_y = (num_components > 1) ? full_field_arr[num_components * t_full_c + 1] : 0.0;
362 PetscScalar aos_c0_z = (num_components > 2) ? full_field_arr[num_components * t_full_c + 2] : 0.0;
363
364 // If source were SoA: full[comp * n_full + t]
365 PetscScalar soa_c0 = full_field_arr[0 * n_full + t_full_c];
366 PetscScalar soa_c1 = full_field_arr[0 * n_full + t_full_c1];
367 PetscScalar soa_c0_y = (num_components > 1) ? full_field_arr[1 * n_full + t_full_c] : 0.0;
368 PetscScalar soa_c0_z = (num_components > 2) ? full_field_arr[2 * n_full + t_full_c] : 0.0;
369
370 double dAOS = fabs((double)(aos_c1 - aos_c0)); // spatial delta along +i
371 double dSOA = fabs((double)(soa_c1 - soa_c0));
372
373 const char *verdict = (dSOA < dAOS) ? "LIKELY SoA (component-major)" : "LIKELY AoS (interleaved)";
374
375 LOG_ALLOW(GLOBAL, LOG_INFO,
376 "Layout probe @ (i=%" PetscInt_FMT ", j=%" PetscInt_FMT ", k=%" PetscInt_FMT ") vs +i neighbor:\n"
377 " AoS-candidate: Ux(center)=%.6e Ux(+i)=%.6e |Δ|=%.6e ; Uy(center)=%.6e Uz(center)=%.6e\n"
378 " SoA-candidate: Ux(center)=%.6e Ux(+i)=%.6e |Δ|=%.6e ; Uy(center)=%.6e Uz(center)=%.6e\n"
379 " Verdict: %s\n",
380 ic, jc, kc,
381 (double)aos_c0, (double)aos_c1, dAOS, (double)aos_c0_y, (double)aos_c0_z,
382 (double)soa_c0, (double)soa_c1, dSOA, (double)soa_c0_y, (double)soa_c0_z,
383 verdict
384 );
385 } else if (num_components == 1) {
386 LOG_ALLOW(GLOBAL, LOG_INFO, "Layout probe: scalar field (dof=1) — no component interleave to detect.\n");
387 }
388 // ======== END diagnostics ========
389 */
390
391 // --- PACK: assumes source is already AoS ---
392 {
393 PetscInt p_out = 0;
394 for (PetscInt k = 0; k < nz; k++) {
395 for (PetscInt j = 0; j < ny; j++) {
396 for (PetscInt i = 0; i < nx; i++) {
397 const PetscInt p_in_start = num_components * (k * (JM * IM) + j * IM + i);
398 for (PetscInt c = 0; c < num_components; c++) {
399 (*out_data)[p_out++] = full_field_arr[p_in_start + c];
400 }
401 }
402 }
403 }
404 }
405
406 // Free temporary gathered array
407 ierr = PetscFree(full_field_arr); CHKERRQ(ierr);
408 } else {
409 // Other ranks: no allocation
410 *out_data = NULL;
411 }
412
413 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Subsampled field data prepared on rank 0 with total %" PetscInt_FMT " points.\n", npoints);
414
415 PetscFunctionReturn(0);
416}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PrepareOutputParticleData()

PetscErrorCode PrepareOutputParticleData ( UserCtx user,
PostProcessParams pps,
VTKMetaData meta,
PetscInt *  p_n_total 
)

Gathers, subsamples, and prepares all particle data for VTK output.

This function is a COLLECTIVE operation. All ranks must enter it. The heavy lifting (memory allocation, subsampling) is performed only on rank 0.

  1. Gathers the full coordinate and field data from the distributed DMSwarm to rank 0.
  2. On rank 0, subsamples the data based on pps->particle_output_freq.
  3. On rank 0, populates the VTKMetaData struct with the new, smaller, subsampled data arrays.
Parameters
[in]userThe UserCtx containing the DMSwarm.
[in]ppsThe PostProcessParams struct for configuration.
[out]metaA pointer to the VTKMetaData struct to be populated (on rank 0).
[out]p_n_totalOn rank 0, the total number of particles before subsampling.
Returns
PetscErrorCode

Gathers, subsamples, and prepares all particle data for VTK output.

Local to this translation unit.

Definition at line 422 of file vtk_io.c.

423{
424 PetscErrorCode ierr;
425 PetscInt n_total_particles, n_components;
426
427 PetscFunctionBeginUser;
428
429 // Initialize output parameters on all ranks
430 *p_n_total = 0;
431
432 // --- The entire preparation process is a series of collective gathers followed by rank-0 processing ---
433
434 // --- Step 1: Gather Full Coordinates from SOURCE Swarm (Collective) ---
435 // This establishes the ground truth for particle positions and total count.
436 PetscScalar *full_coords_arr = NULL;
437 // NOTE: This assumes SwarmFieldToArrayOnRank0 exists and works for PetscScalar fields.
438 // If not, it may need to be generalized like the logic below.
439 ierr = SwarmFieldToArrayOnRank0(user->swarm, "position", &n_total_particles, &n_components, (void**)&full_coords_arr); CHKERRQ(ierr);
440
441 // --- Step 2: Prepare and Subsample Data (Rank 0 Only) ---
442 if (user->simCtx->rank == 0) {
443 *p_n_total = n_total_particles; // Report original count back to the caller.
444
445 if (n_total_particles == 0) {
446 if (full_coords_arr) {
447 ierr = PetscFree(full_coords_arr); CHKERRQ(ierr); // Cleanup the gathered coords array
448 }
449 PetscFunctionReturn(0); // Nothing to prepare.
450 }
451 if (n_components != 3) {
452 SETERRQ(PETSC_COMM_SELF, 1, "Coordinate field position must have 3 components, but has %" PetscInt_FMT, n_components);
453 }
454
455 // --- Subsampling Calculation ---
456 PetscInt stride = pps->particle_output_freq > 0 ? pps->particle_output_freq : 1;
457 meta->npoints = (n_total_particles > 0) ? (n_total_particles - 1) / stride + 1 : 0;
458
459 LOG_ALLOW(LOCAL, LOG_DEBUG, "Subsampling %" PetscInt_FMT " total particles with stride %" PetscInt_FMT " -> %" PetscInt_FMT " output particles.\n",
460 n_total_particles, stride, meta->npoints);
461
462 // --- Prepare Final Coordinates Array ---
463 ierr = PetscMalloc1(3 * meta->npoints, &meta->coords); CHKERRQ(ierr);
464 for (PetscInt i = 0; i < meta->npoints; i++) {
465 PetscInt source_idx = i * stride;
466 for (int d = 0; d < 3; d++) meta->coords[3 * i + d] = full_coords_arr[3 * source_idx + d];
467 }
468 ierr = PetscFree(full_coords_arr); // Free the temporary full array immediately after use.
469
470
471 // --- Prepare Final Field Arrays (Gather -> Subsample -> Store) ---
472 char *fields_copy, *field_name;
473 ierr = PetscStrallocpy(pps->particle_fields, &fields_copy); CHKERRQ(ierr);
474 field_name = strtok(fields_copy, ",");
475 while (field_name && meta->num_point_data_fields < MAX_POINT_DATA_FIELDS) {
476 TrimWhitespace(field_name);
477 if (!*field_name || strcasecmp(field_name, "position") == 0) {
478 field_name = strtok(NULL, ","); continue;
479 }
480
481 // A. Determine which swarm is the source for this field.
482 DM swarm_to_use = user->swarm; // Default to the main solver swarm.
483 const char* internal_name = field_name; // Default internal name to user-facing name.
484
485 PetscErrorCode check_ierr;
486 ierr = PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL); CHKERRQ(ierr);
487 check_ierr = DMSwarmGetField(user->post_swarm, field_name, NULL, NULL, NULL);
488 ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
489 if (!check_ierr) {
490 swarm_to_use = user->post_swarm;
491 ierr = DMSwarmRestoreField(user->post_swarm, field_name, NULL, NULL, NULL); CHKERRQ(ierr);
492 LOG_ALLOW(LOCAL, LOG_DEBUG, "Field '%s' will be sourced from the post-processing swarm.\n", field_name);
493 } else {
494 LOG_ALLOW(LOCAL, LOG_DEBUG, "Field '%s' will be sourced from the main solver swarm.\n", field_name);
495 if (strcasecmp(field_name, "pid") == 0) internal_name = "DMSwarm_pid";
496 else if (strcasecmp(field_name, "CellID") == 0) internal_name = "DMSwarm_CellID";
497 else if (strcasecmp(field_name, "Migration Status") == 0) internal_name = "DMSwarm_location_status";
498 }
499
500 // B. Gather the full data for this field from its determined source.
501 // NOTE: This assumes a generic gather function exists.
502 // A simplified placeholder is used here; you may need to implement this.
503 void* full_field_arr_void = NULL;
504 PetscInt field_total_particles, field_num_components;
505
506 ierr = SwarmFieldToArrayOnRank0(swarm_to_use, internal_name, &field_total_particles, &field_num_components, &full_field_arr_void); CHKERRQ(ierr);
507
508 if (field_total_particles != n_total_particles) {
509 LOG_ALLOW(LOCAL, LOG_WARNING, "Field '%s' has %" PetscInt_FMT " particles, but expected %" PetscInt_FMT ". Skipping.\n", field_name, field_total_particles, n_total_particles);
510 ierr = PetscFree(full_field_arr_void); CHKERRQ(ierr);
511 field_name = strtok(NULL, ","); continue;
512 }
513
514 // C. Allocate final array (ALWAYS as PetscScalar) and copy/cast subsampled data.
515 VTKFieldInfo* current_field = &meta->point_data_fields[meta->num_point_data_fields];
516 strncpy(current_field->name, field_name, MAX_VTK_FIELD_NAME_LENGTH - 1);
517 current_field->name[MAX_VTK_FIELD_NAME_LENGTH - 1] = '\0';
518 current_field->num_components = field_num_components;
519 ierr = PetscMalloc1(current_field->num_components * meta->npoints, (PetscScalar**)&current_field->data); CHKERRQ(ierr);
520
521 PetscScalar* final_data_arr = (PetscScalar*)current_field->data;
522
523 // D. Perform the subsampling and potential casting
524 if (strcasecmp(internal_name, "DMSwarm_pid") == 0) {
525 PetscInt64* source_arr = (PetscInt64*)full_field_arr_void;
526 for (PetscInt i = 0; i < meta->npoints; i++) {
527 final_data_arr[i] = (PetscScalar)source_arr[i * stride];
528 }
529 } else if (strcasecmp(internal_name, "DMSwarm_CellID") == 0 || strcasecmp(internal_name, "DMSwarm_location_status") == 0) {
530 PetscInt* source_arr = (PetscInt*)full_field_arr_void;
531 for (PetscInt i = 0; i < meta->npoints; i++) {
532 PetscInt source_idx = i * stride;
533 for (PetscInt c = 0; c < current_field->num_components; c++) {
534 final_data_arr[current_field->num_components * i + c] = (PetscScalar)source_arr[current_field->num_components * source_idx + c];
535 }
536 }
537 } else { // Default case: The field is already PetscScalar
538 PetscScalar* source_arr = (PetscScalar*)full_field_arr_void;
539 for (PetscInt i = 0; i < meta->npoints; i++) {
540 PetscInt source_idx = i * stride;
541 for (PetscInt c = 0; c < current_field->num_components; c++) {
542 final_data_arr[current_field->num_components * i + c] = source_arr[current_field->num_components * source_idx + c];
543 }
544 }
545 }
546
547 ierr = PetscFree(full_field_arr_void); // Free the temporary full array.
548 meta->num_point_data_fields++;
549 field_name = strtok(NULL, ",");
550 }
551 ierr = PetscFree(fields_copy); CHKERRQ(ierr);
552
553 // --- Finalize VTK MetaData for PolyData ---
554 if (meta->npoints > 0) {
555 meta->fileType = VTK_POLYDATA;
556 ierr = PetscMalloc1(meta->npoints, &meta->connectivity); CHKERRQ(ierr);
557 ierr = PetscMalloc1(meta->npoints, &meta->offsets); CHKERRQ(ierr);
558 for (PetscInt i = 0; i < meta->npoints; i++) {
559 meta->connectivity[i] = i;
560 meta->offsets[i] = i + 1;
561 }
562 }
563 } // End of rank 0 block
564
565 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
566
567 PetscFunctionReturn(0);
568}
PetscErrorCode SwarmFieldToArrayOnRank0(DM swarm, const char *field_name, PetscInt *n_total_particles, PetscInt *n_components, void **gathered_array)
Gathers any DMSwarm field from all ranks to a single, contiguous array on rank 0.
Definition io.c:1953
void TrimWhitespace(char *str)
Helper function to trim leading/trailing whitespace from a string.
Definition io.c:38
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
#define MAX_POINT_DATA_FIELDS
Defines the maximum number of data fields for VTK point data.
Definition variables.h:547
DM post_swarm
Definition variables.h:886
char name[64]
Definition variables.h:590
PetscInt particle_output_freq
Definition variables.h:571
char particle_fields[1024]
Definition variables.h:569
#define MAX_VTK_FIELD_NAME_LENGTH
Maximum length for VTK field names.
Definition variables.h:548
Here is the call graph for this function:
Here is the caller graph for this function: