Add prescribed kinematics for immersed boundaries - #1847
Conversation
Moving immersed boundaries are driven either by two-way coupling or by analytic velocity
expressions in the case file. The analytic path has two drawbacks for a parameter sweep or
an ensemble: the expressions are compiled into `case.fpp`, so every distinct motion forces
a rebuild of the simulation binary, and the motion is obtained by integrating the supplied
rates, so a restart does not reproduce the same trajectory.
Add `patch_ib(i)%kin_model`, which sets the body state directly from the current time at
each Runge-Kutta stage:
1 hinged flapping: roll about the lab x axis through the hinge and pitch about the body
spanwise axis through the same hinge, with amplitudes, mean pitch, frequency, pitch
phase lead, onset time and a raised-cosine onset ramp as inputs
2 the smoothed linear pitch ramp and hold of the AIAA low-Reynolds-number canonical
cases (Eldredge et al. 2009, Ol et al. 2010), with nominal pitch rate and smoothing
parameter as inputs
Because nothing is integrated, restarts are exact and a single case-optimized binary serves
every parameter value and every ensemble member.
This also corrects the angular velocity used for the surface condition. The existing path
advances the three angles componentwise from `angular_vel` and then uses that same array as
a vector in `omega x r` for the ghost-point velocity. Those two uses only agree when a
single angle changes; with the rotation composed as R = Rx(phi) Ry(theta) the lab-frame
angular velocity is phi' e_x + theta' Rx(phi) e_y. The new routine sets that vector
directly. Measured on a plate rolling 30 degrees while pitching 20 degrees, the difference
reaches 10 percent of the surface velocity.
Verified against closed-form kinematics: angles, centroid, centroid velocity and angular
velocity all match to 1e-15 every step, and the velocity carried by the body cells adjacent
to the fluid matches the rigid-body velocity to 8e-13, against 4e-3 to 1e-1 if the
Euler-angle rates are used as the angular velocity.
Claude-Session: https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
There was a problem hiding this comment.
Warning
Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.
Pull request overview
Adds runtime-evaluated prescribed kinematics for immersed boundary (IB) patches, enabling parameter sweeps without recompiling and making restarts reproduce identical trajectories.
Changes:
- Introduces new
patch_ib(i)%kin_*parameters (toolchain schema, validation, docs) for two prescribed kinematics models. - Updates simulation to evaluate prescribed kinematics per Runge–Kutta stage (and at initialization) instead of integrating rates.
- Broadcasts new IB kinematics parameters via MPI and initializes them in global parameter setup.
Reviewed changes
Copilot reviewed 10 out of 10 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/mfc/params/descriptions.py | Adds human-readable descriptions for new kin_* IB parameters. |
| toolchain/mfc/params/definitions.py | Registers new kin_* parameters in the toolchain registry. |
| toolchain/mfc/case_validator.py | Validates kin_model and required inputs; restricts to moving IB in 3D. |
| src/simulation/m_time_steppers.fpp | Calls prescribed kinematics during RK propagation using a computed stage time. |
| src/simulation/m_mpi_proxy.fpp | Adds MPI broadcasts for the new kinematics parameters. |
| src/simulation/m_ibm.fpp | Evaluates kinematics at initialization and implements the kinematics routines. |
| src/simulation/m_global_parameters.fpp | Initializes kinematics parameters for simulation. |
| src/pre_process/m_global_parameters.fpp | Initializes kinematics parameters for preprocessing. |
| src/common/m_derived_types.fpp | Extends IB patch derived type with kinematics parameters. |
| docs/documentation/case.md | Documents new IB kinematics parameters and models. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| t_stage = mytime + dt | ||
| if (time_stepper == time_stepper_rk3 .and. s == 2) t_stage = mytime + 0.5_wp*dt | ||
|
|
||
| $:GPU_PARALLEL_LOOP(private='[i, gbl_id]', copyin='[s, t_stage]') |
| _ib_attrs[f"kin_hinge({j})"] = (REAL, _ib_tags) | ||
| _ib_attrs[f"kin_offset({j})"] = (REAL, _ib_tags) |
|
|
||
| end function f_log_cosh | ||
|
|
||
| !> Prescribed hinged flapping kinematics (kin_model = 1). Roll phi about the lab x axis through the hinge and pitch theta about |
|
Claude Code Review Head SHA: b9fed79 Files changed:
Findings:
|
Two immersed-boundary examples exercising the prescribed-kinematics option, both of which
can be checked without any reference data:
examples/3D_ibm_flapping_plate a wing in glide begins to flap (kin_model = 1). The
motion has a closed form, so angles, centroid, centroid velocity and the lab-frame
angular velocity can all be compared directly against it, and the glide lift can be
compared with lifting-line theory for the aspect ratio.
examples/3D_ibm_pitchup_plate the AIAA low-Reynolds-number canonical pitch-up
(kin_model = 2), validated against the DNS of Jantzen et al. (2014), Phys. Fluids 26,
053606. Lift reaches 76 percent of the reference, consistent with modelling a finite
square-edged section against their infinitely thin plate.
Both readmes record what agrees and what does not, including the streamwise-force defect
for thin inclined plates reported in #1849, for which the pitch-up case is a reproducer
that needs no reference data: the ratio of drag to lift on a plate held at 45 degrees
should be about 1 and instead collapses after the ramp.
They also document why the section is four cells thick, which a companion 2D grid study
showed to be the point where the lift through the ramp converges; two cells appears to
agree better with the reference only because two errors cancel.
Claude-Session: https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
Every directory under examples/ is enumerated as a golden test, so the two example cases added here arrived without goldens and failed on every lane. Generating them exposed a second problem. The Example suite caps the grid at 25 cells per direction, which puts the pitch-up case's 5-percent-chord plate at a third of a cell: the body occupies no cells, the golden is a uniform field identical at step 0 and step 50, and no code change could ever perturb it. That case is skipped with the reason recorded, and kin_model = 2 is instead covered by a purpose-built test whose plate is four cells thick -- the minimum at which the body has an interior. The flapping example keeps its golden: at the capped resolution its section is still 1.25 cells thick, so the body exists and the field responds to it. Claude-Session: https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
|
Pushed
The Pitch Ramp goldens are a separate matter and are not fixed yet. A candidate of exactly 2.5 against a golden of 1.31 is a cell still at its initial condition in one run and evolved in the other, i.e. solid in one and fluid in the other - a classification flip, not drift. The plate's y faces sit at exactly 0.30 and 0.70, which are cell centres on this 25-cell grid, and The fix I would make next is to move the body off the grid so no face lands on a cell centre or boundary, then regenerate. Holding off until you say whether you want that or would rather drop the Pitch Ramp case, since |
Lines of Code
|
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1847 +/- ##
==========================================
+ Coverage 61.26% 61.41% +0.15%
==========================================
Files 84 84
Lines 22330 22409 +79
Branches 3265 3271 +6
==========================================
+ Hits 13680 13763 +83
+ Misses 6207 6195 -12
- Partials 2443 2451 +8 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
Worked through the three review comments. One is fixed, two do not hold - details so the reasoning is checkable rather than just asserted. Docstrings (
There is no scalar-versus-array mismatch either: both sit inside the same
Deriving the stage time from a stored abscissa table would be more robust against a future integrator, and I would take that change if the stage times were ever wrong today - but they are not, and a table would have to be added and kept in step with |
Why
Moving immersed boundaries are driven either by two-way coupling or by analytic velocity expressions in the case file. For a parameter sweep or an ensemble the analytic path has two drawbacks:
case.fpp, so every distinct motion forces a rebuild of the simulation binary — a three-point Strouhal sweep is three builds;What
patch_ib(i)%kin_modelsets the body state directly from the current time at each Runge-Kutta stage:Nothing is integrated, so restarts are exact and one case-optimized binary serves every parameter value and every ensemble member.
A correctness fix that comes with it
The existing path advances the three angles componentwise from
angular_vel, and then uses that same array as a vector inomega x rfor the ghost-point velocity. Those two uses only agree when a single angle is changing. With the rotation composed asR = Rx(phi) Ry(theta), the lab-frame angular velocity isphi' e_x + theta' Rx(phi) e_y. The new routine sets that vector directly. On a plate rolling 30 degrees while pitching 20 degrees the difference reaches 10 percent of the surface velocity.Testing
Against closed-form kinematics on a 3D plate: angles, centroid, centroid velocity and lab-frame angular velocity all match to 1e-15 at every step. The velocity carried by body cells adjacent to the fluid matches the rigid-body velocity to 8e-13, versus 4e-3 to 1e-1 if the Euler-angle rates are used as the angular velocity vector — so the check distinguishes the two. Also exercised on a 3D pitch-ramp validation case against published DNS.
./mfc.sh precheckpasses.Restricted to 3D by the validator, since the in-plane 2D rotation goes through a different path.
https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG