Functions
	couette (harmonic, rot=None, chi=1.0, dt=20.0, N=32, steps=4000, tol=1e-9)

Variables
	e

	c

	e2

	c2

	True

	dt

	steps

	ea

	ca

	ef

	cf

	rot

Detailed Description

Two-layer plane Couette — validates variable viscosity (Phase 4) INCLUDING the
incremental-rotational pressure scheme (large-dt / steady-Stokes capability) at strong contrast.

Bottom wall fixed, top wall moving at U; two viscosity layers (mu1 lower, mu2 upper, 10x jump).
Analytic steady profile: piecewise linear with interface velocity u_i/U = mu2/(mu1+mu2). The
harmonic face-viscosity mean reproduces it to ~1e-5; the arithmetic mean shows O(1%) error at the
jump (why 'harmonic' matters).

Rotational pressure under variable viscosity: the constant-mu Timmermans term -mu*div(u*) is only
valid for HOMOGENEOUS viscosity (Deteix & Yakoubi, Appl. Math. Lett. 79 (2018) 111-117; arXiv:
1902.05643 for the full shear-rate-projection). The default here ('min') uses the constant
coefficient chi*mu_min - stable at any contrast by domination, exact fallback to the constant-mu
scheme for uniform mu, and it KEEPS the incremental predictor: dt=100 converges in ~200 steps below.
'full' (pointwise mu(i)) is demonstrated to diverge at 10x contrast - mild-contrast use only.

Run:  PYTHONPATH=<build> python two_layer_couette.py