|
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
|
|
This page covers two related things: the geometric contract that makes a periodic axis legal at all (sections 2-4), and the driven-flow handlers that sustain a bulk flow through such an axis (section 5). The page name kept its original geometric anchor for reference stability; its subject is now periodic boundaries generally.
PICurv supports geometric periodicity for Eulerian fields. Periodic axes are derived exclusively from paired PERIODIC boundary conditions before the DMDA is created. The resulting internal periodic flags select PETSc periodic index wrapping; they are not user configuration options.
Particles that cross a periodic face re-enter through the paired face, shifted by the validated seam translation; this is the wrap particle policy (p44_cap_ppol_wrap), the default and only policy on periodic faces.
For every active periodic axis, the two opposite nodal surfaces must match pointwise under one nonzero constant Cartesian translation:
\[ \mathbf{x}_{+}(p) - \mathbf{x}_{-}(p) = \mathbf{T} \]
for every corresponding surface point \(p\). The translation may have components in all Cartesian directions, but it may not vary over the surface. Each periodic axis must contain at least four physical nodes. These checks are performed independently, so mixed periodic/non-periodic domains are supported.
This admits translated curvilinear periods, including a repeating wavy channel, while rejecting rotating, shearing, or otherwise nonmatching paired surfaces.
PETSc DM_BOUNDARY_PERIODIC wraps indices but does not translate coordinates or resolve the ownership conventions of every staggered field. PICurv therefore handles periodic Eulerian data by layout:
Ucont,The immutable field catalog records the conceptual layout and ghost-repair class for each persistent field. The actual periodic axes still come from the runtime DMDA and validated boundary pairs; field identity does not replace or override that state. See Field Identity and Layout Catalog.
Face-center and cell-center coordinate ghosts are wrapped and shifted by the validated translation. Higher-depth QUICK stencil preparation is handled separately from ordinary local ghost refresh.
A periodic synchronization (SynchronizePeriodicCellFields(), SynchronizePeriodicStaggeredFields()) refreshes the local vectors once in full, then, after each direction's transfer, refreshes only the local entries that mirror the two end planes that transfer wrote. Because the local vectors were current before the transfer, the result is bitwise the one a full UpdateLocalGhosts() gives, at a small fraction of its communication. The solver option -picurv_verify_ghost_updates checks every such partial refresh, and every caller's claim that the local vectors are already current, against a full update. It also fills each local vector the boundary pass deliberately leaves out of date (because a later full update rewrites it before any read) with NaN, and requires the pass's fields to be finite afterwards, so a read of a stale vector stops the run. It is slow, and meant for tests.
Configure both faces of an axis as periodic:
Do not add i_periodic, j_periodic, or k_periodic under models.domain; those legacy keys are rejected.
A streamwise-periodic domain has no inlet and no outlet, so nothing sustains a bulk flow against wall drag: left alone, the flow simply decays. A driven handler closes that loop. It measures the volumetric flux, compares it with a target, and feeds the difference back as a body force in the momentum equation.
Two handlers implement this. They differ in exactly one respect - where the target comes from - and are otherwise the same code.
| Handler | Target | target_flux parameter |
|---|---|---|
constant_flux | Prescribed by the user | Required |
initial_flux | Measured from the field the run starts with, then held | Rejected - supplying one is a configuration error |
initial_flux is the right choice when you have seeded a field whose flux is already what you want and you simply want it preserved - a restart from a precursor, or a developed field carried over from another case. constant_flux is the right choice when the flux is the thing you are prescribing, for instance when matching a literature bulk Reynolds number.
PreStep takes two different flux measurements in one sweep of lUcont (MeasureDrivenFluxes in src/BC_Handlers.c), because the controller needs a fast signal and a stable one and they cannot be the same number.
boundaryVelocityCorrection, a tactical trim applied directly to the boundary fluxes.bulkVelocityCorrection, which scales the momentum source. It is also the quantity initial_flux latches at startup.Driving the body force from the noisy signal would make the forcing oscillate; trimming the seam from the slow signal would leave the seam inconsistent. Hence two sensors.
These two corrections are refreshed on deliberately different cadences.
ApplyBoundaryConditions() runs the handler sweep three times per call, and it is itself called once per Jameson RK stage under the Picard solver and once per residual evaluation under Newton-Krylov. "Once per `PreStep`" is therefore very far from "once per timestep".
**bulkVelocityCorrection is computed once per physical timestep** and held for the rest of it. Both momentum solvers depend on this. The Picard shadow-Jacobian estimate treats the body force as a constant forcing with zero velocity Jacobian, and a matrix-free Newton solve needs a source that does not drift between residual evaluations. A source that moved with the trial vector would silently change the operator both solvers think they are solving.
Freezing the correction is necessary but not sufficient on its own. The force actually applied is produced downstream by ComputeDrivenChannelFlowSource() (src/BodyForces.c), which smooths it against the previous step with a 50/50 exponential moving average held in SimCtx. That function runs from ComputeRHS, so it too executes once per RK stage or residual evaluation. Advancing the EMA on every call walked the applied force toward its target within a single timestep - 0.5, then 0.75, then 0.875 of the way - making the force depend on how many residual evaluations preceded it. Both the correction and the EMA are therefore gated on simCtx->step, and only the pair of them makes the source genuinely constant across a timestep.
boundaryVelocityCorrection is re-measured on every pass.** Apply() accumulates the trim into Ucont, and re-measuring is exactly what makes that accumulation self-limiting as the seam converges. Freezing it would make repeated Apply() calls add the same trim over and over.Because the momentum source is frozen at the start of a step, the controller holds the flux to first order in dt: the force applied during a step is the one that would have corrected the flux measured at its beginning. Halving dt halves the lag.
The controller is proportional, so it settles below the target. The applied force is f = G * (U_target - U_b) / dt with gain G = 1.5 * 1.8 = 2.7 (the BDF2 coefficient times the internal -driven_flow_scaling_factor, default 1.8, which no YAML key sets). In steady state that force balances the wall drag, so the realized bulk velocity falls short of the target by f * dt / 2.7. For a laminar channel of half-height h, where f = 3 nu U_b / h^2, this gives U_b = U_target / (1 + 3 nu dt / (2.7 h^2)): 0.056% short for the shipped laminar case, and 10% short at nu = 0.1, dt = 1, where the measurement above confirmed the law to second order in the grid. Read the realized flux and the applied force from <run.analysis.metrics>/driven_flow.csv (bulk_velocity, driving_acceleration) rather than assuming the target; in turbulent runs at typical dt the shortfall is far below the statistical noise.
The controller has two actuators, and this option controls the second one.
The body force is the gentle, global actuator: it pushes on every fluid cell and is what actually sustains the flow. The seam trim is a local override. It adds a correction straight into the face fluxes of the single boundary plane where the domain wraps around:
"Trim" is aviation vocabulary - a small standing correction that holds a quantity on target without touching the main control. Here the quantity is the volumetric flux through the seam plane, and the main control is the body force. The option was originally spelled apply_trim, which said that something was trimmed but not what or why; it is now enforce_seam_flux, and the old spelling is still accepted as a deprecated alias in both case.yml and bcs.run.
Note the asymmetry: ucont accumulates (+=) while uch overwrites (=). The trim is re-measured on every boundary pass, and that re-measurement is what makes the accumulation self-limiting - each pass sees a seam that is already closer to target, so the next correction is smaller. This is why the trim correction is deliberately not frozen per timestep the way the body force is (section 5.2); freezing it would make repeated passes add the same correction over and over.
It is off by default, and that is the right default. Reasons to leave it off:
Reasons to turn it on: you have diagnosed an actual flux mismatch at the seam - inspect Bcs.Uch, which records the trim whether or not it is applied - and the body force alone is not closing it.
Practical workflow: run with it off, look at Uch. If the recorded trim is negligible compared with the bulk flux, the seam is fine and enabling it would only add noise. If it is persistently large, that is a real seam problem, and enabling the trim treats the symptom while you look for the cause.
The two handlers restart differently, and the difference is deliberate.
constant_flux re-reads target_flux from the bcs file on every run, so editing it between segments takes effect. Nothing is carried in the checkpoint.
initial_flux cannot do that: its target is a property of the field the run originally started from. Re-measuring it after a restart would silently retarget the controller at whatever the flux had drifted to. The latched value is therefore written into the checkpoint manifest as
and restored by RestoreDrivenFluxTarget() in src/io.c before the first PreStep of the resumed run. A checkpoint written before this metadata existed has no entry; the controller then falls back to re-measuring and logs a warning saying so.
One consequence worth knowing: the latch does not fire during setup. Boundary handlers are initialized before InitializeEulerianState(), and the setup-time boundary sweep runs before the initial condition has been scattered into the ghosted lUcont the measurement reads. Latching there would record zero. The handler therefore waits for the first PreStep of the first real timestep, and does no work at all until then.
Both faces of the driven axis must carry the same handler, and the runtime rejects any INLET, OUTLET, or FARFIELD face while a driven handler is active - the flux controller and an inlet would be fighting over the same degree of freedom.
target_flux is a volumetric flux, U_b times the cross-sectional area, not a velocity. For initial_flux, drop the params block entirely.
examples/periodic_test/driven_channel/laminar/ - plane channel with an exact closed-form answer. For half-height h, viscosity nu and the converged body force f: u(y) = (f / 2nu) (h^2 - (y-h)^2) U_b = f h^2 / 3nu -> f = 3 nu U_b / h^2 u_tau = sqrt(f h)This verifies the source term and the flux controller with no turbulence model involved.
examples/periodic_test/driven_channel/turbulent_retau180/ and .../turbulent_retau395/ - DNS-resolution plane channels against Moser, Kim & Mansour (1999), Phys. Fluids 11, 943.examples/periodic_test/driven_channel/les_retau180/ - the coarse LES repeat. les_wallmodel_retau1000 is the variant validated against Lee & Moser, with dynamic Smagorinsky and the Werner wall model. The channel is periodic in xi and zeta, so averaging.mode: homogeneous derives those two directions from the boundary pairs and produces a wall-normal coefficient profile. See 5. LES Subgrid Models and LES Turbulence Closure.examples/periodic_test/driven_duct/ - square duct at Re_b = 4410, which sustains turbulence-driven secondary flow of the second kind that a plane channel does not. References: Gavrilakis (1992), JFM 244, 101; Huser & Biringen (1993), JFM 257, 65.Grids for all of these are generated from checked-in .cfg files under config/grids/, not shipped as binary .picgrid. Wall spacing is a single number per walled axis, expressed as a fraction of that axis's length, so y+ can be retargeted for a different Re_tau without regenerating anything by hand.
Initial-condition seeding. The channel_spectral_velocity and duct_spectral_velocity providers generate wall-bounded discrete-curl perturbations; see p33_cap_gen_channel_spectral_velocity and p33_cap_gen_duct_spectral_velocity. The turbulent_channel example and every turbulent and LES case under periodic_test/driven_channel/ and periodic_test/driven_duct/ seed with them; only the laminar channel starts from streamwise_constant. The perturbation grows linearly from the wall, as a developed near-wall field does: on the Re_tau = 180 DNS grid its wall-parallel RMS is about 0.013 U_b at y+ ~ 5 and 0.031 U_b at y+ ~ 13, against 0.13 U_b at the centreline (developed turbulence peaks near 0.17 U_b at y+ ~ 13). Transition from it is still not guaranteed. Neither a seeded perturbation nor a short successful solve establishes developed turbulence. The triply periodic spectral_random_velocity provider remains unsuitable for no-slip walls.
Wall-normal profiles from statistics. field_statistics is BC-agnostic - it accumulates per-cell time moments and works unchanged under periodic boundaries. What does not exist is a spatial reduction: nothing in the postprocessor averages a statistics field over homogeneous directions to produce the U+ vs y+ and RMS profiles a DNS comparison needs. That reduction currently has to be done outside the postprocessor from the window payloads.
Momentum convergence on periodic wall-bounded flow.
Re-characterized at current HEAD on 2026-09-18 (measurement periodic-channel-laminar-picard-2026-09-18); the stall recorded on 2026-08-24 does not reproduce. That observation predated the 2026-08-25 move of momentum convergence onto a residual criterion and the change to the pseudo-CFL controller. The shipped laminar channel (Re = 100, 17 x 33 x 17 cells, 4 ranks, dt = 0.05, initial_flux) now converges every step in 9 to 29 pseudo-iterations against a cap of 50. At Re = 10 the steady state reproduces the exact parabola at second order with both driven handlers: the maximum profile error falls 1.8e-2, 5.2e-3, 1.4e-3 on 9, 17 and 33 wall-normal cells.
At startup, the banner reports BC-derived periodic axes and each validated translation. Use make unit-periodic for focused synchronization and geometry coverage and make smoke-periodic for the real runtime periodic workflow. Both are included in the standard unit and check gates.
For driven flows, tests/smoke/run_driven_periodic_regression.sh asserts that the initial_flux target matches the analytic flux of the initial condition and survives a restart unchanged, alongside the Poisson coarse-solve residual check. The driven controller logs its target, both measured fluxes and both corrections at INFO every step; the realized flux should track the target to the first-order-in-dt lag described in section 5.2.
1.9.8