postprocessor.c File Reference

Phase 2 implementation of the post-processing tool. More...

#include "postprocessor.h"

Include dependency graph for postprocessor.c:

Go to the source code of this file.

Macros
#define	__FUNCT__   "SetupPostProcessSwarm"
#define	__FUNCT__   "EulerianDataProcessingPipeline"
#define	__FUNCT__   "WriteEulerianFile"
#define	__FUNCT__   "ParticleDataProcessingPipeline"
#define	__FUNCT__   "WriteParticleFile"
#define	__FUNCT__   "main"

Functions
PetscErrorCode	SetupPostProcessSwarm (UserCtx *user, PostProcessParams *pps)
	Creates a new, dedicated DMSwarm for post-processing tasks.
PetscErrorCode	EulerianDataProcessingPipeline (UserCtx *user, PostProcessParams *pps)
	Parses the processing pipeline string from the config and executes the requested kernels in sequence.
PetscErrorCode	WriteEulerianFile (UserCtx *user, PostProcessParams *pps, PetscInt ti)
	Writes Eulerian field data to a VTK structured grid file (.vts) for visualization.
PetscErrorCode	ParticleDataProcessingPipeline (UserCtx *user, PostProcessParams *pps)
	Parses and executes the particle pipeline using a robust two-pass approach.
PetscErrorCode	WriteParticleFile (UserCtx *user, PostProcessParams *pps, PetscInt ti)
	Writes particle data to a VTP file using the Prepare-Write-Cleanup pattern.
int	main (int argc, char **argv)

Detailed Description

Phase 2 implementation of the post-processing tool.

This phase introduces a dedicated configuration system and performs a single-step data load to verify the I/O and data structures.

Definition in file postprocessor.c.

Macro Definition Documentation

__FUNCT__ [1/6]

#define __FUNCT__   "SetupPostProcessSwarm"

Definition at line 13 of file postprocessor.c.

__FUNCT__ [2/6]

#define __FUNCT__   "EulerianDataProcessingPipeline"

Definition at line 13 of file postprocessor.c.

__FUNCT__ [3/6]

#define __FUNCT__   "WriteEulerianFile"

Definition at line 13 of file postprocessor.c.

__FUNCT__ [4/6]

#define __FUNCT__   "ParticleDataProcessingPipeline"

Definition at line 13 of file postprocessor.c.

__FUNCT__ [5/6]

#define __FUNCT__   "WriteParticleFile"

Definition at line 13 of file postprocessor.c.

__FUNCT__ [6/6]

#define __FUNCT__   "main"

Definition at line 13 of file postprocessor.c.

Function Documentation

SetupPostProcessSwarm()

PetscErrorCode SetupPostProcessSwarm	(	UserCtx *	user,
		PostProcessParams *	pps
	)

Creates a new, dedicated DMSwarm for post-processing tasks.

This function is called once at startup. It creates an empty DMSwarm and associates it with the same grid DM as the primary swarm and registers all the required fields.

Parameters

user	The UserCtx where user->post_swarm will be created.
pps	The PostProcessParams containing the particle_pipeline string for field registration.

Returns

PetscErrorCode

Definition at line 23 of file postprocessor.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    char *pipeline_copy, *step_token, *step_saveptr;
    PetscBool finalize_needed = PETSC_FALSE;
 
    ierr = DMCreate(PETSC_COMM_WORLD, &user->post_swarm); CHKERRQ(ierr);
    ierr = DMSetType(user->post_swarm, DMSWARM); CHKERRQ(ierr);
    ierr = DMSetDimension(user->post_swarm, 3); CHKERRQ(ierr);
    ierr = DMSwarmSetType(user->post_swarm, DMSWARM_BASIC); CHKERRQ(ierr);
    // Associate it with the same grid as the solver's swarm
     if (user->da) {
      ierr = DMSwarmSetCellDM(user->post_swarm, user->da); CHKERRQ(ierr);
      LOG_ALLOW(LOCAL,LOG_INFO,"Associated DMSwarm with Cell DM (user->da).\n");
    } else {
      // If user->da is essential for your simulation logic with particles, this should be a fatal error.
      LOG_ALLOW(GLOBAL, LOG_WARNING, "user->da (Cell DM for Swarm) is NULL. Cell-based swarm operations might fail.\n");
      // SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "CreateParticleSwarm - user->da (Cell DM) is NULL but required.");
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Created dedicated DMSwarm for post-processing.\n");
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, " --- Setting up Post-Processing Pipeline fields-- \n");
 
    ierr = PetscStrallocpy(pps->particle_pipeline, &pipeline_copy); CHKERRQ(ierr);
    step_token = strtok_r(pipeline_copy, ";", &step_saveptr);
    while (step_token) {
        TrimWhitespace(step_token);
        if (strlen(step_token) == 0) { step_token = strtok_r(NULL, ";", &step_saveptr); continue; }
 
        char *keyword = strtok(step_token, ":");
        char *args_str = strtok(NULL, "");
        TrimWhitespace(keyword);
        PetscInt output_field_dimensions = 1; // Default to scalar output fields
 
        if (strcasecmp(keyword, "ComputeSpecificKE") == 0) {
            if (!args_str) SETERRQ(PETSC_COMM_SELF, 1, "Error (ComputeSpecificKE): Missing arguments.");
            char *inputs_str = strtok(args_str, ">");
            char *output_field = strtok(NULL, ">");
            output_field_dimensions = 1;  // SKE is scalar
            if (!output_field) SETERRQ(PETSC_COMM_SELF, 1, "Error (ComputeSpecificKE): Missing output field in 'in>out' syntax.");
            TrimWhitespace(output_field);
            // Register the output field
            ierr = RegisterSwarmField(user->post_swarm, output_field, output_field_dimensions,PETSC_REAL); CHKERRQ(ierr);
            LOG_ALLOW(GLOBAL, LOG_INFO, "Registered particle field '%s' (ComputeSpecificKE).\n", output_field);
            finalize_needed = PETSC_TRUE;
        } else {
            LOG_ALLOW(GLOBAL, LOG_WARNING, "Warning: Unknown particle transformation keyword '%s'. Skipping.\n", keyword);
        }
        
        // Add other 'else if' blocks here for other kernels that create output fields
 
        step_token = strtok_r(NULL, ";", &step_saveptr);
    } // while step_token
    
    ierr = PetscFree(pipeline_copy); CHKERRQ(ierr);
 
    // --- FINALIZE STEP ---
    if (finalize_needed) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Finalizing new field registrations for particle pipeline.\n");
        ierr = DMSwarmFinalizeFieldRegister(user->post_swarm); CHKERRQ(ierr);
    }
        
    LOG_ALLOW(GLOBAL, LOG_INFO, "Post-Processing DMSwarm setup complete.\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

EulerianDataProcessingPipeline()

PetscErrorCode EulerianDataProcessingPipeline	(	UserCtx *	user,
		PostProcessParams *	pps
	)

Parses the processing pipeline string from the config and executes the requested kernels in sequence.

Parses the processing pipeline string and executes the requested kernels.

This function uses a general-purpose parser to handle a syntax of the form: "Keyword1:in1>out1; Keyword2:in1,in2>out2; Keyword3:arg1;"

It tokenizes the pipeline string and dispatches to the appropriate kernel function from processing_kernels.c with the specified field name arguments.

Parameters

user	The UserCtx containing the data to be transformed.
pps	The PostProcessParams struct containing the pipeline string.

Returns

PetscErrorCode

Definition at line 110 of file postprocessor.c.

{
    PetscErrorCode ierr;
    char *pipeline_copy, *step_token, *step_saveptr;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Data Transformation Pipeline ---\n");
    
    // Do nothing if the pipeline string is empty
    if (pps->process_pipeline[0] == '\0') {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Processing pipeline is empty. No transformations will be run.\n");
        PetscFunctionReturn(0);
    }
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pipeline string: [%s]\n", pps->process_pipeline);
    
    // Make a writable copy for strtok_r, as it modifies the string
    ierr = PetscStrallocpy(pps->process_pipeline, &pipeline_copy); CHKERRQ(ierr);
 
    // --- Outer Loop: Tokenize by Semicolon (;) to get each processing step ---
    step_token = strtok_r(pipeline_copy, ";", &step_saveptr);
    while (step_token) {
        TrimWhitespace(step_token);
        if (strlen(step_token) == 0) {
            step_token = strtok_r(NULL, ";", &step_saveptr);
            continue;
        }
 
        char *keyword = strtok(step_token, ":");
        char *args_str = strtok(NULL, ""); // Get the rest of the string as arguments
 
        if (!keyword) { // Should not happen with TrimWhitespace, but is a safe check
            step_token = strtok_r(NULL, ";", &step_saveptr);
            continue;
        }
        
        TrimWhitespace(keyword);
        if (args_str) TrimWhitespace(args_str);
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "Executing Transformation: '%s' on args: '%s'\n", keyword, args_str ? args_str : "None");
 
        // --- DISPATCHER: Route to the correct kernel based on the keyword ---
        if (strcasecmp(keyword, "CellToNodeAverage") == 0) {
            if (!args_str) SETERRQ(PETSC_COMM_SELF, 1, "CellToNodeAverage requires arguments in 'in_field>out_field' format.");
            char *in_field = strtok(args_str, ">");
            char *out_field = strtok(NULL, ">");
            if (!in_field || !out_field) SETERRQ(PETSC_COMM_SELF, 1, "CellToNodeAverage requires 'in>out' syntax (e.g., P>P_nodal).");
            if(strcmp(in_field,out_field)==0) SETERRQ(PETSC_COMM_SELF, 1, "CellToNodeAverage input and output fields must be different.");
            if(user->simCtx->np == 0 && (strcmp(out_field,"Psi_nodal")==0 || strcmp(in_field,"Psi_nodal")==0)){
                LOG(GLOBAL,LOG_WARNING,"CellToNodeAverage cannot process 'Psi_nodal' when no particles are present in the simulation.\n");
                step_token = strtok_r(NULL, ";", &step_saveptr);
                continue;
            }
            TrimWhitespace(in_field); TrimWhitespace(out_field);
            ierr = ComputeNodalAverage(user, in_field, out_field); CHKERRQ(ierr);
        }
        else if (strcasecmp(keyword, "ComputeQCriterion") == 0) {
            ierr = ComputeQCriterion(user); CHKERRQ(ierr);
        }
        else if (strcasecmp(keyword, "NormalizeRelativeField") == 0) {
            if (!args_str) SETERRQ(PETSC_COMM_SELF, 1, "NormalizePressure requires the pressure field name (e.g., 'P') as an argument.");
            ierr = NormalizeRelativeField(user, args_str); CHKERRQ(ierr);
        }
        // *** Add new kernels here in the future using 'else if' ***
        // else if (strcasecmp(keyword, "ComputeVorticity") == 0) { ... }
        else {
            LOG_ALLOW(GLOBAL, LOG_WARNING, "Unknown transformation keyword '%s'. Skipping.\n", keyword);
        }
 
        step_token = strtok_r(NULL, ";", &step_saveptr);
    }
 
    ierr = PetscFree(pipeline_copy); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Data Transformation Pipeline Complete ---\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteEulerianFile()

PetscErrorCode WriteEulerianFile	(	UserCtx *	user,
		PostProcessParams *	pps,
		PetscInt	ti
	)

Writes Eulerian field data to a VTK structured grid file (.vts) for visualization.

Orchestrates the writing of a combined, multi-field VTK file for a single time step.

This function exports instantaneous Eulerian field data (such as pressure, velocity, Q-criterion, and particle concentration) to a VTK structured grid file for the specified time step. The output includes subsampled interior grid coordinates and point data fields as specified in the PostProcessParams configuration.

The function performs the following steps:

1. Initializes VTK metadata structure
2. Prepares output coordinates (subsampled interior grid)
3. Gathers and prepares requested field data on rank 0
4. Writes the VTK structured grid file with the naming convention: {prefix}_{ti:05d}.vts

Only fields listed in pps->output_fields_instantaneous are written. If the field list is empty, the function returns immediately without creating a file.

Note

Fields containing particle data (e.g., Psi_nodal) are automatically skipped when no particles are present in the simulation (simCtx->np == 0).

Parameters

user	The UserCtx containing simulation data and field vectors to be output.
pps	The PostProcessParams struct containing output configuration (field list, prefix).
ti	The time index/step number used for file naming.

Returns

PetscErrorCode indicating success or failure.

Definition at line 217 of file postprocessor.c.

{
    PetscErrorCode ierr;
    VTKMetaData    meta;
    char           filename[MAX_FILENAME_LENGTH];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (pps->output_fields_instantaneous[0] == '\0') {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "No instantaneous fields requested for output at ti=%" PetscInt_FMT ". Skipping.\n", ti);
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting VTK File Writing for ti = %" PetscInt_FMT " ---\n", ti);
 
    /* 1) Metadata init */
    ierr = PetscMemzero(&meta, sizeof(VTKMetaData)); CHKERRQ(ierr);
    meta.fileType = VTK_STRUCTURED;
 
    /* 2) Coordinates (subsampled interior) */
    ierr = PrepareOutputCoordinates(user, &meta.coords, &meta.mx, &meta.my, &meta.mz, &meta.npoints); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Using coords linearization order: fast=i mid=j slow=k  (sizes: %" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT ")\n",
              meta.mx, meta.my, meta.mz);
 
    /* 3) Fields (rank 0) */
    if (user->simCtx->rank == 0) {
        char *fields_copy, *field_name;
        ierr = PetscStrallocpy(pps->output_fields_instantaneous, &fields_copy); CHKERRQ(ierr);
 
        field_name = strtok(fields_copy, ",");
        while (field_name) {
            TrimWhitespace(field_name);
            if (!*field_name) { field_name = strtok(NULL, ","); continue; }
 
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Preparing field '%s' for output.\n", field_name);
 
            Vec       field_vec = NULL;
            PetscInt  num_components = 0;
 
            if (!strcasecmp(field_name, "P_nodal")) {
                field_vec = user->P_nodal; num_components = 1;
            } else if (!strcasecmp(field_name, "Ucat_nodal")) {
                field_vec = user->Ucat_nodal; num_components = 3;
            } else if (!strcasecmp(field_name, "Qcrit")) {
                field_vec = user->Qcrit;    num_components = 1;
            } else if (!strcasecmp(field_name, "Psi_nodal")){
                if(user->simCtx->np==0){
                    LOG_ALLOW(LOCAL, LOG_WARNING, "Field 'Psi_nodal' requested but no particles are present. Skipping.\n");
                    field_name = strtok(NULL, ",");
                    continue;
                }
                field_vec = user->Psi_nodal; num_components = 1;
            } else {
                LOG_ALLOW(LOCAL, LOG_WARNING, "Field '%s' not recognized. Skipping.\n", field_name);
                field_name = strtok(NULL, ",");
                continue;
            }
 
            if (meta.num_point_data_fields >= MAX_POINT_DATA_FIELDS) {
                LOG_ALLOW(LOCAL, LOG_WARNING, "MAX_POINT_DATA_FIELDS reached. Cannot add '%s'.\n", field_name);
                field_name = strtok(NULL, ",");
                continue;
            }
 
            // --- Add field to metadata ---
            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 = num_components;
 
            /* Build interior AoS from NATURAL gathered Vec */
            ierr = PrepareOutputEulerianFieldData(user, field_vec, num_components, &current_field->data); CHKERRQ(ierr);
            
            /*
            // *** DEBUG: Dump Ucat_nodal details and add scalar companions Ux, Uy, Uz for easier visualization ***
            
            // If this is Ucat_nodal, dump a few tuples and add scalar companions 
            if (!strcasecmp(field_name, "Ucat_nodal")) {
                const PetscInt npts = meta.npoints;
                const PetscScalar *a = (const PetscScalar*)current_field->data;
 
                LOG_ALLOW(GLOBAL, LOG_INFO, "DBG Ucat_nodal: ptr=%p  npoints=%" PetscInt_FMT "  num_components=%" PetscInt_FMT "\n",
                          (void*)a, npts, current_field->num_components);
 
                if (a && current_field->num_components == 3 && npts > 0) {
                    const PetscInt nshow = (npts < 5) ? npts : 5;
                    LOG_ALLOW(GLOBAL, LOG_INFO, "DBG Ucat_nodal: showing first %d of %" PetscInt_FMT " tuples (AoS x,y,z):\n",
                              (int)nshow, npts);
                    for (PetscInt t = 0; t < nshow; ++t) {
                        LOG_ALLOW(GLOBAL, LOG_INFO, "  Ucat_nodal[%3" PetscInt_FMT "] = (%g, %g, %g)\n",
                                  t, (double)a[3*t+0], (double)a[3*t+1], (double)a[3*t+2]);
                    }
                    if (npts > 10) {
                        PetscInt mid = npts / 2;
                        LOG_ALLOW(GLOBAL, LOG_INFO, "  Ucat_nodal[mid=%" PetscInt_FMT "] = (%g, %g, %g)\n",
                                  mid, (double)a[3*mid+0], (double)a[3*mid+1], (double)a[3*mid+2]);
                    }
 
                    // Add scalar companions from the AoS we just created 
                    PetscScalar *Ux=NULL,*Uy=NULL,*Uz=NULL;
                    ierr = PetscMalloc1(meta.npoints, &Ux); CHKERRQ(ierr);
                    ierr = PetscMalloc1(meta.npoints, &Uy); CHKERRQ(ierr);
                    ierr = PetscMalloc1(meta.npoints, &Uz); CHKERRQ(ierr);
                    for (PetscInt i = 0; i < meta.npoints; ++i) {
                        Ux[i] = a[3*i+0]; Uy[i] = a[3*i+1]; Uz[i] = a[3*i+2];
                    }
 
                    if (meta.num_point_data_fields + 3 <= MAX_POINT_DATA_FIELDS) {
                        VTKFieldInfo *fx = &meta.point_data_fields[++meta.num_point_data_fields];
                        strncpy(fx->name, "Ux_debug", MAX_VTK_FIELD_NAME_LENGTH-1);
                        fx->name[MAX_VTK_FIELD_NAME_LENGTH-1] = '\0';
                        fx->num_components = 1; fx->data = Ux;
 
                        VTKFieldInfo *fy = &meta.point_data_fields[++meta.num_point_data_fields];
                        strncpy(fy->name, "Uy_debug", MAX_VTK_FIELD_NAME_LENGTH-1);
                        fy->name[MAX_VTK_FIELD_NAME_LENGTH-1] = '\0';
                        fy->num_components = 1; fy->data = Uy;
 
                        VTKFieldInfo *fz = &meta.point_data_fields[++meta.num_point_data_fields];
                        strncpy(fz->name, "Uz_debug", MAX_VTK_FIELD_NAME_LENGTH-1);
                        fz->name[MAX_VTK_FIELD_NAME_LENGTH-1] = '\0';
                        fz->num_components = 1; fz->data = Uz;
 
                        LOG_ALLOW(GLOBAL, LOG_INFO, "DBG: Added scalar companions Ux_debug, Uy_debug, Uz_debug.\n");
                    } else {
                        LOG_ALLOW(GLOBAL, LOG_WARNING, "DBG: Not enough slots to add Ux/Uy/Uz debug fields.\n");
                        PetscFree(Ux); PetscFree(Uy); PetscFree(Uz);
                    }
 
                    // Mid-plane CSV + AoS vs NATURAL compare (component X) 
                    
                    // Gather NATURAL again (small cost, but isolated and clear) 
                    PetscInt   Ng = 0;
                    double    *nat_d = NULL;
                    DM         dmU = NULL;
                    DMDALocalInfo infU;
                    ierr = VecGetDM(field_vec, &dmU); CHKERRQ(ierr);
                    ierr = DMDAGetLocalInfo(dmU, &infU); CHKERRQ(ierr);
 
                    const PetscInt M=infU.mx, N=infU.my, P=infU.mz;
                    const PetscInt mx = meta.mx, my = meta.my, mz = meta.mz;
                    const PetscInt iInnerMid = mx/2; // interior index [0..mx-1] 
                    const PetscInt iGlob     = iInnerMid;
 
                    ierr = VecToArrayOnRank0(field_vec, &Ng, &nat_d); CHKERRQ(ierr);
 
                    if (nat_d) {
                        const PetscScalar *nar = (const PetscScalar*)nat_d;
                        const char *base = pps->output_prefix;
                        char fn[512], fnc[512];
                        snprintf(fn,  sizeof(fn),  "%s_%05" PetscInt_FMT "_iMid.csv",       base, ti);
                        snprintf(fnc, sizeof(fnc), "%s_%05" PetscInt_FMT "_iMid_compare.csv", base, ti);
 
                        FILE *fp  = fopen(fn,  "w");
                        FILE *fpc = fopen(fnc, "w");
                        if (fp)  fprintf(fp,  "jInner,kInner,Ux,Uy,Uz\n");
                        if (fpc) fprintf(fpc, "jInner,kInner,Ux_AoS,Ux_NAT,abs_diff\n");
 
                        double maxAbsDiff = 0.0, sumAbs = 0.0;
                        PetscInt count = 0;
 
                        for (PetscInt kInner = 0; kInner < mz; ++kInner) {
                            const PetscInt k = kInner;
                            for (PetscInt jInner = 0; jInner < my; ++jInner) {
                                const PetscInt j = jInner;
 
                                // AoS tuple index 
                                const PetscInt t = iInnerMid + mx * (jInner + my * kInner);
                                const PetscScalar ux = a[3*t+0], uy = a[3*t+1], uz = a[3*t+2];
 
                                if (fp) fprintf(fp, "%d,%d,%.15e,%.15e,%.15e\n",
                                                (int)jInner,(int)kInner,(double)ux,(double)uy,(double)uz);
 
                                // NATURAL base for (iGlob,j,k) 
                                const PetscInt baseNat = 3 * (((k)*N + j)*M + iGlob);
                                const PetscScalar uxN = nar[baseNat + 0];
 
                                const double diff = fabs((double)ux - (double)uxN);
                                if (diff > maxAbsDiff) maxAbsDiff = diff;
                                sumAbs += diff; ++count;
 
                                if (fpc) fprintf(fpc, "%d,%d,%.15e,%.15e,%.15e\n",
                                                    (int)jInner,(int)kInner,(double)ux,(double)uxN,diff);
                            } // for jInner
                        } // for kInner
                        if (fp)  fclose(fp);
                        if (fpc) fclose(fpc);
 
                        if (count > 0) {
                            const double meanAbs = sumAbs / (double)count;
                            LOG_ALLOW(GLOBAL, LOG_INFO,
                                        "PETSc-Vec vs AoS (i-mid, Ux): max|Δ|=%.6e, mean|Δ|=%.6e  -> CSV: %s\n",
                                        maxAbsDiff, meanAbs, fnc);
                            LOG_ALLOW(GLOBAL, LOG_INFO, "Wrote i-mid plane CSV: %s\n", fn);
                        } // if count>0
 
                        ierr = PetscFree(nat_d); CHKERRQ(ierr);
                        
                        
                    } // if nat_d
 
                } // if a && num_components==3 && npts>0
            } // if Ucat_nodal        
            // --- END DEBUG BLOCK (Ucat_nodal) ---
            */
            
            meta.num_point_data_fields++; /* count the main field we just added */
            field_name = strtok(NULL, ",");
        }
 
        ierr = PetscFree(fields_copy); CHKERRQ(ierr);
 
        // --- DEBUG: Add sanity fields i_idx, j_idx, k_idx and x_pos, y_pos, z_pos ---
        // These are the logical indices and physical coordinates of each point in the subsampled grid.
        // They can be used to verify the grid structure and orientation in visualization tools.
        // They are added as scalar fields with names "i_idx", "j_idx", "k_idx" and "x_pos", "y_pos", "z_pos".
        // Note: these are only added if there is room in the MAX_POINT_DATA_FIELDS limit.
        // They are allocated and owned here, and will be freed below.
        // They are in the same linearization order as meta.coords (AoS x,y,z by point).
        /*
        // Append sanity fields i/j/k indices and coordinates 
 
        // Build i/j/k and x/y/z (length = npoints); these match the same linearization as coords 
        const PetscInt n = meta.npoints;
        PetscScalar *i_idx=NULL,*j_idx=NULL,*k_idx=NULL,*x_pos=NULL,*y_pos=NULL,*z_pos=NULL;
 
        ierr = PetscMalloc1(n, &i_idx); CHKERRQ(ierr);
        ierr = PetscMalloc1(n, &j_idx); CHKERRQ(ierr);
        ierr = PetscMalloc1(n, &k_idx); CHKERRQ(ierr);
        ierr = PetscMalloc1(n, &x_pos); CHKERRQ(ierr);
        ierr = PetscMalloc1(n, &y_pos); CHKERRQ(ierr);
        ierr = PetscMalloc1(n, &z_pos); CHKERRQ(ierr);
 
        // coords is length 3*n, AoS: (x,y,z) by point 
        const PetscScalar *c = (const PetscScalar*)meta.coords;
        for (PetscInt k = 0; k < meta.mz; ++k) {
            for (PetscInt j = 0; j < meta.my; ++j) {
                for (PetscInt i = 0; i < meta.mx; ++i) {
                    const PetscInt t = i + meta.mx * (j + meta.my * k);
                    i_idx[t] = (PetscScalar)i;
                    j_idx[t] = (PetscScalar)j;
                    k_idx[t] = (PetscScalar)k;
                    x_pos[t] = c[3*t+0];
                    y_pos[t] = c[3*t+1];
                    z_pos[t] = c[3*t+2];
                }
            }
        }
 
        const char *nf[6] = {"i_idx","j_idx","k_idx","x_pos","y_pos","z_pos"};
        PetscScalar *arrs[6] = {i_idx,j_idx,k_idx,x_pos,y_pos,z_pos};
        for (int s=0; s<6; ++s) {
            if (meta.num_point_data_fields < MAX_POINT_DATA_FIELDS) {
                VTKFieldInfo *f = &meta.point_data_fields[meta.num_point_data_fields++];
                strncpy(f->name, nf[s], MAX_VTK_FIELD_NAME_LENGTH-1);
                f->name[MAX_VTK_FIELD_NAME_LENGTH-1] = '\0';
                f->num_components = 1;
                f->data = arrs[s];
            } else {
                LOG_ALLOW(GLOBAL, LOG_WARNING, "Sanity field '%s' dropped: MAX_POINT_DATA_FIELDS reached.\n", nf[s]);
                PetscFree(arrs[s]);
            }
        }
        LOG_ALLOW(GLOBAL, LOG_INFO, "DBG: Added sanity fields i_idx/j_idx/k_idx and x_pos/y_pos/z_pos.\n");
        */
        // --- END DEBUG BLOCK (sanity fields) ---
 
        /* Field summary */
        LOG_ALLOW(GLOBAL, LOG_INFO, "PointData fields to write: %d\n", (int)meta.num_point_data_fields);
        for (PetscInt ii=0; ii<meta.num_point_data_fields; ++ii) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "  # %2" PetscInt_FMT "  Field Name = %s  Components = %d\n",
                      ii, meta.point_data_fields[ii].name, (int)meta.point_data_fields[ii].num_components);
        }
    } // if rank 0
 
    /* 4) Write the VTS */
    snprintf(filename, sizeof(filename), "%s_%05" PetscInt_FMT ".vts", pps->output_prefix, ti);
    ierr = CreateVTKFileFromMetadata(filename, &meta, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Eulerian File Writing for ti = %" PetscInt_FMT " Complete ---\n", ti);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ParticleDataProcessingPipeline()

PetscErrorCode ParticleDataProcessingPipeline	(	UserCtx *	user,
		PostProcessParams *	pps
	)

Parses and executes the particle pipeline using a robust two-pass approach.

This function ensures correctness and efficiency by separating field registration from kernel execution.

PASS 1 (Registration): The pipeline string is parsed to identify all new fields that will be created. These fields are registered with the DMSwarm.

Finalize: After Pass 1, DMSwarmFinalizeFieldRegister is called exactly once if any new fields were added, preparing the swarm's memory layout.

PASS 2 (Execution): The pipeline string is parsed again, and this time the actual compute kernels are executed, filling the now-valid fields.

Parameters

user	The UserCtx containing the DMSwarm.
pps	The PostProcessParams struct containing the particle_pipeline string.

Returns

PetscErrorCode

Definition at line 524 of file postprocessor.c.

{
    PetscErrorCode ierr;
    char *pipeline_copy, *step_token, *step_saveptr;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // --- Timestep Setup: Synchronize post_swarm size ---
    PetscInt n_global_source;
    ierr = DMSwarmGetSize(user->swarm, &n_global_source); CHKERRQ(ierr);
 
    // Resize post_swarm to match source swarm
    ierr = ResizeSwarmGlobally(user->post_swarm, n_global_source); CHKERRQ(ierr);
 
    if (pps->particle_pipeline[0] == '\0') {
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Particle Data Transformation Pipeline ---\n");
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Particle Pipeline string: [%s]\n", pps->particle_pipeline);
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Executing compute kernels...\n");
    ierr = PetscStrallocpy(pps->particle_pipeline, &pipeline_copy); CHKERRQ(ierr);
    step_token = strtok_r(pipeline_copy, ";", &step_saveptr);
    while (step_token) {
        TrimWhitespace(step_token);
        if (strlen(step_token) == 0) { step_token = strtok_r(NULL, ";", &step_saveptr); continue; }
 
        char *keyword = strtok(step_token, ":");
        char *args_str = strtok(NULL, "");
        TrimWhitespace(keyword);
        if (args_str) TrimWhitespace(args_str);
        
        LOG_ALLOW(GLOBAL, LOG_INFO, "Executing Particle Transformation: '%s' on args: '%s'\n", keyword, args_str ? args_str : "None");
        
        if (strcasecmp(keyword, "ComputeSpecificKE") == 0) {
            char *velocity_field = strtok(args_str, ">");
            char *ske_field = strtok(NULL, ">");
            TrimWhitespace(velocity_field); TrimWhitespace(ske_field);
            
            ierr = ComputeSpecificKE(user, velocity_field, ske_field); CHKERRQ(ierr);
        }
        else {
             LOG_ALLOW(GLOBAL, LOG_WARNING, "Unknown particle transformation keyword '%s'. Skipping.\n", keyword);
        }
 
        step_token = strtok_r(NULL, ";", &step_saveptr);
    }
    ierr = PetscFree(pipeline_copy); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Particle Data Transformation Pipeline Complete ---\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteParticleFile()

PetscErrorCode WriteParticleFile	(	UserCtx *	user,
		PostProcessParams *	pps,
		PetscInt	ti
	)

Writes particle data to a VTP file using the Prepare-Write-Cleanup pattern.

Definition at line 587 of file postprocessor.c.

{
    PetscErrorCode ierr;
    VTKMetaData    part_meta;
    char           filename[MAX_FILENAME_LENGTH];
    PetscInt       n_total_particles_before_subsample;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // These checks can be done on all ranks
    if (!pps->outputParticles || pps->particle_fields[0] == '\0') {
        PetscFunctionReturn(0);
    }
    PetscInt n_global;
    ierr = DMSwarmGetSize(user->swarm, &n_global); CHKERRQ(ierr);
    if (n_global == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Swarm is empty for ti=%" PetscInt_FMT ". Skipping particle file write.\n", ti);
        PetscFunctionReturn(0);
    }
 
    ierr = PetscMemzero(&part_meta, sizeof(VTKMetaData)); CHKERRQ(ierr);
 
    // --- 1. PREPARE (Collective Call) ---
    ierr = PrepareOutputParticleData(user, pps, &part_meta, &n_total_particles_before_subsample); CHKERRQ(ierr);
 
    // --- 2. WRITE and CLEANUP (Rank 0 only) ---
    if (user->simCtx->rank == 0) {
        if (part_meta.npoints > 0) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting VTP Particle File Writing for ti = %" PetscInt_FMT " (writing %" PetscInt_FMT " of %" PetscInt_FMT " particles) ---\n",
                      ti, part_meta.npoints, n_total_particles_before_subsample);
 
            /* Field summary */
            LOG_ALLOW(GLOBAL, LOG_INFO, "Particle Data fields to write: %d\n", (int)part_meta.num_point_data_fields);
            for (PetscInt ii=0; ii<part_meta.num_point_data_fields; ++ii) {
                LOG_ALLOW(GLOBAL, LOG_INFO, "  # %2" PetscInt_FMT "  Field Name = %s  Components = %d\n",
                          ii, part_meta.point_data_fields[ii].name, (int)part_meta.point_data_fields[ii].num_components);
            }
 
            snprintf(filename, sizeof(filename), "%s_%05" PetscInt_FMT ".vtp", pps->particle_output_prefix, ti);
            ierr = CreateVTKFileFromMetadata(filename, &part_meta, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
        } else {
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "No particles to write at ti=%" PetscInt_FMT " after subsampling. Skipping.\n", ti);
        }
 
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Particle File Writing for ti = %" PetscInt_FMT " Complete ---\n", ti);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 642 of file postprocessor.c.

{
    PetscErrorCode    ierr;
    SimCtx            *simCtx = NULL;
 
    // === I. INITIALIZE PETSC & MPI ===========================================
    ierr = PetscInitialize(&argc, &argv, (char *)0, "Unified Post-Processing Tool"); CHKERRQ(ierr);
 
    // === II. CONFIGURE SIMULATION & POST-PROCESSING CONTEXTS =================
    ierr = CreateSimulationContext(argc, argv, &simCtx); CHKERRQ(ierr);
    // === IIB. SET EXECUTION MODE (SOLVER vs POST-PROCESSOR) =====
    simCtx->exec_mode = EXEC_MODE_POSTPROCESSOR;
    // == IIC. CONFIGURE SIMULATION ENVIRONMENT & DIRECTORIES =====
  ierr = SetupSimulationEnvironment(simCtx); CHKERRQ(ierr);
    // === III. SETUP GRID & DATA STRUCTURES ===================================
    ierr = SetupGridAndSolvers(simCtx); CHKERRQ(ierr);
    // === IV. SETUP DOMAIN DECOMPOSITION INFORMATION =========================
    ierr = SetupDomainRankInfo(simCtx); CHKERRQ(ierr);
    // === V. SETUP BOUNDARY CONDITIONS ====================================
    ierr = SetupBoundaryConditions(simCtx); CHKERRQ(ierr);
    // === VI. SETUP USER CONTEXT & DATA STRUCTURES ============================
    // Get the finest-level user context, as this is where we'll load data
    UserCtx *user = simCtx->usermg.mgctx[simCtx->usermg.mglevels-1].user;
    PostProcessParams *pps = simCtx->pps;
 
    // === VI. PARTICLE INITIALIZATION (if needed) ============================
    if(pps->outputParticles) {
        if(simCtx->np > 0){
            ierr = InitializeParticleSwarm(simCtx); CHKERRQ(ierr);
            // Create a post-processing specific DMSwarm
            ierr = SetupPostProcessSwarm(user,pps); CHKERRQ(ierr);
        }else{
            SETERRQ(PETSC_COMM_SELF,1,"Particle output requested but np=0. Please set np>0 to enable particle output.");
        }
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "=============================================================\n");
 
 
    // === VII. MAIN POST-PROCESSING LOOP ======================================
    for (PetscInt ti = pps->startTime; ti <= pps->endTime; ti += pps->timeStep) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "--- Processing Time Step %" PetscInt_FMT " ---\n", ti);
 
        // 1. Load Data (UpdateLocalGhosts is called inside the kernels)
        ierr = ReadSimulationFields(user, ti); CHKERRQ(ierr);
        
        // 2. Transform Data
        ierr = EulerianDataProcessingPipeline(user, pps); CHKERRQ(ierr);
 
        // 3. Write Output
        ierr = WriteEulerianFile(user, pps, ti); CHKERRQ(ierr);
 
        if(pps->outputParticles) {
            // 1. Resize swarm based on particle count in this timestep's file
            ierr = PreCheckAndResizeSwarm(user, ti, pps->particleExt); CHKERRQ(ierr);
 
            // 2. Load particle data into the correctly sized swarm
            ierr = ReadAllSwarmFields(user, ti); CHKERRQ(ierr);
            
            // 3. Transform particle data
            ierr = ParticleDataProcessingPipeline(user, pps); CHKERRQ(ierr);
            
            // 4. Write particle output
           ierr = WriteParticleFile(user, pps, ti); CHKERRQ(ierr);            
        }
 
        if(simCtx->rank == 0){
            PetscInt StepsToRun = pps->endTime - pps->startTime;
            PetscReal currentTime = (PetscReal)ti*simCtx->dt;
            PrintProgressBar(ti-1,pps->startTime,StepsToRun,currentTime);
            if(get_log_level()>LOG_ERROR)PetscPrintf(PETSC_COMM_SELF,"\n");
        }
    }
 
    // After the loop, print the 100% complete bar on rank 0 and add a newline
    // to ensure subsequent terminal output starts on a fresh line.
    if (simCtx->rank == 0) {
        PetscInt endTime = pps->endTime-1; // needs to be verified.
        PetscInt StepsToRun = pps->endTime - pps->startTime;
        PetscReal endTimeValue = (PetscReal)pps->endTime*simCtx->dt;
        PrintProgressBar(endTime, pps->startTime, StepsToRun, endTimeValue);
        PetscPrintf(PETSC_COMM_SELF, "\n");
        fflush(stdout);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "=============================================================\n");
    LOG_ALLOW(GLOBAL, LOG_INFO, "Post-processing finished successfully.\n");
    
 
    // === VIII. FINALIZE =========================================================
   // ierr = FinalizeSimulation(simCtx); CHKERRQ(ierr);
    ierr = ProfilingFinalize(simCtx); CHKERRQ(ierr);
    ierr = PetscFinalize();
    return ierr;
}

Here is the call graph for this function: