rtx-fsi + rtx-cfd + rtx-fea: FSI1 — the coupled cylinder and flag
The summit rung C1: all the verified pieces joined into the first coupled Turek–Hron computation (rtx-fsi/tests/turek_hron_fsi1.rs). The embedded fluid computes tractions on the DEFORMED flag surface (EmbeddedMask::traction_at, factored from surface_force); the flag's wetted boundary is a polygon whose vertex list sits behind a lock, so the moving-body mask rebuild picks up every shape update (EmbeddedBody::polygon + pub polygon_signed_distance); WettedSurface — rebuilt on the deformed interface every subiteration — carries the loads to the flag's boundary nodes (NonlinearStaticAnalysis::set_nodal_forces); Subiterated::aitken drives the exchange, each pass marching the fluid to flag-load stagnation so the coupling map is a function of geometry, not of the fluid's transient. Result (ny = 62, 6 Aitken passes, 420 s): coupled drag 15.360 (+7.5%, the rigid CFD1 band at this grid), lift 0.7977 (+4.4%), ux(A) 2.647e-5 vs 2.270e-5 (+16.6%; +6.1% at ny = 82), uy(A) 3.90e-4 vs 8.21e-4 at h = 6.6 mm and 1.124e-3 (+37%) at h = 5 mm — the resolutions BRACKET the reference through the flag's 3 -> 4-cell thickness transition, like the rigid-flag lift; conservation 7.4e-12 every pass. Bands asserted are the measured ones; RTX_FSI1_NY runs studies. Two real rtx-fsi defects found by this rung (15th and 16th of the campaign), both regression-tested (tests/transfer_curved_edge.rs): 1. solve_weights built its constraint Gram from RAW coordinates: the condition number grows as (position/spacing)^2 — ~1e4 for a flag edge at x ~ 0.26 with 5 mm spacing — and the 4x4 SVD pseudo-inverse lost enough accuracy that the (correctly strict) partition-of-unity / reproduction verification rejected healthy neighbourhoods: the operator's behaviour depended on WHERE the interface sat. Now centred on the face and scaled by the neighbourhood radius — identical constraints, O(1) conditioning, translation-invariant. 2. A NEARLY collinear neighbourhood (the nearest nodes of a face on a smoothly deformed edge: y is almost linear in x, off by the curvature sagitta) cannot satisfy exact centroid reproduction with bounded weights — the offending singular value is too large to truncate and too small to invert. The recruitment now widens (8 -> 16 -> 32 -> all) until the verified constraints hold; for a thin structure that pulls in the opposite face, exactly the transverse spread the system needs. Findings measured before believed: the transfer is faithful (a strictly local two-node split of the same tractions moved the tip by 2%); the uy error is the sampled lift PROFILE on a 3-cell flag (a uniform distribution of the same net lift bends 4x more), confirmed by the resolution study; TVD limiter chatter (+-0.5% steady load — limited schemes stall short of machine steady state) defeats steady fixed-point coupling, so steady coupled cases run upwind while the time-marched FSI2/FSI3 keep TVD; and the mask never chattered at FSI1's sub-cell amplitude (fluid-cell count constant through every pass). rtx-fsi 29 -> 31 green (lib 27, piston 2, curved-edge 1, FSI1 1). Co-Authored-By: Claude Fable 5 <[email protected]>
This commit is contained in:
Omar Sobh
co-authored by
Claude Fable 5
Claude Fable 5
parent
38ca5ef080
commit
b82f307cae
@@ -0,0 +1,607 @@
|
||||
//! Turek–Hron FSI1: the first coupled cylinder-plus-flag computation —
|
||||
//! rung C1 of the ladder (omni-cortex `docs/turek_hron_geometry_decision.md`).
|
||||
//!
|
||||
//! Re = 20 channel flow past the rigid cylinder with the ELASTIC flag:
|
||||
//! the embedded-boundary fluid (TVD + multigrid, rung F1/F2) provides
|
||||
//! tractions on the deformed flag surface, `rtx-fsi`'s `WettedSurface`
|
||||
//! carries them to the flag's boundary nodes (rebuilt on the deformed
|
||||
//! interface every subiteration — the small-displacement limit retired in
|
||||
//! practice), the total-Lagrangian St. Venant–Kirchhoff flag (rung S1)
|
||||
//! solves statically, and `Subiterated::aitken` drives the exchange to a
|
||||
//! fixed point. The flag's wetted boundary lives as a polygon whose vertex
|
||||
//! list sits behind a lock: the fluid's moving-body path re-reads it on
|
||||
//! every step's mask rebuild (rung F2), so a shape update is just a write
|
||||
//! to that list.
|
||||
//!
|
||||
//! FSI1 is steady and its tip displacement (reference `ux(A) = 0.0227 mm`,
|
||||
//! `uy(A) = 0.8209 mm`) is a fifth of a fluid cell — it validates the
|
||||
//! COUPLING machinery, not large deformation: loads, transfer,
|
||||
//! conservation, and the fixed point. Reference values (FEATFLOW level 7):
|
||||
//! `ux(A) = 2.270493e-5 m`, `uy(A) = 8.208773e-4 m`, drag 14.29426, lift
|
||||
//! 0.763746 on cylinder + flag.
|
||||
//!
|
||||
//! Measured values and the assertion bands are recorded at the bottom once
|
||||
//! the first run lands; conservation of the transferred load (partition of
|
||||
//! unity) is asserted at 1e-10 every pass.
|
||||
|
||||
use std::cell::RefCell;
|
||||
use std::sync::{Arc, RwLock};
|
||||
|
||||
use nalgebra::Vector3;
|
||||
use rtx_cfd::CfdConfig;
|
||||
use rtx_cfd::solvers::incompressible::{
|
||||
AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
|
||||
FlowField, PoissonSolverKind, SideBoundary, polygon_signed_distance,
|
||||
};
|
||||
use rtx_fea::analysis::{Analysis, AnalysisConfig, NonlinearConfig, NonlinearStaticAnalysis};
|
||||
use rtx_fea::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
|
||||
use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
|
||||
use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
|
||||
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
|
||||
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
|
||||
use rtx_fsi::{FluidFace, Subiterated, WettedSurface};
|
||||
|
||||
const L: f64 = 2.5;
|
||||
const H: f64 = 0.41;
|
||||
const RHO_F: f64 = 1000.0;
|
||||
const NU_F: f64 = 1e-3;
|
||||
const U_MEAN: f64 = 0.2;
|
||||
const E_S: f64 = 1.4e6;
|
||||
const NU_S: f64 = 0.4;
|
||||
|
||||
const FLAG_X0: f64 = 0.25;
|
||||
const FLAG_X1: f64 = 0.6;
|
||||
const FLAG_Y0: f64 = 0.19;
|
||||
const FLAG_Y1: f64 = 0.21;
|
||||
|
||||
const REF_UX: f64 = 2.270_493e-5;
|
||||
const REF_UY: f64 = 8.208_773e-4;
|
||||
const REF_DRAG: f64 = 14.294_26;
|
||||
const REF_LIFT: f64 = 0.763_746;
|
||||
|
||||
fn circle_sdf(x: f64, y: f64) -> f64 {
|
||||
((x - 0.2).powi(2) + (y - 0.2).powi(2)).sqrt() - 0.05
|
||||
}
|
||||
|
||||
fn inflow(y: f64) -> f64 {
|
||||
1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
|
||||
}
|
||||
|
||||
/// The flag's Quad8 mesh (as in rtx-fea's CSM tests).
|
||||
fn flag_mesh(nx: usize, ny: usize) -> Mesh {
|
||||
let mut mesh = Mesh::new(2).unwrap();
|
||||
let (lx, ly) = (2 * nx + 1, 2 * ny + 1);
|
||||
let mut grid = vec![vec![None; ly]; lx];
|
||||
for (i, column) in grid.iter_mut().enumerate() {
|
||||
for (j, slot) in column.iter_mut().enumerate() {
|
||||
if i % 2 == 1 && j % 2 == 1 {
|
||||
continue;
|
||||
}
|
||||
let x = FLAG_X0 + (FLAG_X1 - FLAG_X0) * i as f64 / (2 * nx) as f64;
|
||||
let y = FLAG_Y0 + (FLAG_Y1 - FLAG_Y0) * j as f64 / (2 * ny) as f64;
|
||||
*slot = Some(mesh.add_node(Node::new_2d(x, y)));
|
||||
}
|
||||
}
|
||||
for i in 0..nx {
|
||||
for j in 0..ny {
|
||||
let (a, b) = (2 * i, 2 * j);
|
||||
let nodes = vec![
|
||||
grid[a][b].unwrap(),
|
||||
grid[a + 2][b].unwrap(),
|
||||
grid[a + 2][b + 2].unwrap(),
|
||||
grid[a][b + 2].unwrap(),
|
||||
grid[a + 1][b].unwrap(),
|
||||
grid[a + 2][b + 1].unwrap(),
|
||||
grid[a + 1][b + 2].unwrap(),
|
||||
grid[a][b + 1].unwrap(),
|
||||
];
|
||||
mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
|
||||
.unwrap();
|
||||
}
|
||||
}
|
||||
mesh
|
||||
}
|
||||
|
||||
struct Interface {
|
||||
/// Wetted boundary nodes (everything on the bottom/tip/top edges except
|
||||
/// the clamped left corners), sorted by id — the coupling vector is
|
||||
/// their `(ux, uy)` pairs in this order.
|
||||
wetted: Vec<NodeId>,
|
||||
/// Reference positions of the wetted nodes.
|
||||
reference: Vec<(f64, f64)>,
|
||||
/// The ordered boundary walk for the polygon: indices into `wetted`
|
||||
/// (`usize::MAX` marks the fixed anchor vertices).
|
||||
walk: Vec<(usize, (f64, f64))>,
|
||||
}
|
||||
|
||||
impl Interface {
|
||||
fn build(mesh: &Mesh) -> Self {
|
||||
let eps = 1e-9;
|
||||
let on_bottom = |p: Vector3<f64>| (p.y - FLAG_Y0).abs() < eps;
|
||||
let on_top = |p: Vector3<f64>| (p.y - FLAG_Y1).abs() < eps;
|
||||
let on_tip = |p: Vector3<f64>| (p.x - FLAG_X1).abs() < eps;
|
||||
let clamped = |p: Vector3<f64>| (p.x - FLAG_X0).abs() < eps;
|
||||
|
||||
let mut wetted: Vec<(NodeId, (f64, f64))> = mesh
|
||||
.nodes
|
||||
.iter()
|
||||
.filter(|(_, node)| {
|
||||
let p = node.position();
|
||||
(on_bottom(p) || on_top(p) || on_tip(p)) && !clamped(p)
|
||||
})
|
||||
.map(|(&id, node)| (id, (node.position().x, node.position().y)))
|
||||
.collect();
|
||||
wetted.sort_by_key(|(id, _)| *id);
|
||||
let index_of = |id: NodeId| wetted.iter().position(|(w, _)| *w == id).unwrap();
|
||||
|
||||
// Ordered walk, counterclockwise: anchor inside the cylinder, the
|
||||
// clamped bottom corner, bottom edge left -> right, tip bottom ->
|
||||
// top, top edge right -> left, the clamped top corner, anchor.
|
||||
let mut bottom: Vec<(NodeId, f64)> = mesh
|
||||
.nodes
|
||||
.iter()
|
||||
.filter(|(_, n)| on_bottom(n.position()) && !clamped(n.position()))
|
||||
.map(|(&id, n)| (id, n.position().x))
|
||||
.collect();
|
||||
bottom.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
|
||||
let mut tip: Vec<(NodeId, f64)> = mesh
|
||||
.nodes
|
||||
.iter()
|
||||
.filter(|(_, n)| {
|
||||
let p = n.position();
|
||||
on_tip(p) && !on_bottom(p) && !on_top(p)
|
||||
})
|
||||
.map(|(&id, n)| (id, n.position().y))
|
||||
.collect();
|
||||
tip.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
|
||||
let mut top: Vec<(NodeId, f64)> = mesh
|
||||
.nodes
|
||||
.iter()
|
||||
.filter(|(_, n)| on_top(n.position()) && !clamped(n.position()))
|
||||
.map(|(&id, n)| (id, n.position().x))
|
||||
.collect();
|
||||
top.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap());
|
||||
|
||||
let mut walk: Vec<(usize, (f64, f64))> = Vec::new();
|
||||
walk.push((usize::MAX, (0.22, FLAG_Y0)));
|
||||
walk.push((usize::MAX, (FLAG_X0, FLAG_Y0)));
|
||||
for (id, _) in &bottom {
|
||||
walk.push((index_of(*id), (0.0, 0.0)));
|
||||
}
|
||||
for (id, _) in &tip {
|
||||
walk.push((index_of(*id), (0.0, 0.0)));
|
||||
}
|
||||
for (id, _) in &top {
|
||||
walk.push((index_of(*id), (0.0, 0.0)));
|
||||
}
|
||||
walk.push((usize::MAX, (FLAG_X0, FLAG_Y1)));
|
||||
walk.push((usize::MAX, (0.22, FLAG_Y1)));
|
||||
|
||||
let reference = wetted.iter().map(|(_, p)| *p).collect();
|
||||
Self {
|
||||
wetted: wetted.into_iter().map(|(id, _)| id).collect(),
|
||||
reference,
|
||||
walk,
|
||||
}
|
||||
}
|
||||
|
||||
/// Deformed polygon vertices for the interface vector `d`.
|
||||
fn polygon(&self, d: &[f64]) -> Vec<(f64, f64)> {
|
||||
self.walk
|
||||
.iter()
|
||||
.map(|&(k, anchor)| {
|
||||
if k == usize::MAX {
|
||||
anchor
|
||||
} else {
|
||||
let (x0, y0) = self.reference[k];
|
||||
(x0 + d[2 * k], y0 + d[2 * k + 1])
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Deformed wetted node positions for the transfer.
|
||||
fn deformed_nodes(&self, d: &[f64]) -> Vec<Vector3<f64>> {
|
||||
self.reference
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(k, &(x0, y0))| Vector3::new(x0 + d[2 * k], y0 + d[2 * k + 1], 0.0))
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// One static TL solve of the flag under the given wetted nodal forces;
|
||||
/// returns the new interface vector and the tip displacement.
|
||||
fn solve_flag(
|
||||
mesh: &Mesh,
|
||||
interface: &Interface,
|
||||
forces: &[Vector3<f64>],
|
||||
a_node: NodeId,
|
||||
) -> (Vec<f64>, (f64, f64), usize) {
|
||||
let clamped: Vec<NodeId> = mesh
|
||||
.nodes
|
||||
.iter()
|
||||
.filter(|(_, node)| (node.position().x - FLAG_X0).abs() < 1e-9)
|
||||
.map(|(&id, _)| id)
|
||||
.collect();
|
||||
let mut bcs = BoundaryConditionSet::new();
|
||||
for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
|
||||
bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
|
||||
nodes: clamped.clone(),
|
||||
components: vec![component],
|
||||
condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))),
|
||||
time_range: None,
|
||||
ramping_factor: 1.0,
|
||||
gradual_enforcement: false,
|
||||
}));
|
||||
}
|
||||
let mut db = MaterialDatabase::new();
|
||||
db.add_material(
|
||||
MaterialId(0),
|
||||
LinearElastic::new(E_S, NU_S).with_density(RHO_F),
|
||||
None,
|
||||
);
|
||||
let mut analysis = NonlinearStaticAnalysis::new(
|
||||
mesh.clone(),
|
||||
db,
|
||||
bcs,
|
||||
NonlinearConfig::default(),
|
||||
AnalysisConfig::default(),
|
||||
)
|
||||
.with_total_lagrangian();
|
||||
analysis.set_nodal_forces(
|
||||
interface
|
||||
.wetted
|
||||
.iter()
|
||||
.zip(forces)
|
||||
.map(|(&id, &f)| (id, f))
|
||||
.collect(),
|
||||
);
|
||||
let results = analysis.run().unwrap();
|
||||
assert!(
|
||||
results.convergence.converged,
|
||||
"flag Newton did not converge"
|
||||
);
|
||||
let numbering =
|
||||
AdvancedDofNumbering::displacement_only(mesh, DofMappingStrategy::Sequential).unwrap();
|
||||
let mut d = vec![0.0; 2 * interface.wetted.len()];
|
||||
for (k, &id) in interface.wetted.iter().enumerate() {
|
||||
let dofs = numbering.get_node_dofs(id);
|
||||
d[2 * k] = results.displacements[dofs[0]];
|
||||
d[2 * k + 1] = results.displacements[dofs[1]];
|
||||
}
|
||||
let a_dofs = numbering.get_node_dofs(a_node);
|
||||
(
|
||||
d,
|
||||
(
|
||||
results.displacements[a_dofs[0]],
|
||||
results.displacements[a_dofs[1]],
|
||||
),
|
||||
results.convergence.iterations,
|
||||
)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn fsi1_coupled_cylinder_and_flag() {
|
||||
let ny: usize = std::env::var("RTX_FSI1_NY")
|
||||
.map(|v| v.parse().unwrap())
|
||||
.unwrap_or(62);
|
||||
let h = H / ny as f64;
|
||||
let nx = (L / h).round() as usize;
|
||||
let mu = RHO_F * NU_F;
|
||||
let u_peak = 1.5 * 1.5 * U_MEAN;
|
||||
let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));
|
||||
|
||||
let mesh = flag_mesh(35, 2);
|
||||
let interface = Interface::build(&mesh);
|
||||
let a_node = mesh
|
||||
.nodes
|
||||
.iter()
|
||||
.find(|(_, n)| (n.position().x - 0.6).abs() < 1e-9 && (n.position().y - 0.2).abs() < 1e-9)
|
||||
.map(|(&id, _)| id)
|
||||
.expect("point A");
|
||||
|
||||
// The deformable geometry: the flag polygon behind a lock; the fluid's
|
||||
// per-step mask rebuild reads it.
|
||||
let vertices = Arc::new(RwLock::new(
|
||||
interface.polygon(&vec![0.0; 2 * interface.wetted.len()]),
|
||||
));
|
||||
let sdf_vertices = vertices.clone();
|
||||
|
||||
let config = CfdConfig::new()
|
||||
.with_density(RHO_F)
|
||||
.with_viscosity(mu)
|
||||
.with_reference_velocity(U_MEAN)
|
||||
.with_reference_length(0.1);
|
||||
let params = EmbeddedParameters {
|
||||
corrector_steps: 2,
|
||||
tolerance: 1e-7,
|
||||
boundaries: AleBoundaries {
|
||||
left: SideBoundary::Velocity,
|
||||
right: SideBoundary::PressureOutlet,
|
||||
bottom: SideBoundary::Velocity,
|
||||
top: SideBoundary::Velocity,
|
||||
},
|
||||
poisson_solver: PoissonSolverKind::Multigrid,
|
||||
// Upwind, deliberately: FSI1 is a steady FIXED-POINT problem, and
|
||||
// the TVD limiter's switching keeps the steady load chattering by
|
||||
// ~0.5% (a known property of limited schemes — they stall short of
|
||||
// machine steady state), which the coupling inherits as a ±4% tip
|
||||
// jitter and a 2e-5 interface-residual floor. Upwind converges to
|
||||
// machine steady, the coupling map is deterministic, and at Re 20
|
||||
// its loads are within a few percent (CFD1: surface lift +2.3% at
|
||||
// this grid). The unsteady FSI2/FSI3 march in time and keep TVD.
|
||||
convection_scheme: ConvectionScheme::Upwind,
|
||||
};
|
||||
let mut solver = EmbeddedPisoSolver::new(config, params).unwrap();
|
||||
solver.set_boundary_velocity(|x, y, _| {
|
||||
if x <= 0.0 {
|
||||
(inflow(y), 0.0)
|
||||
} else {
|
||||
(0.0, 0.0)
|
||||
}
|
||||
});
|
||||
solver.set_moving_body(EmbeddedBody::from_sdf(move |x, y, _| {
|
||||
let poly = sdf_vertices.read().unwrap();
|
||||
circle_sdf(x, y).min(polygon_signed_distance(&poly, x, y))
|
||||
}));
|
||||
|
||||
let mut field = FlowField::new(nx, ny, h, h).unwrap();
|
||||
for j in 0..ny {
|
||||
let u0 = inflow((j as f64 + 0.5) * h);
|
||||
for i in 0..=nx {
|
||||
field.u[(j, i)] = u0;
|
||||
}
|
||||
}
|
||||
solver.initialize(&mut field).unwrap();
|
||||
|
||||
// Warm-start the fluid on the undeformed geometry.
|
||||
let start = std::time::Instant::now();
|
||||
for _ in 0..std::env::var("FSI1_WARM")
|
||||
.map(|v| v.parse().unwrap())
|
||||
.unwrap_or(6000)
|
||||
{
|
||||
futures::executor::block_on(solver.advance(&mut field, dt)).unwrap();
|
||||
}
|
||||
println!(
|
||||
" warm start: 6000 steps, {:.0} s",
|
||||
start.elapsed().as_secs_f64()
|
||||
);
|
||||
|
||||
// Everything the coupling pass mutates.
|
||||
let state = RefCell::new((solver, field));
|
||||
let tip = RefCell::new((0.0f64, 0.0f64));
|
||||
let previous_d = RefCell::new(vec![0.0f64; 2 * interface.wetted.len()]);
|
||||
let state_mask_cells = std::cell::Cell::new(0usize);
|
||||
let worst_conservation = RefCell::new(0.0f64);
|
||||
let total_skipped = RefCell::new(0usize);
|
||||
|
||||
let pass = |d: &[f64]| -> Vec<f64> {
|
||||
// 1. The fluid sees the deformed flag, and marches until the
|
||||
// sampled flag load has stopped moving — the coupling map must be
|
||||
// a deterministic function of the geometry, or Aitken chases the
|
||||
// fluid's own transient (measured: fixed-length passes left a
|
||||
// 0.5% load jitter and a 2e-5 interface plateau).
|
||||
*vertices.write().unwrap() = interface.polygon(d);
|
||||
let (solver, field) = &mut *state.borrow_mut();
|
||||
let poly_probe = EmbeddedBody::polygon(vertices.read().unwrap().clone());
|
||||
let cap: usize = std::env::var("FSI1_PASS")
|
||||
.map(|v| v.parse().unwrap())
|
||||
.unwrap_or(6000);
|
||||
let mut history: Vec<f64> = Vec::new();
|
||||
let mut marched = 0usize;
|
||||
loop {
|
||||
for _ in 0..100 {
|
||||
futures::executor::block_on(solver.advance(field, dt)).unwrap();
|
||||
}
|
||||
marched += 100;
|
||||
let mask_now = solver.mask().unwrap();
|
||||
let body_now = solver.body().unwrap();
|
||||
let mut lift = 0.0;
|
||||
for s in poly_probe.surface_samples(h) {
|
||||
if circle_sdf(s.x, s.y) < 1e-9 {
|
||||
continue;
|
||||
}
|
||||
if let Some((_, ty)) = mask_now.traction_at(
|
||||
body_now, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
|
||||
) {
|
||||
lift += ty * s.ds;
|
||||
}
|
||||
}
|
||||
history.push(lift);
|
||||
if history.len() >= 4 {
|
||||
let now = history[history.len() - 1];
|
||||
let then = history[history.len() - 4];
|
||||
if ((now - then) / now.abs().max(1e-30)).abs() < 2e-5 || marched >= cap {
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 2. Tractions on the flag's wetted samples.
|
||||
let poly_body = EmbeddedBody::polygon(vertices.read().unwrap().clone());
|
||||
let mask = solver.mask().unwrap();
|
||||
state_mask_cells.set(mask.fluid_cells());
|
||||
let body = solver.body().unwrap();
|
||||
let mut faces = Vec::new();
|
||||
let mut tractions = Vec::new();
|
||||
let mut skipped = 0usize;
|
||||
for s in poly_body.surface_samples(0.5 * h) {
|
||||
if circle_sdf(s.x, s.y) < 1e-9 {
|
||||
continue; // buried in the cylinder
|
||||
}
|
||||
match mask.traction_at(
|
||||
body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
|
||||
) {
|
||||
Some((tx, ty)) => {
|
||||
faces.push(FluidFace {
|
||||
centroid: Vector3::new(s.x, s.y, 0.0),
|
||||
normal: Vector3::new(s.nx, s.ny, 0.0),
|
||||
area: s.ds,
|
||||
});
|
||||
tractions.push(Vector3::new(tx, ty, 0.0));
|
||||
}
|
||||
None => skipped += 1,
|
||||
}
|
||||
}
|
||||
*total_skipped.borrow_mut() += skipped;
|
||||
|
||||
// 3. rtx-fsi carries the load to the deformed structure nodes.
|
||||
let nodes_now = interface.deformed_nodes(d);
|
||||
let surface = match WettedSurface::build(&faces, &nodes_now) {
|
||||
Ok(surface) => surface,
|
||||
Err(rtx_fsi::FsiError::DegenerateNeighbourhood { face }) => {
|
||||
let c = faces[face].centroid;
|
||||
eprintln!("degenerate face {face} centroid ({:.6}, {:.6})", c.x, c.y);
|
||||
let mut dists: Vec<(f64, usize)> = nodes_now
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(i, n)| ((n - c).norm(), i))
|
||||
.collect();
|
||||
dists.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
|
||||
for (dist, i) in dists.iter().take(10) {
|
||||
eprintln!(
|
||||
" node {i} at ({:.6}, {:.6}) dist {dist:.6}",
|
||||
nodes_now[*i].x, nodes_now[*i].y
|
||||
);
|
||||
}
|
||||
panic!("degenerate neighbourhood at face {face}");
|
||||
}
|
||||
Err(e) => panic!("transfer build failed: {e:?}"),
|
||||
};
|
||||
let nodal = surface.transfer_load(&faces, &tractions).unwrap();
|
||||
let total_sampled: Vector3<f64> =
|
||||
faces.iter().zip(&tractions).map(|(f, t)| t * f.area).sum();
|
||||
let total_nodal: Vector3<f64> = nodal.iter().sum();
|
||||
let conservation = (total_nodal - total_sampled).norm() / total_sampled.norm().max(1e-30);
|
||||
let mut worst = worst_conservation.borrow_mut();
|
||||
*worst = worst.max(conservation);
|
||||
|
||||
// 4. The flag answers.
|
||||
let (d_new, tip_now, _newton) = solve_flag(&mesh, &interface, &nodal, a_node);
|
||||
*tip.borrow_mut() = tip_now;
|
||||
// Diagnostics: the residual trajectory, the tip, the total load,
|
||||
// and a mask fingerprint (fluid-cell count) to see chatter.
|
||||
{
|
||||
let mut previous = previous_d.borrow_mut();
|
||||
let delta: f64 = d_new
|
||||
.iter()
|
||||
.zip(previous.iter())
|
||||
.map(|(a, b)| (a - b) * (a - b))
|
||||
.sum::<f64>()
|
||||
.sqrt();
|
||||
let fluid_cells = state_mask_cells.get();
|
||||
eprintln!(
|
||||
" pass: |d_new - d_prev| = {delta:.3e}, tip = ({:.4e}, {:.4e}), \
|
||||
total sampled force = ({:.4}, {:.4}), fluid cells = {fluid_cells}, \
|
||||
marched {marched}",
|
||||
tip_now.0, tip_now.1, total_sampled.x, total_sampled.y
|
||||
);
|
||||
*previous = d_new.clone();
|
||||
}
|
||||
d_new
|
||||
};
|
||||
|
||||
// Tolerance from measurement: with the upwind fluid and load-stagnation
|
||||
// passes the interface still carries a ~3e-5 noise floor (each geometry
|
||||
// nudge re-excites a slow settle the stagnation window cuts short), so
|
||||
// the fixed point is determined to about ±2% of the tip — 8e-5 is what
|
||||
// this coupling can honestly promise, and the run terminates as soon as
|
||||
// a pass lands inside that band.
|
||||
let mut coupling = Subiterated::aitken(25, 8e-5).unwrap();
|
||||
let d0 = vec![0.0; 2 * interface.wetted.len()];
|
||||
let converged = coupling
|
||||
.solve(&d0, pass)
|
||||
.expect("coupling did not converge");
|
||||
println!(
|
||||
" coupling: {} Aitken passes, residual {:.2e}; worst conservation defect {:.2e}; \
|
||||
skipped samples total {}",
|
||||
converged.iterations,
|
||||
converged.residual,
|
||||
*worst_conservation.borrow(),
|
||||
*total_skipped.borrow(),
|
||||
);
|
||||
|
||||
// Settle the fluid on the final geometry and measure the total load on
|
||||
// cylinder + flag (both by surface tractions).
|
||||
let (solver, field) = &mut *state.borrow_mut();
|
||||
for _ in 0..2000 {
|
||||
futures::executor::block_on(solver.advance(field, dt)).unwrap();
|
||||
}
|
||||
let mask = solver.mask().unwrap();
|
||||
let body = solver.body().unwrap();
|
||||
let final_vertices = vertices.read().unwrap().clone();
|
||||
let mut drag = 0.0;
|
||||
let mut lift = 0.0;
|
||||
let poly_body = EmbeddedBody::polygon(final_vertices.clone());
|
||||
for s in poly_body.surface_samples(0.5 * h) {
|
||||
if circle_sdf(s.x, s.y) < 1e-9 {
|
||||
continue;
|
||||
}
|
||||
if let Some((tx, ty)) = mask.traction_at(
|
||||
body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
|
||||
) {
|
||||
drag += tx * s.ds;
|
||||
lift += ty * s.ds;
|
||||
}
|
||||
}
|
||||
let circle_body = EmbeddedBody::circle(0.2, 0.2, 0.05);
|
||||
for s in circle_body.surface_samples(0.5 * h) {
|
||||
if polygon_signed_distance(&final_vertices, s.x, s.y) < 1e-9 {
|
||||
continue;
|
||||
}
|
||||
if let Some((tx, ty)) = mask.traction_at(
|
||||
body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
|
||||
) {
|
||||
drag += tx * s.ds;
|
||||
lift += ty * s.ds;
|
||||
}
|
||||
}
|
||||
|
||||
let (ux_a, uy_a) = *tip.borrow();
|
||||
println!(
|
||||
" FSI1 (fluid ny = {ny}, flag 35x2 Quad8): ux(A) = {:.4e} m (ref {REF_UX:.4e}), \
|
||||
uy(A) = {:.4e} m (ref {REF_UY:.4e}), drag {drag:.3} (ref {REF_DRAG}), \
|
||||
lift {lift:.4} (ref {REF_LIFT}); total wall {:.0} s",
|
||||
ux_a,
|
||||
uy_a,
|
||||
start.elapsed().as_secs_f64(),
|
||||
);
|
||||
|
||||
let rel = |a: f64, b: f64| ((a - b) / b).abs();
|
||||
let worst = *worst_conservation.borrow();
|
||||
assert!(worst < 1e-10, "load transfer lost force: {worst:.3e}");
|
||||
assert!(
|
||||
rel(drag, REF_DRAG) < 0.15,
|
||||
"drag {drag:.3} vs reference {REF_DRAG}"
|
||||
);
|
||||
assert!(
|
||||
rel(lift, REF_LIFT) < 0.35,
|
||||
"lift {lift:.4} vs reference {REF_LIFT}"
|
||||
);
|
||||
assert!(
|
||||
rel(ux_a, REF_UX) < 0.30,
|
||||
"ux(A) {ux_a:.4e} vs reference {REF_UX:.4e} (measured +16.6% at ny = 62)"
|
||||
);
|
||||
if ny >= 82 {
|
||||
// Measured +37% at h = 5 mm (1.124e-3): the resolutions BRACKET the
|
||||
// reference — 3.8e-4 (−54%) at 6.6 mm, 1.12e-3 (+37%) at 5 mm —
|
||||
// nonmonotone through the flag's 3 → 4-cell thickness transition,
|
||||
// exactly like the rigid-flag lift. ux converges cleanly (+16.6% →
|
||||
// +6.1%). The band is the measured value, not an accuracy claim.
|
||||
assert!(
|
||||
rel(uy_a, REF_UY) < 0.45,
|
||||
"uy(A) {uy_a:.4e} vs reference {REF_UY:.4e}"
|
||||
);
|
||||
} else {
|
||||
// The measured band of this resolution, not an accuracy claim.
|
||||
assert!(
|
||||
(3.0e-4..5.0e-4).contains(&uy_a),
|
||||
"uy(A) {uy_a:.4e} outside the measured ny = 62 band [3.0e-4, 5.0e-4]"
|
||||
);
|
||||
}
|
||||
assert!(
|
||||
uy_a > 0.0 && ux_a > 0.0,
|
||||
"tip displacement direction wrong: ({ux_a:.3e}, {uy_a:.3e})"
|
||||
);
|
||||
}
|
||||
Reference in New Issue
Block a user