postprocessor.h File Reference

#include "io.h"
#include "variables.h"
#include "logging.h"
#include "ParticleSwarm.h"
#include "interpolation.h"
#include "grid.h"
#include "setup.h"
#include "Metric.h"
#include "postprocessing_kernels.h"
#include "vtk_io.h"
#include "particle_statistics.h"

Include dependency graph for postprocessor.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
PetscErrorCode	SetupPostProcessSwarm (UserCtx *user, PostProcessParams *pps)
	Creates a new, dedicated DMSwarm for post-processing tasks.
PetscErrorCode	WriteEulerianFile (UserCtx *user, PostProcessParams *pps, PetscInt ti)
	Orchestrates the writing of a combined, multi-field VTK file for a single time step.
PetscErrorCode	EulerianDataProcessingPipeline (UserCtx *user, PostProcessParams *pps)
	Parses the processing pipeline string and executes the requested kernels.
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.
PetscErrorCode	GlobalStatisticsPipeline (UserCtx *user, PostProcessParams *pps, PetscInt ti)
	Executes the global statistics pipeline, computing aggregate reductions over all particles.
PetscErrorCode	FieldStatisticsPipeline (UserCtx *user, PostProcessParams *pps, PetscInt ti)
	Derives and writes field statistics for the windows a recipe requests.

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

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

Local to this translation unit.

Definition at line 21 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 *input_field = strtok(args_str, ">");
            char *output_field = strtok(NULL, ">");
            output_field_dimensions = 1;  // SKE is scalar
            if (!input_field) SETERRQ(PETSC_COMM_SELF, 1, "Error (ComputeSpecificKE): Missing input field in 'in>out' syntax.");
            if (!output_field) SETERRQ(PETSC_COMM_SELF, 1, "Error (ComputeSpecificKE): Missing output field in 'in>out' syntax.");
            TrimWhitespace(input_field);
            TrimWhitespace(output_field);
            if (strlen(input_field) == 0) SETERRQ(PETSC_COMM_SELF, 1, "Error (ComputeSpecificKE): Empty input field name.");
            if (strlen(output_field) == 0) SETERRQ(PETSC_COMM_SELF, 1, "Error (ComputeSpecificKE): Empty output field name.");
            // 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 input='%s').\n", output_field, input_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 registered particle fields for the post-processing swarm.\n");
    } else {
        LOG_ALLOW(GLOBAL, LOG_INFO, "No custom particle fields requested; finalizing an empty post-processing swarm for safe use.\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:

WriteEulerianFile()

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

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

This function is the primary driver for generating output. It performs these steps:

1. Prepares the subsampled coordinate array required for the legacy grid format.
2. Parses the user-requested list of fields from the configuration.
3. For each field, prepares a corresponding subsampled data array.
4. Assembles all prepared arrays into a single VTKMetaData struct.
5. Calls the low-level VTK writer to generate the final .vts file.
6. Frees all temporary memory allocated during the preparation phase.

Parameters

user	The UserCtx for the finest grid level.
pps	The post-processing configuration struct.
ti	The current time step index.

Returns

PetscErrorCode

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

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

See also

WriteEulerianFile()

Definition at line 201 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);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting VTK File Writing for ti = %" PetscInt_FMT " ---\n", ti);
 
    /* 1) Metadata init */
    /* 2) Metadata and coordinates, through the shared assembly the statistics
     *    stage also uses, so both producers build a file the same way. */
    ierr = BeginStructuredVTKOutput(user, &meta); 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) Field preparation is collective because each append gathers a
     * distributed Vec. Rank zero alone owns the packed output buffers, but
     * every rank must resolve the same field list and enter each append. */
    {
        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;
            }
 
            // --- Add field to metadata ---
            ierr = AppendStructuredVTKField(user, &meta, field_name, field_vec, num_components); 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) ---
            */
            
            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) ---
 
        if (user->simCtx->rank == 0) {
            /* 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);
            }
        }
    }
 
    /* 4) Write the VTS */
    ierr = PetscSNPrintf(filename, sizeof(filename), "%s_%05" PetscInt_FMT ".vts", pps->output_prefix, ti); CHKERRQ(ierr);
    ierr = FinishStructuredVTKOutput(&meta, filename); 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:

EulerianDataProcessingPipeline()

PetscErrorCode EulerianDataProcessingPipeline	(	UserCtx *	user,
		PostProcessParams *	pps
	)

Parses the processing pipeline string and executes the requested kernels.

Parameters

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

Returns

PetscErrorCode

Parses the processing pipeline string and executes the requested kernels.

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

See also

EulerianDataProcessingPipeline()

Definition at line 107 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");
        PROFILE_FUNCTION_END;
        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, "DimensionalizeAllLoadedFields") == 0) {
            /* Emitted by global_operations.dimensionalize. It had no dispatch branch,
             * so the option was accepted, serialized, and then silently skipped. */
            ierr = DimensionalizeAllLoadedFields(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:

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

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

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

See also

ParticleDataProcessingPipeline()

Definition at line 626 of file postprocessor.c.

{
    PetscErrorCode ierr;
    char *pipeline_copy, *step_token, *step_saveptr;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    if (pps->particle_pipeline[0] == '\0') {
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    // --- Timestep Setup: Synchronize post_swarm size ---
    PetscInt n_local_source;
    ierr = DMSwarmGetLocalSize(user->swarm, &n_local_source); CHKERRQ(ierr);
 
    // Derived entries use the same local index as their source particle.
    ierr = DMSwarmSetLocalSizes(user->post_swarm, n_local_source, -1); CHKERRQ(ierr);
 
    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) {
            if (!args_str) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ComputeSpecificKE requires 'input_field>output_field' arguments.");
            char *velocity_field = strtok(args_str, ">");
            char *ske_field = strtok(NULL, ">");
            if (!velocity_field || !ske_field) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "ComputeSpecificKE requires 'input_field>output_field' arguments.");
            }
            TrimWhitespace(velocity_field); TrimWhitespace(ske_field);
            if (strlen(velocity_field) == 0 || strlen(ske_field) == 0) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "ComputeSpecificKE does not allow empty input/output field names.");
            }
            
            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.

Parameters

user	Primary UserCtx input for the operation.
pps	Post-processing configuration for the operation.
ti	Source timestep written into the VTP filename and metadata.

Returns

PetscErrorCode 0 on success.

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

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

See also

WriteParticleFile()

Definition at line 749 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') {
        PROFILE_FUNCTION_END;
        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);
        PROFILE_FUNCTION_END;
        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);
            }
 
            ierr = PetscSNPrintf(filename, sizeof(filename), "%s_%05" PetscInt_FMT ".vtp", pps->particle_output_prefix, ti); CHKERRQ(ierr);
            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);
        }
 
        for (PetscInt field = 0; field < part_meta.num_point_data_fields; ++field) {
            ierr = PetscFree(part_meta.point_data_fields[field].data); CHKERRQ(ierr);
        }
        ierr = PetscFree(part_meta.coords); CHKERRQ(ierr);
        ierr = PetscFree(part_meta.connectivity); CHKERRQ(ierr);
        ierr = PetscFree(part_meta.offsets); CHKERRQ(ierr);
    }
 
    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:

GlobalStatisticsPipeline()

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

Executes the global statistics pipeline, computing aggregate reductions over all particles.

Parses the semicolon-delimited pps->statistics_pipeline string and dispatches to the appropriate kernel (e.g. ComputeParticleMSD). Each kernel appends one row to its own CSV file and logs a summary via LOG_INFO. All MPI reductions happen inside each kernel. This pipeline is independent of the per-particle VTK pipeline; it produces no .vtp output.

Parameters

user	The UserCtx containing the primary DMSwarm (user->swarm).
pps	The PostProcessParams containing statistics_pipeline and statistics_output_prefix.
ti	Current time-step index (passed through to kernels for time computation).

Returns

PetscErrorCode

Executes the global statistics pipeline, computing aggregate reductions over all particles.

Local to this translation unit.

Definition at line 698 of file postprocessor.c.

{
    PetscErrorCode ierr;
    char *pipeline_copy, *step_token, *step_saveptr;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (pps->statistics_pipeline[0] == '\0') { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    PetscInt n_global;
    ierr = DMSwarmGetSize(user->swarm, &n_global); CHKERRQ(ierr);
    if (n_global == 0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Global Statistics Pipeline ---\n");
 
    ierr = PetscStrallocpy(pps->statistics_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, ":");
        TrimWhitespace(keyword);
 
        if (strcasecmp(keyword, "ComputeMSD") == 0) {
            ierr = ComputeParticleMSD(user, pps->statistics_output_prefix, ti); CHKERRQ(ierr);
        } else {
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "Unknown statistics keyword '%s'. Skipping.\n", keyword);
        }
        /* Additional kernels should add else-if branches here when implemented. */
 
        step_token = strtok_r(NULL, ";", &step_saveptr);
    }
    ierr = PetscFree(pipeline_copy); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Global Statistics Pipeline Complete ---\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FieldStatisticsPipeline()

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

Derives and writes field statistics for the windows a recipe requests.

Runs once per processed step, beside the particle statistics pipeline, because a window's accumulated state is a property of the bundle at that step. A recipe spanning several steps therefore produces a convergence history rather than the same picture repeated; pinning field_statistics_source_step collapses it back to a single bundle.

Each window is written to its own file. The VTK point-data cap is per file, and one window carrying every output already fills most of it.

Does nothing when no window is requested. A window the run does not configure is fatal; a window that has not yet accumulated a sample is skipped with a note, since a recipe covering a whole run legitimately reaches bundles from before it began.

Parameters

[in]	user	Finest-level block context holding accumulators and staging fields.
[in]	pps	Post-processing recipe naming the windows, outputs, and formats.
[in]	ti	Step being processed; names the output and selects the bundle.

Returns

Zero on success, or a PETSc error.

Derives and writes field statistics for the windows a recipe requests.

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

See also

FieldStatisticsPipeline()

Definition at line 485 of file postprocessor.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = NULL;
    char *windows_copy = NULL;
    char *window_name = NULL;
    PetscInt source_step = 0;
    PetscBool want_vtk = PETSC_FALSE, want_csv = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    if (pps->field_statistics_windows[0] == '\0') { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
    simCtx = user->simCtx;
 
    PetscCheck(FieldStatisticsIsActive(simCtx), PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
               "Field-statistics post-processing was requested, but the run's control configures "
               "no statistics window.");
 
    {
        char formats[MAX_FIELD_LIST_LENGTH];
        char *token = NULL;
 
        ierr = PetscStrncpy(formats, pps->field_statistics_formats, sizeof(formats)); CHKERRQ(ierr);
        token = strtok(formats, ",");
        while (token) {
            TrimWhitespace(token);
            if (!strcasecmp(token, "vtk")) want_vtk = PETSC_TRUE;
            else if (!strcasecmp(token, "csv")) want_csv = PETSC_TRUE;
            else LOG_ALLOW(GLOBAL, LOG_WARNING,
                           "Unknown field-statistics format '%s'. Known formats are vtk and csv.\n", token);
            token = strtok(NULL, ",");
        }
    }
    if (!want_vtk && !want_csv) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    /* An explicit source step pins every processed step to one bundle; otherwise each
     * step derives from its own, which is what turns a multi-step recipe into a
     * convergence history rather than the same picture repeated. */
    source_step = (pps->field_statistics_source_step >= 0) ? pps->field_statistics_source_step : ti;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Starting Field Statistics Pipeline (step %" PetscInt_FMT ") ---\n",
              source_step);
    simCtx->fieldStatisticsContinue = PETSC_TRUE;
    ierr = RestoreFieldStatisticsState(simCtx, source_step); CHKERRQ(ierr);
 
    ierr = PetscStrallocpy(pps->field_statistics_windows, &windows_copy); CHKERRQ(ierr);
    window_name = strtok(windows_copy, ",");
    while (window_name) {
        PetscInt window_index = -1;
        const PicurvWindow *window = NULL;
 
        TrimWhitespace(window_name);
        if (!*window_name) { window_name = strtok(NULL, ","); continue; }
 
        for (PetscInt w = 0; w < simCtx->fieldStatisticsWindowCount; ++w) {
            PetscBool matches = PETSC_FALSE;
 
            ierr = PetscStrcmp(simCtx->fieldStatisticsWindows[w].definition.name,
                               window_name, &matches); CHKERRQ(ierr);
            if (matches) { window_index = w; break; }
        }
        if (window_index < 0) {
            ierr = PetscFree(windows_copy); CHKERRQ(ierr);
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                    "Field-statistics post-processing requested window '%s', which this run does "
                    "not configure.", window_name);
        }
        window = &simCtx->fieldStatisticsWindows[window_index];
 
        /* A window that has not started by this step is not an error: a recipe
         * spanning a whole run legitimately reaches bundles from before it began. */
        if (window->sample_count <= 0) {
            LOG_ALLOW(GLOBAL, LOG_INFO,
                      "Statistics window '%s' had accumulated no sample by step %" PetscInt_FMT
                      "; nothing to derive yet.\n", window->definition.name, source_step);
            window_name = strtok(NULL, ",");
            continue;
        }
 
        if (want_vtk) {
            VTKMetaData meta;
            PetscInt derived_count = 0;
            char filename[PETSC_MAX_PATH_LEN];
 
            ierr = PicurvWindowDerivedCount(&window->definition,
                                            &user->fieldStatisticsStorage[window_index],
                                            pps->field_statistics_outputs, &derived_count); CHKERRQ(ierr);
            /* An output kind resolves against what the window accumulated, so asking
             * for stresses from a means-only window yields nothing. Report it the way
             * an unrecognized Eulerian output field is reported, rather than writing a
             * quietly short file. */
            if (derived_count == 0) {
                LOG_ALLOW(LOCAL, LOG_WARNING,
                          "Outputs '%s' produce no field for window '%s'; it accumulates none of "
                          "the state they need. Skipping.\n",
                          pps->field_statistics_outputs, window->definition.name);
            }
            /* One file per window. The point-data cap is per file, and a single
             * window with every output already fills most of it. */
            ierr = BeginStructuredVTKOutput(user, &meta); CHKERRQ(ierr);
            for (PetscInt index = 0; index < derived_count; ++index) {
                char name[MAX_VTK_FIELD_NAME_LENGTH];
                Vec nodal = NULL;
                PetscInt components = 0;
 
                ierr = ComputeWindowStatisticNodal(user, window_index,
                                                   pps->field_statistics_outputs, index,
                                                   name, sizeof(name), &nodal, &components); CHKERRQ(ierr);
                /* The staging vector is reused for the next field, which is safe
                 * because appending copies the values out. */
                ierr = AppendStructuredVTKField(user, &meta, name, nodal, components); CHKERRQ(ierr);
            }
            ierr = PetscSNPrintf(filename, sizeof(filename), "%s_statistics_%s_%05" PetscInt_FMT ".vts",
                                 pps->output_prefix, window->definition.name, ti); CHKERRQ(ierr);
            LOG_ALLOW(GLOBAL, LOG_INFO, "Wrote %d derived field(s) for window '%s' to %s\n",
                      (int)meta.num_point_data_fields, window->definition.name, filename);
            ierr = FinishStructuredVTKOutput(&meta, filename); CHKERRQ(ierr);
        }
 
        if (want_csv) {
            ierr = ComputeWindowStatisticsSummary(user, window_index, pps->output_prefix,
                                                  ti); CHKERRQ(ierr);
        }
 
        window_name = strtok(NULL, ",");
    }
    ierr = PetscFree(windows_copy); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Field Statistics Pipeline Complete ---\n");
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: