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 <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 124 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 36 of file particle_statistics.c.

37{
38 char *endptr = NULL;
39 long parsed_step;
40
41 if (!line) return PETSC_FALSE;
42 parsed_step = strtol(line, &endptr, 10);
43 if (endptr == line) return PETSC_FALSE;
44 while (*endptr == ' ' || *endptr == '\t') endptr++;
45 if (*endptr != ',') return PETSC_FALSE;
46 return (PetscBool)(parsed_step == (long)ti ? PETSC_TRUE : PETSC_FALSE);
47}
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 56 of file particle_statistics.c.

57{
58 char tmp_path[PETSC_MAX_PATH_LEN];
59 FILE *src = NULL;
60 FILE *dst = NULL;
61 char line[4096];
62 PetscBool wrote_header = PETSC_FALSE;
63 PetscBool inserted_row = PETSC_FALSE;
64
65 PetscFunctionBeginUser;
66 PetscCall(PetscSNPrintf(tmp_path, sizeof(tmp_path), "%s.tmp", csv_path));
67
68 dst = fopen(tmp_path, "w");
69 if (!dst) {
70 LOG_ALLOW(GLOBAL, LOG_WARNING,
71 "ComputeParticleMSD: could not open '%s' for writing.\n", tmp_path);
72 PetscFunctionReturn(0);
73 }
74
75 src = fopen(csv_path, "r");
76 if (src) {
77 while (fgets(line, sizeof(line), src)) {
78 if (!wrote_header) {
79 if (strncmp(line, "step,", 5) == 0) {
80 fputs(line, dst);
81 wrote_header = PETSC_TRUE;
82 continue;
83 }
84 fputs(PARTICLE_MSD_CSV_HEADER, dst);
85 wrote_header = PETSC_TRUE;
86 }
87
88 if (MSDCSVLineMatchesStep(line, ti)) {
89 if (!inserted_row) {
90 fputs(row_line, dst);
91 inserted_row = PETSC_TRUE;
92 }
93 continue;
94 }
95 fputs(line, dst);
96 }
97 fclose(src);
98 src = NULL;
99 }
100
101 if (!wrote_header) {
102 fputs(PARTICLE_MSD_CSV_HEADER, dst);
103 }
104 if (!inserted_row) {
105 fputs(row_line, dst);
106 }
107
108 if (fclose(dst) != 0) {
109 remove(tmp_path);
110 LOG_ALLOW(GLOBAL, LOG_WARNING,
111 "ComputeParticleMSD: could not finalize '%s'.\n", tmp_path);
112 PetscFunctionReturn(0);
113 }
114
115 if (rename(tmp_path, csv_path) != 0) {
116 remove(tmp_path);
117 LOG_ALLOW(GLOBAL, LOG_WARNING,
118 "ComputeParticleMSD: could not replace '%s'.\n", csv_path);
119 }
120 PetscFunctionReturn(0);
121}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
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:199
LOG_WARNING
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:29
PARTICLE_MSD_CSV_HEADER
static const char * PARTICLE_MSD_CSV_HEADER
Definition particle_statistics.c:24
MSDCSVLineMatchesStep
static PetscBool MSDCSVLineMatchesStep(const char *line, PetscInt ti)
Return whether a CSV row belongs to the requested timestep.
Definition particle_statistics.c:36
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 129 of file particle_statistics.c.

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

Generated by doxygen 1.9.8