diff --git a/WRITING-PLAN.md b/WRITING-PLAN.md index f23c95f..0559be0 100644 --- a/WRITING-PLAN.md +++ b/WRITING-PLAN.md @@ -227,14 +227,49 @@ it is the one a reader arrives with. On a box, "no flow through this wall" is a component of the velocity and you constrain it. On an annulus, a sphere, a boundary with topography, or any mesh that has been moved, it is not a component of anything — and that is the whole -difficulty. The note is about what you can do instead, and what each choice -costs: - -- **Penalty and Nitsche.** What they enforce, and that they leak — order `1e-3` - where a strong constraint holds to machine precision. -- **Rotating the degrees of freedom.** A per-node rotation `Q` and a strong - `v_n = 0`. Exact, and correct on curved, tilted and deformed boundaries - because the normal is taken per node. +difficulty. **Three approaches, in that order** — each is the previous one's +answer, which is the spine of the note: + +- **Direct penalty.** Add a term to the weak form that punishes `v.n != 0`. + One line, works anywhere, and never quite holds: the leak is set by the + penalty parameter, and driving it down to close the leak conditions the + operator badly. You are trading one error for another. +- **Nitsche.** The consistent version of the same idea. Carrying the boundary + traction terms as well as the penalty makes the discrete problem consistent + for any stabilisation above a threshold, rather than only in the limit — so + it converges at the optimal order without the conditioning price. It is still + a weak imposition and still leaks, order `1e-3` in what we measure. +- **Rotating the degrees of freedom.** Stop asking for the constraint and + impose it: a per-node rotation `Q` and a strong `v_n = 0`. Exact to machine + precision, and correct on curved, tilted and deformed boundaries because the + normal is taken per node. + + **This is the classical answer, not a new one** — it goes back to the early + finite-element texts, and Engelman, Sani & Gresho were already reviewing the + alternatives in 1982. The note should present it as the textbook method + recovered, and then explain why it is nonetheless the least used of the three. + + **The reason is structural, not numerical.** Rotating the degrees of freedom + leaves the discrete vector in a *mixed basis*: interior nodes hold + `(v_x, v_y)`, constrained boundary nodes hold `(v_n, v_t)`, and every piece + of machinery downstream has to know which is which. That is a solver-wide + obligation, and it is where the cost actually lands. Ours, concretely: + + - the multigrid prolongation has to be rotated too, which is why the rotated + path cannot use the DM-coupled hierarchy at all and needs custom-P + transfers; + - the rotated solve builds its own KSP under a per-solve prefix, so + `stokes.petsc_options` does not reach it — a trap that has cost us time + more than once; + - the Schur block and the preconditioner both had to be revisited for the + rotated operator. + + None of that is an argument against the method. It is an argument for + knowing what you are taking on, and it is the honest reason a weakly imposed + condition survives in codes that could do this instead. + +The leak numbers are the argument for the ordering, and they should be measured +in the note rather than asserted. - **Which normal, which is subtler than it looks.** A node-averaged normal weighted by facet measure matches the straight-facet integral the assembler actually evaluates; an analytic normal is exact for the *geometry* and @@ -253,6 +288,111 @@ costs: condition that has to evolve in time — a Dirichlet-to-traction ramp — still wants Nitsche. +**Where this note came from, and it should say so.** The rotated boundary +conditions were not built to tidy up free slip on an annulus — they were built +because the free surface needs an accurate surface traction, and that is the +one quantity a weakly imposed constraint gets wrong. Leading with that gives +the note a reason to exist beyond completeness, and it ties it to S1/S2, which +should be written near it (see the free-surface section). + +**The hard part of this note is that the three usually agree.** Solve a +convection model with any of them and the velocity field is the same to +plotting accuracy; a `1e-3` leak in `v.n` is invisible in anything that +consumes the velocity, which is most of what a model does. A note that compares +three methods on a problem where they agree has no argument, and a reader who +suspects the comparison was staged is right to. + +So the note is organised around **where the difference is actually visible**, +and the clearest case is **surface stress**. When the wall-normal traction is +the answer rather than a by-product — dynamic topography, plate-boundary force +balance, anything compared against a geoid or a gravity field — the three stop +agreeing: + +- Under a penalty or Nitsche condition the constraint is approximate, so the + traction recovered from it inherits the approximation. You are differentiating + a field that was never made to satisfy the condition exactly. +- Under the rotated constraint the reaction **is** `sigma_nn`. It is not + recovered, post-processed or split off — it is the multiplier the solve + already computed, available through `boundary_normal_traction` / + `dynamic_topography`. + +That contrast is the note's worked example and it should be measured, not +described: same model, three boundary treatments, compare the surface traction +against a case with a known answer. + +Secondary discriminators, worth a paragraph each rather than a section: a +boundary that is genuinely curved or has been deformed, where the leak is not +merely small but geometrically inconsistent; composition with transverse +isotropy, which the rotated constraint survives and Nitsche does not; and +conditioning as the penalty parameter is driven down. + +And say plainly, early, that for a model which only consumes the velocity +field, the simplest thing that works is the right choice. The note is more +useful if it tells the reader when they can stop reading. + +**This one needs the mathematics written out**, unlike the measurement-led +notes. The three approaches differ in their weak forms, and the differences are +the argument — a reader cannot be asked to take "consistent only in the limit" +on trust. What has to appear: + +- The Stokes weak form and where the boundary term `int_G (sigma.n).w` comes + from, because every method below is a statement about that term. +- That free slip is *two* conditions — `v.n = 0` and zero tangential traction — + and the second is natural, which is why it is the one people forget. +- **Direct penalty**: add `(gamma/h) int_G (v.n)(w.n)`. The discrete problem is + a perturbed problem, and the perturbation is what the leak is. +- **Nitsche**: the consistency term and the adjoint-consistency term alongside + the penalty, and that `gamma` has a threshold set by an inverse inequality + rather than being a free dial. This is the part that most needs writing out, + because "add two more terms and it becomes consistent" is not believable + without seeing them. +- **Rotated**: the per-node `Q`, solving in `(v_n, v_t)`, constraining `v_n` + strongly, and the reaction falling out as `sigma_nn`. +- **The normal on a faceted boundary**: the assembled constraint is an integral + over straight facets, so the node normal consistent with it is the one + weighted by facet measure. This is where our own #560 landed, and it is worth + deriving rather than asserting. + +**Starting bibliography.** Verified references, not a reading list yet: + +- Engelman, Sani & Gresho, *The implementation of normal and/or tangential + boundary conditions in finite element codes for incompressible fluid flow*, + Int. J. Numer. Methods Fluids **2** (1982) 225-238. The classic statement of + the rotated-degrees-of-freedom approach; reviews the alternatives and uses + global mass conservation to choose between them. Our rotated BCs are this + idea, and the note should say so rather than presenting it as new. +- Behr, *On the application of slip boundary condition on curved boundaries*, + Int. J. Numer. Methods Fluids **45** (2004) 43-51. Directly the "which + normal" question on a discretised curved boundary. ⚠️ Bibliographic details + confirmed, contents NOT yet read — whether it reaches the same + measure-weighted normal we did is exactly what to check, and if it does, #560 + was a rediscovery and should be described as one. +- Nitsche, *Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei + Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind*, + Abh. Math. Semin. Univ. Hamburg **36** (1971) 9-15, + `10.1007/BF02995904`. The original. + +Still to find: a modern treatment of Nitsche for *slip* specifically (as +opposed to no-slip), and whatever the geodynamics codes cite for free slip on a +spherical shell. + +**Formulating in another coordinate system is the same idea, globally.** +Solving in spherical or cylindrical components makes the wall-normal direction +a coordinate direction again, so the constraint returns to being "hold one +component" — which is rotating the degrees of freedom, imposed once for the +whole domain instead of node by node. Worth saying explicitly, because it +explains why it is attractive: applied globally there is no mixed basis and +none of the structural cost above. + +It does not generalise, and that is the whole point. It works exactly when the +boundary lies along a coordinate surface — a sphere, an annulus, a cylinder — +and does nothing for topography, a deformed mesh, or a tilted internal surface. +The per-node rotation is what you are left with once the geometry stops +cooperating, and paying its structural price is what buys the generality. + +Treat it as the third approach's special case rather than a fourth approach, +and do not develop it further than that. + Curved boundaries under *refinement* are G1's, not this note's: the snapping callback that keeps a refined boundary on the true surface is already written up there, and this note should link rather than repeat it. @@ -322,7 +462,21 @@ should be acknowledged as such. ## Free surface Its own development, and it needs discussion before it is written. Two notes, -and the split matters: +and the split matters. + +**Write these near R1.** The free surface is *why* the rotated boundary +conditions were implemented: the surface evolves under the traction it carries, +so the wall-normal stress stops being a diagnostic and becomes the thing that +drives the model. That is the strongest possible case for a constraint whose +reaction is `sigma_nn` exactly rather than recovered from an approximately +satisfied condition, and it is the motivation R1 should lead with rather than +arriving at. + +The dependency runs one way — R1 is the machinery, S1 and S2 are what it was +built for — so R1 either goes first or they go out together. **M1 and G1 are +the cautionary example**: they were meant to publish together, went four days +apart, and now owe a v2 for a cross-link that could have been in v1. Decide +which of the two patterns this pair follows before drafting, not after. ### S1. The algorithm