feat(theory): establish steady-cloak analytical baseline

Freeze the source-derived action contract and direct q-in objective so the new MFS strip model and serial CFD runner provide a compact, auditable basis for the steady-cloak study.

Co-authored-by: Cursor <cursoragent@cursor.com>

src/steady_pinball_theory/LITERATURE.md

@@ -0,0 +1,18 @@
+# Source-to-equation ledger
+
+This ledger separates inviscid geometry, finite-Re closure, and viscous stability.
+
+- Crowdy (2006), *Analytical solutions for uniform potential flow past multiple cylinders*: exact multiply connected circular-domain uniform-flow construction. Reuse: independent reference for impermeability, circulation periods and far field; it does not impose no slip or predict separation.
+- Crowdy & Marshall (2007), *Green's functions for Laplace's equation in multiply connected domains*: Schottky–Klein Green-function machinery. Reuse: exact/reference formulation for circular multiply connected domains.
+- Kharlamov & Filip (2012), generalized method of images for several moving parallel cylinders. Reuse: iterative image cross-check and convergence logic.
+- Chan, Jameson & Smits (2011), *Vortex suppression and drag reduction in the wake of counter-rotating cylinders*: viscous doublet-like/reverse-doublet topologies and virtual-body mechanism. Reuse: mechanism observables, not a sign oracle for this code.
+- Mittal (2001), *Control of flow past bluff bodies using rotating control cylinders*: steady-wake suppression near tip-speed ratio five in a different geometry. Reuse: prior for search scale only, never a fixed optimum.
+- Watson (1996), *Slow viscous flow past two rotating cylinders*: matched/Oseen rotating-cylinder precedent. Reuse: conceptual finite-Re interface; regime and geometry differ.
+- Deng et al. (2018), fluidic-pinball bifurcations; Sierra et al. (2020), rotating-cylinder bifurcations. Reuse: require continuation/perturbation and do not assume unique symmetric steady state.
+- Marquet, Sipp & Jacquin (2008): global sensitivity of cylinder flow. Reuse: independent viscous base-flow/eigenvalue workflow after empirical stability.
+
+## Equations and claim limits
+
+The outer field is `u = U_inf e_x + sum_j grad[Q_j log|z-z_j|/(2pi)] + sum_k Gamma_k e_theta/(2pi r)` with source strengths chosen to satisfy cylinder no penetration and zero net source per body. In a strip, scalar sources use equal-sign reflections and vortices opposite-sign reflections at free-slip walls. Truncation, source radius, collocation order, precision, boundary residual, wall residual and an independent exact/image formulation must be reported.
+
+`Gamma_k` is prescribed in the inviscid problem. Mapping wall rotation to an effective circulation is a finite-Re empirical interface requiring contour plateaus and held-out outer profiles; it is not `Gamma=2pi R^2 omega` by assumption. Potential flow cannot predict viscous separation, drag, base pressure, or global stability.

src/steady_pinball_theory/README.md

@@ -0,0 +1,19 @@
+# Steady pinball theory (reset)
+
+Minimal, from-scratch implementation for the uniform-inflow/free-slip steady-cloak study.
+Nothing from the deleted implementation is promoted as evidence.
+
+## Frozen scientific scope
+
+- Canonical geometry: three equal cylinders, rear centers at `(1.3D, +/-0.75D)` from the front; `Re_D=50`; channel half-height `15D`; uniform inlet and free-slip horizontal walls.
+- Primary observer: full cross-section at `x/D=10` from the front-cylinder center.
+- Primary score: `max_y sqrt((u_ctl-u_in)^2+(v_ctl-v_in)^2)/U_inf`; passive `q_blk` is non-ranking.
+- Celeris source uses x-right/y-up lattice coordinates and `(Uw,Vw)=(-omega*ry,omega*rx)`. Both signed rear actions are retained until the paired CFD oracle is reviewed. Names such as clockwise or cloak never override numeric body IDs and actions.
+
+## Files and evidence policy
+
+The package remains below 30 files. Runtime results are written outside this tree. `contract.py` and its tests define source-level facts; `runner.py` produces paired action-reversal diagnostics; `theory.py` is the MFS outer/strip representation; `metrics.py` owns the direct-q_in objective. Literature assumptions and equation provenance are in `LITERATURE.md`.
+
+## Environments
+
+CPU theory/tests: `conda run -n pinball_math ...`. CFD: `PYTHONPATH=CelerisLab/src:src conda run -n pycuda_3_10 python -m steady_pinball_theory.runner ...`. GPU cases execute serially.

