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)
	Computes node-centered data by averaging 8 surrounding cell-centered values, exactly replicating the legacy code's indexing and boundary handling.
PetscErrorCode	ComputeQCriterion (UserCtx *user)
	Computes the Q-Criterion, a scalar value identifying vortex cores.
PetscErrorCode	NormalizeRelativeField (UserCtx *user, const char *relative_field_name)
	Normalizes a relative field by subtracting a reference value.
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.

Function Documentation

ComputeNodalAverage()

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

Computes node-centered data by averaging 8 surrounding cell-centered values, exactly replicating the legacy code's indexing and boundary handling.

This kernel uses a stencil that averages the 8 cells from the corner (i,j,k) to (i+1, j+1, k+1) and stores the result at the output node (i,j,k). This matches the legacy code's behavior. It operates on the full range of output nodes necessary for the subsampled grid, preventing zero-padding at the boundaries.

Parameters

user	The UserCtx, providing access to DMs and Vecs.
in_field_name	The string name of the source Vec (e.g., "P", "Ucat").
out_field_name	The string name of the destination Vec (e.g., "P_nodal").

Returns

PetscErrorCode

Definition at line 166 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, a scalar value identifying vortex cores.

This function is self-contained. It ensures all its required input data (Ucat and grid metrics) have up-to-date ghost values before proceeding with the calculation. The result is stored in the global vector user->Qcrit.

Parameters

user	The UserCtx, providing access to all necessary data.

Returns

PetscErrorCode

Definition at line 264 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 field by subtracting a reference value.

This kernel finds the relative field at a specific grid point (i,j,k) and subtracts this value from the entire field. The reference point is configurable via command-line options (-ip, -jp, -kp). The operation is performed in-place on the provided relative field vector.

Parameters

user	The UserCtx, providing access to DMs and Vecs.
relative_field_name	The string name of the relative field Vec to normalize (e.g., "P").

Returns

PetscErrorCode

Definition at line 400 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

Definition at line 20 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

Definition at line 118 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

Definition at line 493 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, (const 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, (const 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: