vtk_io.c

Go to the documentation of this file.

#include "vtk_io.h"
#include "particle_field_catalog.h"
 
//================================================================================
//                 STATIC (PRIVATE) HELPER FUNCTION PROTOTYPES
//================================================================================
 
static PetscErrorCode WriteVTKAppendedBlock(FILE *fp, const void *data, PetscInt num_elements, size_t element_size);
static PetscErrorCode WriteVTSXMLHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset);
static PetscErrorCode WriteVTPXMLHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset);
static PetscErrorCode WriteVTKFileHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset);
static PetscErrorCode WriteVTKFileFooter(FILE *fp, const VTKMetaData *meta);
 
//================================================================================
//                 IMPLEMENTATION OF PRIVATE HELPERS
//================================================================================
 
/**
 * @brief Write one binary data block in VTK appended-data format.
 */
static PetscErrorCode WriteVTKAppendedBlock(FILE *fp, const void *data, PetscInt num_elements, size_t element_size) {
    uint32_t block_size = num_elements * (uint32_t)element_size;
    if (fwrite(&block_size, sizeof(uint32_t), 1, fp) != 1) return PETSC_ERR_FILE_WRITE;
    if (fwrite(data, element_size, num_elements, fp) != (size_t)num_elements) return PETSC_ERR_FILE_WRITE;
    return 0;
}
 
/**
 * @brief Write the XML header for a VTK structured-grid (`.vts`) file.
 */
static PetscErrorCode WriteVTSXMLHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset)
{
    const char *byte_order = "LittleEndian";
    const char *precision  = "Float64";
 
    const PetscInt header = (PetscInt)sizeof(uint32_t);
 
    fprintf(fp, "<?xml version=\"1.0\"?>\n");
    fprintf(fp, "<VTKFile type=\"StructuredGrid\" version=\"1.0\" byte_order=\"%s\">\n", byte_order);
    fprintf(fp, "  <StructuredGrid WholeExtent=\"0 %" PetscInt_FMT " 0 %" PetscInt_FMT " 0 %" PetscInt_FMT "\">\n", meta->mx - 1, meta->my - 1, meta->mz - 1);
    fprintf(fp, "    <Piece Extent=\"0 %" PetscInt_FMT " 0 %" PetscInt_FMT " 0 %" PetscInt_FMT "\">\n", meta->mx - 1, meta->my - 1, meta->mz - 1);
 
    fprintf(fp, "      <Points>\n");
    fprintf(fp, "        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",
            precision, *boffset);
    *boffset += header + 3 * meta->npoints * (PetscInt)sizeof(PetscScalar);
    fprintf(fp, "      </Points>\n");
 
    if (meta->num_point_data_fields > 0) {
        fprintf(fp, "      <PointData>\n");
        for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
            const VTKFieldInfo* field = &meta->point_data_fields[i];
            fprintf(fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",
                    precision, field->name, field->num_components, *boffset);
            *boffset += header + field->num_components * meta->npoints * (PetscInt)sizeof(PetscScalar);
        }
        fprintf(fp, "      </PointData>\n");
    }
 
    fprintf(fp, "    </Piece>\n");
    fprintf(fp, "  </StructuredGrid>\n");
    fprintf(fp, "  <AppendedData encoding=\"raw\">\n_");
    
    return 0;
}
 
/**
 * @brief Write the XML header for a VTK polydata (`.vtp`) file.
 */
static PetscErrorCode WriteVTPXMLHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset)
{
    const char *byte_order = "LittleEndian";
    const char *precision     = "Float64";
    const char *int_type_str  = (sizeof(PetscInt) == 8) ? "Int64" : "Int32";
    const PetscInt header = (PetscInt)sizeof(uint32_t);
 
    fprintf(fp, "<?xml version=\"1.0\"?>\n");
    fprintf(fp, "<VTKFile type=\"PolyData\" version=\"1.0\" byte_order=\"%s\">\n", byte_order);
    fprintf(fp, "  <PolyData>\n");
    fprintf(fp, "    <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfVerts=\"%" PetscInt_FMT "\" NumberOfLines=\"0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">\n",
            meta->npoints, meta->npoints);
 
    fprintf(fp, "      <Points>\n");
    fprintf(fp, "        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",
            precision, *boffset);
    *boffset += header + 3 * meta->npoints * (PetscInt)sizeof(PetscScalar);
    fprintf(fp, "      </Points>\n");
    
    if (meta->num_point_data_fields > 0) {
        fprintf(fp, "      <PointData>\n");
        for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
            const VTKFieldInfo* field = &meta->point_data_fields[i];
            fprintf(fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",
                    precision, field->name, field->num_components, *boffset);
            *boffset += header + field->num_components * meta->npoints * (PetscInt)sizeof(PetscScalar);
        }
        fprintf(fp, "      </PointData>\n");
    }
 
    fprintf(fp, "      <Verts>\n");
    fprintf(fp, "        <DataArray type=\"%s\" Name=\"connectivity\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",
            int_type_str, *boffset);
    *boffset += (uint32_t)sizeof(uint32_t) + meta->npoints * (PetscInt)sizeof(PetscInt);
    fprintf(fp, "        <DataArray type=\"%s\" Name=\"offsets\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",
            int_type_str, *boffset);
    *boffset += (uint32_t)sizeof(uint32_t) + meta->npoints * (PetscInt)sizeof(PetscInt);
    fprintf(fp, "      </Verts>\n");
 
    fprintf(fp, "    </Piece>\n");
    fprintf(fp, "  </PolyData>\n");
    fprintf(fp, "  <AppendedData encoding=\"raw\">\n_");
 
    return 0;
}
 
/**
 * @brief Open a VTK XML document and write its file-level header.
 */
static PetscErrorCode WriteVTKFileHeader(FILE *fp, const VTKMetaData *meta, PetscInt *boffset)
{
    if (meta->fileType == VTK_STRUCTURED) {
        return WriteVTSXMLHeader(fp, meta, boffset);
    } else if (meta->fileType == VTK_POLYDATA) {
        return WriteVTPXMLHeader(fp, meta, boffset);
    }
    return PETSC_ERR_ARG_WRONG;
}
 
/**
 * @brief Close a VTK XML document after all appended data have been written.
 */
static PetscErrorCode WriteVTKFileFooter(FILE *fp, const VTKMetaData *meta)
{
    (void)meta;
    fprintf(fp, "\n  </AppendedData>\n");
    fprintf(fp, "</VTKFile>\n");
    return 0;
}
 
 
 
/**
 * @brief Implementation of \ref CreateVTKFileFromMetadata().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/vtk_io.h`.
 * @see CreateVTKFileFromMetadata()
 */
 
PetscErrorCode CreateVTKFileFromMetadata(const char *filename, const VTKMetaData *meta, MPI_Comm comm)
{
    PetscMPIInt rank;
    MPI_Comm_rank(comm, &rank);
    PetscErrorCode ierr = 0;
 
    if (!rank) {
        char tmp_filename[PETSC_MAX_PATH_LEN];
        FILE *fp = NULL;
 
        ierr = PetscSNPrintf(tmp_filename, sizeof(tmp_filename), "%s.tmp", filename);CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 writing combined VTK file '%s'.\n", filename);
        fp = fopen(tmp_filename, "wb");
        if (!fp) {
            LOG_ALLOW(GLOBAL, LOG_ERROR, "fopen failed for %s.\n", tmp_filename);
            return PETSC_ERR_FILE_OPEN;
        }
 
        PetscInt boffset = 0;
 
        ierr = WriteVTKFileHeader(fp, meta, &boffset);
        if (ierr) {
            fclose(fp);
            remove(tmp_filename);
            return ierr;
        }
 
        if (meta->coords) {
            ierr = WriteVTKAppendedBlock(fp, meta->coords, 3 * meta->npoints, sizeof(PetscScalar));
            if (ierr) {
                fclose(fp);
                remove(tmp_filename);
                return ierr;
            }
        }
 
        for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
            const VTKFieldInfo *field = &meta->point_data_fields[i];
            if (field->data) {
                ierr = WriteVTKAppendedBlock(fp, field->data, field->num_components * meta->npoints, sizeof(PetscScalar));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
        }
        if (meta->fileType == VTK_POLYDATA) {
            if (meta->connectivity) {
                ierr = WriteVTKAppendedBlock(fp, meta->connectivity, meta->npoints, sizeof(PetscInt));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
            if (meta->offsets) {
                ierr = WriteVTKAppendedBlock(fp, meta->offsets, meta->npoints, sizeof(PetscInt));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
        }
        ierr = WriteVTKFileFooter(fp, meta);
        if (ierr) {
            fclose(fp);
            remove(tmp_filename);
            return ierr;
        }
        if (fclose(fp) != 0) {
            remove(tmp_filename);
            return PETSC_ERR_FILE_WRITE;
        }
        if (rename(tmp_filename, filename) != 0) {
            remove(tmp_filename);
            return PETSC_ERR_FILE_WRITE;
        }
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 finished writing VTK file '%s'.\n", filename);
    }
    return 0;
}
 
//================================================================================
//                 PUBLIC COORDINATE PREPARATION FUNCTION
//================================================================================
 
/**
 * @brief Internal helper implementation: `PrepareOutputCoordinates()`.
 * @details Local to this translation unit.
 */
PetscErrorCode PrepareOutputCoordinates(UserCtx* user, PetscScalar** out_coords, PetscInt* out_nx, PetscInt* out_ny, PetscInt* out_nz, PetscInt* out_npoints)
{
    PetscErrorCode ierr;
    Vec            coords_global;
    PetscInt       IM, JM, KM;
    PetscInt       N_full_coords;
    PetscScalar    *full_coords_arr = NULL;
 
    PetscFunctionBeginUser;
    ierr = DMDAGetInfo(user->da, NULL, &IM, &JM, &KM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
 
    // Set the dimensions of the output grid
    *out_nx = IM - 1;
    *out_ny = JM - 1;
    *out_nz = KM - 1;
    *out_npoints = (*out_nx) * (*out_ny) * (*out_nz);
 
    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);
 
    // --- Step 1: Gather the full coordinate vector to rank 0 ---
    ierr = DMGetCoordinates(user->da, &coords_global); CHKERRQ(ierr);
    // Reuse of your existing, proven utility function's logic
    ierr = VecToArrayOnRank0(coords_global, &N_full_coords, &full_coords_arr); CHKERRQ(ierr);
 
    // --- Step 2: On rank 0, subsample the gathered array ---
    if (user->simCtx->rank == 0) {
        // We expect N_full_coords to be IM * JM * KM * 3
        if (N_full_coords != IM * JM * KM * 3) {
            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);
        }
 
        // Allocate the smaller output C array
        ierr = PetscMalloc1(3 * (*out_npoints), out_coords); CHKERRQ(ierr);
 
        // Loop over the smaller grid dimensions and copy data
        PetscInt p_out = 0; // Index for the small output array
        for (PetscInt k = 0; k < *out_nz; k++) {
            for (PetscInt j = 0; j < *out_ny; j++) {
                for (PetscInt i = 0; i < *out_nx; i++) {
                    // Calculate the index in the full, 1D source array
                    PetscInt p_in = 3 * (k * (JM * IM) + j * IM + i);
                    
                    (*out_coords)[p_out++] = full_coords_arr[p_in + 0]; // x
                    (*out_coords)[p_out++] = full_coords_arr[p_in + 1]; // y
                    (*out_coords)[p_out++] = full_coords_arr[p_in + 2]; // z
                }
            }
        }
        // Free the temporary full array that was allocated by VecToArrayOnRank0
        ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
    } else {
        // On other ranks, the output pointer is NULL
        *out_coords = NULL;
    }
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Subsampled coordinates prepared on rank 0 with total %" PetscInt_FMT " points.\n", *out_npoints);
 
    PetscFunctionReturn(0);
}
 
 
/**
 * @brief Implementation of \ref BeginStructuredVTKOutput().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/vtk_io.h`.
 * @see BeginStructuredVTKOutput()
 */
PetscErrorCode BeginStructuredVTKOutput(UserCtx *user, VTKMetaData *meta)
{
    PetscFunctionBeginUser;
    PetscCheck(user != NULL && meta != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context and metadata are required.");
    PetscCall(PetscMemzero(meta, sizeof(*meta)));
    meta->fileType = VTK_STRUCTURED;
    PetscCall(PrepareOutputCoordinates(user, &meta->coords, &meta->mx, &meta->my, &meta->mz,
                                       &meta->npoints));
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref AppendStructuredVTKField().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/vtk_io.h`.
 * @see AppendStructuredVTKField()
 */
PetscErrorCode AppendStructuredVTKField(UserCtx *user, VTKMetaData *meta, const char *name,
                                        Vec field_vec, PetscInt components)
{
    VTKFieldInfo *entry = NULL;
 
    PetscFunctionBeginUser;
    PetscCheck(user != NULL && meta != NULL && name != NULL && field_vec != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Context, metadata, name, and vector are required.");
    if (meta->num_point_data_fields >= MAX_POINT_DATA_FIELDS) {
        LOG_ALLOW(GLOBAL, LOG_WARNING,
                  "Reached the %d point-data field limit; '%s' is not written.\n",
                  MAX_POINT_DATA_FIELDS, name);
        PetscFunctionReturn(0);
    }
    entry = &meta->point_data_fields[meta->num_point_data_fields];
    PetscCall(PetscStrncpy(entry->name, name, MAX_VTK_FIELD_NAME_LENGTH));
    entry->num_components = components;
    PetscCall(PrepareOutputEulerianFieldData(user, field_vec, components, &entry->data));
    meta->num_point_data_fields++;
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref FinishStructuredVTKOutput().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/vtk_io.h`.
 * @see FinishStructuredVTKOutput()
 */
PetscErrorCode FinishStructuredVTKOutput(VTKMetaData *meta, const char *filename)
{
    PetscFunctionBeginUser;
    PetscCheck(meta != NULL && filename != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Metadata and filename are required.");
    PetscCall(CreateVTKFileFromMetadata(filename, meta, PETSC_COMM_WORLD));
    for (PetscInt index = 0; index < meta->num_point_data_fields; ++index) {
        if (meta->point_data_fields[index].data) {
            PetscCall(PetscFree(meta->point_data_fields[index].data));
        }
    }
    if (meta->coords) PetscCall(PetscFree(meta->coords));
    meta->num_point_data_fields = 0;
    PetscFunctionReturn(0);
}
 
/**
 * @brief Implementation of \ref PrepareOutputEulerianFieldData().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/vtk_io.h`.
 * @see PrepareOutputEulerianFieldData()
 */
 
PetscErrorCode PrepareOutputEulerianFieldData(UserCtx *user, Vec field_vec, PetscInt num_components, PetscScalar** out_data)
{
    PetscErrorCode ierr;
    PetscInt       IM, JM, KM, nx, ny, nz, npoints;
    PetscInt       N_full_field = 0;
    PetscScalar   *full_field_arr = NULL;
 
    PetscFunctionBeginUser;
    ierr = DMDAGetInfo(user->da, NULL, &IM, &JM, &KM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    nx = IM - 1; ny = JM - 1; nz = KM - 1;
    npoints = nx * ny * nz;
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "Preparing subsampled field data (dof=%" PetscInt_FMT ") for a %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT " grid.\n",
              num_components, nx, ny, nz);
 
    // --- Step 1: Gather the full field vector to rank 0 ---
    ierr = VecToArrayOnRank0(field_vec, &N_full_field, &full_field_arr); CHKERRQ(ierr);
 
    // --- Step 2: On rank 0, subsample the gathered array ---
    if (user->simCtx->rank == 0) {
        // Sanity check the size of the gathered array
        if (N_full_field != IM * JM * KM * num_components) {
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
                    "Gathered field array has wrong size. Expected %" PetscInt_FMT ", got %" PetscInt_FMT,
                    IM * JM * KM * num_components, N_full_field);
        }
 
        // Allocate the smaller output C array (AoS layout expected by your writer)
        ierr = PetscMalloc1(num_components * npoints, out_data); CHKERRQ(ierr);
 
        /*
        // ======== Layout diagnostics (rank 0 only, read-only) ========
        if (num_components > 1 && full_field_arr) {
            const PetscInt n_full = IM * JM * KM;
 
            // Choose a safe interior probe point and its +i neighbor
            const PetscInt ic = (IM > 4) ? IM/2 : 1;
            const PetscInt jc = (JM > 4) ? JM/2 : 1;
            const PetscInt kc = (KM > 4) ? KM/2 : 1;
            const PetscInt ic1 = (ic + 1 < IM-1) ? (ic + 1) : ic; // stay interior
 
            const PetscInt t_full_c  = (kc * JM + jc) * IM + ic;
            const PetscInt t_full_c1 = (kc * JM + jc) * IM + ic1;
 
            // If source were AoS:  full[num_components * t + comp]
            PetscScalar aos_c0   = full_field_arr[num_components * t_full_c  + 0];
            PetscScalar aos_c1   = full_field_arr[num_components * t_full_c1 + 0];
            PetscScalar aos_c0_y = (num_components > 1) ? full_field_arr[num_components * t_full_c + 1] : 0.0;
            PetscScalar aos_c0_z = (num_components > 2) ? full_field_arr[num_components * t_full_c + 2] : 0.0;
 
            // If source were SoA:  full[comp * n_full + t]
            PetscScalar soa_c0   = full_field_arr[0 * n_full + t_full_c];
            PetscScalar soa_c1   = full_field_arr[0 * n_full + t_full_c1];
            PetscScalar soa_c0_y = (num_components > 1) ? full_field_arr[1 * n_full + t_full_c] : 0.0;
            PetscScalar soa_c0_z = (num_components > 2) ? full_field_arr[2 * n_full + t_full_c] : 0.0;
 
            double dAOS = fabs((double)(aos_c1 - aos_c0)); // spatial delta along +i
            double dSOA = fabs((double)(soa_c1 - soa_c0));
 
            const char *verdict = (dSOA < dAOS) ? "LIKELY SoA (component-major)" : "LIKELY AoS (interleaved)";
 
            LOG_ALLOW(GLOBAL, LOG_INFO,
                "Layout probe @ (i=%" PetscInt_FMT ", j=%" PetscInt_FMT ", k=%" PetscInt_FMT ") vs +i neighbor:\n"
                "  AoS-candidate:  Ux(center)=%.6e  Ux(+i)=%.6e  |Δ|=%.6e ; Uy(center)=%.6e  Uz(center)=%.6e\n"
                "  SoA-candidate:  Ux(center)=%.6e  Ux(+i)=%.6e  |Δ|=%.6e ; Uy(center)=%.6e  Uz(center)=%.6e\n"
                "  Verdict: %s\n",
                ic, jc, kc,
                (double)aos_c0,  (double)aos_c1,  dAOS, (double)aos_c0_y, (double)aos_c0_z,
                (double)soa_c0,  (double)soa_c1,  dSOA, (double)soa_c0_y, (double)soa_c0_z,
                verdict
            );
        } else if (num_components == 1) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "Layout probe: scalar field (dof=1) — no component interleave to detect.\n");
        }
        // ======== END diagnostics ========
        */
 
        // --- PACK: assumes source is already AoS ---
        {
            PetscInt p_out = 0;
            for (PetscInt k = 0; k < nz; k++) {
                for (PetscInt j = 0; j < ny; j++) {
                    for (PetscInt i = 0; i < nx; i++) {
                        const PetscInt p_in_start = num_components * (k * (JM * IM) + j * IM + i);
                        for (PetscInt c = 0; c < num_components; c++) {
                            (*out_data)[p_out++] = full_field_arr[p_in_start + c];
                        }
                    }
                }
            }
        }
 
        // Free temporary gathered array
        ierr = PetscFree(full_field_arr); CHKERRQ(ierr);
    } else {
        // Other ranks: no allocation
        *out_data = NULL;
    }
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Subsampled field data prepared on rank 0 with total %" PetscInt_FMT " points.\n", npoints);
 
    PetscFunctionReturn(0);
}
 
/**
 * @brief Internal helper implementation: `PrepareOutputParticleData()`.
 * @details Local to this translation unit.
 */
PetscErrorCode PrepareOutputParticleData(UserCtx* user, PostProcessParams* pps, VTKMetaData* meta, PetscInt* p_n_total)
{
    PetscErrorCode ierr;
    PetscInt n_total_particles, n_components;
    PetscInt stride = pps->particle_output_freq > 0 ? pps->particle_output_freq : 1;
 
    PetscFunctionBeginUser;
 
    // Initialize output parameters on all ranks
    *p_n_total = 0;
 
    // --- The entire preparation process is a series of collective gathers followed by rank-0 processing ---
 
    // --- Step 1: Gather Full Coordinates from SOURCE Swarm (Collective) ---
    // This establishes the ground truth for particle positions and total count.
    PetscReal *full_coords_arr = NULL;
    PetscDataType coordinate_type;
    // NOTE: This assumes SwarmFieldToArrayOnRank0 exists and works for PetscScalar fields.
    // If not, it may need to be generalized like the logic below.
    ierr = SwarmFieldToArrayOnRank0(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),
                                    &n_total_particles, &n_components, &coordinate_type,
                                    (void**)&full_coords_arr); CHKERRQ(ierr);
    PetscCheck(coordinate_type == PETSC_REAL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Particle position field must use PETSC_REAL storage, not %s.",
               PetscDataTypes[coordinate_type]);
 
    *p_n_total = n_total_particles;
    if (n_total_particles == 0) {
        ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
        PetscFunctionReturn(0);
    }
    PetscCheck(n_components == 3,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Coordinate field position must have 3 components, but has %" PetscInt_FMT,
               n_components);
 
    // --- Step 2: Prepare and Subsample Coordinates (Rank 0 Only) ---
    if (user->simCtx->rank == 0) {
        // --- Subsampling Calculation ---
        meta->npoints = (n_total_particles > 0) ? (n_total_particles - 1) / stride + 1 : 0;
 
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Subsampling %" PetscInt_FMT " total particles with stride %" PetscInt_FMT " -> %" PetscInt_FMT " output particles.\n",
                  n_total_particles, stride, meta->npoints);
 
        // --- Prepare Final Coordinates Array ---
        ierr = PetscMalloc1(3 * meta->npoints, &meta->coords); CHKERRQ(ierr);
        for (PetscInt i = 0; i < meta->npoints; i++) {
            PetscInt source_idx = i * stride;
            for (int d = 0; d < 3; d++) meta->coords[3 * i + d] = full_coords_arr[3 * source_idx + d];
        }
    }
    ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
 
    // --- Step 3: Gather Fields (Collective), then Subsample and Store (Rank 0 Only) ---
    char *fields_copy, *field_name;
    PetscInt num_fields = 0;
    ierr = PetscStrallocpy(pps->particle_fields, &fields_copy); CHKERRQ(ierr);
    field_name = strtok(fields_copy, ",");
    while (field_name && num_fields < MAX_POINT_DATA_FIELDS) {
        TrimWhitespace(field_name);
        if (!*field_name || strcasecmp(field_name, "position") == 0) {
            field_name = strtok(NULL, ","); continue;
        }
 
        // A. Determine which swarm is the source for this field on every rank.
        DM swarm_to_use = user->swarm; // Default to the main solver swarm.
        const char* internal_name = field_name; // Default internal name to user-facing name.
 
        PetscErrorCode check_ierr;
        ierr = PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL); CHKERRQ(ierr);
        check_ierr = DMSwarmGetField(user->post_swarm, field_name, NULL, NULL, NULL);
        ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
        if (!check_ierr) {
            swarm_to_use = user->post_swarm;
            ierr = DMSwarmRestoreField(user->post_swarm, field_name, NULL, NULL, NULL); CHKERRQ(ierr);
            if (user->simCtx->rank == 0) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Field '%s' will be sourced from the post-processing swarm.\n", field_name);
            }
        } else {
            if (user->simCtx->rank == 0) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Field '%s' will be sourced from the main solver swarm.\n", field_name);
            }
            if (strcasecmp(field_name, "pid") == 0) internal_name = "DMSwarm_pid";
            else if (strcasecmp(field_name, "CellID") == 0) internal_name = "DMSwarm_CellID";
            else if (strcasecmp(field_name, "Migration Status") == 0) internal_name = "DMSwarm_location_status";
        }
 
        // B. Gather the complete field on rank 0 before applying the output stride.
        // SwarmFieldToArrayOnRank0 handles both scalar and vector field layouts.
        void* full_field_arr_void = NULL;
        PetscInt field_total_particles, field_num_components;
        PetscDataType field_data_type;
 
        ierr = SwarmFieldToArrayOnRank0(swarm_to_use, internal_name,
                                        &field_total_particles, &field_num_components,
                                        &field_data_type, &full_field_arr_void); CHKERRQ(ierr);
 
        if (field_total_particles != n_total_particles) {
            if (user->simCtx->rank == 0) {
                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);
            }
            ierr = PetscFree(full_field_arr_void); CHKERRQ(ierr);
            field_name = strtok(NULL, ","); continue;
        }
 
        if (user->simCtx->rank == 0) {
            // C. Allocate final array (ALWAYS as PetscScalar) and copy/cast subsampled data.
            VTKFieldInfo* current_field = &meta->point_data_fields[meta->num_point_data_fields];
            strncpy(current_field->name, field_name, MAX_VTK_FIELD_NAME_LENGTH - 1);
            current_field->name[MAX_VTK_FIELD_NAME_LENGTH - 1] = '\0';
            current_field->num_components = field_num_components;
            ierr = PetscMalloc1(current_field->num_components * meta->npoints, (PetscScalar**)&current_field->data); CHKERRQ(ierr);
 
            PetscScalar* final_data_arr = (PetscScalar*)current_field->data;
 
            // D. Perform the subsampling and potential casting
            if (field_data_type == PETSC_INT64) {
                PetscInt64* source_arr = (PetscInt64*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    final_data_arr[i] = (PetscScalar)source_arr[i * stride];
                }
            } else if (field_data_type == PETSC_INT) {
                PetscInt* source_arr = (PetscInt*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    PetscInt source_idx = i * stride;
                    for (PetscInt c = 0; c < current_field->num_components; c++) {
                        final_data_arr[current_field->num_components * i + c] = (PetscScalar)source_arr[current_field->num_components * source_idx + c];
                    }
                }
            } else if (field_data_type == PETSC_REAL) {
                PetscReal* source_arr = (PetscReal*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    PetscInt source_idx = i * stride;
                    for (PetscInt c = 0; c < current_field->num_components; c++) {
                        final_data_arr[current_field->num_components * i + c] = source_arr[current_field->num_components * source_idx + c];
                    }
                }
#if defined(PETSC_USE_COMPLEX)
            } else if (field_data_type == PETSC_SCALAR) {
                PetscScalar* source_arr = (PetscScalar*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    PetscInt source_idx = i * stride;
                    for (PetscInt c = 0; c < current_field->num_components; c++) {
                        final_data_arr[current_field->num_components * i + c] = source_arr[current_field->num_components * source_idx + c];
                    }
                }
#endif
            } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
                           "VTK particle output does not support field '%s' with PETSc type %s.",
                           field_name, PetscDataTypes[field_data_type]);
 
            meta->num_point_data_fields++;
        }
        ierr = PetscFree(full_field_arr_void); CHKERRQ(ierr);
        num_fields++;
        field_name = strtok(NULL, ",");
    }
    ierr = PetscFree(fields_copy); CHKERRQ(ierr);
 
    // --- Step 4: Finalize VTK MetaData for PolyData (Rank 0 Only) ---
    if (user->simCtx->rank == 0) {
        if (meta->npoints > 0) {
            meta->fileType = VTK_POLYDATA;
            ierr = PetscMalloc1(meta->npoints, &meta->connectivity); CHKERRQ(ierr);
            ierr = PetscMalloc1(meta->npoints, &meta->offsets); CHKERRQ(ierr);
            for (PetscInt i = 0; i < meta->npoints; i++) {
                meta->connectivity[i] = i;
                meta->offsets[i] = i + 1;
            }
        }
    }
 
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}