Omar Sobh
f9d38bee55
CI / Build (macos-latest) (push) Waiting to run
CI / Test (macos-latest) (push) Blocked by required conditions
CI / Test (ubuntu-latest) (push) Blocked by required conditions
CI / Python Bindings (maturin) (macos-latest) (push) Blocked by required conditions
CI / Python Bindings (maturin) (ubuntu-latest) (push) Blocked by required conditions
CI / WASM Build + Size Check (push) Blocked by required conditions
CI / Distributed Training Tests (push) Blocked by required conditions
CI / CI Success (push) Blocked by required conditions
CI / Build (ubuntu-latest) (push) Failing after 6s
CI / Format Check (push) Failing after 7s
Documentation / Build API Documentation (push) Failing after 5s
CI / Clippy Check (push) Failing after 7s
Performance Benchmarks / Run Benchmarks (push) Failing after 8s
CI / Build CPU-Only (Explicit) (push) Failing after 5s
Documentation / Build User Guide (push) Successful in 5s
610 lines
23 KiB
Rust
610 lines
23 KiB
Rust
//! 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,
|
||
poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
|
||
poisson_smoother: rtx_cfd::solvers::incompressible::MgSmoother::Lexicographic,
|
||
// 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})"
|
||
);
|
||
}
|