particle_statistics.h File Reference

Global statistics kernels for the Statistics Pipeline. More...

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

Include dependency graph for particle_statistics.h:

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

Go to the source code of this file.

Functions
PetscErrorCode	ComputeParticleMSD (UserCtx *user, const char *stats_prefix, PetscInt ti)
	Computes the mean-squared displacement (MSD) of a particle cloud.

Detailed Description

Global statistics kernels for the Statistics Pipeline.

These kernels compute global aggregate quantities via MPI reduction — they do NOT produce per-particle VTK data. Each kernel appends one CSV row per call and logs a one-line summary via LOG_INFO.

Output file convention: {stats_prefix}_{kernel_name}.csv All MPI_Allreduce operations are internal; file I/O is rank-0-only.

To add a new statistic:

1. Implement PetscErrorCode ComputeXxx(UserCtx*, const char*, PetscInt) in src/particle_statistics.c
2. Declare it below.
3. Add one else-if in GlobalStatisticsPipeline() in src/postprocessor.c No other files need to change.

Definition in file particle_statistics.h.

Function Documentation

ComputeParticleMSD()

PetscErrorCode ComputeParticleMSD	(	UserCtx *	user,
		const char *	stats_prefix,
		PetscInt	ti
	)

Computes the mean-squared displacement (MSD) of a particle cloud.

Reference point r0 = (simCtx->psrc_x, psrc_y, psrc_z). Diffusivity D = 1 / (Re * Sc), time t = ti * dt. Computes MSD_x, MSD_y, MSD_z (isotropy check), MSD_total, r_rms, centre-of-mass drift, and fractions inside 1σ/2σ/3σ theoretical shells. Appends one row to {stats_prefix}_msd.csv; writes header on first call.

Parameters

user	The UserCtx containing the DMSwarm (user->swarm).
stats_prefix	Base filename prefix (e.g. "brownian_stats").
ti	Current time-step index.

Returns

PetscErrorCode

Computes the mean-squared displacement (MSD) of a particle cloud.

Local to this translation unit.

Definition at line 130 of file particle_statistics.c.

{
    PetscErrorCode ierr;
    SimCtx        *simCtx = user->simCtx;
    PetscMPIInt    rank;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    /* ------------------------------------------------------------------ *
     * Physics parameters                                                   *
     * ------------------------------------------------------------------ */
    const PetscReal Re  = (simCtx->ren > 0.0) ? simCtx->ren : 1.0;
    const PetscReal Sc  = (simCtx->schmidt_number > 0.0) ? simCtx->schmidt_number : 1.0;
    const PetscReal D   = 1.0 / (Re * Sc);
    const PetscReal t   = (PetscReal)ti * simCtx->dt;
    const PetscReal x0  = simCtx->psrc_x;
    const PetscReal y0  = simCtx->psrc_y;
    const PetscReal z0  = simCtx->psrc_z;
 
    /* ------------------------------------------------------------------ *
     * Pass 1: local accumulation                                           *
     * ------------------------------------------------------------------ */
    PetscInt       n_local;
    const PetscReal (*pos_arr)[3];
 
    ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
    ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
 
    PetscReal local_sq_x = 0.0, local_sq_y = 0.0, local_sq_z = 0.0;
    PetscReal local_sx   = 0.0, local_sy   = 0.0, local_sz   = 0.0;
 
    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;
        local_sq_x += dx * dx;
        local_sq_y += dy * dy;
        local_sq_z += dz * dz;
        local_sx   += dx;
        local_sy   += dy;
        local_sz   += dz;
    }
 
    ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
 
    /* ------------------------------------------------------------------ *
     * MPI reduction — 7 doubles + count                                   *
     * ------------------------------------------------------------------ */
    PetscReal local_buf[7] = { local_sq_x, local_sq_y, local_sq_z,
                                local_sx,   local_sy,   local_sz,
                                (PetscReal)n_local };
    PetscReal global_buf[7] = {0};
    ierr = MPI_Allreduce(local_buf, global_buf, 7, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    const PetscReal N_total  = global_buf[6];
    if (N_total < 1.0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); } /* no particles */
 
    const PetscReal MSD_x    = global_buf[0] / N_total;
    const PetscReal MSD_y    = global_buf[1] / N_total;
    const PetscReal MSD_z    = global_buf[2] / N_total;
    const PetscReal com_x    = global_buf[3] / N_total;
    const PetscReal com_y    = global_buf[4] / N_total;
    const PetscReal com_z    = global_buf[5] / N_total;
 
    const PetscReal MSD_total    = MSD_x + MSD_y + MSD_z;
    const PetscReal r_rms_meas   = PetscSqrtReal(MSD_total);
    const PetscReal r_theory     = (t > 1e-300) ? PetscSqrtReal(6.0 * D * t) : 0.0;
    const PetscReal rel_err_pct  = (r_theory > 1e-12)
                                   ? PetscAbsReal(r_rms_meas - r_theory) / r_theory * 100.0
                                   : 0.0;
 
    /* ------------------------------------------------------------------ *
     * Pass 2: fraction inside σ-shells                                    *
     * ------------------------------------------------------------------ */
    ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
 
    PetscInt local_n1 = 0, local_n2 = 0, local_n3 = 0;
    const PetscReal r1 = r_theory;
    const PetscReal r2 = 2.0 * r_theory;
    const PetscReal r3 = 3.0 * r_theory;
 
    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;
        const PetscReal r  = PetscSqrtReal(dx*dx + dy*dy + dz*dz);
        if (r < r1) local_n1++;
        if (r < r2) local_n2++;
        if (r < r3) local_n3++;
    }
 
    ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
 
    PetscReal local_counts[3] = { (PetscReal)local_n1,
                                   (PetscReal)local_n2,
                                   (PetscReal)local_n3 };
    PetscReal global_counts[3] = {0};
    ierr = MPI_Allreduce(local_counts, global_counts, 3, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
 
    const PetscReal frac_1s = global_counts[0] / N_total * 100.0;
    const PetscReal frac_2s = global_counts[1] / N_total * 100.0;
    const PetscReal frac_3s = global_counts[2] / N_total * 100.0;
 
    /* ------------------------------------------------------------------ *
     * Output: CSV + LOG_INFO (rank 0 only)                                *
     * ------------------------------------------------------------------ */
    if (rank == 0) {
        char csv_path[PETSC_MAX_PATH_LEN];
        char row_line[1024];
        PetscSNPrintf(csv_path, sizeof(csv_path), "%s_msd.csv", stats_prefix);
        PetscSNPrintf(
            row_line,
            sizeof(row_line),
            "%d,%.6e,%.0f,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.4f,%.6e,%.6e,%.6e,%.2f,%.2f,%.2f\n",
            (int)ti, t, N_total,
            MSD_x, MSD_y, MSD_z, MSD_total,
            r_rms_meas, r_theory, rel_err_pct,
            com_x, com_y, com_z,
            frac_1s, frac_2s, frac_3s
        );
        PetscCall(RewriteParticleMSDCSV(csv_path, ti, row_line));
 
        LOG_ALLOW(GLOBAL, LOG_INFO,
                  "[MSD ti=%d t=%.4f] total=%.4e theory=%.4e err=%.2f%% | "
                  "fracs: %.1f/%.1f/%.1f%% | COM: (%.2e,%.2e,%.2e)\n",
                  (int)ti, t, MSD_total, r_theory * r_theory, rel_err_pct,
                  frac_1s, frac_2s, frac_3s,
                  com_x, com_y, com_z);
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: