postprocessor.c

Go to the documentation of this file.

/**
 * @file postprocessor.c
 * @brief Offline post-processing tool driving the derived-output pipelines.
 *
 * Reads committed checkpoint bundles over a step range and dispatches the
 * Eulerian, Lagrangian, particle-statistics, and field-statistics pipelines that
 * a recipe selected, writing one output family per enabled pipeline.
 */
 
#include "postprocessor.h" // Use our new header
#include "statistics_accumulator.h"
#include "statistics_window.h"
 
 
#undef __FUNCT__
#define __FUNCT__ "SetupPostProcessSwarm"
/**
 * @brief Internal helper implementation: `SetupPostProcessSwarm()`.
 * @details Local to this translation unit.
 */
PetscErrorCode SetupPostProcessSwarm(UserCtx* user, PostProcessParams* pps)
{
    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);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "EulerianDataProcessingPipeline"
/**
 * @brief Implementation of \ref EulerianDataProcessingPipeline().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/postprocessor.h`.
 * @see EulerianDataProcessingPipeline()
 */
PetscErrorCode EulerianDataProcessingPipeline(UserCtx* user, PostProcessParams* pps)
{
    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);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "WriteEulerianFile"
/**
 * @brief Implementation of \ref WriteEulerianFile().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/postprocessor.h`.
 * @see WriteEulerianFile()
 */
PetscErrorCode WriteEulerianFile(UserCtx* user, PostProcessParams* pps, PetscInt ti)
{
    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);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "FieldStatisticsPipeline"
/**
 * @brief Implementation of \ref FieldStatisticsPipeline().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/postprocessor.h`.
 * @see FieldStatisticsPipeline()
 */
PetscErrorCode FieldStatisticsPipeline(UserCtx *user, PostProcessParams *pps, PetscInt ti)
{
    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);
}
 
#undef __FUNCT__
#define __FUNCT__ "ParticleDataProcessingPipeline"
/**
 * @brief Implementation of \ref ParticleDataProcessingPipeline().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/postprocessor.h`.
 * @see ParticleDataProcessingPipeline()
 */
PetscErrorCode ParticleDataProcessingPipeline(UserCtx* user, PostProcessParams* pps)
{
    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);
}
 
#undef __FUNCT__
#define __FUNCT__ "GlobalStatisticsPipeline"
/**
 * @brief Internal helper implementation: `GlobalStatisticsPipeline()`.
 * @details Local to this translation unit.
 */
PetscErrorCode GlobalStatisticsPipeline(UserCtx *user, PostProcessParams *pps, PetscInt ti)
{
    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);
}
 
#undef __FUNCT__
#define __FUNCT__ "WriteParticleFile"
/**
 * @brief Implementation of \ref WriteParticleFile().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/postprocessor.h`.
 * @see WriteParticleFile()
 */
PetscErrorCode WriteParticleFile(UserCtx* user, PostProcessParams* pps, PetscInt ti)
{
    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);
}
 
#undef __FUNCT__
#define __FUNCT__ "main"
#ifndef PICURV_POSTPROCESSOR_NO_MAIN
/**
 * @brief Entry point for the postprocessor executable.
 * @details Initializes PETSc, loads post-processing inputs, executes the
 *          requested pipelines, and finalizes runtime resources before exit.
 */
int main(int argc, char **argv)
{
    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);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Postprocessor MPI processes: %d\n", (int)simCtx->size); 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. CAPABILITY DISPATCH ============================================
    // Each stage declares what it needs rather than inferring it from another
    // stage's configuration. Field statistics are Eulerian and must not require a
    // swarm: a turbulence run normally carries no particles at all.
    PetscBool needs_particle_stage = (pps->outputParticles || pps->particle_pipeline[0] != '\0' || pps->statistics_pipeline[0] != '\0') ? PETSC_TRUE : PETSC_FALSE;
    if(needs_particle_stage) {
        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 post-processing requested (particle output or particle statistics pipeline) "
                    "but np=0. Please set np>0 during solver run to enable particle post-processing.");
        }
    }
 
    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(needs_particle_stage) {
            // 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 (optional)
            if (pps->outputParticles) {
                ierr = WriteParticleFile(user, pps, ti); CHKERRQ(ierr);
            }
 
            // 5. Global statistical reductions (MSD, etc.) → CSV files
            ierr = GlobalStatisticsPipeline(user, pps, ti); CHKERRQ(ierr);
        }
 
        // 4. Accumulated Eulerian window statistics → derived fields and history.
        //    Eulerian and independent of the swarm, so it sits outside the particle
        //    block: a turbulence run normally carries no particles at all.
        ierr = FieldStatisticsPipeline(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");
        }
        ierr = RuntimeMemoryLogSample(simCtx, ti, "Post", "-"); CHKERRQ(ierr);
    }
 
    // 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 = RuntimeMemoryLogSample(simCtx, pps->endTime, "Final", "Complete"); CHKERRQ(ierr);
    ierr = ProfilingFinalize(simCtx); CHKERRQ(ierr);
    ierr = FinalizeSimulation(simCtx); CHKERRQ(ierr);
    ierr = PetscFinalize();
    return ierr;
}
#endif