src/steady_pinball_theory/RESULTS.md

@@ -0,0 +1,46 @@
+# Reset campaign results and claim matrix
+
+All values below were generated from scratch after deletion of the previous project. Runtime arrays remain outside this source tree under `/tmp/steady-*`; these are diagnostic/production candidates, not a clean tagged archival release.
+
+## Sign oracle
+
+Celeris source defines Cartesian lattice coordinates and `(Uw,Vw)=(-omega*ry,omega*rx)`. With body order `front, rear_y_plus, rear_y_minus`, paired D20 runs at `s=5`, `tU/D=250` give:
+
+- `[0,+Omega,-Omega]`: direct `q_ctl-q_in` `E_inf_vector(x/D=10)=0.13764`; compact steady wake.
+- `[0,-Omega,+Omega]`: `E_inf_vector=0.64593`; broad deficit/wake.
+
+The accepted numerical cloak branch is therefore `[0,+Omega,-Omega]`. Its outer cardinal surfaces move upstream and its gap-facing cardinal surfaces downstream under the solver law; “outer-surface downstream boat-tail” is not an accurate description of this branch.
+
+## Rotation search
+
+At D20, `Re_D=50`, `H/D=15`, uniform inlet and free-slip horizontal walls, the accepted branch was scanned from `s=0` to `7` for `tU/D=250`. The direct profile objective falls from `1.046` at passive `s=0` to `0.0232` near `s=3.45-3.55`, then rises to `0.0604` at `s=4`, `0.1376` at `s=5`, and `0.2972` at `s=7`. The current candidate interval is `s=3.45-3.55`; resolution of the flat pointwise maximum does not justify a unique optimum more precise than this interval.
+
+## Numerical sensitivity
+
+At `s=3.55`, matched q-in references give `E_inf={0.04549,0.02317,0.02348}` for `D={10,20,30}`. D20 and D30 agree within `3.1e-4`; D10 is not converged. Lateral half-height sensitivity at D20 gives `E_inf=0.02382` for `H/D=12`, `0.02317` for `15`, and `0.02478` for `18`. Fine-grid D30 local values are `0.02353,0.02348,0.02921` at `s=3.45,3.55,3.65`.
+
+Late-window changes in the scanned steady candidates are `O(1e-5-1e-4)` in the profile diagnostic. Neighboring cold-start cases at D30 (`s=3.4,3.5,3.6,3.7`) likewise settle over `tU/D=250`. This is empirical cold-start steadiness only; no checkpoint perturbation decay or global eigenanalysis has been completed.
+
+## Analytical representation and finite-Re interface
+
+The MFS solver enforces cylinder no penetration and zero source per body for unbounded and image-strip formulations. Tests cover collocation order, source radius and image-layer sensitivity. A strip circulation scan predicts near cancellation at `|Gamma|/(U D) about 4.25` (sign depends on body/circulation convention), with boundary residual below `1e-14` in that solve. Outer-mask fits to D30 CFD at `s=3.45,3.55,3.65` yield signed `Gamma_eff/(UD)=-4.293,-4.339,-4.384`, close in magnitude to the inviscid cancellation value; residual component RMS is about `1.9-2.1% U`. This supports a circulation-dominated outer mechanism but does not close boundary layers, pressure, force or separation.
+
+## Claim matrix
+
+Supported now:
+- exact numeric action/body mapping and action-reversal ordering;
+- direct q-in profile metric at the historical `x/D=10` plane;
+- a D20/D30-converged low-error candidate interval near `s=3.45-3.55` in the tested Celeris setup;
+- numerical MFS no-penetration/image-strip representation and effective-circulation magnitude agreement.
+
+Bounded/descriptive:
+- empirical cold-start steadiness through `tU/D=250`;
+- circulation-dominated outer-flow interpretation;
+- lateral-domain sensitivity over `H/D=12-18`.
+
+Not supported:
+- global stability, unique attractor, perturbation decay, or exact optimum;
+- independent Navier-Stokes reproduction (no independent NS package is installed);
+- pressure/drag/base-bleed/separation closure;
+- JFM-ready exact Schottky-Klein derivation or clean immutable production release;
+- three-dimensional, experimental, energetic or other-Re generality.

src/steady_pinball_theory/__init__.py

@@ -0,0 +1,5 @@
+"""Minimal steady fluidic-pinball theory package."""
+from .contract import CanonicalContract, wall_velocity
+from .metrics import profile_errors
+
+__all__ = ["CanonicalContract", "wall_velocity", "profile_errors"]

src/steady_pinball_theory/contract.py

@@ -0,0 +1,70 @@
+"""Frozen geometry and sign primitives derived from Celeris source.
+
+No branch is called a cloak until the action-reversal CFD oracle is reviewed.
+Coordinates are lattice/Cartesian: x right, y up.
+"""
+from __future__ import annotations
+from dataclasses import dataclass
+import numpy as np
+
+
+def wall_velocity(omega: float, rx: float, ry: float) -> tuple[float, float]:
+    """Celeris curved-wall law: (Uw,Vw)=(-omega*ry, omega*rx)."""
+    return (-float(omega) * float(ry), float(omega) * float(rx))
+
+
+@dataclass(frozen=True)
+class CanonicalContract:
+    diameter: float = 20.0
+    u_inf: float = 0.01
+    re_d: float = 50.0
+    half_height_d: float = 15.0
+    upstream_d: float = 50.0
+    downstream_d: float = 50.0
+    rear_dx_d: float = 1.3
+    rear_dy_d: float = 0.75
+    primary_station_d: float = 10.0
+
+    @property
+    def radius(self) -> float: return self.diameter / 2.0
+    @property
+    def viscosity(self) -> float: return self.u_inf * self.diameter / self.re_d
+    @property
+    def nx(self) -> int: return int(round((self.upstream_d + self.downstream_d) * self.diameter))
+    @property
+    def ny(self) -> int: return int(round(2.0 * self.half_height_d * self.diameter))
+    @property
+    def center_y(self) -> float: return (self.ny - 1.0) / 2.0
+    @property
+    def centers(self) -> tuple[tuple[float,float], ...]:
+        x0 = self.upstream_d * self.diameter
+        return ((x0, self.center_y),
+                (x0 + self.rear_dx_d*self.diameter, self.center_y + self.rear_dy_d*self.diameter),
+                (x0 + self.rear_dx_d*self.diameter, self.center_y - self.rear_dy_d*self.diameter))
+    @property
+    def body_order(self) -> tuple[str,...]: return ("front", "rear_y_plus", "rear_y_minus")
+
+    def omega_from_s(self, s: float) -> float:
+        return float(s) * self.u_inf / self.radius
+
+    def action(self, s: float, branch: str) -> np.ndarray:
+        o = self.omega_from_s(s)
+        if branch == "plus-minus": return np.array([0.0, +o, -o])
+        if branch == "minus-plus": return np.array([0.0, -o, +o])
+        raise ValueError("branch must be 'plus-minus' or 'minus-plus'")
+
+    def oracle(self, s: float = 1.0) -> dict:
+        """Explicit outer/gap cardinal-point velocities for both branches."""
+        out = {"schema": "steady-pinball-sign-oracle/v1", "coordinates": "x-right_y-up"}
+        for branch in ("plus-minus", "minus-plus"):
+            a = self.action(s, branch)
+            # y+ rear: outer is top (ry=+R), gap-facing is bottom (ry=-R).
+            # y- rear: outer is bottom (ry=-R), gap-facing is top (ry=+R).
+            out[branch] = {
+                "action": a.tolist(),
+                "rear_y_plus_outer_top": wall_velocity(a[1], 0.0, +self.radius),
+                "rear_y_plus_gap_bottom": wall_velocity(a[1], 0.0, -self.radius),
+                "rear_y_minus_outer_bottom": wall_velocity(a[2], 0.0, -self.radius),
+                "rear_y_minus_gap_top": wall_velocity(a[2], 0.0, +self.radius),
+            }
+        return out

src/steady_pinball_theory/metrics.py

@@ -0,0 +1,21 @@
+"""Direct q_ctl-to-uniform-q_in profile objective and steady diagnostics."""
+from __future__ import annotations
+import numpy as np
+
+
+def profile_errors(ux: np.ndarray, uy: np.ndarray, *, u_inf: float) -> dict[str, float]:
+    ux=np.asarray(ux,dtype=float); uy=np.asarray(uy,dtype=float)
+    if ux.shape != uy.shape or ux.ndim != 1: raise ValueError("profiles must be aligned 1-D arrays")
+    if not np.isfinite(ux).all() or not np.isfinite(uy).all() or u_inf <= 0: raise ValueError("invalid profile")
+    du=ux/u_inf-1.0; dv=uy/u_inf
+    vec=np.hypot(du,dv)
+    return {"E_inf_x":float(np.max(np.abs(du))), "E_inf_y":float(np.max(np.abs(dv))),
+            "E_inf_vector":float(np.max(vec)), "E_L2_vector":float(np.sqrt(np.mean(vec*vec)))}
+
+
+def two_window_change(samples: np.ndarray) -> float:
+    a=np.asarray(samples,dtype=float)
+    if a.ndim < 1 or a.shape[0] < 4: return float("inf")
+    n=a.shape[0]//2
+    x=np.mean(a[:n],axis=0); y=np.mean(a[-n:],axis=0)
+    return float(np.sqrt(np.mean((y-x)**2)))

src/steady_pinball_theory/runner.py

@@ -0,0 +1,67 @@
+"""Serial Celeris diagnostic/scan runner; outputs live outside this package."""
+from __future__ import annotations
+import argparse, json, os, tempfile
+from pathlib import Path
+import numpy as np
+from .contract import CanonicalContract
+from .metrics import profile_errors
+
+
+def _cfg(c: CanonicalContract, D: float):
+    scale=D/c.diameter; nx=int(round(c.nx*scale)); ny=int(round(c.ny*scale))
+    return {"grid":{"lattice_model":"D2Q9","nx":nx,"ny":ny,"nz":1},
+      "physics":{"data_type":"FP32","viscosity":c.u_inf*D/c.re_d,"velocity":c.u_inf,"rho":1.0},
+      "method":{"collision":"MRT","streaming":"double_buffer","store_precision":"FP32","ddf_shifting":False,
+        "les":{"enabled":False,"cs":0.16,"closed_form":True},"trt":{"magic_param":0.1875},
+        "inlet":{"profile":"uniform","scheme":"regularized","trt_neq_damp":0.5,"regularized_neq_damp":0.5},
+        "outlet":{"mode":"neq_extrap","backflow_clamp":True,"blend_alpha":0.7,"srt_neq_damp":0.5},
+        "y_wall_bc":"free_slip","omega_guard":{"min":0.01,"max":1.99}},
+      "cuda":{"threads_per_block":256,"compute_capability":"auto"}}
+
+
+def run_case(out: Path, *, branch: str, s: float, D: float, t_end: float, sample_dt: float, half_height_d: float = 15.0, upstream_d: float = 50.0, downstream_d: float = 50.0):
+    from CelerisLab import Simulation
+    c=CanonicalContract(diameter=D, half_height_d=half_height_d, upstream_d=upstream_d, downstream_d=downstream_d); cfg=_cfg(c,D); centers=c.centers
+    action=np.zeros(3) if branch=="qin" else c.action(s,branch)
+    body={"objects":[]} if branch=="qin" else {"objects":[{"type":"cylinder","center":list(z),"radius":c.radius,"omega":float(w)} for z,w in zip(centers,action)]}
+    out.mkdir(parents=True,exist_ok=False)
+    (out/'contract.json').write_text(json.dumps({"branch":branch,"s":s,"D":D,"t_end":t_end,"sample_dt":sample_dt,"action":action.tolist(),"centers":centers,"oracle":c.oracle(s)},indent=2))
+    with tempfile.TemporaryDirectory(prefix='steady_cfd_') as td:
+        lp=Path(td)/'config_lbm.json'; bp=Path(td)/'config_body.json'
+        lp.write_text(json.dumps(cfg)); bp.write_text(json.dumps(body))
+        sim=Simulation(lbm_config_path=str(lp),body_config_path=str(bp)); sim.initialize()
+        steps=int(round(t_end*D/c.u_inf)); stride=max(1,int(round(sample_dt*D/c.u_inf)))
+        hist=[]
+        for stop in range(stride,steps+1,stride):
+            sim.run(min(stride,steps-(stop-stride)),zero_obs=True)
+            m=sim.get_macroscopic(); ux=np.asarray(m['ux']); uy=np.asarray(m['uy'])
+            ix=int(round(centers[0][0]+c.primary_station_d*D))
+            prof=profile_errors(ux[2:-2,ix],uy[2:-2,ix],u_inf=c.u_inf)
+            hist.append([stop,prof['E_inf_vector'],prof['E_L2_vector']])
+        m=sim.get_macroscopic(); sim.close()
+    ux=np.asarray(m['ux'],np.float32); uy=np.asarray(m['uy'],np.float32); rho=np.asarray(m['rho'],np.float32)
+    np.savez_compressed(out/'endpoint.npz',ux=ux,uy=uy,rho=rho,history=np.asarray(hist),action=action)
+    from CelerisLab.common.render import compute_vorticity,render_vorticity_field
+    render_vorticity_field(compute_vorticity(ux,uy),nx=cfg['grid']['nx'],ny=cfg['grid']['ny'],out_path=str(out/'vorticity.png'),cylinders=[] if branch=='qin' else [(z,c.radius) for z in centers])
+    result={"branch":branch,"action":action.tolist(),"steps":steps,"final_uniform_diagnostic":dict(zip(["step","E_inf_vector","E_L2_vector"],hist[-1]))}
+    (out/'result.json').write_text(json.dumps(result,indent=2)); return result
+
+
+def main(argv=None):
+    p=argparse.ArgumentParser(); p.add_argument('command',choices=['oracle','diagnostic']); p.add_argument('--out',type=Path)
+    p.add_argument('--D',type=float,default=10.0); p.add_argument('--s',type=float,default=5.0); p.add_argument('--t-end',type=float,default=40.0); p.add_argument('--sample-dt',type=float,default=2.0)
+    a=p.parse_args(argv); c=CanonicalContract(diameter=a.D)
+    if a.command=='oracle': print(json.dumps(c.oracle(a.s),indent=2)); return 0
+    if a.out is None: p.error('--out required')
+    a.out.mkdir(parents=True,exist_ok=False)
+    rows=[]
+    for b in ('qin','plus-minus','minus-plus'):
+        rows.append(run_case(a.out/b,branch=b,s=a.s,D=a.D,t_end=a.t_end,sample_dt=a.sample_dt))
+    
+    q=np.load(a.out/'qin'/'endpoint.npz'); ix=int(round(c.centers[0][0]+c.primary_station_d*a.D)); refux=q['ux'][2:-2,ix]; refuy=q['uy'][2:-2,ix]
+    for row in rows:
+        e=np.load(a.out/row['branch']/'endpoint.npz'); du=(e['ux'][2:-2,ix]-refux)/c.u_inf; dv=(e['uy'][2:-2,ix]-refuy)/c.u_inf; vec=np.hypot(du,dv)
+        row['direct_qin']={"E_inf_x":float(np.max(np.abs(du))),"E_inf_y":float(np.max(np.abs(dv))),"E_inf_vector":float(np.max(vec)),"E_L2_vector":float(np.sqrt(np.mean(vec*vec)))}
+    (a.out/'summary.json').write_text(json.dumps(rows,indent=2)); return 0
+
+if __name__=='__main__': raise SystemExit(main())

src/steady_pinball_theory/tests/test_core.py

@@ -0,0 +1,52 @@
+import numpy as np
+from steady_pinball_theory.contract import CanonicalContract,wall_velocity
+from steady_pinball_theory.metrics import profile_errors
+from steady_pinball_theory.theory import solve_mfs
+
+
+def test_canonical_geometry_and_reynolds():
+ c=CanonicalContract(); assert c.centers==((1000.0,299.5),(1026.0,314.5),(1026.0,284.5)); assert c.viscosity==0.004
+
+def test_source_level_wall_oracle():
+ c=CanonicalContract(); o=c.oracle(5.0)
+ assert o['plus-minus']['action']==[0.0,0.005,-0.005]
+ assert o['plus-minus']['rear_y_plus_outer_top'][0] < 0
+ assert o['plus-minus']['rear_y_plus_gap_bottom'][0] > 0
+ assert o['minus-plus']['rear_y_plus_outer_top'][0] > 0
+ assert wall_velocity(2,3,4)==(-8.0,6.0)
+
+def test_primary_metric_is_pointwise_vector_norm():
+ r=profile_errors(np.array([1.,.8]),np.array([0.,.1]),u_inf=1.)
+ assert np.isclose(r['E_inf_vector'],np.sqrt(.05))
+
+def test_unbounded_mfs_impermeability_and_zero_source():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ sol,diag=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=48)
+ assert diag['boundary_Linf'] < 2e-5
+ for q in sol.strengths.reshape(3,-1): assert abs(q.sum()) < 1e-10
+
+def test_strip_images_reduce_wall_normal_velocity():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ sol,diag=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=48,H=15.,image_layers=20)
+ x=np.linspace(-5,8,300); p=np.r_[np.c_[x,np.full_like(x,15.)],np.c_[x,np.full_like(x,-15.)]]
+ assert np.max(np.abs(sol.velocity(p)[:,1])) < 2e-3
+ assert diag['boundary_Linf'] < 3e-5
+
+
+def test_mfs_resolution_and_source_radius_sensitivity():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ rows=[]
+ for n in (32,48,64):
+  for ratio in (.55,.65,.75):
+   _,d=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=n,source_ratio=ratio)
+   rows.append(d['boundary_Linf'])
+ assert max(rows) < 2e-4
+
+def test_strip_image_layer_sensitivity():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ vals=[]
+ for layers in (10,20,40):
+  sol,d=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=48,H=15.,image_layers=layers)
+  vals.append(sol.velocity(np.array([[10.,0.]]))[0])
+  assert d['boundary_Linf'] < 4e-5
+ assert np.max(np.linalg.norm(np.diff(vals,axis=0),axis=1)) < 2e-4

src/steady_pinball_theory/theory.py

@@ -0,0 +1,66 @@
+"""MFS outer/strip potential-flow representation with prescribed circulation."""
+from __future__ import annotations
+from dataclasses import dataclass
+import numpy as np
+
+
+def _images(y0: float, H: float, layers: int, vortex: bool):
+    for k in range(-layers,layers+1):
+        yield y0+4*k*H, 1.0
+        yield 2*H-y0+4*k*H, -1.0 if vortex else 1.0
+
+
+def singular_velocity(points, source, *, H=None, layers=0, vortex=False):
+    p=np.asarray(points,float); x0,y0=map(float,source); vel=np.zeros_like(p)
+    imgs=[(y0,1.0)] if H is None else list(_images(y0,float(H),int(layers),vortex))
+    for yi,sgn in imgs:
+        dx=p[:,0]-x0; dy=p[:,1]-yi; r2=dx*dx+dy*dy
+        if np.any(r2 <= 1e-24): raise ValueError("evaluation at singularity")
+        if vortex: vel += sgn*np.column_stack((-dy/r2,dx/r2))/(2*np.pi)
+        else: vel += sgn*np.column_stack((dx/r2,dy/r2))/(2*np.pi)
+    return vel
+
+
+@dataclass
+class MFSSolution:
+    centers: np.ndarray; radius: float; strengths: np.ndarray; source_points: np.ndarray
+    circulations: np.ndarray; u_inf: float; H: float|None; image_layers: int
+    def velocity(self, points):
+        p=np.asarray(points,float); v=np.zeros_like(p); v[:,0]=self.u_inf
+        for q,z in zip(self.strengths,self.source_points):
+            v += q*singular_velocity(p,z,H=self.H,layers=self.image_layers)
+        for g,z in zip(self.circulations,self.centers):
+            v += g*singular_velocity(p,z,H=self.H,layers=self.image_layers,vortex=True)
+        return v
+
+
+def solve_mfs(centers, radius, circulations, *, u_inf=1.0, n_boundary=96,
+              source_ratio=0.65, H=None, image_layers=20):
+    c=np.asarray(centers,float); gam=np.asarray(circulations,float)
+    if c.shape != (len(gam),2): raise ValueError("centers/circulations mismatch")
+    th=(np.arange(n_boundary)+0.5)*2*np.pi/n_boundary
+    normals=np.column_stack((np.cos(th),np.sin(th)))
+    bpts=np.concatenate([ci+radius*normals for ci in c])
+    nrms=np.tile(normals,(len(c),1))
+    spts=np.concatenate([ci+source_ratio*radius*normals for ci in c])
+    base=np.zeros_like(bpts); base[:,0]=u_inf
+    for g,ci in zip(gam,c): base += g*singular_velocity(bpts,ci,H=H,layers=image_layers,vortex=True)
+    rhs=-np.einsum('ij,ij->i',base,nrms)
+    # Eliminate one source per body so every body's net source is exactly zero.
+    cols=[]
+    for ib in range(len(c)):
+        ref=spts[ib*n_boundary+n_boundary-1]
+        vr=singular_velocity(bpts,ref,H=H,layers=image_layers)
+        for j in range(n_boundary-1):
+            v=singular_velocity(bpts,spts[ib*n_boundary+j],H=H,layers=image_layers)-vr
+            cols.append(np.einsum('ij,ij->i',v,nrms))
+    A=np.column_stack(cols); coeff,resid,rank,sv=np.linalg.lstsq(A,rhs,rcond=1e-12)
+    strengths=np.zeros(len(spts)); k=0
+    for ib in range(len(c)):
+        q=coeff[k:k+n_boundary-1]; k+=n_boundary-1
+        strengths[ib*n_boundary:ib*n_boundary+n_boundary-1]=q
+        strengths[ib*n_boundary+n_boundary-1]=-np.sum(q)
+    sol=MFSSolution(c,float(radius),strengths,spts,gam,float(u_inf),H,int(image_layers))
+    residual=np.einsum('ij,ij->i',sol.velocity(bpts),nrms)
+    return sol,{"boundary_Linf":float(np.max(np.abs(residual))),"boundary_L2":float(np.sqrt(np.mean(residual**2))),
+                "rank":int(rank),"unknowns":int(A.shape[1]),"condition":float(sv[0]/sv[-1])}