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	ExtendToLayoutBoundary (UserCtx *user, Vec global, PetscInt components)
	Populates the layout boundary of a field that was written on the interior only.
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	DimensionalizeAllLoadedFields (UserCtx *user)
	Orchestrates the dimensionalization of all relevant fields loaded from a file.
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 236 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; }
    /* 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, "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);
 
    // --- 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:

ExtendToLayoutBoundary()

PetscErrorCode ExtendToLayoutBoundary	(	UserCtx *	user,
		Vec	global,
		PetscInt	components
	)

Populates the layout boundary of a field that was written on the interior only.

This is not a boundary condition. It enforces no physics and satisfies no equation. It replicates the layout convention that boundary-condition application already gives solver state fields, so that stencil kernels — ComputeNodalAverage() above all — read defined values instead of the structural zeros an interior-only producer leaves behind.

Fields written by PicurvWindowDerive() and ComputeQCriterion() cover only the physical interior, so their layout boundary holds zeros that mean "never written" rather than "zero". A node on the domain edge averages four such entries with four real cells and comes out at roughly half its true value. Solver state fields never need this: Ucat, P, Nu_t, and CS all carry boundary values written when their boundary conditions were applied.

Only periodic directions are handled. There the correct value exists and is exact — the low dummy plane repeats the last physical plane and the high dummy plane repeats the first. On a non-periodic face the correct value depends on both the quantity and the boundary type, and no single convention serves a staging buffer that carries stresses, pressure, and eddy viscosity in turn; nothing is written there rather than something invented. The design for the non-periodic case is recorded in Field Statistics Planned Extensions.

Operates on the global vector, so it must run before the caller's UpdateLocalGhosts(), which then carries the written values into every local ghost region. Rank interfaces are that scatter's responsibility and are never touched here.

Multi-block interface boundaries are out of scope.

Parameters

[in]	user	Block context supplying the DMDA layout and periodicity.
[in,out]	global	Global vector whose layout boundary is populated.
[in]	components	Degrees of freedom carried: 1 or 3.

Returns

Zero on success, or PETSC_ERR_ARG_NULL for a null argument, or PETSC_ERR_ARG_OUTOFRANGE for an unsupported component count.

Populates the layout boundary of a field that was written on the interior only.

Full API contract is documented with the header declaration in include/postprocessing_kernels.h.

See also

ExtendToLayoutBoundary()

Definition at line 142 of file postprocessing_kernels.c.

{
    SimCtx        *simCtx = NULL;
    DM             dm = NULL;
    Vec            local = NULL;
    PetscInt       periodic[3];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && global != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context and field are required.");
    PetscCheck(components == 1 || components == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "ExtendToLayoutBoundary handles 1- or 3-component fields; got %d.", (int)components);
    simCtx = user->simCtx;
    dm = (components == 1) ? user->da : user->fda;
 
    periodic[0] = simCtx->i_periodic;
    periodic[1] = simCtx->j_periodic;
    periodic[2] = simCtx->k_periodic;
    if (!periodic[0] && !periodic[1] && !periodic[2]) {
        /* Nothing is defined on a non-periodic layout boundary, so nothing is written.
         * The boundary node keeps the value the producing kernel left there, and the
         * spatial reductions that matter scientifically never read it. */
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    PetscCall(DMGetLocalVector(dm, &local));
    /* One pass per periodic direction, each preceded by its own scatter. The order
     * matters: after the i pass the i layout boundary carries data, so the j pass
     * reading across it fills edges correctly, and the k pass then fills corners. */
    for (PetscInt dir = 0; dir < 3; ++dir) {
        DMDALocalInfo info;
        const PetscReal ****source = NULL;
        PetscReal      ****target = NULL;
        PetscInt        extent = 0;
 
        if (!periodic[dir]) continue;
 
        PetscCall(DMGlobalToLocalBegin(dm, global, INSERT_VALUES, local));
        PetscCall(DMGlobalToLocalEnd(dm, global, INSERT_VALUES, local));
        PetscCall(DMDAGetLocalInfo(dm, &info));
        extent = (dir == 0) ? info.mx : ((dir == 1) ? info.my : info.mz);
 
        PetscCall(DMDAVecGetArrayDOFRead(dm, local, &source));
        PetscCall(DMDAVecGetArrayDOF(dm, global, &target));
        for (PetscInt k = info.zs; k < info.zs + info.zm; ++k) {
            for (PetscInt j = info.ys; j < info.ys + info.ym; ++j) {
                for (PetscInt i = info.xs; i < info.xs + info.xm; ++i) {
                    const PetscInt index = (dir == 0) ? i : ((dir == 1) ? j : k);
                    PetscInt       ss = i, sj = j, sk = k;
 
                    /* The layout wraps one plane inside the DMDA extent: the low dummy
                     * plane repeats the last physical plane, and the high dummy plane
                     * repeats the first. PETSc's periodic ghosts supply both sources
                     * even when they live on another rank. */
                    if (index == 0) {
                        if (dir == 0) ss = extent - 2;
                        else if (dir == 1) sj = extent - 2;
                        else sk = extent - 2;
                    } else if (index == extent - 1) {
                        if (dir == 0) ss = extent + 1;
                        else if (dir == 1) sj = extent + 1;
                        else sk = extent + 1;
                    } else {
                        continue;
                    }
                    for (PetscInt c = 0; c < components; ++c) {
                        target[k][j][i][c] = source[sk][sj][ss][c];
                    }
                }
            }
        }
        PetscCall(DMDAVecRestoreArrayDOF(dm, global, &target));
        PetscCall(DMDAVecRestoreArrayDOFRead(dm, local, &source));
    }
    PetscCall(DMRestoreLocalVector(dm, &local));
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "-> KERNEL: Extended %d-component field across the periodic layout boundary "
              "(i=%d, j=%d, k=%d).\n", (int)components,
              (int)periodic[0], (int)periodic[1], (int)periodic[2]);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

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 481 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);
 
    /* The loop above skips the layout boundary, so extend 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 only that case is written. */
    ierr = ExtendToLayoutBoundary(user, user->Qcrit, 1); 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 616 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.

This function acts as a dispatcher. It takes the string name of a field, identifies the corresponding PETSc Vec object and the correct physical scaling factor (e.g., U_ref for velocity, P_ref for pressure), and then performs an in-place VecScale operation. It correctly handles the different physical dimensions of Cartesian velocity vs. contravariant volume flux.

Parameters

[in,out]	user	The UserCtx containing the PETSc Vecs to be modified.
[in]	field_name	The case-insensitive string name of the field to dimensionalize (e.g., "Ucat", "P", "Ucont", "Coordinates", "ParticlePosition", "ParticleVelocity").

Returns

PetscErrorCode

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 = user->simCtx;
    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
 
    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.");
 
    // --- 1. Identify the target Vec and the correct scaling factor ---
    if (strcasecmp(field_name, "Ucat") == 0) {
        target_vec = user->Ucat;
        scale_factor = simCtx->scaling.U_ref;
        strcpy(field_type, "Cartesian Velocity (L/T)");
    } else if (strcasecmp(field_name, "Ucont") == 0) {
        target_vec = user->Ucont;
        scale_factor = simCtx->scaling.U_ref * simCtx->scaling.L_ref * simCtx->scaling.L_ref;
        strcpy(field_type, "Contravariant Volume Flux (L^3/T)");
    } else if (strcasecmp(field_name, "P") == 0) {
        target_vec = user->P;
        scale_factor = simCtx->scaling.P_ref;
        strcpy(field_type, "Pressure (M L^-1 T^-2)");
    } else if (strcasecmp(field_name, "Coordinates") == 0) {
        ierr = DMGetCoordinates(user->da, &target_vec); CHKERRQ(ierr);
        scale_factor = simCtx->scaling.L_ref;
        strcpy(field_type, "Grid Coordinates (L)");
    } else if (strcasecmp(field_name, "ParticlePosition") == 0) {
        is_swarm_field = PETSC_TRUE;
        swarm_field_name = "position";
        scale_factor = simCtx->scaling.L_ref;
        strcpy(field_type, "Particle Position (L)");
    } else if (strcasecmp(field_name, "ParticleVelocity") == 0) {
        is_swarm_field = PETSC_TRUE;
        swarm_field_name = "velocity";
        scale_factor = simCtx->scaling.U_ref;
        strcpy(field_type, "Particle Velocity (L/T)");
    } else {
        LOG(GLOBAL, LOG_WARNING, "DimensionalizeField: Unknown or unhandled field_name '%s'. Field will not be scaled.\n", field_name);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
    
    // --- 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:

DimensionalizeAllLoadedFields()

PetscErrorCode DimensionalizeAllLoadedFields

(

UserCtx *

user

)

Orchestrates the dimensionalization of all relevant fields loaded from a file.

This function is intended to be called in the post-processor immediately after all solver output has been read into memory. It calls DimensionalizeField() for each of the core physical quantities to convert the entire loaded state from non-dimensional to dimensional units, preparing it for analysis and visualization.

Parameters

[in,out]

user

The UserCtx containing all the fields to be dimensionalized.

Returns

PetscErrorCode

Orchestrates the dimensionalization of all relevant fields loaded from a file.

Local to this translation unit.

Definition at line 102 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG(GLOBAL, LOG_INFO, "--- Converting all loaded fields to dimensional units ---\n");
 
    // Scale the grid itself first
    ierr = DimensionalizeField(user, "Coordinates"); CHKERRQ(ierr);
 
    // Scale primary fluid fields
    ierr = DimensionalizeField(user, "Ucat"); CHKERRQ(ierr);
    ierr = DimensionalizeField(user, "Ucont"); CHKERRQ(ierr);
    ierr = DimensionalizeField(user, "P"); CHKERRQ(ierr);
 
    // If particles are present, scale their fields
    if (simCtx->np > 0 && user->swarm) {
        ierr = DimensionalizeField(user, "ParticlePosition"); CHKERRQ(ierr);
        ierr = DimensionalizeField(user, "ParticleVelocity"); CHKERRQ(ierr);
    }
 
    LOG(GLOBAL, LOG_INFO, "--- Field dimensionalization complete ---\n");
 
    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 689 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 731 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.

Derived statistics are left non-dimensional even when global_operations.dimensionalize is set. A Reynolds stress scales as velocity squared and a co-moment as a product of two different scales, none of which the existing per-field scaling table expresses; silently applying a velocity scale would be wrong rather than merely incomplete.

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 343 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. Extend it before the
     * ghost update, which then carries the written values outward. */
    if (derived.components == 1) {
        ierr = ExtendToLayoutBoundary(user, user->PostScalar, 1); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_POST_SCALAR); CHKERRQ(ierr);
        ierr = ComputeNodalAverage(user, "PostScalar", "PostScalarNodal"); CHKERRQ(ierr);
        *out_vec = user->PostScalarNodal;
    } else {
        ierr = ExtendToLayoutBoundary(user, user->PostVector, 3); 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 400 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;
    }
 
    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, (double)window->represented_time,
                (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: