PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions | Variables
particle_statistics.c File Reference

Global statistical reduction kernels for the Statistics Pipeline. More...

#include "particle_statistics.h"
#include "particle_field_catalog.h"
#include <math.h>
#include <mpi.h>
#include <stdlib.h>
#include <string.h>
Include dependency graph for particle_statistics.c:

Go to the source code of this file.

Macros

#define __FUNCT__   "ComputeParticleMSD"
 

Functions

static PetscBool MSDCSVLineMatchesStep (const char *line, PetscInt ti)
 Return whether a CSV row belongs to the requested timestep.
 
static PetscErrorCode RewriteParticleMSDCSV (const char *csv_path, PetscInt ti, const char *row_line)
 Rewrite the MSD CSV so the requested timestep appears exactly once.
 
PetscErrorCode ComputeParticleMSD (UserCtx *user, const char *stats_prefix, PetscInt ti)
 Internal helper implementation: ComputeParticleMSD().
 

Variables

static const char * PARTICLE_MSD_CSV_HEADER
 

Detailed Description

Global statistical reduction kernels for the Statistics Pipeline.

Each kernel in this file:

To add a new statistic:

  1. Implement PetscErrorCode ComputeXxx(UserCtx*, const char*, PetscInt) here.
  2. Declare it in include/particle_statistics.h.
  3. Add one else-if in GlobalStatisticsPipeline() in src/postprocessor.c.

Definition in file particle_statistics.c.

Macro Definition Documentation

◆ __FUNCT__

#define __FUNCT__   "ComputeParticleMSD"

Definition at line 125 of file particle_statistics.c.

Function Documentation

◆ MSDCSVLineMatchesStep()

static PetscBool MSDCSVLineMatchesStep ( const char *  line,
PetscInt  ti 
)
static

Return whether a CSV row belongs to the requested timestep.

Parameters
[in]lineCSV row text to inspect.
[in]tiRequested timestep id.
Returns
PETSC_TRUE when the row begins with the requested step id.

Definition at line 37 of file particle_statistics.c.

38{
39 char *endptr = NULL;
40 long parsed_step;
41
42 if (!line) return PETSC_FALSE;
43 parsed_step = strtol(line, &endptr, 10);
44 if (endptr == line) return PETSC_FALSE;
45 while (*endptr == ' ' || *endptr == '\t') endptr++;
46 if (*endptr != ',') return PETSC_FALSE;
47 return (PetscBool)(parsed_step == (long)ti ? PETSC_TRUE : PETSC_FALSE);
48}
Here is the caller graph for this function:

◆ RewriteParticleMSDCSV()

static PetscErrorCode RewriteParticleMSDCSV ( const char *  csv_path,
PetscInt  ti,
const char *  row_line 
)
static

Rewrite the MSD CSV so the requested timestep appears exactly once.

Parameters
[in]csv_pathTarget CSV path.
[in]tiTimestep being written.
[in]row_lineFully formatted CSV row for the timestep.
Returns
Petsc error code.

Definition at line 57 of file particle_statistics.c.

58{
59 char tmp_path[PETSC_MAX_PATH_LEN];
60 FILE *src = NULL;
61 FILE *dst = NULL;
62 char line[4096];
63 PetscBool wrote_header = PETSC_FALSE;
64 PetscBool inserted_row = PETSC_FALSE;
65
66 PetscFunctionBeginUser;
67 PetscCall(PetscSNPrintf(tmp_path, sizeof(tmp_path), "%s.tmp", csv_path));
68
69 dst = fopen(tmp_path, "w");
70 if (!dst) {
71 LOG_ALLOW(GLOBAL, LOG_WARNING,
72 "ComputeParticleMSD: could not open '%s' for writing.\n", tmp_path);
73 PetscFunctionReturn(0);
74 }
75
76 src = fopen(csv_path, "r");
77 if (src) {
78 while (fgets(line, sizeof(line), src)) {
79 if (!wrote_header) {
80 if (strncmp(line, "step,", 5) == 0) {
81 fputs(line, dst);
82 wrote_header = PETSC_TRUE;
83 continue;
84 }
85 fputs(PARTICLE_MSD_CSV_HEADER, dst);
86 wrote_header = PETSC_TRUE;
87 }
88
89 if (MSDCSVLineMatchesStep(line, ti)) {
90 if (!inserted_row) {
91 fputs(row_line, dst);
92 inserted_row = PETSC_TRUE;
93 }
94 continue;
95 }
96 fputs(line, dst);
97 }
98 fclose(src);
99 src = NULL;
100 }
101
102 if (!wrote_header) {
103 fputs(PARTICLE_MSD_CSV_HEADER, dst);
104 }
105 if (!inserted_row) {
106 fputs(row_line, dst);
107 }
108
109 if (fclose(dst) != 0) {
110 remove(tmp_path);
111 LOG_ALLOW(GLOBAL, LOG_WARNING,
112 "ComputeParticleMSD: could not finalize '%s'.\n", tmp_path);
113 PetscFunctionReturn(0);
114 }
115
116 if (rename(tmp_path, csv_path) != 0) {
117 remove(tmp_path);
118 LOG_ALLOW(GLOBAL, LOG_WARNING,
119 "ComputeParticleMSD: could not replace '%s'.\n", csv_path);
120 }
121 PetscFunctionReturn(0);
122}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
PARTICLE_MSD_CSV_HEADER
static const char * PARTICLE_MSD_CSV_HEADER
Definition particle_statistics.c:25
MSDCSVLineMatchesStep
static PetscBool MSDCSVLineMatchesStep(const char *line, PetscInt ti)
Return whether a CSV row belongs to the requested timestep.
Definition particle_statistics.c:37
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeParticleMSD()

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

Internal helper implementation: ComputeParticleMSD().

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

Local to this translation unit.

Definition at line 130 of file particle_statistics.c.

131{
132 PetscErrorCode ierr;
133 SimCtx *simCtx = user->simCtx;
134 PetscMPIInt rank;
135
136 PetscFunctionBeginUser;
137 PROFILE_FUNCTION_BEGIN;
138
139 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
140
141 /* ------------------------------------------------------------------ *
142 * Physics parameters *
143 * ------------------------------------------------------------------ */
144 const PetscReal Re = (simCtx->ren > 0.0) ? simCtx->ren : 1.0;
145 const PetscReal Sc = (simCtx->schmidt_number > 0.0) ? simCtx->schmidt_number : 1.0;
146 const PetscReal D = 1.0 / (Re * Sc);
147 const PetscReal t = (PetscReal)ti * simCtx->dt;
148 const PetscReal x0 = simCtx->psrc_x;
149 const PetscReal y0 = simCtx->psrc_y;
150 const PetscReal z0 = simCtx->psrc_z;
151
152 /* ------------------------------------------------------------------ *
153 * Pass 1: local accumulation *
154 * ------------------------------------------------------------------ */
155 PetscInt n_local;
156 const PetscReal (*pos_arr)[3];
157
158 ierr = DMSwarmGetLocalSize(user->swarm, &n_local); CHKERRQ(ierr);
159 ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
160
161 PetscReal local_sq_x = 0.0, local_sq_y = 0.0, local_sq_z = 0.0;
162 PetscReal local_sx = 0.0, local_sy = 0.0, local_sz = 0.0;
163
164 for (PetscInt p = 0; p < n_local; p++) {
165 const PetscReal dx = pos_arr[p][0] - x0;
166 const PetscReal dy = pos_arr[p][1] - y0;
167 const PetscReal dz = pos_arr[p][2] - z0;
168 local_sq_x += dx * dx;
169 local_sq_y += dy * dy;
170 local_sq_z += dz * dz;
171 local_sx += dx;
172 local_sy += dy;
173 local_sz += dz;
174 }
175
176 ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
177
178 /* ------------------------------------------------------------------ *
179 * MPI reduction — 7 doubles + count *
180 * ------------------------------------------------------------------ */
181 PetscReal local_buf[7] = { local_sq_x, local_sq_y, local_sq_z,
182 local_sx, local_sy, local_sz,
183 (PetscReal)n_local };
184 PetscReal global_buf[7] = {0};
185 ierr = MPI_Allreduce(local_buf, global_buf, 7, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
186
187 const PetscReal N_total = global_buf[6];
188 if (N_total < 1.0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); } /* no particles */
189
190 const PetscReal MSD_x = global_buf[0] / N_total;
191 const PetscReal MSD_y = global_buf[1] / N_total;
192 const PetscReal MSD_z = global_buf[2] / N_total;
193 const PetscReal com_x = global_buf[3] / N_total;
194 const PetscReal com_y = global_buf[4] / N_total;
195 const PetscReal com_z = global_buf[5] / N_total;
196
197 const PetscReal MSD_total = MSD_x + MSD_y + MSD_z;
198 const PetscReal r_rms_meas = PetscSqrtReal(MSD_total);
199 const PetscReal r_theory = (t > 1e-300) ? PetscSqrtReal(6.0 * D * t) : 0.0;
200 const PetscReal rel_err_pct = (r_theory > 1e-12)
201 ? PetscAbsReal(r_rms_meas - r_theory) / r_theory * 100.0
202 : 0.0;
203
204 /* ------------------------------------------------------------------ *
205 * Pass 2: fraction inside σ-shells *
206 * ------------------------------------------------------------------ */
207 ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
208
209 PetscInt local_n1 = 0, local_n2 = 0, local_n3 = 0;
210 const PetscReal r1 = r_theory;
211 const PetscReal r2 = 2.0 * r_theory;
212 const PetscReal r3 = 3.0 * r_theory;
213
214 for (PetscInt p = 0; p < n_local; p++) {
215 const PetscReal dx = pos_arr[p][0] - x0;
216 const PetscReal dy = pos_arr[p][1] - y0;
217 const PetscReal dz = pos_arr[p][2] - z0;
218 const PetscReal r = PetscSqrtReal(dx*dx + dy*dy + dz*dz);
219 if (r < r1) local_n1++;
220 if (r < r2) local_n2++;
221 if (r < r3) local_n3++;
222 }
223
224 ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void**)&pos_arr); CHKERRQ(ierr);
225
226 PetscReal local_counts[3] = { (PetscReal)local_n1,
227 (PetscReal)local_n2,
228 (PetscReal)local_n3 };
229 PetscReal global_counts[3] = {0};
230 ierr = MPI_Allreduce(local_counts, global_counts, 3, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD); CHKERRMPI(ierr);
231
232 const PetscReal frac_1s = global_counts[0] / N_total * 100.0;
233 const PetscReal frac_2s = global_counts[1] / N_total * 100.0;
234 const PetscReal frac_3s = global_counts[2] / N_total * 100.0;
235
236 /* ------------------------------------------------------------------ *
237 * Output: CSV + LOG_INFO (rank 0 only) *
238 * ------------------------------------------------------------------ */
239 if (rank == 0) {
240 char csv_path[PETSC_MAX_PATH_LEN];
241 char row_line[1024];
242 PetscSNPrintf(csv_path, sizeof(csv_path), "%s_msd.csv", stats_prefix);
243 PetscSNPrintf(
244 row_line,
245 sizeof(row_line),
246 "%d,%.6e,%.0f,%.6e,%.6e,%.6e,%.6e,%.6e,%.6e,%.4f,%.6e,%.6e,%.6e,%.2f,%.2f,%.2f\n",
247 (int)ti, t, N_total,
248 MSD_x, MSD_y, MSD_z, MSD_total,
249 r_rms_meas, r_theory, rel_err_pct,
250 com_x, com_y, com_z,
251 frac_1s, frac_2s, frac_3s
252 );
253 PetscCall(RewriteParticleMSDCSV(csv_path, ti, row_line));
254
255 LOG_ALLOW(GLOBAL, LOG_INFO,
256 "[MSD ti=%d t=%.4f] total=%.4e theory=%.4e err=%.2f%% | "
257 "fracs: %.1f/%.1f/%.1f%% | COM: (%.2e,%.2e,%.2e)\n",
258 (int)ti, t, MSD_total, r_theory * r_theory, rel_err_pct,
259 frac_1s, frac_2s, frac_3s,
260 com_x, com_y, com_z);
261 }
262
263 PROFILE_FUNCTION_END;
264 PetscFunctionReturn(0);
265}
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:859
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:31
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:850
ParticleFieldName
const char * ParticleFieldName(ParticleFieldId field_id)
Return the canonical PETSc DMSwarm name for an ID.
Definition particle_field_catalog.c:111
PARTICLE_FIELD_ID_POSITION
@ PARTICLE_FIELD_ID_POSITION
Definition particle_field_catalog.h:25
RewriteParticleMSDCSV
static PetscErrorCode RewriteParticleMSDCSV(const char *csv_path, PetscInt ti, const char *row_line)
Rewrite the MSD CSV so the requested timestep appears exactly once.
Definition particle_statistics.c:57
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:787
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
SimCtx::ren
PetscReal ren
Definition variables.h:744
SimCtx::dt
PetscReal dt
Definition variables.h:710
SimCtx
The master context for the entire simulation.
Definition variables.h:695
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ PARTICLE_MSD_CSV_HEADER

const char* PARTICLE_MSD_CSV_HEADER
static
Initial value:
=
"step,t,N,MSD_x,MSD_y,MSD_z,MSD_total,"
"r_rms_meas,r_rms_theory,rel_err_pct,"
"com_x,com_y,com_z,"
"frac_1sigma_pct,frac_2sigma_pct,frac_3sigma_pct\n"

Definition at line 25 of file particle_statistics.c.

Generated by doxygen 1.9.8