postprocessing_kernels.h File Reference

#include "variables.h"
#include "logging.h"
#include "io.h"

Include dependency graph for postprocessing_kernels.h:

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

Go to the source code of this file.

Functions
PetscErrorCode	ComputeNodalAverage (UserCtx *user, const char *in_field_name, const char *out_field_name)
	Interpolates a cell-centered field to nodal locations using local stencil averaging.
PetscErrorCode	ComputeQCriterion (UserCtx *user)
	Computes the Q-criterion diagnostic from the local velocity-gradient tensor.
PetscErrorCode	NormalizeRelativeField (UserCtx *user, const char *relative_field_name)
	Normalizes pressure using the value at the configured logical grid point.
PetscErrorCode	DimensionalizeField (UserCtx *user, const char *field_name)
	Scales a specified field from non-dimensional to dimensional units in-place.
PetscErrorCode	ComputeSpecificKE (UserCtx *user, const char *velocity_field, const char *ske_field)
	Computes the specific kinetic energy (KE per unit mass) for each particle.
PetscErrorCode	ComputeDisplacement (UserCtx *user, const char *disp_field)
	Computes the displacement magnitude |r_i - r_0| for each particle (per-particle VTK kernel).
PetscErrorCode	ComputeWindowStatisticNodal (UserCtx *user, PetscInt window_index, const char *outputs, PetscInt output_index, char *out_name, size_t name_size, Vec *out_vec, PetscInt *out_components)
	Derives one accumulated statistic and converts it to nodal values.
PetscErrorCode	ComputeWindowStatisticsSummary (UserCtx *user, PetscInt window_index, const char *output_prefix, PetscInt ti)
	Appends one convergence row for an accumulated window to its CSV history.

Function Documentation

ComputeNodalAverage()

PetscErrorCode ComputeNodalAverage	(	UserCtx *	user,
		const char *	in_field_name,
		const char *	out_field_name
	)

Interpolates a cell-centered field to nodal locations using local stencil averaging.

The kernel reads the input field by name, computes nodal values, and stores the output in the named destination field. Both fields must already exist in the current UserCtx.

Parameters

[in,out]	user	Block-level context that owns the source and destination vectors.
[in]	in_field_name	Name of the input field to sample.
[in]	out_field_name	Name of the output field to populate.

Returns

PetscErrorCode 0 on success.

Interpolates a cell-centered field to nodal locations using local stencil averaging.

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

See also

ComputeNodalAverage()

Definition at line 99 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    FieldId        in_field_id;
    Vec            in_vec_local = NULL, out_vec_global = NULL;
    DM             dm_in = NULL, dm_out = NULL;
    PetscInt       dof = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeNodalAverage on '%s' -> '%s'.\n", in_field_name, out_field_name);
 
    // --- 1. Map string names to PETSc objects ---
    if (strcasecmp(in_field_name, "P") == 0)             { in_vec_local = user->lP;         dm_in = user->da;   dof = 1; }
    else if (strcasecmp(in_field_name, "Ucat") == 0)    { in_vec_local = user->lUcat;      dm_in = user->fda;  dof = 3; }
    else if (strcasecmp(in_field_name, "Psi") == 0)     { in_vec_local = user->lPsi;       dm_in = user->da;   dof = 1; }
    else if (strcasecmp(in_field_name, "Qcrit") == 0)   { in_vec_local = user->lQcrit;     dm_in = user->da;   dof = 1; }
    /* The staging pair carries derived statistics, which are config-counted and so
     * cannot be named by a compile-time member of their own. */
    else if (strcasecmp(in_field_name, "PostScalar") == 0) { in_vec_local = user->lPostScalar; dm_in = user->da;  dof = 1; }
    else if (strcasecmp(in_field_name, "PostVector") == 0) { in_vec_local = user->lPostVector; dm_in = user->fda; dof = 3; }
    // ... (add other fields as needed) ...
    else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Unknown input field name for nodal averaging: %s", in_field_name);
 
    if (strcasecmp(out_field_name, "P_nodal") == 0)      { out_vec_global = user->P_nodal;    dm_out = user->da; }
    else if (strcasecmp(out_field_name, "Ucat_nodal") == 0) { out_vec_global = user->Ucat_nodal; dm_out = user->fda; }
    else if (strcasecmp(out_field_name, "Psi_nodal") == 0)   { out_vec_global = user->Psi_nodal;  dm_out = user->da; }
    else if (strcasecmp(out_field_name, "Qcrit_nodal") == 0) { out_vec_global = user->Qcrit_nodal; dm_out = user->da; }
    else if (strcasecmp(out_field_name, "PostScalarNodal") == 0) { out_vec_global = user->PostScalarNodal; dm_out = user->da; }
    else if (strcasecmp(out_field_name, "PostVectorNodal") == 0) { out_vec_global = user->PostVectorNodal; dm_out = user->fda; }
    // ... (add other fields as needed) ...
    else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Unknown output field name for nodal averaging: %s", out_field_name);
 
    // --- 2. Ensure Input Data Ghosts are Up-to-Date ---
    ierr = FieldIdFromName(in_field_name, &in_field_id); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, in_field_id); CHKERRQ(ierr);
    /* The boundary node average reads the dummy cells, so show them before averaging. */
    if (get_log_level() == LOG_VERBOSE && is_function_allowed(__FUNCT__)) {
        ierr = LOG_FIELD_ANATOMY(user, in_field_id, "PreNodalAverage"); CHKERRQ(ierr);
    }
 
    // --- 3. Get DMDA info and array pointers ---
    DMDALocalInfo info;
    ierr = DMDAGetLocalInfo(dm_out, &info); CHKERRQ(ierr);
    /* Every owned output point is valid except the global high layout plane.
     * A rank interface is not a boundary: its +1 source value is in the halo. */
    const PetscInt i_end = PetscMin(info.xs + info.xm, info.mx - 1);
    const PetscInt j_end = PetscMin(info.ys + info.ym, info.my - 1);
    const PetscInt k_end = PetscMin(info.zs + info.zm, info.mz - 1);
 
    if (dof == 1) { // --- Scalar Field Averaging ---
        const PetscReal ***l_in_arr;
        PetscReal       ***g_out_arr;
        ierr = DMDAVecGetArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
 
        // Loop over the output NODE locations. The loop bounds match the required
        // size of the final subsampled grid.
        for (PetscInt k = info.zs; k < k_end; k++) {
            for (PetscInt j = info.ys; j < j_end; j++) {
                for (PetscInt i = info.xs; i < i_end; i++) {
                    g_out_arr[k][j][i] = 0.125 * (l_in_arr[k][j][i]     + l_in_arr[k][j][i+1] +
                                                  l_in_arr[k][j+1][i]   + l_in_arr[k][j+1][i+1] +
                                                  l_in_arr[k+1][j][i]   + l_in_arr[k+1][j][i+1] +
                                                  l_in_arr[k+1][j+1][i] + l_in_arr[k+1][j+1][i+1]);
                }
            }
        }
        ierr = DMDAVecRestoreArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
 
    } else if (dof == 3) { // --- Vector Field Averaging ---
        const Cmpnts ***l_in_arr;
        Cmpnts       ***g_out_arr;
        ierr = DMDAVecGetArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
 
        for (PetscInt k = info.zs; k < k_end; k++) {
            for (PetscInt j = info.ys; j < j_end; j++) {
                for (PetscInt i = info.xs; i < i_end; i++) {
                    g_out_arr[k][j][i].x = 0.125 * (l_in_arr[k][j][i].x + l_in_arr[k][j][i+1].x +
                                                    l_in_arr[k][j+1][i].x + l_in_arr[k][j+1][i+1].x +
                                                    l_in_arr[k+1][j][i].x + l_in_arr[k+1][j][i+1].x +
                                                    l_in_arr[k+1][j+1][i].x + l_in_arr[k+1][j+1][i+1].x);
 
                    g_out_arr[k][j][i].y = 0.125 * (l_in_arr[k][j][i].y + l_in_arr[k][j][i+1].y +
                                                    l_in_arr[k][j+1][i].y + l_in_arr[k][j+1][i+1].y +
                                                    l_in_arr[k+1][j][i].y + l_in_arr[k+1][j][i+1].y +
                                                    l_in_arr[k+1][j+1][i].y + l_in_arr[k+1][j+1][i+1].y);
 
                    g_out_arr[k][j][i].z = 0.125 * (l_in_arr[k][j][i].z + l_in_arr[k][j][i+1].z +
                                                    l_in_arr[k][j+1][i].z + l_in_arr[k][j+1][i+1].z +
                                                    l_in_arr[k+1][j][i].z + l_in_arr[k+1][j][i+1].z +
                                                    l_in_arr[k+1][j+1][i].z + l_in_arr[k+1][j+1][i+1].z);
                }
            }
        }
        ierr = DMDAVecRestoreArrayRead(dm_in,in_vec_local, (void*)&l_in_arr); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(dm_out,out_vec_global, (void*)&g_out_arr); CHKERRQ(ierr);
    }
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeQCriterion()

PetscErrorCode ComputeQCriterion

(

UserCtx *

user

)

Computes the Q-criterion diagnostic from the local velocity-gradient tensor.

This kernel evaluates rotational versus strain-rate dominance and writes the result into the configured Q-criterion output vector for visualization and flow feature identification.

Parameters

[in,out]

user

Block-level context containing velocity fields and target output storage.

Returns

PetscErrorCode 0 on success.

Computes the Q-criterion diagnostic from the local velocity-gradient tensor.

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

See also

ComputeQCriterion()

Definition at line 366 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    const Cmpnts   ***lucat, ***lcsi, ***leta, ***lzet;
    const PetscReal***laj, ***lnvert;
    PetscReal      ***gq;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeQCriterion.\n");
 
    // --- 1. Ensure all required ghost values are up-to-date ---
    ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT);  CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_CSI);   CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_ETA);   CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_ZET);   CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_AJ);    CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, FIELD_ID_NVERT); CHKERRQ(ierr);
 
    // --- 2. Get DMDA info and array pointers ---
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    
    ierr = DMDAVecGetArrayRead(user->fda, user->lUcat,  (void*)&lucat);   CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi,   (void*)&lcsi);    CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta,   (void*)&leta);    CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet,   (void*)&lzet);    CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lAj,    (void*)&laj);     CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da,  user->lNvert, (void*)&lnvert);  CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da,  user->Qcrit, (void*)&gq);       CHKERRQ(ierr);
 
    // --- 3. Define Loop Bounds for INTERIOR Cells ---
    PetscInt i_start = (info.xs == 0) ? 1 : info.xs;
    PetscInt i_end   = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
    PetscInt j_start = (info.ys == 0) ? 1 : info.ys;
    PetscInt j_end   = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
    PetscInt k_start = (info.zs == 0) ? 1 : info.zs;
    PetscInt k_end   = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
 
    // --- 4. Main Computation Loop ---
    for (PetscInt k = k_start; k < k_end; k++) {
        for (PetscInt j = j_start; j < j_end; j++) {
            for (PetscInt i = i_start; i < i_end; i++) {
                
                // Calculate velocity derivatives in computational space (central differences)
                PetscReal uc = 0.5 * (lucat[k][j][i+1].x - lucat[k][j][i-1].x);
                PetscReal vc = 0.5 * (lucat[k][j][i+1].y - lucat[k][j][i-1].y);
                PetscReal wc = 0.5 * (lucat[k][j][i+1].z - lucat[k][j][i-1].z);
 
                PetscReal ue = 0.5 * (lucat[k][j+1][i].x - lucat[k][j-1][i].x);
                PetscReal ve = 0.5 * (lucat[k][j+1][i].y - lucat[k][j-1][i].y);
                PetscReal we = 0.5 * (lucat[k][j+1][i].z - lucat[k][j-1][i].z);
 
                PetscReal uz = 0.5 * (lucat[k+1][j][i].x - lucat[k-1][j][i].x);
                PetscReal vz = 0.5 * (lucat[k+1][j][i].y - lucat[k-1][j][i].y);
                PetscReal wz = 0.5 * (lucat[k+1][j][i].z - lucat[k-1][j][i].z);
 
                // Average metrics to the cell center
                PetscReal csi1 = 0.5 * (lcsi[k][j][i].x + lcsi[k][j][i-1].x) * laj[k][j][i];
                PetscReal csi2 = 0.5 * (lcsi[k][j][i].y + lcsi[k][j][i-1].y) * laj[k][j][i];
                PetscReal csi3 = 0.5 * (lcsi[k][j][i].z + lcsi[k][j][i-1].z) * laj[k][j][i];
                
                PetscReal eta1 = 0.5 * (leta[k][j][i].x + leta[k][j-1][i].x) * laj[k][j][i];
                PetscReal eta2 = 0.5 * (leta[k][j][i].y + leta[k][j-1][i].y) * laj[k][j][i];
                PetscReal eta3 = 0.5 * (leta[k][j][i].z + leta[k][j-1][i].z) * laj[k][j][i];
 
                PetscReal zet1 = 0.5 * (lzet[k][j][i].x + lzet[k-1][j][i].x) * laj[k][j][i];
                PetscReal zet2 = 0.5 * (lzet[k][j][i].y + lzet[k-1][j][i].y) * laj[k][j][i];
                PetscReal zet3 = 0.5 * (lzet[k][j][i].z + lzet[k-1][j][i].z) * laj[k][j][i];
 
                // Calculate velocity gradient tensor components d_ij = du_i/dx_j
                PetscReal d11 = uc * csi1 + ue * eta1 + uz * zet1;
                PetscReal d12 = uc * csi2 + ue * eta2 + uz * zet2;
                PetscReal d13 = uc * csi3 + ue * eta3 + uz * zet3;
 
                PetscReal d21 = vc * csi1 + ve * eta1 + vz * zet1;
                PetscReal d22 = vc * csi2 + ve * eta2 + vz * zet2;
                PetscReal d23 = vc * csi3 + ve * eta3 + vz * zet3;
                
                PetscReal d31 = wc * csi1 + we * eta1 + wz * zet1;
                PetscReal d32 = wc * csi2 + we * eta2 + wz * zet2;
                PetscReal d33 = wc * csi3 + we * eta3 + wz * zet3;
 
                // Strain-Rate Tensor S_ij = 0.5 * (d_ij + d_ji)
                PetscReal s11 = d11;
                PetscReal s12 = 0.5 * (d12 + d21);
                PetscReal s13 = 0.5 * (d13 + d31);
                PetscReal s22 = d22;
                PetscReal s23 = 0.5 * (d23 + d32);
                PetscReal s33 = d33;
 
                // Vorticity Tensor Omega_ij = 0.5 * (d_ij - d_ji)
                PetscReal w12 = 0.5 * (d12 - d21);
                PetscReal w13 = 0.5 * (d13 - d31);
                PetscReal w23 = 0.5 * (d23 - d32);
 
                // Squared norms of the tensors
                PetscReal s_norm_sq = s11*s11 + s22*s22 + s33*s33 + 2.0*(s12*s12 + s13*s13 + s23*s23);
                PetscReal w_norm_sq = 2.0 * (w12*w12 + w13*w13 + w23*w23);
                
                gq[k][j][i] = 0.5 * (w_norm_sq - s_norm_sq);
 
                if (lnvert[k][j][i] > 0.1) {
                    gq[k][j][i] = 0.0;
                }
            }
        }
    }
 
    // --- 5. Restore arrays ---
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lUcat,  (void*)&lucat);   CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi,   (void*)&lcsi);    CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta,   (void*)&leta);    CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet,   (void*)&lzet);    CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lAj,    (void*)&laj);     CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da,  user->lNvert, (void*)&lnvert);  CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da,  user->Qcrit, (void*)&gq);       CHKERRQ(ierr);
 
    /* With dimensionalize on, Ucat already carries U_ref but the metrics are still
     * nondimensional, so the gradients above are U_ref/L_ref short of physical by one
     * factor of L_ref each: Q came out U_ref^2 Q* instead of (U_ref/L_ref)^2 Q*. */
    if (user->simCtx->pps && user->simCtx->pps->dimensionalize) {
        PetscReal length_scale = 1.0;
        if (!PicurvFieldReferenceScale(user->simCtx, "Coordinates", &length_scale, NULL, 0) &&
            length_scale > 0.0) {
            ierr = VecScale(user->Qcrit, 1.0 / (length_scale * length_scale)); CHKERRQ(ierr);
        }
    }
 
    /* The loop above skips the layout boundary, so fill it before anything reads
     * Qcrit with a stencil. Q is not a moment: it does not vanish at a wall, so only
     * the periodic case is defined, and the cell synchronizer writes only that case. */
    {
        const FieldId qcrit[] = {FIELD_ID_QCRIT};
        ierr = SynchronizePeriodicCellFields(user, 1, qcrit); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

NormalizeRelativeField()

PetscErrorCode NormalizeRelativeField	(	UserCtx *	user,
		const char *	relative_field_name
	)

Normalizes pressure using the value at the configured logical grid point.

The owning rank reads the reference value, shares it collectively, and every rank subtracts it in-place from its distributed portion of the field.

Parameters

[in,out]	user	Block-level context containing pressure and reference configuration.
[in]	relative_field_name	Name of the field to normalize.

Returns

PetscErrorCode 0 on success.

Normalizes pressure using the value at the configured logical grid point.

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

See also

NormalizeRelativeField()

Definition at line 515 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    Vec            P_vec = NULL;
    DMDALocalInfo  info;
    PetscInt       ip=1, jp=1, kp=1; // Default reference point
    PetscReal      p_ref = 0.0;
    PetscReal      p_ref_local = 0.0;
    PetscInt       found_local = 0, found_global = 0;
    PostProcessParams *pps = user->simCtx->pps;
 
    // Fetch the logical reference point from pps.
    ip = pps->reference[0];
    jp = pps->reference[1];
    kp = pps->reference[2];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running NormalizeRelativeField on '%s'.\n", relative_field_name);
 
    // --- 1. Map string argument to the PETSc Vec ---
    if (strcasecmp(relative_field_name, "P") == 0) {
        P_vec = user->P;
    } else {
        SETERRQ(PETSC_COMM_SELF, 1, "NormalizeRelativeField only supports the primary 'P' field , not '%s' currently.", relative_field_name);
    }
 
    // --- 2. Read the logical reference point from whichever rank owns it ---
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    PetscCheck(ip >= 0 && ip < info.mx && jp >= 0 && jp < info.my &&
               kp >= 0 && kp < info.mz,
               PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE,
               "Reference point (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT
               ") lies outside the %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT
               " pressure layout.", ip, jp, kp, info.mx, info.my, info.mz);
    if (ip >= info.xs && ip < info.xs + info.xm &&
        jp >= info.ys && jp < info.ys + info.ym &&
        kp >= info.zs && kp < info.zs + info.zm) {
        const PetscReal ***pressure = NULL;
 
        ierr = DMDAVecGetArrayRead(user->da, P_vec, &pressure); CHKERRQ(ierr);
        p_ref_local = pressure[kp][jp][ip];
        ierr = DMDAVecRestoreArrayRead(user->da, P_vec, &pressure); CHKERRQ(ierr);
        found_local = 1;
    }
    ierr = MPI_Allreduce(&p_ref_local, &p_ref, 1, MPIU_REAL, MPI_SUM,
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&found_local, &found_global, 1, MPIU_INT, MPI_SUM,
                         PETSC_COMM_WORLD); CHKERRQ(ierr);
    PetscCheck(found_global == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB,
               "Reference pressure point must have exactly one owner; found %" PetscInt_FMT ".",
               found_global);
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "%s reference point (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT
              ") has value %g.\n", relative_field_name, ip, jp, kp, (double)p_ref);
 
    // --- 3. Perform the normalization (in-place shift) on the full distributed vector ---
    ierr = VecShift(P_vec, -p_ref); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "%s field normalized by subtracting %g.\n", relative_field_name, p_ref);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

DimensionalizeField()

PetscErrorCode DimensionalizeField	(	UserCtx *	user,
		const char *	field_name
	)

Scales a specified field from non-dimensional to dimensional units in-place.

Resolves the name in the Eulerian field catalog, then the particle field catalog, and scales the field's storage (its global Vec, the DM coordinates, or its DMSwarm field) by the reference scale of the dimension recorded on its catalog entry.

Parameters

[in,out]	user	The UserCtx whose storage is modified.
[in]	field_name	Case-insensitive catalogued name or alias, e.g. "Ucat", "Nu_t", "Coordinates", "ParticlePosition".

Returns

Zero on success; PETSC_ERR_ARG_UNKNOWN_TYPE for a name in neither catalog, and PETSC_ERR_ARG_WRONGSTATE for a field without a fixed dimension or without storage in this run.

Scales a specified field from non-dimensional to dimensional units in-place.

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

See also

DimensionalizeField()

Definition at line 15 of file postprocessing_kernels.c.

{
    PetscErrorCode  ierr;
    SimCtx          *simCtx = NULL;
    Vec             target_vec = NULL;
    PetscReal       scale_factor = 1.0;
    char            field_type[64] = "Unknown";
    PetscBool       is_swarm_field = PETSC_FALSE; // Flag for special swarm handling
    const char      *swarm_field_name = NULL;    // Name of the field within the swarm
    FieldId         field_id = FIELD_ID_INVALID;
    ParticleFieldId particle_id = PARTICLE_FIELD_ID_INVALID;
    PetscBool       found = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx is NULL.");
    if (!field_name) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "field_name is NULL.");
    simCtx = user->simCtx;
 
    // --- 1. Resolve the field in either catalog; its entry records its dimension ---
    ierr = FieldTryIdFromName(field_name, &field_id, &found); CHKERRQ(ierr);
    if (found && field_id == FIELD_ID_COORDINATES) {
        ierr = DMGetCoordinates(user->da, &target_vec); CHKERRQ(ierr);
    } else if (found) {
        FieldView view;
 
        ierr = FieldGetView(user, field_id, &view); CHKERRQ(ierr);
        target_vec = view.global_vec;
    } else {
        const ParticleFieldDescriptor *descriptor = NULL;
 
        ierr = ParticleFieldTryIdFromName(field_name, &particle_id, &found); CHKERRQ(ierr);
        PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE,
                   "DimensionalizeField: '%s' is in neither field catalog.", field_name);
        ierr = ParticleFieldGetDescriptor(particle_id, &descriptor); CHKERRQ(ierr);
        is_swarm_field = PETSC_TRUE;
        swarm_field_name = descriptor->canonical_name;
    }
    ierr = PicurvFieldReferenceScale(simCtx, field_name, &scale_factor,
                                     field_type, sizeof(field_type)); CHKERRQ(ierr);
 
    // --- 2. Check for trivial scaling ---
    if (PetscAbsReal(scale_factor - 1.0) < PETSC_MACHINE_EPSILON) {
        LOG(GLOBAL, LOG_DEBUG, "DimensionalizeField: Scaling factor for '%s' is 1.0. Skipping operation.\n", field_name);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    // --- 3. Perform the in-place scaling operation ---
    LOG(GLOBAL, LOG_INFO, "Scaling '%s' field (%s) by factor %.4e.\n", field_name, field_type, scale_factor);
 
    if (is_swarm_field) {
        // Special handling for DMSwarm fields
        ierr = DMSwarmCreateGlobalVectorFromField(user->swarm, swarm_field_name, &target_vec); CHKERRQ(ierr);
        ierr = VecScale(target_vec, scale_factor); CHKERRQ(ierr);
        ierr = DMSwarmDestroyGlobalVectorFromField(user->swarm, swarm_field_name, &target_vec); CHKERRQ(ierr);
    } else {
        // Standard handling for PETSc Vecs
        if (target_vec) {
            ierr = VecScale(target_vec, scale_factor); CHKERRQ(ierr);
        } else {
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "Target vector for field '%s' was not found or is NULL.", field_name);
        }
    }
    
    // --- 4. Post-scaling updates for special cases ---
    if (strcasecmp(field_name, "Coordinates") == 0) {
        ierr = UpdateLocalGhosts(user, FIELD_ID_COORDINATES); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeSpecificKE()

PetscErrorCode ComputeSpecificKE	(	UserCtx *	user,
		const char *	velocity_field,
		const char *	ske_field
	)

Computes the specific kinetic energy (KE per unit mass) for each particle.

This kernel calculates SKE = 0.5 * |velocity|^2. It requires that the velocity field exists and will populate the specific kinetic energy field. The output field must be registered before this kernel is called.

Parameters

user	The UserCtx containing the DMSwarm.
velocity_field	The name of the input vector field for particle velocity.
ske_field	The name of the output scalar field to store specific KE.

Returns

PetscErrorCode

Computes the specific kinetic energy (KE per unit mass) for each particle.

Local to this translation unit.

Definition at line 588 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    PetscInt n_local;
    const PetscScalar (*vel_arr)[3]; // Access velocity as array of 3-component vectors
    PetscScalar *ske_arr;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeSpecificKE ('%s' -> '%s').\n", velocity_field, ske_field);
 
    // Get local data arrays from the DMSwarm
    ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
    if (n_local == 0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    // Get read-only access to velocity and write access to the output field
    ierr = DMSwarmGetField(user->swarm, velocity_field, NULL, NULL, (void**)&vel_arr); CHKERRQ(ierr);
    ierr = DMSwarmGetField(user->post_swarm, ske_field, NULL, NULL, (void**)&ske_arr); CHKERRQ(ierr);
 
    // Main computation loop
    for (PetscInt p = 0; p < n_local; p++) {
        const PetscScalar u = vel_arr[p][0];
        const PetscScalar v = vel_arr[p][1];
        const PetscScalar w = vel_arr[p][2];
        const PetscScalar vel_sq = u*u + v*v + w*w;
        ske_arr[p] = 0.5 * vel_sq;
    }
 
    // Restore arrays
    ierr = DMSwarmRestoreField(user->swarm, velocity_field, NULL, NULL, (void**)&vel_arr); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(user->post_swarm, ske_field, NULL, NULL, (void**)&ske_arr); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ComputeDisplacement()

PetscErrorCode ComputeDisplacement	(	UserCtx *	user,
		const char *	disp_field
	)

Computes the displacement magnitude |r_i - r_0| for each particle (per-particle VTK kernel).

Reference point r_0 = (simCtx->psrc_x, psrc_y, psrc_z). Writes the scalar displacement to post_swarm[disp_field]. This is a visualisation kernel only — use ComputeParticleMSD from particle_statistics.h for quantitative global statistics.

Parameters

user	The UserCtx containing the DMSwarms.
disp_field	Name of the output scalar field in post_swarm.

Returns

PetscErrorCode

Computes the displacement magnitude |r_i - r_0| for each particle (per-particle VTK kernel).

Local to this translation unit.

Definition at line 630 of file postprocessing_kernels.c.

{
    PetscErrorCode  ierr;
    PetscInt        n_local;
    const PetscReal (*pos_arr)[3];
    PetscScalar    *disp_out;
    SimCtx         *simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    LOG_ALLOW(GLOBAL, LOG_INFO, "-> KERNEL: Running ComputeDisplacement (-> '%s').\n", disp_field);
 
    ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
    if (n_local == 0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    const PetscReal x0 = simCtx->psrc_x;
    const PetscReal y0 = simCtx->psrc_y;
    const PetscReal z0 = simCtx->psrc_z;
 
    ierr = DMSwarmGetField(user->swarm,      ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
    ierr = DMSwarmGetField(user->post_swarm, disp_field, NULL, NULL, (void**)&disp_out);       CHKERRQ(ierr);
 
    for (PetscInt p = 0; p < n_local; p++) {
        const PetscReal dx = pos_arr[p][0] - x0;
        const PetscReal dy = pos_arr[p][1] - y0;
        const PetscReal dz = pos_arr[p][2] - z0;
        disp_out[p] = PetscSqrtReal(dx*dx + dy*dy + dz*dz);
    }
 
    ierr = DMSwarmRestoreField(user->swarm,      ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(user->post_swarm, disp_field, NULL, NULL, (void**)&disp_out);       CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeWindowStatisticNodal()

PetscErrorCode ComputeWindowStatisticNodal	(	UserCtx *	user,
		PetscInt	window_index,
		const char *	outputs,
		PetscInt	output_index,
		char *	out_name,
		size_t	name_size,
		Vec *	out_vec,
		PetscInt *	out_components
	)

Derives one accumulated statistic and converts it to nodal values.

The counterpart of ComputeNodalAverage() for accumulated window state: it takes an index into a window's requested output set, normalizes the centered state behind it, and leaves the result in the shared post-processing staging field ready for output. The staging vectors are reused between calls, so a caller must consume the result before requesting the next one.

With global_operations.dimensionalize set, the result carries physical units: the source field's reference scale raised to the power the derived kind carries: a mean and an RMS are linear in it, a Reynolds stress and a turbulent kinetic energy quadratic. A co-moment flux relates two possibly different fields, so its factor is the product of their scales rather than one squared. That is what the per-field scaling table alone cannot express, and why one blanket velocity factor would have been wrong for three of the five kinds.

Parameters

[in]	user	Block context holding the accumulators and staging fields.
[in]	window_index	Window whose state is derived.
[in]	outputs	Comma-separated output kinds the recipe requested.
[in]	output_index	Index into that output set.
[out]	out_name	Name of the derived field, window qualified.
[in]	name_size	Capacity of out_name.
[out]	out_vec	Nodal vector holding the result; borrowed, not owned.
[out]	out_components	Component count of the result.

Returns

Zero on success, or a PETSc error.

Derives one accumulated statistic and converts it to nodal values.

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

See also

ComputeWindowStatisticNodal()

Definition at line 212 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = NULL;
    PicurvDerivedField derived;
    const PicurvWindowDefinition *definition = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && out_name != NULL && out_vec != NULL && out_components != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Context and outputs are required.");
    simCtx = user->simCtx;
    PetscCheck(FieldStatisticsIsActive(simCtx) && user->fieldStatisticsStorage != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "No accumulated window state exists to derive.");
    definition = &simCtx->fieldStatisticsWindows[window_index].definition;
 
    /* Derive into the cell-centred staging pair, then reach the nodal path by
     * catalogued name. A per-window accumulator has no compile-time offset, so
     * staging is what lets the shared ghost and nodal kernels address it at all. */
    ierr = PicurvWindowDerive(user, definition, &user->fieldStatisticsStorage[window_index],
                              outputs, output_index, user->PostScalar, user->PostVector,
                              &derived); CHKERRQ(ierr);
    /* The derivation writes the physical interior only, so the layout boundary would
     * otherwise reach the nodal average as structural zeros. On periodic axes the
     * shared cell synchronizer fills the dummy planes from the wrapped physical
     * planes; the ghost update then carries the written values outward. */
    if (derived.components == 1) {
        const FieldId staged[] = {FIELD_ID_POST_SCALAR};
        ierr = SynchronizePeriodicCellFields(user, 1, staged); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_POST_SCALAR); CHKERRQ(ierr);
        ierr = ComputeNodalAverage(user, "PostScalar", "PostScalarNodal"); CHKERRQ(ierr);
        *out_vec = user->PostScalarNodal;
    } else {
        const FieldId staged[] = {FIELD_ID_POST_VECTOR};
        ierr = SynchronizePeriodicCellFields(user, 1, staged); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_POST_VECTOR); CHKERRQ(ierr);
        ierr = ComputeNodalAverage(user, "PostVector", "PostVectorNodal"); CHKERRQ(ierr);
        *out_vec = user->PostVectorNodal;
    }
    *out_components = derived.components;
    ierr = PetscStrncpy(out_name, derived.name, name_size); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "-> KERNEL: Derived '%s' (%d component(s)).\n",
              out_name, (int)derived.components);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeWindowStatisticsSummary()

PetscErrorCode ComputeWindowStatisticsSummary	(	UserCtx *	user,
		PetscInt	window_index,
		const char *	output_prefix,
		PetscInt	ti
	)

Appends one convergence row for an accumulated window to its CSV history.

The counterpart of ComputeParticleMSD() for Eulerian window state: it reduces the window to the few numbers that answer whether it has run long enough — sample count, total weight, represented time, the per-point valid-fraction range, and the mean turbulent kinetic energy — and appends them as one row per processed step. No single field snapshot can answer that question, which is why the history exists beside the field output rather than instead of it.

The mean is taken over the fluid cells the window actually sampled. Cells outside the target domain, and cells a moving mask excluded, hold zeros that mean "never measured"; averaging over them would scale the result down by the fraction of the grid the window never covered.

Parameters

[in]	user	Block context holding the accumulators and staging fields.
[in]	window_index	Window to summarize.
[in]	output_prefix	Output path prefix; the window name and .csv are appended.
[in]	ti	Step being processed, recorded as the row's key.

Returns

Zero on success, or a PETSc error.

Appends one convergence row for an accumulated window to its CSV history.

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

See also

ComputeWindowStatisticsSummary()

Definition at line 272 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = NULL;
    const PicurvWindow *window = NULL;
    const PicurvWindowStorage *storage = NULL;
    PetscReal lowest = 1.0, highest = 0.0;
    PetscReal mean_tke = 0.0;
    PetscBool has_tke = PETSC_FALSE;
    PetscInt derived_count = 0;
    char path[PETSC_MAX_PATH_LEN];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && output_prefix != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context and output prefix are required.");
    simCtx = user->simCtx;
    PetscCheck(FieldStatisticsIsActive(simCtx) && user->fieldStatisticsStorage != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "No accumulated window state exists to summarize.");
    window = &simCtx->fieldStatisticsWindows[window_index];
    storage = &user->fieldStatisticsStorage[window_index];
 
    ierr = PicurvWindowValidFractionRange(user, &window->definition, storage,
                                          window->sample_count, &lowest, &highest); CHKERRQ(ierr);
 
    /* A single domain number is what makes the row plottable against step. It is
     * produced by the same derivation the field output uses, so the scalar and the
     * field can never disagree about what the window holds. */
    ierr = PicurvWindowDerivedCount(&window->definition, storage, "tke", &derived_count); CHKERRQ(ierr);
    if (derived_count > 0) {
        PicurvDerivedField derived;
 
        ierr = PicurvWindowDerive(user, &window->definition, storage, "tke", 0,
                                  user->PostScalar, user->PostVector, &derived); CHKERRQ(ierr);
        /* Averaged over the fluid cells the window actually sampled. A whole-vector
         * mean would divide by boundary and dummy entries the derivation never
         * writes, scaling the answer down by the fraction it never covered. */
        ierr = PicurvWindowSpatialMean(user, &window->definition, storage,
                                       user->PostScalar, &mean_tke); CHKERRQ(ierr);
        has_tke = PETSC_TRUE;
    }
 
    /* The window's clock is solver time. A dimensionalized summary reports its derived
     * statistics physically, so the times beside them are reported in seconds too; a
     * time-weighted total weight is a time, a sample-weighted one is a count. */
    PetscReal time_scale = 1.0;
    if (simCtx->pps && simCtx->pps->dimensionalize) {
        const FieldDimension time_dimension = FIELD_DIM_TIME;
 
        ierr = FieldDimensionReferenceScale(&simCtx->scaling, time_dimension, &time_scale); CHKERRQ(ierr);
    }
    const PetscReal weight_scale =
        (window->definition.weighting == PICURV_WEIGHTING_PHYSICAL_TIME) ? time_scale : 1.0;
 
    if (simCtx->rank == 0) {
        FILE *csv = NULL;
        PetscBool exists = PETSC_FALSE;
 
        ierr = PetscSNPrintf(path, sizeof(path), "%s_statistics_%s.csv",
                             output_prefix, window->definition.name); CHKERRQ(ierr);
        ierr = PetscTestFile(path, 'r', &exists); CHKERRQ(ierr);
        csv = fopen(path, exists ? "a" : "w");
        PetscCheck(csv != NULL, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
                   "Unable to open statistics summary '%s'.", path);
        if (!exists) {
            fprintf(csv, "step,state,samples,total_weight,represented_time,"
                         "valid_fraction_min,valid_fraction_max,mean_tke\n");
        }
        fprintf(csv, "%" PetscInt_FMT ",%s,%d,%.10e,%.10e,%.6f,%.6f,",
                ti, PicurvWindowStateName(window->state), window->sample_count,
                (double)(window->total_weight * weight_scale),
                (double)(window->represented_time * time_scale),
                (double)lowest, (double)highest);
        if (has_tke) fprintf(csv, "%.10e\n", (double)mean_tke);
        else fprintf(csv, "\n");
        PetscCheck(fclose(csv) == 0, PETSC_COMM_SELF, PETSC_ERR_FILE_WRITE,
                   "Unable to close statistics summary '%s'.", path);
    }
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "-> KERNEL: Summarized window '%s' at step %" PetscInt_FMT ".\n",
              window->definition.name, ti);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: