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 a relative scalar field using the configured reference pressure scale.
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).

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 143 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    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; }
    // ... (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; }
    // ... (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 = UpdateLocalGhosts(user, in_field_name); CHKERRQ(ierr);
 
    // --- 3. Get DMDA info and array pointers ---
    DMDALocalInfo info;
    ierr = DMDAGetLocalInfo(dm_out, &info); CHKERRQ(ierr);
 
    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 < info.zs + info.zm - 1; k++) {
            for (PetscInt j = info.ys; j < info.ys + info.ym - 1; j++) {
                for (PetscInt i = info.xs; i < info.xs + info.xm - 1; 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 < info.zs + info.zm - 1; k++) {
            for (PetscInt j = info.ys; j < info.ys + info.ym - 1; j++) {
                for (PetscInt i = info.xs; i < info.xs + info.xm - 1; 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 237 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, "Ucat");  CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Csi");   CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Eta");   CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Zet");   CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Aj");    CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "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);
 
    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 a relative scalar field using the configured reference pressure scale.

This is primarily used for pressure-normalized outputs in post-processing. The operation is in-place on the selected field.

Parameters

[in,out]	user	Block-level context containing scaling information.
[in]	relative_field_name	Name of the field to normalize.

Returns

PetscErrorCode 0 on success.

Normalizes a relative scalar field using the configured reference pressure scale.

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

See also

NormalizeRelativeField()

Definition at line 367 of file postprocessing_kernels.c.

{
    PetscErrorCode ierr;
    Vec            P_vec = NULL;
    PetscMPIInt    rank;
    PetscInt       ip=1, jp=1, kp=1; // Default reference point
    PetscReal      p_ref = 0.0;
    PetscInt       ref_point_global_idx[1];
    PetscScalar    ref_value_local[1];
    IS             is_from, is_to;
    VecScatter     scatter_ctx;
    Vec            ref_vec_seq;
    PostProcessParams *pps = user->simCtx->pps;
 
    // Fetch referenc 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. Get reference point from options and calculate its global DA index ---
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // Convert the (i,j,k) logical grid coordinates to the global 1D index used by the DMDA vector
    ref_point_global_idx[0] = kp * (user->IM * user->JM) + jp * user->IM + ip;
 
    // --- 3. Robustly Scatter the single reference value to rank 0 ---
    // Create an Index Set (IS) for the source point (from the global vector)
    ierr = ISCreateGeneral(PETSC_COMM_WORLD, 1, ref_point_global_idx, PETSC_COPY_VALUES, &is_from); CHKERRQ(ierr);
 
    // Create a sequential vector on rank 0 to hold the result
    ierr = VecCreateSeq(PETSC_COMM_SELF, 1, &ref_vec_seq); CHKERRQ(ierr);
    
    // Create an Index Set for the destination point (index 0 of the new sequential vector)
    PetscInt dest_idx[1] = {0};
    ierr = ISCreateGeneral(PETSC_COMM_SELF, 1, dest_idx, PETSC_COPY_VALUES, &is_to); CHKERRQ(ierr);
 
    // Create the scatter context and perform the scatter
    ierr = VecScatterCreate(P_vec, is_from, ref_vec_seq, is_to, &scatter_ctx); CHKERRQ(ierr);
    ierr = VecScatterBegin(scatter_ctx, P_vec, ref_vec_seq, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecScatterEnd(scatter_ctx, P_vec, ref_vec_seq, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
 
    // On rank 0, get the value. On other ranks, this will do nothing.
    if (rank == 0) {
        ierr = VecGetValues(ref_vec_seq, 1, dest_idx, ref_value_local); CHKERRQ(ierr);
        p_ref = ref_value_local[0];
        LOG_ALLOW(LOCAL, LOG_DEBUG, "%s reference point (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ") has value %g.\n", relative_field_name, jp, kp, ip, p_ref);
    }
    
    // --- 4. Broadcast the reference value from rank 0 to all other processes ---
    ierr = MPI_Bcast(&p_ref, 1, MPIU_REAL, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    // --- 5. 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);
    
    // --- 6. Cleanup ---
    ierr = ISDestroy(&is_from); CHKERRQ(ierr);
    ierr = ISDestroy(&is_to); CHKERRQ(ierr);
    ierr = VecScatterDestroy(&scatter_ctx); CHKERRQ(ierr);
    ierr = VecDestroy(&ref_vec_seq); CHKERRQ(ierr);
 
    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 12 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) {
        Vec l_coords;
        ierr = DMGetCoordinates(user->da, &target_vec); CHKERRQ(ierr); // Re-get the global vector
        ierr = DMGetCoordinatesLocal(user->da, &l_coords); CHKERRQ(ierr);
        ierr = DMGlobalToLocalBegin(user->fda, target_vec, INSERT_VALUES, l_coords); CHKERRQ(ierr);
        ierr = DMGlobalToLocalEnd(user->fda, target_vec, INSERT_VALUES, l_coords); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

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 103 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:

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 452 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) 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 494 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) PetscFunctionReturn(0);
 
    const PetscReal x0 = simCtx->psrc_x;
    const PetscReal y0 = simCtx->psrc_y;
    const PetscReal z0 = simCtx->psrc_z;
 
    ierr = DMSwarmGetField(user->swarm,      "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,      "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 caller graph for this function: