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.
Computes the mean-squared displacement (MSD) of a particle cloud.
Two-pass algorithm: Pass 1 — accumulate per-component squared displacements and centre-of-mass sums. MPI reduction — allreduce 7 values (6 sums + count) to get global MSD and COM. Compute scalars — MSD_total, r_rms_meas, r_rms_theory, rel_err_pct. Pass 2 — count particles inside 1σ/2σ/3σ theoretical shells; allreduce counts. Output (rank 0) — append one row to {stats_prefix}_msd.csv; LOG_INFO summary.
Physics: D = 1/(Re·Sc), t = ti·dt, r_theory = sqrt(6·D·t). Expected fractions inside r_theory / 2·r_theory / 3·r_theory ≈ 61%/97%/99.9% for a 3D Maxwell-Boltzmann distribution.
39{
40 PetscErrorCode ierr;
42 PetscMPIInt rank;
43
44 PetscFunctionBeginUser;
46
47 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
48
49
50
51
52 const PetscReal Re = (simCtx->
ren > 0.0) ? simCtx->
ren : 1.0;
54 const PetscReal D = 1.0 / (Re * Sc);
55 const PetscReal t = (PetscReal)ti * simCtx->
dt;
56 const PetscReal x0 = simCtx->psrc_x;
57 const PetscReal y0 = simCtx->psrc_y;
58 const PetscReal z0 = simCtx->psrc_z;
59
60
61
62
63 PetscInt n_local;
64 const PetscReal (*pos_arr)[3];
65
66 ierr = DMSwarmGetLocalSize(user->
swarm, &n_local); CHKERRQ(ierr);
67 ierr = DMSwarmGetField(user->
swarm,
"position", NULL, NULL, (
const void**)&pos_arr); CHKERRQ(ierr);
68
69 PetscReal local_sq_x = 0.0, local_sq_y = 0.0, local_sq_z = 0.0;
70 PetscReal local_sx = 0.0, local_sy = 0.0, local_sz = 0.0;
71
72 for (PetscInt p = 0; p < n_local; p++) {
73 const PetscReal dx = pos_arr[p][0] - x0;
74 const PetscReal dy = pos_arr[p][1] - y0;
75 const PetscReal dz = pos_arr[p][2] - z0;
76 local_sq_x += dx * dx;
77 local_sq_y += dy * dy;
78 local_sq_z += dz * dz;
79 local_sx += dx;
80 local_sy += dy;
81 local_sz += dz;
82 }
83
84 ierr = DMSwarmRestoreField(user->
swarm,
"position", NULL, NULL, (
const void**)&pos_arr); CHKERRQ(ierr);
85
86
87
88
89 PetscReal local_buf[7] = { local_sq_x, local_sq_y, local_sq_z,
90 local_sx, local_sy, local_sz,
91 (PetscReal)n_local };
92 PetscReal global_buf[7] = {0};
93 MPI_Allreduce(local_buf, global_buf, 7, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
94
95 const PetscReal N_total = global_buf[6];
96 if (N_total < 1.0) PetscFunctionReturn(0);
97
98 const PetscReal MSD_x = global_buf[0] / N_total;
99 const PetscReal MSD_y = global_buf[1] / N_total;
100 const PetscReal MSD_z = global_buf[2] / N_total;
101 const PetscReal com_x = global_buf[3] / N_total;
102 const PetscReal com_y = global_buf[4] / N_total;
103 const PetscReal com_z = global_buf[5] / N_total;
104
105 const PetscReal MSD_total = MSD_x + MSD_y + MSD_z;
106 const PetscReal r_rms_meas = PetscSqrtReal(MSD_total);
107 const PetscReal r_theory = (t > 1e-300) ? PetscSqrtReal(6.0 * D * t) : 0.0;
108 const PetscReal rel_err_pct = (r_theory > 1e-12)
109 ? PetscAbsReal(r_rms_meas - r_theory) / r_theory * 100.0
110 : 0.0;
111
112
113
114
115 ierr = DMSwarmGetField(user->
swarm,
"position", NULL, NULL, (
const void**)&pos_arr); CHKERRQ(ierr);
116
117 PetscInt local_n1 = 0, local_n2 = 0, local_n3 = 0;
118 const PetscReal r1 = r_theory;
119 const PetscReal r2 = 2.0 * r_theory;
120 const PetscReal r3 = 3.0 * r_theory;
121
122 for (PetscInt p = 0; p < n_local; p++) {
123 const PetscReal dx = pos_arr[p][0] - x0;
124 const PetscReal dy = pos_arr[p][1] - y0;
125 const PetscReal dz = pos_arr[p][2] - z0;
126 const PetscReal r = PetscSqrtReal(dx*dx + dy*dy + dz*dz);
127 if (r < r1) local_n1++;
128 if (r < r2) local_n2++;
129 if (r < r3) local_n3++;
130 }
131
132 ierr = DMSwarmRestoreField(user->
swarm,
"position", NULL, NULL, (
const void**)&pos_arr); CHKERRQ(ierr);
133
134 PetscReal local_counts[3] = { (PetscReal)local_n1,
135 (PetscReal)local_n2,
136 (PetscReal)local_n3 };
137 PetscReal global_counts[3] = {0};
138 MPI_Allreduce(local_counts, global_counts, 3, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
139
140 const PetscReal frac_1s = global_counts[0] / N_total * 100.0;
141 const PetscReal frac_2s = global_counts[1] / N_total * 100.0;
142 const PetscReal frac_3s = global_counts[2] / N_total * 100.0;
143
144
145
146
147 if (rank == 0) {
148 char csv_path[PETSC_MAX_PATH_LEN];
149 PetscSNPrintf(csv_path, sizeof(csv_path), "%s_msd.csv", stats_prefix);
150
151 FILE *f = fopen(csv_path, "a");
152 if (!f) {
154 "ComputeParticleMSD: could not open '%s' for writing.\n", csv_path);
155 } else {
156
157 if (ftell(f) == 0) {
158 fprintf(f, "step,t,N,MSD_x,MSD_y,MSD_z,MSD_total,"
159 "r_rms_meas,r_rms_theory,rel_err_pct,"
160 "com_x,com_y,com_z,"
161 "frac_1sigma_pct,frac_2sigma_pct,frac_3sigma_pct\n");
162 }
163 fprintf(f, "%d,%.6e,%.0f,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.4f,%.6e,%.6e,%.6e,%.2f,%.2f,%.2f\n",
164 (int)ti, t, N_total,
165 MSD_x, MSD_y, MSD_z, MSD_total,
166 r_rms_meas, r_theory, rel_err_pct,
167 com_x, com_y, com_z,
168 frac_1s, frac_2s, frac_3s);
169 fclose(f);
170 }
171
173 "[MSD ti=%d t=%.4f] total=%.4e theory=%.4e err=%.2f%% | "
174 "fracs: %.1f/%.1f/%.1f%% | COM: (%.2e,%.2e,%.2e)\n",
175 (int)ti, t, MSD_total, r_theory * r_theory, rel_err_pct,
176 frac_1s, frac_2s, frac_3s,
177 com_x, com_y, com_z);
178 }
179
181 PetscFunctionReturn(0);
182}
#define GLOBAL
Scope for global logging across all processes.
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
@ LOG_INFO
Informational messages about program execution.
@ LOG_WARNING
Non-critical issues that warrant attention.
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
SimCtx * simCtx
Back-pointer to the master simulation context.
The master context for the entire simulation.
1.9.8