|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
|
Loading...
Searching...
No Matches
|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
|
Field Statistics describes what the statistics pipeline does today: pointwise weighted centered first and second moments, accumulated per named window, carried across a restart, and derived into Reynolds stresses, RMS, turbulent kinetic energy, and fluxes.
This page records what it deliberately does not do, why each extension was left out, and what it depends on. Each is opt-in and builds on the abstractions that already exist, so none of them requires reworking what is there. Configuration written against the current contract keeps working through all of them.
One rule decides where each extension belongs, and it is worth stating before the list because it removes more work than it adds.
Averaging commutes with linear operations and not with anything else.
A spatial reduction over accumulated pointwise state is linear, so it is exact after the fact. Reducing per-point accumulators over a region reproduces the region-wide centered moment about the region mean, including the spatial variance of the local means, through the same weighted combination PicurvMomentStateMerge already implements:
\[ M_{2,AB} = M_{2,A} + M_{2,B} + \frac{W_A W_B}{W_A + W_B}\,(\mu_B - \mu_A)^2 . \]
That identity extends to higher orders, so whatever order is accumulated pointwise can be reduced to any spatial subset at that order or below without approximation. Bins, profiles, regions, and clipping are therefore post-processing operations, not accumulator kinds. They need no online machinery, no additional checkpoint state, and no monitor-side configuration.
What genuinely must happen online is everything the identity does not reach:
Memory is the one remaining argument for accumulating a reduced quantity online rather than reducing later, and it is stratified: negligible for moments, decisive for per-point histograms, prohibitive for retained histories.
A physical caveat applies to spatial reduction however it is computed: merging across points is only meaningful where those points are statistically equivalent. Over an inhomogeneous region the result is a well-defined total variance that mixes in the spatial variation of the mean, which is rarely what was wanted. Homogeneous-direction averaging is the case this serves.
| Extension | Depends on | Unlocks |
|---|---|---|
| 4. Spatial Reduction In Post-Processing | nothing | profiles, regions, bins, clipping |
| 3. A View Bindable To An Explicit Vector Pair | nothing | nodal statistics output, accumulator field logging |
| 9. Effective Sample Size Reporting | nothing | uncertainty reporting under unequal weights |
| 6. Raising The Accumulated Moment Order | nothing | skewness, flatness, and their spatial reductions |
| 8. Further Products | online face targets, for the face case | face-flux and single-component products |
| 5. Online Targets For The Non-Commuting Cases | 3. A View Bindable To An Explicit Vector Pair for surface output | surfaces, probes, conditional sampling |
| 10. One Reducer Behind Flow Observables | 4. Spatial Reduction In Post-Processing | one reducer behind flow observables |
| 11. Dimensional Derived Statistics | nothing | dimensional derived statistics |
| 12. Derived Statistics At Non-Periodic Boundaries | nothing | meaningful statistics at walls, inlets, and outlets |
| 13. Offline Line And Plane Spectra | selected Cartesian planes/lines | ensembles, masks, and expanded plotting |
| 14. Temporal Spectra From Bounded Probe Histories | bounded probe history | frequency spectra where nothing is homogeneous |
| 15. Online Spatial Spectral Accumulation | 14. Temporal Spectra From Bounded Probe Histories for the transpose | spatial spectra without a snapshot series |
Spatial reduction is listed first because it is the largest capability for the least new code: the kernel it needs is already written and tested, and it has no online counterpart to design.
What is missing. UpdateLocalGhosts resolves its vectors through a field view built on compile-time UserCtx offsets recorded in the catalog, and the nodal averaging kernel maps fixed field names onto fixed members. Per-window accumulators are counted from configuration and therefore have no compile-time offset, so neither routine accepts one.
Current workaround. Derived results are staged through two catalogued scratch fields: a derived quantity is written there, reached by catalogued name, converted to nodal values, and copied out before the next quantity reuses the buffer. This works with the offset constraint rather than against it, and follows the precedent of the corner-staging buffers, which are also workspace rather than simulation state.
What the extension is. A FieldView that can be bound to an explicitly supplied global and local vector pair, carrying its own degree of freedom, layout, and synchronization class rather than reading them from a compile-time table. After that the whole ghost path applies to an accumulator unchanged.
Why it is worth doing. Two surfaces depend on it at once. Nodal statistics output stops needing a staging round trip, and accumulator vectors become addressable by the field-anatomy and min/max loggers, which today can only report window-level scalars. A statistics-specific logging path would solve half of that and add a second mechanism doing what the first already does.
Profiles, regions, bins, and clipping are derived from the pointwise state a window already holds, by the reduction described in 1. What Must Be Online, And What Must Not. This is where the majority of the remaining scientific capability lives, and it needs no new online state.
What it requires:
i and k at fixed j, or a named region. This sits in post.yml because the choice is made after the run, not committed at solve time.PicurvWindowSpatialMean is the shape to follow; it already performs the domain-wide case of exactly this reduction, and it already excludes points the mask never sampled..vts.The practical gain beyond cost is that regions are chosen after the solve. Because the pointwise state is checkpointed, a wake window or a wall-normal profile that nobody anticipated can be extracted from a completed run without re-solving, and revised as often as the analysis demands.
A target must be resolved online only where 1. What Must Be Online, And What Must Not says the reduction cannot follow:
SpatialTargetPlan exists as an abstraction with only the pointwise identity mapping implemented, so these extend it rather than retrofit it, and the identity hash already reserves its own group for the mask so a later mask key extends that group instead of renumbering the groups after it and invalidating every saved window.
The user-facing surface for these is a target: block with kind and mask as optional keys defaulting to current behaviour, so no configuration written today breaks. They are not exposed now because a single-valued key reserves syntax without earning it.
The second moment is what is accumulated today, which fixes the ceiling on what any later reduction can produce: order cannot be raised after the fact. Raising it to the third and fourth moments would make skewness and flatness available, and — because the merge identity extends to higher orders — would make them available for every spatial reduction at the same time, not only for the whole domain.
What this needs: additional per-point state, the higher-order terms of the weighted combination, the corresponding checkpoint payloads, and moment entries in the window definition so the requested order is part of the identity hash. The per-point cost is one additional field per order per accumulated quantity, which is modest next to what a per-point histogram would cost.
This is the extension that most enlarges what post-processing can do, because it raises the ceiling rather than adding a mechanism.
Per-point probability density functions and histograms are deferred rather than designed. They would reduce spatially without difficulty — histograms over shared bin edges are additive, so a region histogram is the sum of its points' — but the per-point cost is one field per bin per quantity, which is one to two orders of magnitude above the moment path. That cost, not the mathematics, is what defers them.
Should they arrive, bin edges must be fixed and global rather than per-point adaptive, since additivity is what makes the spatial reduction exact.
Retained time histories remain outside this pipeline at any point. Whatever retains history must stay bounded, for reduced probes or regions only; turning the statistics manager into a second full-field output system is not on this list. A run that needs the states themselves selects ordinary instantaneous field output.
Explicit Ucont face-flux products. Ucont is component-staggered, so pairs involving it are rejected on layout compatibility today. Face-flux statistics need a face target and a product definition that respects which component lives on which face, not a relaxation of the compatibility check.
Single-component products. Requesting one component of a vector's self-product rather than all six, for a run that needs only a shear stress and cannot afford the rest.
Vector-vector cross products. These need a full nine-component carrier, which nothing currently allocates. The symmetric six-component form exists precisely because a self-product is symmetric; a cross product between two different vectors is not.
Under physical-time weighting the weights are unequal, so a raw sample count overstates how much independent information a window holds. The Kish effective sample size \(W^2/W_2\) is the right quantity to report instead, and the kernel that computes it already exists and is tested.
What is missing is plumbing rather than mathematics: a window-level sum of squared weights, which means a new checkpoint metadata field and a restore path for it. Per-point squared weight is already accumulated and checkpointed; only the window-level scalar is absent.
ComputeCurrentFlowObservables computes its own domain reductions for the runtime convergence monitor. The reduction driver in 4. Spatial Reduction In Post-Processing performs the same traversal, so it could serve that monitor too and remove one of the two implementations.
This is permitted only after tests demonstrate identical results, and the replacement must correctly follow field layout, periodicity, masks, blocks, and MPI ownership. It is listed last among the functional extensions because the existing implementation is correct and the gain is structural.
Function-level logging, logging allow-lists, PETSc monitors, and performance profiling are not subsumed into scientific statistics at any point.
Derived statistics are written in nondimensional form even when global_operations.dimensionalize is set. A Reynolds stress scales as velocity squared, and a co-moment as the product of two different scales; the existing per-field scaling table expresses neither.
The extension is a scaling rule per derived quantity rather than per field, so a product knows the scales of both its factors. Applying the velocity scale as things stand would be wrong rather than merely incomplete, which is why the current behaviour is to leave derived output nondimensional and say so.
What exists today. SynchronizePeriodicCellFields populates the layout boundary of a derived statistics field on periodic faces only, where the value is exact. See p58_output_sec. On a wall, inlet, or outlet nothing is written, so the outermost node layer of the .vts keeps whatever the interior-only derivation left there.
Why nothing is written rather than something approximate. The correct ghost value is a property of the quantity crossed with the boundary type, and the staging buffer carries stresses, pressure, and eddy viscosity through the same vector on successive calls, so the convention cannot be attached to the buffer. Two candidate defaults were considered and both are wrong somewhere:
At a boundary node on a face the eight-cell stencil is four dummy plus four interior, so zero-gradient yields the mean of the adjacent interior cells and mirror yields exactly zero. They are not close, and picking either one globally is wrong for some quantity.
Split by moment order. Second moments need no per-field classification: a fluctuation about the window mean vanishes at any steady Dirichlet boundary whatever the field is, and follows the flow out at an outlet.
| Face | First moment | Second moment and co-moment |
|---|---|---|
PERIODIC (geometric, constant_flux) | wrap (implemented) | wrap (implemented) |
WALL (noslip) | per field class below | mirror |
INLET (constant_velocity, parabolic, prescribed_flow) | the inlet value | mirror |
OUTLET (conservation) | zero-gradient | zero-gradient |
First moments at a wall need to know whether the field vanishes there:
| Field | At a no-slip wall | Rule |
|---|---|---|
Ucat | vanishes | mirror |
Nu_t | vanishes at a resolved wall | mirror |
P | Neumann | zero-gradient |
CS | no imposed value; uniform in constant-coefficient mode | zero-gradient |
CS is a model coefficient rather than a transported quantity. The requirement that eddy viscosity vanish at a wall lands on Nu_t through \(\nu_t = (C_s \Delta)^2 |S|\), not on the coefficient, and in constant-coefficient mode CS is a uniform field that zero-gradient reproduces exactly while mirror would corrupt.
Nu_t with a wall model.** Mirror is correct only for a resolved no-slip wall. A wall function deliberately supplies a non-zero effective viscosity at the first cell, and mirroring would zero it. If wall_function.enabled is set and Nu_t is accumulated, refuse; do not silently mirror.ResolveStatisticsField accepts any catalogued name, so a window may accumulate Diffusivity or a metric field. A fixed table would apply a wrong default to the first one outside it. Refuse, naming the field and the missing classification, so the next person edits a table instead of debugging a boundary.Both refusals are narrow: periodic faces never consult the table, and second moments never consult the field classification, so a fully periodic case reaches neither however many fields it accumulates.
"Mean at the boundary equals the boundary value" and "fluctuations vanish at a Dirichlet
boundary" both hold only if that boundary is steady over the window. True for noslip, constant_velocity, and parabolic. Not guaranteed for prescribed_flow, which reads a profile that could in principle vary in time. Treating prescribed_flow as steady is a documented assumption rather than a checked one.
Note also that the rule must be applied to the statistics field, never copied as a value from the source field. At a wall the source's ghost is a function of the instantaneous interior cell, so copying it would place an instantaneous value in a mean field. The inlet is the one exception, because a steady Dirichlet value is constant in time and copying it is therefore exact.
The boundary node layer is one node thick and the convergence-history CSV already carries the authoritative domain mean, computed with never-sampled points excluded. The defect this would close is that an unwritten boundary looks like data; leaving it unwritten and documented removes the trap at a fraction of the cost. The work becomes worthwhile when wall-bounded statistics are read quantitatively from the field — which is the same class of case that needs 13. Offline Line And Plane Spectra.
Selected physical plane and line spectra are now experimental implementations in generators/spectra.gen; see p10_cap_spec_plane_spectrum and p10_cap_spec_line_spectrum for their current contracts. The remaining proposals below concern broader reductions, mask handling, plotting, and scale. They are not implemented options. Current tasks require Cartesian unmasked samples and explicit axes/physical-cell indices; they do not infer a reduction from topology.
The obvious design — infer the task from how many directions are periodic — does not survive contact with what people actually want. Two homogeneous directions admit at least four distinct deliverables, and they answer different questions:
The most common channel deliverable is not "the plane spectrum" at all: it is a plane transform reduced two different ways. So the map from topology to task is one-to-many, and the topology cannot pick for the user.
Periodicity is also necessary rather than sufficient. The pairwise consistency check in validate_and_prepare_boundary_conditions makes "this axis is periodic" a reliable, already-validated property, but homogeneity is the real requirement. An immersed body inside a periodic box breaks it while every boundary condition still reads PERIODIC, and the Nvert mask that reveals it is a runtime property that moves with the body, so no configuration-time inference can see it.
What may be inferred. The direction may, where it is unambiguous: a line_spectrum on a case with exactly one periodic axis has only one sensible direction, and defaulting it removes real tedium without removing a scientific decision. The task may not, for the reasons above.
That distinction matters beyond convenience. The value of the precondition gate is that it refuses: ask for something undefined and the reason is named before any compute is spent. Inferring the task inverts that — when the inference is right a configuration line is saved, and when it is wrong there is no error at all, only a plausible curve. A pipeline built on 1. What Must Be Online, And What Must Not should not trade a loud failure for a silent one to save a line of YAML.
What is worth doing instead is making the refusal constructive: a case whose spanwise axis is periodic could be told that line_spectrum along that axis would be valid, rather than only being told what is wrong.
The transform is the easy part. What genuinely changes:
shell_spectrum yields \(E(k)\) per step. A line spectrum yields \(E(k_1)\) per retained index, a plane spectrum \(E(k_x,k_z)\) per station. The long-format CSV generalizes with an index column; the plot does not, because "a spectrum at every wall-normal station" is a family or a contour rather than a curve. --plot-spectrum needs a station selector.shell_spectrum sidesteps immersed boundaries by refusing them outright. A line spectrum must test each transform line for blanked points and drop or refuse that line individually.generators/spectra.gen reads whole fields in Python. A 64-cubed box is under seven megabytes; a 512-cubed channel is a few gigabytes per snapshot. For channel-scale cases this, rather than sample count, is what eventually forces 15. Online Spatial Spectral Accumulation.Test case. examples/flat_channel, which is homogeneous in the streamwise and spanwise directions. Not examples/bent_channel: a 90-degree bend develops streamwise and is bounded on all four sides, so it has no homogeneous direction and belongs to 14. Temporal Spectra From Bounded Probe Histories.
What exists today. post.yml -> spectra measures spatial spectra offline from the instantaneous velocity in each committed checkpoint, documented in 8. spectra. The one implemented task, shell_spectrum, requires a triply periodic uniform Cartesian box.
What is missing, and why it is not a gap in the above. A spatial spectrum needs a direction that is uniformly spaced, periodic, and statistically homogeneous. A curved duct, a wall-bounded channel section, or any immersed geometry supplies none of the three: the streamwise direction develops, and the cross-stream directions are bounded. For these the meaningful quantity is a frequency spectrum \(E(f)\) at a point, which needs no homogeneity at all and is what resolves shedding and Strouhal content.
A frequency spectrum needs the signal itself, retained in time. That is the one thing 7. Distributions, Deferred rules out at full-field scale and explicitly leaves open at reduced scale: whatever retains history must stay bounded, for reduced probes or regions only. A handful of probes sampled every step is kilobytes; the same information as full-field snapshots is terabytes.
Why the snapshot path cannot substitute. A temporal transform of the existing checkpoint series is limited to \(f_{\max} = 1/(2\,\Delta t\,\texttt{data\_output\_frequency})\). Reaching a useful band means writing full three-dimensional fields every step or two, which is precisely the cost this extension avoids. The snapshot series can reach a low-frequency peak; it cannot reach an inertial range.
What it requires:
temporal_spectrum task in the spectra recipe, with Welch averaging over segments and a window function, since a finite record is not periodic in time.This is the extension that unlocks spectra for every geometry the spatial tasks refuse, and it is the smaller of the two remaining spectral pieces.
What is missing. Spatial spectra are measured offline from written snapshots, so the number of independent samples is bounded by how many full fields a run can afford to write, and the snapshot cadence has to be chosen before the run.
What the extension is. Accumulate \(\langle |\hat{u}(k)|^2 \rangle\) online: transform along the homogeneous direction each accepted state and add the modulus squared into a per-wavenumber accumulator.
Why it must be online rather than derived later. The transform is linear, so a time-averaged \(\hat{u}\) is recoverable from a time-averaged field. The modulus squared is not: \(\langle |\hat{u}|^2 \rangle \neq |\langle \hat{u} \rangle|^2\). This is the same rule 1. What Must Be Online, And What Must Not states, applied in Fourier space, and it is also why no reduction of the existing pointwise accumulators can yield a spectrum — a spectrum is the transform of a two-point correlation, and pointwise state holds no covariance between distinct points.
Why it is nonetheless second in priority. It is an efficiency change, not a capability change: it produces spectra the snapshot path already produces, for runs where writing the snapshots is the binding constraint. 14. Temporal Spectra From Bounded Probe Histories produces spectra that are otherwise unavailable at all.
What it requires:
MatCreateFFT needs an FFTW-enabled build, which the current configuration does not have. A single-direction transform is cheap enough that gathering the pencil and transforming directly is a viable alternative.Storage is negligible — wavenumber by wall-normal index by component — which is orders of magnitude below the per-point histograms 7. Distributions, Deferred defers.
These are recorded so they are not proposed again as gaps.
su0/su1/su2/sp importer. The old averaging state was removed rather than retained as a second workflow, and no reader, writer, or converter for it is planned.Nu_t and CS can be averaged like any other field, but modelled k is never presented as turbulent kinetic energy derived from resolved fluctuations.
1.9.8