diff --git a/.gitignore b/.gitignore index 461c280..1a0f02b 100644 --- a/.gitignore +++ b/.gitignore @@ -54,3 +54,8 @@ articles/*/*.pdf # macOS .DS_Store + +# Underworld writes a mesh cache beside the script that asked for a mesh. It is +# regenerated by running the example, and one note's examples committed 650 kB +# of it before this line existed. +.meshes/ diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 1465a59..fb9d833 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -167,6 +167,27 @@ something works, not how hard it was to make work. repaired afterwards. - **Say which quantity, on which mesh, against what.** "Converges" and "agrees well" are not measurements. +- **Write display maths as a `math` directive, not as `$$`.** Both build, but + the fenced form survives an editor that is not a MyST editor — Typora, which + several of these notes are written in, rewrites a `$$ ... $$ (label)` block + and moves the label onto a line of its own, which silently breaks it. A fence + is left alone. + + ```` + ```{math} + :label: eq-free-slip + \mathbf{u}\cdot\hat{\mathbf{n}} = 0 + ``` + ```` + + With `$$`, two things are then load-bearing and neither is reported by the + build. The label belongs on the CLOSING delimiter — `$$ (eq-free-slip)` — and + on its own line it renders as literal text while every `{eq}` reference to it + fails the Typst build. And a `$$` block that opens on the line straight after + a paragraph is a lazy continuation: fine for one-line arithmetic, but a line + inside it beginning `-`, `*`, `>` or `#` is claimed by markdown, the maths + never closes, and the page publishes raw LaTeX with half an equation in a + bullet. `pixi run test-unit` guards both. ## Things that are load-bearing diff --git a/WRITING-PLAN.md b/WRITING-PLAN.md index f23c95f..1c32c3b 100644 --- a/WRITING-PLAN.md +++ b/WRITING-PLAN.md @@ -220,9 +220,9 @@ beginning of a solvers paper, but neither is being written to fit one. ### R1. Boundary conditions on non-planar boundaries -Status: not started. **One post**, not a series. Rescoped 2026-08-17: rotated -boundary conditions are the answer, but the *question* is the better frame, and -it is the one a reader arrives with. +Status: drafted as UWTN 2026-016, in review on PR #26. **One post**, not a +series. Rescoped 2026-08-17: rotated boundary conditions are the answer, but the +*question* is the better frame, and 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 @@ -451,6 +451,6 @@ Listed as candidates, not commitments. | 6 | F5 The comparison | Last, with Thyagarajulu's benchmarks as its evidence | | 7 | S1, S2 free surface | After the discussion about how to split it | | — | C1 Launching from any repository | Standalone; the capability is live and undocumented | -| — | R1 BCs on non-planar boundaries | Standalone; write whenever it suits. Much of the evidence exists — see the rescoped entry | +| — | R1 BCs on non-planar boundaries | Drafted as UWTN 2026-016, in review on PR #26 | | ✓ | G1 Setting up FMG | Published 2026-08-17 (UWTN 2026-014) | | 8 | #1 Release announcement | Written last; links to everything | diff --git a/articles/boundary-conditions-on-non-planar-boundaries/boundary-conditions-on-non-planar-boundaries.md b/articles/boundary-conditions-on-non-planar-boundaries/boundary-conditions-on-non-planar-boundaries.md new file mode 100644 index 0000000..bbc4c8f --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/boundary-conditions-on-non-planar-boundaries.md @@ -0,0 +1,733 @@ +--- +title: Free-Slip boundary conditions on curved boundaries +description: >- + "No flow through this wall" is a single velocity component on a box and is + not a component of anything on a sphere, a deformed mesh, or a surface with + topography. Three ways to impose it — a direct penalty, Nitsche, and rotating + the degrees of freedom — what each costs, and the one measurement that tells + them apart. +date: 2026-08-18 +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliations: + - Australian National University +license: CC-BY-4.0 +banner: figures/banner.jpg +keywords: + - Underworld Code + - Tricks of the Trade + - development +exports: + - format: typst + logo: ../../static/uwtn-logo.png + series: "Underworld Technical Notes" + origin_url: https://www.underworldcode.org/boundary-conditions-on-non-planar-boundaries/ + template: ../../templates/pdf + output: boundary-conditions-on-non-planar-boundaries.pdf + article_id: UWTN 2026-016 + article_version: 1.0.0 + software_version: underworld3 development @ 8b7c8b9e +--- +
Photo by Alexander Londoño / Unsplash
+ +Free-slip boundary conditions are used to simplify the physical behaviour at a domain boundary. +It may be a free-surface where the boundary deforms slightly in response to the internal flow, or it may +be an interface where the boundary layer thickness is so small that it cannot be resolved at the same time +time as the interior flow. The simplifying assumption: ignore the changes in shape, ignore the thin boundary layer, +treat the surface as impenetrable, and the tangential stresses as vanishingly small. + +$$ +\mathbf{u}\cdot\hat{\mathbf{n}} = 0 +\qquad\text{and}\qquad +\hat{\mathbf{t}}\cdot\boldsymbol{\sigma}\cdot\hat{\mathbf{n}} = 0. +$$ (eq-free-slip) + +In the weak form, boundary tractions appear as surface integrals. Multiplying the momentum +balance by a test function $\mathbf{w}$ and integrating by parts gives + +$$ +\int_\Omega \boldsymbol{\sigma} : \nabla\mathbf{w} \; \mathrm{d}V +- \int_{\partial\Omega} (\boldsymbol{\sigma}\cdot\hat{\mathbf{n}})\cdot\mathbf{w} + \; \mathrm{d}S + = \int_\Omega \mathbf{f}\cdot\mathbf{w} \; \mathrm{d}V . +$$ + +Drop the surface integral and you have imposed zero traction in all +directions — free *everything* (a free surface), not free slip. Constrain the surface-normal degrees of freedom +and the surface integral only addresses the tangential traction terms. + +On a Cartesian box, the first expression in {eq}`eq-free-slip` constrains a single velocity component. If you hold $u_x$ fixed on a vertical wall, the solver removes a row of unknowns, and there is nothing further to +discuss. On a sphere, an annulus, a mesh that does not align with the coordinates, or a surface with +topography, $\mathbf{u}\cdot\hat{\mathbf{n}}$ is not a single component of the unknown — it constrains a combination of unknowns at +each point, and leaves other combinations free. That has the potential to make a simplifying +assumption complicated to implement. + +Let's assume, for a moment, we confine ourselves to simple domains such as an annulus, or a spherical shell, +which are commonly used for planetary modelling. For each of these cases, there are +coordinate systems, and well known forms of the differential operators that do restore the boundary condition +to being a constraint in a single direction. Admittedly this requires reformulating all the equations, but for +a symbolic-first code such as underworld, this is quite straightforward. +This is the strategy used by CITCOMS [Zhong et al, 2008](https://doi.org/10.1029/2008GC002048). But not every domain boundary has a convenient +coordinate system to follow. Even accounting for slight ellipticity introduces + significant complexity in all the differential +operators; anything more complicated will not have a useful +coordinate reformulation. + +The second condition in {eq}`eq-free-slip` is also worth noting. +We don't generally think about this when we constrain a degree of freedom, +the other one/s, left unconstrained are *natural* to the problem. +They fall out as traction-free surface conditions +automatically in a finite element weak form. +If we cannot simply eliminate one degree of freedom at each point, +how **do** we satisfy all the parts of {eq}`eq-free-slip` ? + +## Four possibilities + +We outline four possible approaches (all of which you can try out in Underworld3). +They fall into two pairs: two +impose the constraint **weakly**, by adding a term to the momentum equation and +letting the solution satisfy the condition to within the accuracy of that term: a +direct penalty, and Nitsche's method, which differ in whether the term is +consistent. Two impose it **exactly**: by construction, changing the basis so +that the constraint is a component that can be struck out, or by a Lagrange +multiplier, adding an equation that enforces it. The weak pair have a parameter +to select that may need to be tuned for each problem and a floor + they cannot go below. The exact pair are not tuneable, and they return the +boundary traction as a side-effect of the solution. So does the direct penalty; +Nitsche is the one that does not. + +### 1. A direct penalty + +This could not be more simple, conceptually. We are working in a variational +environment, so we just add into our equation system, a term that punishes any flow through the boundary: + +$$ +\dots + \kappa\int_{\partial\Omega} +(\mathbf{u}\cdot\hat{\mathbf{n}})(\mathbf{w}\cdot\hat{\mathbf{n}}) +\; \mathrm{d}S . +$$ + +$\kappa$ is a single scalar. It has to absorb the scale of the problem itself, which is why +the value that works is a property of the model rather than a default. + +One line, no new machinery, and it works on any geometry. +What we are solving is a mildly perturbed problem and it is perturbed by +exactly the amount the constraint cannot be satisfied: the discrete solution +sits where the penalty term balances the boundary +traction and this leaves $\mathbf{u}\cdot\hat{\mathbf{n}}$ small but not zero. + +Making the residual $\mathbf{u}\cdot\hat{\mathbf{n}}$ smaller means pushing harder, and pushing harder +degrades the condition-number of the operator. The error is traded against the conditioning, +and (discussed below), this trade-off eventually stops returning any benefit. + +Underworld codes this term it as a boundary traction opposing normal flow, using the +surface normal at the quadrature points ($\Gamma$) provided by PETSc: + +```python +G = mesh.Gamma +penalty = 10000 +stokes.add_natural_bc(penalty * G.dot(v.sym) * G, "Upper") +``` + +**The topography** is the penalty term itself. That term is the traction the +condition holds the wall with, so on the boundary + +$$ +\sigma_{nn} = -\kappa\,(\mathbf{u}\cdot\hat{\mathbf{n}}), +\qquad +h = -\frac{\sigma_{nn} - \overline{\sigma_{nn}}}{\Delta\rho\,g} , +$$ + +which is arithmetic on values the solve already returned — nothing is +differentiated and nothing is solved. + +```python +n = mesh.boundary_normal("Upper") +sigma_nn = -penalty * n.dot(v.sym) +``` + +That is a saving in cost and not in accuracy. The leak and the traction are the +same quantity scaled by $\kappa$, so a coefficient too small to hold the +boundary reports a traction that is short in the same proportion. + +### Nitsche's method + +The reason the penalty is only accurate in the limit is that it is not +*consistent*: substituting the true solution does not fully satisfy the equation, +because the true solution is subject to a separate boundary traction the +penalty form ignores. Nitsche's method [@Nitsche_1971] restores +consistency by carrying that traction explicitly: + +$$ +\dots +- \int_{\partial\Omega} (\hat{\mathbf{n}}\cdot\boldsymbol{\sigma}(\mathbf{u}) + \cdot \hat{\mathbf{n}})(\mathbf{w}\cdot\hat{\mathbf{n}}) \; \mathrm{d}S +- \int_{\partial\Omega} (\hat{\mathbf{n}}\cdot\boldsymbol{\sigma}(\mathbf{w}) + \cdot \hat{\mathbf{n}})(\mathbf{u}\cdot\hat{\mathbf{n}}) \; \mathrm{d}S ++ \frac{\gamma}{h}\int_{\partial\Omega} + (\mathbf{u}\cdot\hat{\mathbf{n}})(\mathbf{w}\cdot\hat{\mathbf{n}}) + \; \mathrm{d}S . +$$ + +The first of the three is the consistency term: it is the boundary traction the +integration by parts produced, and including this makes the true +solution satisfy the discrete equations exactly. The second is its transpose, +which keeps the form symmetric and buys optimal convergence in $L^2$. The third +is the penalty again, and it is still needed — but now for *stability* rather +than for accuracy, and $\gamma$ has a threshold set by an inverse inequality +rather than being an unspecified free parameter. + +This is a real improvement and it is still done through a weak imposition. The constraint +holds to the accuracy of the discretisation, not to the accuracy of the +arithmetic — measured below, it leaks a few parts in a thousand on a typical mesh, + and the leak falls with increasing mesh resolution. + +**The topography** has to be recovered from the solved fields, which is where +Nitsche parts company with the direct penalty. Its boundary term is the penalty +part *less* the consistency terms, and those are written in +$\boldsymbol{\sigma}(\mathbf{u})$, so the traction cannot be read off without +differentiating the velocity. Recover $\sigma_{nn}$ that way and divide by +$\Delta\rho\,g$. + +### A constraint equation, with a multiplier + +The two strategies above add a *term* to the weak form of the equation. +This approach adds an *equation*. + +Carry a scalar field $\lambda$ on the boundary and require, as a row of the +system in its own right, + +$$ +\int_{\partial\Omega} (\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n)\, q +\; \mathrm{d}S = 0 \quad \text{for all } q , +$$ + +where $\tilde{u}_n$ is the prescribed wall-normal velocity — zero for free slip, +and a datum if the wall is being driven — and $\lambda$ enters the momentum row +as the traction $\lambda\hat{\mathbf{n}}$ that holds the constraint. It is a +Lagrange multiplier, and the system becomes a larger saddle point: velocity, +pressure, and now $\lambda$. + +$\lambda$ has units of stress. At convergence it *is* $\sigma_{nn}$ on that +boundary, so dividing by $\Delta\rho\,g$ (density contrast $\times$ gravity) +is the dynamic topography. + +The constraint row is exact, so unlike a penalty there is no parameter whose +size decides how well it holds. Two practical things do have to be dealt with. + +- $\lambda$ is carried as a full-domain field but only its boundary trace means + anything, so the interior degrees of freedom are constrained out of the global + system in the section before the solve sees it. They are not solved for and + are not stored in the $[p, \lambda]$ block. +- The $[p, \lambda]$ Schur complement is poorly conditioned on its own, so an + augmented-Lagrangian term + $r(\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n)\,\hat{\mathbf{n}}$ is added to the momentum row. It does not change what the constraint enforces, + because the $\lambda$ row still carries the exact constraint. + +```python +stokes = uw.systems.Stokes_Constrained(mesh, velocityField=v, pressureField=p) +lam = stokes.add_constraint_bc(0.0, "Upper") +stokes.solve() +``` + +The `Stokes_Constrained` solver carries three fields in PETSc — velocity, pressure and +$\lambda$ — and splits them two ways: velocity against the pair, with the volume +constraint (incompressibility) and the boundary constraint grouped into a single Schur +block. The boundary constraint applies to surface nodes only. + +:::{note} The boundary traction is an unknown of the solve, not a post-processing step +At convergence the momentum row's boundary term *is* the normal traction, so there is +nothing to recover from the velocity field afterwards. The solver returns it through +`traction` and, divided by $\Delta\rho g$, through `topography`. +::: + +That boundary term is the whole boundary load, +$\lambda + r(\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n)$, not the +multiplier alone. The second part vanishes only where the constraint row is +satisfied exactly; discretely it is satisfied to the solver's tolerance, and $r$ +multiplies that residual back into the traction. With a viscosity-weighted $r$ +and a lateral viscosity contrast it can be the largest term of the result, +so the two parts are not separable in practice. + +### Rotating the degrees of freedom + +In this approach, we stop *"asking for"* the constraint and just impose it. +At each constrained node, we change +the coordinate basis in which the velocity unknowns are expressed, from the global +Cartesian frame to the local $(\hat{\mathbf{n}}, \hat{\mathbf{t}})$ frame. In +that basis "no flow through the boundary" is again a single component, and it +is removed the same way it would be on a box. + +Collect the per-node rotations into a block-diagonal $Q$, equal to the identity +at every node that is not constrained. The rotated system is + +$$ +\hat{A} = Q^{T} A Q, \qquad \hat{\mathbf{b}} = Q^{T}\mathbf{b}, +\qquad \mathbf{u} = Q\hat{\mathbf{u}} , +$$ + +and the wall-normal row of $\hat{A}$ is struck out. The constraint then holds to +machine precision, because it is not being solved for at all. + +**This is the classical strategy**. It is in the early +finite-element literature, and Engelman, Sani and Gresho [@Engelman_1982] +were already reviewing +the alternatives and choosing between them on grounds of global mass +conservation in 1982. What is worth explaining is not the idea but why, given +that it is exact and the others are not, it is the least used of the three. + +**The topography** is the reaction of that struck row — the force the constraint +had to supply — de-smeared by the boundary mass to turn an integrated nodal load +into a pointwise stress, + +$$ +\sigma_{nn} = -M_\Gamma^{-1}\left.(A\mathbf{u} - \mathbf{b})\right|_\Gamma, +\qquad +h = -\frac{\sigma_{nn} - \overline{\sigma_{nn}}}{\Delta\rho\,g} , +$$ + +which is the consistent boundary flux of Zhong, Gurnis and Hulbert +[@Zhong_1993]. In Underworld3, the solver's `boundary_normal_traction()` and `dynamic_topography()` return +these. Nothing is differentiated and nothing is solved: in two dimensions +$M_\Gamma$ is lumped and the de-smear is a division. + +## Trade-offs + +Rotating the degrees of freedom leaves the discrete problem in a **mixed +basis**. Interior nodes hold $(u_x, u_y)$; constrained nodes hold +$(u_n, u_t)$. Nothing about that is difficult in itself, but everything +downstream has to agree about which nodes are which. + +```{figure} figures/rotated-basis.svg +:alt: Two panels. On the left, a meshed domain bounded above by a free surface that rises on the left and falls on the right with an inflection between, so that the outward normal points in a different direction at every surface node. Surface nodes are drawn as filled circles each carrying its own rotated pair of arrows labelled n and t; interior nodes are open circles, with one carrying the unrotated x and y arrows shared by all of them. On the right, a block diagram. A red block labelled "Velocity solve, rotated" contains the rotated operator and right-hand side, and encloses a smaller block labelled "Multigrid" listing three rows: prolongation becomes Q-transpose P, coarse operators inherit Q through RAP, and the coarse solve uses SVD for the rigid rotations. A separate green block beside it, labelled "Fieldsplit / Schur solve", carries the pressure and constraints and is marked as never seeing a rotated vector. A single arrow labelled v equals Q v-hat leaves the velocity block at its boundary and branches, one branch entering the Schur block and the other leaving for output, advection and the surface update. + +Where the rotation lives. The obligation is contained: the velocity solve is +rotated and carries its multigrid with it, while the Schur complement and the +pressure solve beside it never handle a rotated vector, because the pressure +block carries no boundary condition of this kind. One un-rotation sits on the +boundary between them and feeds both. +``` + +Four objects carry $Q$: the operator, the right-hand side, the solution on the +way out, and the multigrid prolongation. The coarse operators inherit it +through the Galerkin triple product rather than being rotated separately, and +the coarse solve should be an SVD, because a Galerkin-coarsened rotated +operator inherits any rigid-rotation null space of the constrained problem +(the exact constraint makes the null space of the sphere and the annulus a dominant +feature of the solve). + +We do not know the cost of the addtional complexity on the solver and setup times, +or on the accuracy of the solution but this can be measured and will differ from problem +to problem. + +## Choice of the surface normal + +In a discrete representation of a curved surface, the normal can be defined in various +ways. The boundary of a +discretised domain is a set of straight facets, and the assembled constraint is +an integral over those facets. +The node normal consistent with that integral is +the average of the adjacent facet normals **weighted by facet measure** — not +the normal of the smooth surface the mesh approximates, and not the facet +normal on its own.**This is the consistent normal of Engelman, Sani and Gresho** +[@Engelman_1982]. They derived this result in 1982 from global conservation of mass. + +The analytic normal is exact for the geometry and therefore inconsistent +with the discretisation: the solver is not solving +on the sphere (or annulus), it is solving on the polyhedral approximation to the sphere. + +Using the **facet** normal is worse than inconsistent, and this is the one place +where the wrong choice does real damage. In 2D, Imposing +$\mathbf{u}\cdot\hat{\mathbf{n}} = 0$ facet by facet asks a node shared by two +facets to satisfy two different constraints, and two independent constraints on a +two-component velocity provide no freedom. Push the penalty higher, and the vertex +velocities go to zero: the flow is being asked to stay inside a polygon rather +than a circle, and the discrete limit is a different problem from the smooth one. +Refining the mesh does not approach the smooth answer, because it is not +converging to it. + +On an annulus with a free slip boundary, the direct-penalty approach locks at high penalty values ($\sim 10^6$) if facet normals are used in the +constraint equation. In the figure below, the node-normal approach does solve and reproduces the analytic solution (described in detail in the next section) + +```{figure} figures/locking.png +:alt: Three annulus solutions side by side on one colour scale from zero to 5.0e-3, blue for slow speeds and red for fast, with the triangular mesh drawn over each. The left panel is the exact solution: two deep red patches of fast flow sit against the outer boundary on the left and right of the annulus, with a blue slow ring inside them. The middle panel is the same problem solved with a direct penalty against the facet normal: the red patches at the outer boundary are gone and the whole outer half is blue, the peak speed having fallen from 5.0e-3 to 3.8e-3, while a pale ring survives near the inner boundary. The right panel is the same penalty against the measure-weighted node normal and is indistinguishable from the exact panel, with a peak speed of 5.0e-3. + +The same problem, the same coefficient, the same colour scale. Against the facet +normal the flow along the outer boundary is suppressed — the peak speed falls by +a quarter and the two fast lobes at the boundary are gone. Against the +measure-weighted node normal it is the exact solution. +``` + +Everything that follows uses the node normal, which is what `add_nitsche_bc` and +`add_rotated_freeslip_bc` take by default and what `mesh.boundary_normal` returns. +The facet normal does not appear again. + +## When the choice of constraint matters + +Solve a convection model with any of these approaches, and the +velocity field is the same to plotting accuracy. +Generally speaking, a leak of order $10^{-3}$ or $10^{-4}$ in $\mathbf{u}\cdot\hat{\mathbf{n}}$ is +within the expected accuracy of the solution on the mesh and the main +driver of which method to choose should be solver efficiency (wall time). + +The difference in the methods appears when the wall-normal traction is a +required output of the model: dynamic topography, geoid, gravity, +or a plate-boundary force balance require accurate integration of boundary stresses. +Here the choice becomes more subtle, +and the methods have quite different accuracies, and different efficiencies. + +### The benchmark + +Kramer, Davies and Wilson [@Kramer_2021] give exact Stokes solutions in a +cylindrical annulus, and their `assess` package publishes the radial stress as +well as the velocity, which is what makes it an oracle for this question rather +than only for the flow. Underworld wraps it as `uw.analytic.CylindricalStokes`. + +The case used throughout is the smooth one: a density anomaly +$(r/r_o)^k \cos n\theta$ with $n = 2$ and $k = 3$, viscosity 1, free slip on both +radii. On the outer boundary the exact radial stress is a single harmonic, + +$$ +\sigma_{rr}(r_o, \theta) = 0.1506696\,\cos 2\theta , +$$ + +fitted to a residual of $10^{-16}$, so the whole of the surface stress is that one +amplitude and the error in it is one number. The treatment under test is on the +outer radius; the inner carries the exact analytic velocity as a Dirichlet +condition, so it is the only free-slip condition in the model. + +Two things are measured on every solve: + +- **the surface permeability** — the largest $\mathbf{u}\cdot\hat{\mathbf{n}}$ on the outer + boundary, against the true radial direction, divided by the flow speed. The + fraction of the flow going through an impermeable boundary. +- **the boundary stress error** — the relative error in that harmonic amplitude, recovered + from the solved fields by projection, which is the route every method has available. + +Penalty at $\kappa = 10^4$, Nitsche at $\gamma = 10$. + +| cell size | penalty | Nitsche | multiplier | rotated | +|---|---|---|---|---| +| 0.150 | 3.1 × 10⁻³ / 2.5 × 10⁻² | 1.0 × 10⁻² / 5.8 × 10⁻² | 2.3 × 10⁻⁴ / 2.4 × 10⁻² | 5.3 × 10⁻¹¹ / 2.4 × 10⁻² | +| 0.100 | 3.0 × 10⁻³ / 1.1 × 10⁻² | 2.4 × 10⁻³ / 2.4 × 10⁻² | 1.0 × 10⁻⁴ / 1.0 × 10⁻² | 5.8 × 10⁻¹¹ / 1.0 × 10⁻² | +| 0.075 | 3.0 × 10⁻³ / 7.2 × 10⁻³ | 1.2 × 10⁻³ / 1.5 × 10⁻² | 8.6 × 10⁻⁵ / 6.3 × 10⁻³ | 1.2 × 10⁻¹⁰ / 6.2 × 10⁻³ | +| 0.050 | 3.0 × 10⁻³ / 3.6 × 10⁻³ | 2.7 × 10⁻⁴ / 6.3 × 10⁻³ | 6.7 × 10⁻⁵ / 2.7 × 10⁻³ | 1.2 × 10⁻¹⁰ / 2.7 × 10⁻³ | + +Reading the leak first, Nitsche leaks parts in a thousand and improves with the mesh +— the rate consistency buys. The multiplier is an order of magnitude better and +improves faster. The rotated constraint does not move: it sits at the solver's +floor at every resolution, because the mesh has nothing to do with it. The +penalty does not improve either, and for the opposite reason — its leak is set by +the penalty coefficient rather than by the discretisation. + +Now read the stress beside it, and the ranking is not the same. **Every treatment +that imposes the constraint properly lands on the same stress error at a given +mesh**: 6.3 × 10⁻³ for the multiplier and 6.2 × 10⁻³ for the rotated constraint at +cell 0.075, where their leaks differ by nine orders of magnitude. What sets that +number is the recovery — a projection of a stress differentiated out of a +piecewise-quadratic velocity — and not the boundary condition underneath it. A +constraint held to $10^{-10}$ buys nothing over one held to $10^{-4}$ if the +answer is then recovered the same way. + +Nitsche is the exception, at twice the error of the others on the coarser meshes. +Its $\gamma = 10$ is enough for the leak and not for the stress: at +$\gamma = 100$ the leak improves by a factor of nearly forty and the stress by a +factor of two, onto the same floor as everything else, after which more $\gamma$ +buys nothing. + +### Traction extraction v. Stress recovery + +Three of the four treatments do not have to recover anything. The two exact ones +carry the traction as a constraint reaction or as an unknown, and the direct +penalty carries it as the term it adds. Against the same exact answer: + +| cell size | rotated, reaction | multiplier, traction | penalty, $\kappa(\mathbf{u}\cdot\hat{\mathbf{n}})$ | recovered by projection | +|---|---|---|---|---| +| 0.150 | 6.8 × 10⁻³ | 8.6 × 10⁻³ | 9.0 × 10⁻³ | 2.4 × 10⁻² | +| 0.100 | 3.3 × 10⁻³ | 8.5 × 10⁻⁴ | 1.2 × 10⁻³ | 1.0 × 10⁻² | +| 0.075 | 2.1 × 10⁻³ | 1.7 × 10⁻³ | 2.3 × 10⁻³ | 6.3 × 10⁻³ | +| 0.050 | 1.1 × 10⁻³ | 1.4 × 10⁻³ | 2.2 × 10⁻³ | 2.7 × 10⁻³ | + +Better than the projection on the same solve at every resolution, using the +expressions given with each method above — by an order of magnitude at best, and +by very little where the penalty's coefficient sets its floor. None of the +three falls smoothly with $h$: part of what they report is the constraint +residual, and how far a particular solve drove that is not a function of the +mesh. For the penalty that is the whole story below cell 0.10 — its leak is +3 × 10⁻³ at every resolution here, set by $\kappa$, and the traction it reports +cannot be better than the constraint it holds. + +### Influence of penalty parameters + +The two weak methods look alike in the comparison above, but this is for a fixed, tuned penalty parameter. + +| $\kappa$ (penalty) | leak | | $\gamma$ (Nitsche) | leak | +|---|---|---|---|---| +| 10² | 2.6 × 10⁻¹ | | 1 | diverged | +| 10³ | 2.6 × 10⁻² | | 10 | 1.7 × 10⁻³ | +| 10⁴ | 2.6 × 10⁻³ | | 100 | 2.7 × 10⁻⁴ | +| 10⁵ | 3.0 × 10⁻⁴ | | 1000 | 3.0 × 10⁻⁵ | +| 10⁶ | 4.5 × 10⁻⁵ | | 10⁴ and above | diverged | + +Nitsche is bounded at both ends. Below $\gamma \sim 1$ the form is no longer +coercive and no amount of solver tuning recovers a solution; from $\gamma \sim 10^4$ in this +problem, the line search stops converging. The virtue of $\gamma \sim 10$ is that it +sits within that window on any mesh, because $\gamma$ is dimensionless +and the term it scales already carries $\mu / h$. + +The penalty coefficient scales differently because it directly penalises the +value of the velocity across the boundary. It should therefore scale with the +characteristic velocity which is best estimated from the magnitude of the forcing +terms and the resisting viscosity. + +Written against the node normal, the penalty simply trades: a decade of +coefficient for a decade of leak, all the way to $10^6$, with no wall in this +problem. What it does not do is converge with resolution — the leak is bought with the parameter +rather than with the mesh resolution — so the coefficient has to be re-chosen whenever the +forcing or the viscosity changes. + +### Multiplier and CBF equivalence + +The two expressions given above for computing topography from the boundary reaction are exactly equivalent. +Write the momentum row's boundary term out and the identity is immediate: the assembled +load is $M_\Gamma\,(\lambda + r(\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n))$, +and at convergence it +balances the volume residual restricted to the boundary, which is precisely the +nodal load the consistent boundary flux back-calculation reads +[@Zhong_1993]. So + +$$ +\lambda + r(\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n) + = -M_\Gamma^{-1} \left. (A\mathbf{u} - \mathbf{b}) \right|_\Gamma , +$$ + +which is the rotated constraint's reaction, de-smeared with the same boundary +mass. The multiplier is not a second, independent estimate of the surface stress: +it is the same computation, arrived at by carrying the traction as an unknown +instead of reading it out of the residual afterwards. + +### Lateral viscosity contrast case + +No published solution has both a curved boundary and a laterally varying viscosity, +so the case where weak constraints are most often reported to give trouble is a +separate test with a simple geometry. SolCx is a good example: unit box, free +slip on all four walls, viscosity 1 to the left of $x = 0.5$ and $\eta_B$ to the +right. `uw.analytic.SolCx` publishes the exact dynamic topography on the top +wall. Three walls carry the ordinary component condition and the treatment under +test is on the top wall alone. + +On a box every treatment reduces to holding one velocity component. +What it can say is whether a treatment holds the traction +it was given when the viscosity beside it jumps. However, on a box the rotated constraint's per-node rotation is the identity, so the +table below exercises none of the rotation machinery. We check that separately on an +equivalent problem: the domain, the gravity vector and the exact solution all turned by +45$^\circ$, with every wall then carrying the rotated constraint, because a component +condition cannot express $\mathbf{u}\cdot\hat{\mathbf{n}} = 0$ on a tilted wall. Turned, +the constraint still holds to machine precision and the velocity error is unchanged at +8.8 × 10⁻⁶; imposing the un-turned condition on those same walls instead lets 71% of the +flow through the boundary. + +In the table below, we show the relative $l_2$ error of the surface topography along the top wall, mean removed, +at 32 × 32 elements. Each entry is the whole wall and then the wall with two +elements trimmed from each end. + +| $\eta_B/\eta_A$ | component Dirichlet | penalty, $10^4$ | multiplier | rotated | +|---|---|---|---|---| +| 10 | 0.048 / 0.054 | 0.045 / 0.051 | 0.048 / 0.054 | 0.048 / 0.054 | +| 10² | 0.072 / 0.081 | 0.056 / 0.060 | 0.072 / 0.081 | 0.072 / 0.081 | +| 10³ | 0.075 / 0.084 | 0.234 / 0.230 | 0.075 / 0.084 | 0.075 / 0.084 | +| 10⁴ | 0.076 / 0.085 | 0.698 / 0.697 | 0.076 / 0.085 | 0.076 / 0.085 | +| 10⁶ | 0.076 / 0.085 | 0.992 / 1.000 | 0.075 / 0.084 | 0.076 / 0.085 | + +```{figure} figures/topography.png +:alt: Two line plots of surface topography along the top wall from x=0 to x=1, mean removed, at viscosity contrasts of 100 and a million. In both, the exact answer is a thick grey curve falling from +0.29 at the left, flattening near +0.21, dropping sharply at the viscosity step at x=0.5 and continuing down to -0.38 at the right. At a contrast of 100 every curve lies on the grey one. At a contrast of a million they separate: the component Dirichlet, the rotated reaction and the traction lambda + r(u.n - u_n) still lie on the exact curve, while the multiplier field lambda alone is a nearly flat line near zero reaching only 0.04, and the penalty at 1e4 is a second nearly flat line near zero. Nitsche does not solve at either contrast and is absent. + +Surface topography along the top wall, against the exact answer. At a contrast of +100 nothing distinguishes the treatments. At $10^6$ the multiplier field +$\lambda$ carries almost none of the traction on its own — the augmentation +holds the rest — while $\lambda + r(\mathbf{u}\cdot\hat{\mathbf{n}} - +\tilde{u}_n)$, which is what `traction()` returns, lies on the exact curve. The +penalty has failed by this contrast: its coefficient is a bare number and cannot +be large against $10^6$ and moderate against 1 at the same time. +``` + +**The three exact treatments agree to three figures at every contrast**, whole +wall and trimmed alike. That is the result to take from this half, and it took +the multiplier reporting the whole traction rather than $\lambda$ alone, and the +rotated constraint holding the corner where it meets the side walls. + +**Read the first column as the reference.** The component Dirichlet condition is +exact and has no parameter, and its velocity error is 8.8 × 10⁻⁶ at a contrast of +$10^6$. It still reads 0.085. That number is the recovery's error, not a boundary +condition's: on the stiff half the recovered $\sigma_{zz}$ is a difference between +the pressure and $2\eta\,\partial_z u_z$ with $\eta = 10^6$, so a relative velocity +error of $10^{-5}$ appears in the stress. + +**A bare penalty coefficient cannot serve both halves.** At $10^4$ it is the best +column in the table at low contrast (the constraint is weak enough not to fight +the recovery) but by $10^6$ it is meaningless: 0.992, which is to say the recovered +topography carries none of the signal. Reading $\kappa(\mathbf{u}\cdot\hat{\mathbf{n}})$ +instead of recovering the stress gives the same numbers to three figures, here and +at every coefficient tried, which is the point made above from the other side: the +term is the traction, so it inherits the error in the constraint rather than +curing it. Scaling the coefficient by the local +viscosity is the obviously right thing to want, but the solver does not converge here at any +magnitude we tried, from $\eta$ to $10^3\eta$. + +**Nitsche has no column in this table.** In our implementation it is unreliable on a +boundary that mixes essential patches with Nitsche patches, which is what this test +asks for: the top wall weak, the other three held strongly. We could not reach a +converged solution for this example at any penalty. Imposed weakly on all +four walls does converge, but that is a different problem. + +### Solver timing + +Seconds on the annulus, uniform viscosity, one core, direct solver: the solve, +and then the surface traction by whatever route that treatment has. Median of +three timed repeats after an untimed warm-up, run sequentially. Two sizes, +because below about ten thousand nodes the four are separated by less than the +run-to-run spread and there is nothing to read. + +| velocity nodes | penalty | Nitsche | multiplier | rotated | +|---|---|---|---|---| +| 28 338 | 0.17 / — | 0.18 / 0.267 | 0.26 / — | 0.16 / 0.010 | +| 71 424 | 0.45 / — | 0.47 / 0.657 | 0.67 / — | 0.43 / 0.016 | + +The dash indicates that the multiplier and the penalty do not require any additional *solver* — +$\lambda$ is a finite element field in its own right, so its nodal values are the +traction, pointwise, and $\lambda + r(\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n)$ is an +expression evaluated where it is wanted. The penalty traction is recovered similarly as the penalty scaling +the leakage velocity. + +**The rotated constraint's reaction does need one step**. +The reaction is an *integrated* nodal load, +$\int_\Gamma \sigma_{nn}\,\phi_i\,\mathrm{d}S$, which is $M_\Gamma$ times the +pointwise traction. Turning it into a pointwise value means undoing that boundary +mass. On a 2-D trace, and on 3-D P1 triangles, the lumped mass is diagonal and +undoing it is a division — the 10 to 16 ms above. On **3-D P2 triangles it is a +true solve**: the lumped row sums vanish at the vertices, so the consistent +trace mass has to be assembled and solved. That is the one place the CBF +route pays for being a +back-calculation. + +**The multiplier's solve costs about 50% more**, consistently — 0.67 s against +0.43 s at 71 000 nodes. That is the extra field and the larger saddle point. +Rotating the degrees of freedom costs nothing measurable against the weak forms: +the rotation is a sparse orthogonal transform on a boundary's worth of rows. + +**Nitsche has to recover surface stress by differentiating the solution**, and +the projection that does it costs *more than the Stokes solve did* — 0.63 s +against 0.45 s — so asking a Nitsche model for its surface stress roughly +doubles the timestep. The direct penalty need not pay that: its own term is the +traction, read as arithmetic on the boundary nodes. + +The obvious question is whether that is the method's cost or the recovery's. A +global $L^2$ projection to get values on a thousand boundary nodes is plainly +more work than the job requires, and the consistent boundary flux is an available alternative, +reading the assembled residual rather than differentiating anything. However, we find that **it does +not work for a weakly imposed condition**, and the reason is the same one that +makes it unavailable to the multiplier: + +| cell 0.0125 | projection | CBF back-calculation | +|---|---|---| +| penalty, node normal | 0.651 s, error 1.1 × 10⁻³ | 0.224 s, **error 1.00** | +| Nitsche | 0.666 s, error 3.6 × 10⁻⁴ | 0.230 s, **error 1.00** | +| rotated | 0.602 s, error 1.6 × 10⁻⁴ | 0.206 s, error 1.6 × 10⁻⁴ | + +An error of 1.00 is the metric reporting that nothing was recovered. A reaction +exists in the residual only where a row has been *constrained*; a weak condition +supplies its traction as a term inside the row it acts on, so the residual there +is balanced at convergence and there is nothing left to read. What rescues a weak +condition is not the residual but its own term: the multiplier supplies $\lambda$, +which is the traction as a field, and the penalty supplies +$\kappa(\mathbf{u}\cdot\hat{\mathbf{n}})$, which is the traction as boundary +arithmetic. Nitsche's term is written in terms of +$\boldsymbol{\sigma}(\mathbf{u})$, so reading it back still requires differentiating +the answer. + +That is the structural statement the timings are really making, and it follows +the two pairs exactly: + +| how the constraint is imposed | the traction is | to read it | +|---|---|---| +| weakly, by a penalty term | the term itself | evaluate $\kappa(\mathbf{u}\cdot\hat{\mathbf{n}})$ | +| weakly, by Nitsche | inside a term written in $\boldsymbol{\sigma}(\mathbf{u})$ | differentiate the solution | +| exactly, by construction (rotated) | the constraint reaction | de-smear the nodal load | +| exactly, by a multiplier | an unknown of the system | read the field | + +Only the second row pays. Its cost is negotiable — a recovery restricted to the +boundary would be cheaper than a global projection — but the differentiation is +not. + +:::{note} What these timings are not +Two-dimensional, one core, direct solver. What they measure is the difference +between a recovery *solve* and boundary arithmetic, which is structural and +survives scaling. They say nothing about a large parallel spherical shell, where +the rotated velocity block's multigrid and the multiplier's larger Schur +complement are the terms that matter and neither is exercised here. +::: + +### Which one to use + +For a model that consumes the velocity and nothing else, all four are the same to +plotting accuracy — provided the constraint is written against the node normal. +That proviso is the only one that can spoil the velocity, and it costs one line. + +When the wall-normal traction is needed: + +- **Rotated free slip is the default.** It holds the constraint to machine + precision rather than to the discretisation, its reaction is the most accurate + surface stress measured here, and that reaction is nearly free. The price is + structural — a mixed basis that the multigrid has to carry — and it is paid + once, inside the solver, rather than by the person setting up the model. +- **The multiplier is its equal on accuracy** and returns the same object by a + different route; take it when you want the traction as an unknown of the + system, or when the constrained problem's conditioning suits you better. It + costs about 50% more to solve. +- **Nitsche is the one to reach for when the boundary condition must change + during the model** — a wall that begins as a prescribed velocity and relaxes to + a prescribed traction is a Nitsche problem, because a hard constraint cannot + morph. Budget for tuning $\gamma$ against the stress and not against the leak, + and expect the window to move with the viscosity contrast. If the viscosity + contrast is large with jumps or strong gradients along the boundary, be very careful + if you choose Nitsche. +- **A direct penalty is fine for a velocity-only model** and needs the node + normal, a coefficient chosen per problem, and a check on something physical + before the answer is believed. It will also hand back the traction it holds + the wall with, at no cost, and that reading is worth no more than the + coefficient behind it. It has the advantage that this is pure, direct penalty + on the weak form and can be used for many things beyond simply boundary conditions. + Good for a first pass on a very general idea. + + +## Using it + +```python +import underworld3 as uw + +mesh = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, cellSize=0.05) +stokes = uw.systems.Stokes(mesh) + +# Value first: 0 is free slip. A non-zero scalar or expression prescribes the +# wall-normal datum u.n = u_n strongly instead. +stokes.add_rotated_freeslip_bc(0.0, "Upper") +stokes.add_rotated_freeslip_bc(0.0, "Lower") + +stokes.solve() + +# The constraint reaction, which is the boundary normal traction. +sigma_nn = stokes.boundary_normal_traction("Upper") +``` + +Leave the normal to Underworld unless the constraint has to follow the true +surface rather than the mesh. Passing an analytic normal — `X / |X|` on a +sphere — is exact for the geometry and keeps a consistency error against the +faceted assembly, which is usually not what you want. + +Reach for Nitsche when the boundary condition has to **change during the +model**. A hard constraint cannot morph: a wall that begins as a prescribed +velocity and relaxes to a prescribed traction is a Nitsche problem, because the +rotated constraint is either imposed or it is not. + +
Comments
Discussion of these notes happens in GitHub Discussions, so it stays with the source and is searchable alongside it.
Read the discussionStart one
diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/README-stress-test.md b/articles/boundary-conditions-on-non-planar-boundaries/examples/README-stress-test.md new file mode 100644 index 0000000..87730a0 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/README-stress-test.md @@ -0,0 +1,49 @@ +# The examples, and what each one measures + +Every table and figure in the note is produced by one of these. Run them from +this directory with an Underworld3 environment on the path. + +| script | what it produces | +|---|---| +| `leak.py` | the constraint tables — how much flow each treatment lets through, under refinement (`sweep`) and against its own parameter (`params`) | +| `stress.py` | the surface-stress tables against the exact `uw.analytic.CylindricalStokes` answer: `sweep`, `params`, `locking`, `control` | +| `solcx.py` | the lateral-viscosity half, against `uw.analytic.SolCx`: `sweep`, `contrast`, `params`, `control` | +| `generate-locking-figure.py` | `figures/locking.png` and `figures/banner.png` | +| `generate-rotated-basis.py` | the data behind `rotated-basis.typ`, which draws `figures/rotated-basis.svg` | + +## The traps, all paid for once + +- **`v.array` is `(N, 1, dim)`.** It broadcasts silently against `(N, dim)` + normals and returns ~1e-16 projections for a 1e-2 velocity, with no error. + `np.squeeze` it. Every number in the first leak run was this artefact. +- **Give every metric a negative control and run it.** Two free-slip circles + leave the rigid rotation unconstrained, and that nullspace is purely + tangential — so a radial leak metric read 2e-14 on a solve that had diverged + with `|u| = 2.7e5`. The inner boundary is now held. +- **Check `snes.getConvergedReason() > 0` before tabulating anything.** Diverged + runs leave plausible numbers in the array; two nearly reached the note. +- **Vertex against edge midpoint.** On a curved boundary, vertex values of + `sigma_nn` carry the O(h) facet error and midpoints are superconvergent + (underworld3#414). `stress.py` splits them and the note says which it uses. +- **Select a boundary trace by the mesh LABEL, not by a radius band.** A band + that narrows with the mesh admits a different node set at each resolution, and + the earlier version of this comparison drifted by four nodes between methods + because of it. +- **Enclosed domain**: the multiplier, the pressure and the traction are each + determined only up to a constant. Compare deviations, which is what topography + is anyway. +- **A `Piecewise` viscosity inside a boundary penalty term does not solve** on + SolCx, at any magnitude tried. + +## What is still open + +- The spherical case. `uw.analytic.Zhong2008` publishes + `.response().surface_topography = 0.4191904156575601` for the default + degree-2 case (load r = 0.775, r_inner = 0.55, isoviscous), which is the right + oracle for dynamic topography proper. It is a 3-D shell, and that is the cost + step this note stopped short of. +- Behr (2004) reports non-physical recirculation at curved walls even with the + consistent normal. We have not looked for it here. +- underworld3 issues #607 (the multiplier misses the augmented-Lagrangian share) + and #608 (the rotated constraint at a corner shared with a component + condition) both came out of these runs and are open. diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-locking-figure.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-locking-figure.py new file mode 100644 index 0000000..e2f51b3 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-locking-figure.py @@ -0,0 +1,106 @@ +"""The figure for "A constraint that is satisfied, and wrong". + +Three solves of the same annulus problem, rendered with the same colour scale: +the exact solution, a direct penalty written against the facet normal, and the +same penalty written against the measure-weighted node normal. The coefficient +is 1e6 in both penalty panels, so the only difference between them is the +normal. + +Writes: + figures/locking.png the three-panel figure + figures/banner.png a wide two-panel crop for the article banner + + python3 generate-locking-figure.py + +Run against underworld3 `bugfix/multiplier-traction` (PR #617); the +constrained solver's `traction()` is the fix this note prompted. +""" +import pathlib + +import numpy as np +import pyvista as pv + +pv.OFF_SCREEN = True + +import underworld3 as uw +import underworld3.visualisation as vis + +import stress as S + +CELL = 0.075 +PENALTY = 1.0e6 +HERE = pathlib.Path(__file__).resolve().parent +FIGURES = HERE.parent / "figures" + + +def speed_mesh(mode): + """(pyvista mesh carrying |u|, mesh edges, max speed) for one treatment.""" + mesh, stokes, v, exact = S.build(mode, cell=CELL, penalty=PENALTY) + stokes.solve() + assert S.converged(stokes), "%s did not converge" % mode + pv_v = vis.meshVariable_to_pv_mesh_object(v) + u = np.squeeze(np.asarray(v.array)) + pv_v.point_data["speed"] = np.linalg.norm(u, axis=1) + edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() + return pv_v, edges, exact, v + + +def exact_speed_mesh(reference): + """The same object built from the exact velocity, for the first panel.""" + mesh, stokes, v, exact = reference + field = uw.discretisation.MeshVariable("Uplot", v.mesh, v.mesh.dim, degree=2) + field.data[:] = exact.evaluate("velocity", field.coords) + pv_v = vis.meshVariable_to_pv_mesh_object(field) + pv_v.point_data["speed"] = np.linalg.norm( + np.squeeze(np.asarray(field.array)), axis=1) + return pv_v + + +def panel(plotter, index, pv_mesh, edges, title, clim, zoom=1.3): + plotter.subplot(0, index) + plotter.set_background("white") + plotter.add_mesh(pv_mesh, scalars="speed", cmap="RdBu_r", clim=clim, + show_edges=False, lighting=False, show_scalar_bar=False) + plotter.add_mesh(edges, color="black", line_width=0.4, lighting=False) + plotter.add_text(title, position="upper_left", font_size=12, color="black") + plotter.view_xy() + plotter.camera.zoom(zoom) + + +def main(): + facet_pv, facet_edges, exact, facet_v = speed_mesh("penalty") + node_pv, node_edges, _exact, node_v = speed_mesh("penalty_node") + + # The exact field on the node-normal run's mesh, which is the same mesh. + truth = uw.discretisation.MeshVariable("Utruth", node_v.mesh, node_v.mesh.dim, + degree=2) + truth.data[:] = _exact.evaluate("velocity", truth.coords) + truth_pv = vis.meshVariable_to_pv_mesh_object(truth) + truth_pv.point_data["speed"] = np.linalg.norm( + np.squeeze(np.asarray(truth.array)), axis=1) + + top = float(truth_pv.point_data["speed"].max()) + clim = (0.0, top) + print("colour scale 0 to %.4e" % top) + for name, pv_mesh in (("exact", truth_pv), ("facet", facet_pv), ("node", node_pv)): + print("%-6s max speed %.4e" % (name, float(pv_mesh.point_data["speed"].max()))) + + FIGURES.mkdir(exist_ok=True) + + plotter = pv.Plotter(off_screen=True, shape=(1, 3), window_size=(1650, 620)) + panel(plotter, 0, truth_pv, node_edges, "exact", clim) + panel(plotter, 1, facet_pv, facet_edges, "penalty, facet normal", clim) + panel(plotter, 2, node_pv, node_edges, "penalty, node normal", clim) + plotter.screenshot(str(FIGURES / "locking.png")) + plotter.close() + + banner = pv.Plotter(off_screen=True, shape=(1, 2), window_size=(1600, 560)) + panel(banner, 0, facet_pv, facet_edges, "", clim, zoom=1.9) + panel(banner, 1, truth_pv, node_edges, "", clim, zoom=1.9) + banner.screenshot(str(FIGURES / "banner.png")) + banner.close() + print("wrote", FIGURES / "locking.png", "and", FIGURES / "banner.png") + + +if __name__ == "__main__": + main() diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-rotated-basis.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-rotated-basis.py new file mode 100644 index 0000000..bfa3a74 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-rotated-basis.py @@ -0,0 +1,102 @@ +"""Geometry for the rotated-boundary-conditions figure. + +The cetz skill's rule: geometry computation happens in Python and Typst just +draws. Doing the triangulation here rather than in the figure is what stops +nodes being left out of the mesh -- an earlier version connected nodes by a +distance threshold and silently missed several. + +Emits the schema in the skill's `underworld-bridge.md`: + + {"vertices": [[x, y], ...], + "triangles": [[i, j, k], ...], + "surface": [i, ...], indices of the constrained nodes + "frames": [{"p": [x, y], "n": [nx, ny], "t": [tx, ty]}, ...], + "curve": [[x, y], ...]} the surface, finely sampled + +Run: python3 generate-rotated-basis.py +""" +import json +import pathlib + +import numpy as np +from scipy.spatial import Delaunay + +OUT = pathlib.Path(__file__).with_name("rotated-basis-data.json") + +X0, X1 = -3.5, 3.5 +BASE = -2.9 +NX = 9 # columns of nodes +NY = 4 # rows, surface included + + +def surface_y(x): + """A deformed surface: rises on the left, falls on the right, with an + inflection between, so the normal swings through a wide range and the + curvature changes sign. No global rotation straightens this out, which is + the reason the figure exists.""" + x = np.asarray(x, dtype=float) + return (0.95 * np.exp(-(((x + 1.75) / 1.30) ** 2)) + - 0.80 * np.exp(-(((x - 1.70) / 1.15) ** 2)) + + 0.75) + + +def surface_slope(x, eps=1.0e-4): + return (surface_y(x + eps) - surface_y(x - eps)) / (2 * eps) + + +# Nodes on a grid warped to sit under the surface. Columns are staggered on +# alternate rows so the Delaunay triangulation comes out as triangles rather +# than as near-degenerate right angles on a perfect lattice. +xs = np.linspace(X0, X1, NX) +pts = [] +surface_idx = [] +for row in range(NY): + frac = row / (NY - 1) # 0 at the base, 1 at the surface + offset = 0.0 if row % 2 == 0 else 0.5 * (xs[1] - xs[0]) + cols = xs + offset + if row == NY - 1: + cols = xs # surface row unstaggered + for x in cols: + if x < X0 - 1e-9 or x > X1 + 1e-9: + continue + top = float(surface_y(x)) + y = BASE + (top - BASE) * frac + if row == NY - 1: + surface_idx.append(len(pts)) + pts.append([float(x), float(y)]) + +pts = np.array(pts) +tri = Delaunay(pts) + +# Drop the slivers Delaunay leaves along a non-convex top edge: any triangle +# whose centroid sits above the surface is outside the domain. +keep = [] +for simplex in tri.simplices: + c = pts[simplex].mean(axis=0) + if c[1] <= float(surface_y(c[0])) + 1.0e-9: + keep.append([int(i) for i in simplex]) + +frames = [] +for i in surface_idx: + x = pts[i][0] + m = float(surface_slope(x)) + n = np.array([-m, 1.0]) + n /= np.linalg.norm(n) + t = np.array([1.0, m]) + t /= np.linalg.norm(t) + frames.append({"p": [float(pts[i][0]), float(pts[i][1])], + "n": [float(n[0]), float(n[1])], + "t": [float(t[0]), float(t[1])]}) + +curve_x = np.linspace(X0, X1, 121) +data = { + "vertices": [[float(a), float(b)] for a, b in pts], + "triangles": keep, + "surface": [int(i) for i in surface_idx], + "frames": frames, + "curve": [[float(a), float(b)] for a, b in zip(curve_x, surface_y(curve_x))], +} +OUT.write_text(json.dumps(data, indent=1)) +print("wrote %s: %d vertices, %d triangles, %d surface nodes" + % (OUT.name, len(data["vertices"]), len(data["triangles"]), + len(data["surface"]))) diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-topography-figure.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-topography-figure.py new file mode 100644 index 0000000..e815024 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/generate-topography-figure.py @@ -0,0 +1,234 @@ +"""Figure: the surface topography each treatment predicts, against the exact one. + +SolCx publishes the exact dynamic topography on the top wall +(`uw.analytic.SolCx.topography_top`, which is -sigma_zz), so the comparison can +be drawn rather than tabulated. Two panels, the same five curves in each, at a +viscosity contrast of 100 and of a million. + +The topography is mean-removed in every case: the box is enclosed, so the +pressure and with it the level of sigma_zz is fixed only up to a constant, and +the deviation is the part that is determined and the part topography is built +from. + +Three routes are drawn, and they are not the same measurement: + + * the recovered traction, projected out of the solved velocity and pressure, + which is the only route Nitsche has; + * the constraint reaction -- `boundary_normal_traction` for the rotated + constraint and the multiplier field for the constraint method -- which the + solve returns as an unknown. Its sign convention is the traction holding the + wall, which is the topography's sign directly; + * the penalty's own term, kappa (u.n), which is the traction that condition + holds the wall with. Drawn rather than the projection off the same solve + because it is what the note recommends reading; the two agree to about one + per cent, which is invisible here. + +Colour carries the treatment and the exact answer is a thick grey line behind +everything, so no curve is identified by colour alone against the reference. + + python3 generate-topography-figure.py + +Writes figures/topography.png. Run against underworld3 `development` at commit +`8b7c8b9e`. +""" +import json +import pathlib +import sys + +import numpy as np + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +import solcx as C + +HERE = pathlib.Path(__file__).resolve().parent +OUT = HERE.parent / "figures" / "topography.png" +# The solves take about ten minutes and the figure is redrawn far more often +# than it is recomputed, so the curves are cached beside the script the way +# `rotated-basis-data.json` is. Delete it, or pass --recompute, to re-solve. +CACHE = HERE / "topography-data.json" + +RES = 32 +CONTRASTS = (1.0e2, 1.0e6) +MODES = ("dirichlet", "penalty", "nitsche", "constraint", "rotated") + +# Validated categorical palette, as in the multigrid note's figure. +COLOUR = { + "dirichlet": "#0b0b0b", + "penalty": "#eb6834", + "nitsche": "#c13ec1", + "constraint": "#2a78d6", + "constraint+r": "#2a78d6", + "rotated": "#1e9e6a", +} +LABEL = { + "dirichlet": "component Dirichlet", + "penalty": "penalty, $10^4$", + "nitsche": r"Nitsche, $\gamma = 10$", + "constraint": r"multiplier field $\lambda$", + "constraint+r": r"traction $\lambda + r(\mathbf{u}\cdot\hat{\mathbf{n}} - \tilde{u}_n)$", + "rotated": "rotated (reaction)", +} +STYLE = {"dirichlet": (0, (4, 2)), "penalty": "-", "nitsche": "-", + "constraint": "-", "constraint+r": (0, (1, 1.6)), "rotated": "-"} +INK = "#0b0b0b" +INK_MUTED = "#52514e" +GRID = "#e4e3df" +EXACT = "#b9b7b2" + + +def profile(mode, eta_B): + """{name: (x, topography)} along the top wall, mean removed, and the exact one. + + The reaction where there is one, the recovered traction otherwise: that is + what a user of each treatment would actually have. The constraint method + returns two curves -- the multiplier as the API returns it, and the + multiplier plus the augmented-Lagrangian share r(u.n - g), which is the + other half of the traction the momentum row carries (underworld3#607). + """ + mesh, stokes, v, exact = C.build(mode, res=RES, eta_B=eta_B) + stokes.solve() + if not C.converged(stokes): + print("%-10s eta_B %.0e diverged" % (mode, eta_B), flush=True) + return None + read = C.reaction_traction(stokes, mode, v=v) + if read is None: + coords, values = C.recovered_traction(mesh, stokes) + curves = {mode: -np.asarray(values)} # h = -sigma_zz + else: + coords, values = read + curves = {mode: np.asarray(values)} + if mode == "constraint": + # BOTH curves come from the solver, and neither is assembled here. That + # matters: this script used to add r(u.n) to what `reaction_traction` + # returned, which was right while that returned the bare multiplier and + # became a DOUBLE COUNT the moment it returned `traction()` instead -- + # the corrected curve drew at twice its augmentation share and left the + # panel. Two copies of one expression, one of them stale. There is now + # one copy, and it lives in the solver. + # multiplier() -> lambda, the field + # traction() -> lambda + r(u.n - u~_n), the whole boundary load + curves["constraint+r"] = curves[mode] # traction(), as read + coords_bare, bare = C.trace(stokes, 2, stokes.multiplier("Top")) + assert np.allclose(coords_bare, coords), "the two traces disagree" + curves[mode] = np.asarray(bare) # the multiplier alone + + order = np.argsort(coords[:, 0]) + x = coords[order, 0] + truth = exact.topography_top(coords)[order] + truth = truth - truth.mean() + out = {} + for name, values in curves.items(): + got = values[order] - values.mean() + # The sign convention is checked rather than assumed, and a curve that + # comes back anti-correlated is drawn AND named rather than flipped. + correlation = float(np.dot(got, truth) + / (np.linalg.norm(got) * np.linalg.norm(truth) + 1e-300)) + print("%-14s eta_B %.0e max|h| %.4f corr %+.3f l2 %.3f" + % (name, eta_B, np.abs(got).max(), correlation, + np.linalg.norm(got - truth) / np.linalg.norm(truth)), flush=True) + out[name] = (x, got) + return out, (x, truth) + + +def panel(ax, results, truth, title, ylim): + x, exact = truth + ax.plot(x, exact, color=EXACT, linewidth=5.0, solid_capstyle="round", + zorder=1, label="exact") + for mode, (xs, got) in results.items(): + ax.plot(xs, got, linestyle=STYLE[mode], color=COLOUR[mode], + linewidth=1.6, zorder=3, label=LABEL[mode]) + ax.axvline(0.5, color=GRID, linewidth=1.0, zorder=0) + ax.text(0.505, ylim[1] * 0.80, "viscosity step", fontsize=8.5, + color=INK_MUTED, ha="left", va="top") + ax.set_xlabel("$x$ along the top wall", fontsize=9.5, color=INK_MUTED) + ax.set_title(title, fontsize=10.5, color=INK, pad=10) + ax.set_xlim(0.0, 1.0) + ax.set_ylim(*ylim) + ax.grid(True, which="major", color=GRID, linewidth=0.8, zorder=0) + ax.set_axisbelow(True) + for side in ("top", "right"): + ax.spines[side].set_visible(False) + for side in ("left", "bottom"): + ax.spines[side].set_color(GRID) + ax.tick_params(colors=INK_MUTED, labelsize=9) + + +def compute(): + data, truth = {}, {} + for eta_B in CONTRASTS: + got = {} + for mode in MODES: + out = profile(mode, eta_B) + if out is None: + continue + curves, truth[eta_B] = out + got.update(curves) + data[eta_B] = got + CACHE.write_text(json.dumps( + {"contrasts": list(CONTRASTS), + "truth": {str(k): [v[0].tolist(), v[1].tolist()] for k, v in truth.items()}, + "data": {str(k): {m: [c[0].tolist(), c[1].tolist()] for m, c in v.items()} + for k, v in data.items()}})) + return data, truth + + +def cached(): + raw = json.loads(CACHE.read_text()) + truth = {float(k): (np.array(v[0]), np.array(v[1])) + for k, v in raw["truth"].items()} + data = {float(k): {m: (np.array(c[0]), np.array(c[1])) for m, c in v.items()} + for k, v in raw["data"].items()} + return data, truth + + +def main(): + if CACHE.exists() and "--recompute" not in sys.argv: + data, truth = cached() + print("drawn from", CACHE.name, "-- pass --recompute to re-solve") + else: + data, truth = compute() + + span = max(np.abs(truth[e][1]).max() for e in CONTRASTS) + ylim = (-1.6 * span, 1.6 * span) + + fig, axes = plt.subplots(1, 2, figsize=(9.6, 4.3), sharey=True) + fig.patch.set_facecolor("white") + for ax, eta_B in zip(axes, CONTRASTS): + panel(ax, data[eta_B], truth[eta_B], + r"$\eta_B/\eta_A = 10^{%d}$" % int(round(np.log10(eta_B))), ylim) + # Name what left the panel rather than leaving a curve to run off it. + for mode, (xs, got) in data[eta_B].items(): + if np.abs(got).max() > ylim[1]: + ax.annotate("%s leaves the panel:\npeaks at %.2f" + % (LABEL[mode], np.abs(got).max()), + xy=(0.03, ylim[0] * 0.72), fontsize=8.5, + color=COLOUR[mode], linespacing=1.35) + missing = [m for m in MODES if m not in data[eta_B]] + if missing: + ax.annotate("does not solve: %s" % ", ".join(LABEL[m] for m in missing), + xy=(0.03, ylim[1] * 0.86), fontsize=8.5, color=INK_MUTED) + axes[0].set_ylabel("surface topography, mean removed", + fontsize=9.5, color=INK_MUTED) + # The legend goes below the panels: at 1e6 every corner of both axes has a + # curve in it, and a legend inside covered the multiplier's collapse. + # Both panels, deduplicated: the penalty only converges in the left one and + # would otherwise be an unlabelled curve. + handles, labels = [], [] + for ax in axes: + for handle, label in zip(*ax.get_legend_handles_labels()): + if label not in labels: + handles.append(handle) + labels.append(label) + fig.legend(handles, labels, frameon=False, fontsize=9, labelcolor=INK_MUTED, + loc="lower center", ncol=3, handlelength=2.6, + bbox_to_anchor=(0.5, -0.01)) + fig.tight_layout(rect=(0, 0.15, 1, 1)) + fig.savefig(OUT, dpi=200, facecolor="white") + print("wrote", OUT) + + +if __name__ == "__main__": + main() diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/leak.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/leak.py new file mode 100644 index 0000000..df43248 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/leak.py @@ -0,0 +1,212 @@ +"""How much flow does each boundary treatment let through? + +The note claims a direct penalty and Nitsche satisfy `u.n = 0` only to the +accuracy of the discretisation, while rotating the degrees of freedom satisfies +it to the accuracy of the arithmetic. This measures that rather than asserting +it. + +The test is an annulus -- a boundary with no preferred direction, which is the +whole point -- driven by a single-wavenumber density anomaly, free slip on both +radii. + +TWO leaks are reported, and the difference between them is the subject of the +note's section on which normal to use: + + * against the FACET normal, which is what the discrete constraint actually + imposes; + * against the TRUE radial direction of the circle the mesh approximates. + +A strongly imposed constraint should be at machine precision against the first +and at the faceting error against the second. A weakly imposed one is limited +by its own penalty long before either. + +Reproduces the table in the note: + + python3 leak.py sweep # the resolution table + python3 leak.py free # the control: outer boundary left natural + +Run against underworld3 `bugfix/multiplier-traction` (PR #617); the +constrained solver's `traction()` is the fix this note prompted. +""" +import sys + +import numpy as np +import sympy + +import underworld3 as uw + +CELL = 0.075 +GAMMA = 10.0 + + +def build(mode, cell=CELL, penalty=1.0e4, gamma=GAMMA): + mesh = uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, cellSize=cell) + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + solver_class = (uw.systems.Stokes_Constrained if mode == "constraint" + else uw.systems.Stokes) + stokes = solver_class(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + + # a degree-4 density anomaly: enough structure that the boundary has to work + theta = sympy.atan2(y, x) + stokes.bodyforce = sympy.Matrix( + [sympy.cos(4 * theta) * x / r, sympy.cos(4 * theta) * y / r]) + + # The INNER boundary is no-slip throughout. Two free-slip circles leave + # the rigid rotation unconstrained, and the resulting nullspace is purely + # TANGENTIAL -- so a radial leak metric reads zero on a diverged solve and + # cannot tell a working constraint from a broken one. Measured: the + # unconstrained control diverged with |u| = 2.7e5 and still reported a + # radial leak of 2e-14. + stokes.add_dirichlet_bc((0.0, 0.0), "Lower") + + for boundary in ("Upper",): + if mode == "free": + continue # the control: outer boundary left natural + if mode == "constraint": + # A multiplier field h coupled into the saddle-point system, so the + # constraint is a ROW of the system rather than a term added to + # one. At convergence h on the boundary is the normal traction, + # which is why this one is also a way of getting the stress. + stokes.add_constraint_bc(0.0, boundary) + elif mode == "rotated": + stokes.add_rotated_freeslip_bc(0.0, boundary) + elif mode == "nitsche": + stokes.add_nitsche_bc(0.0, boundary, gamma=gamma, theta=1) + elif mode == "penalty": + # A DIRECT penalty: a boundary traction opposing normal flow, and + # nothing else -- no consistency term, which is exactly what leaves + # it consistent only in the limit. The documented form, from + # docs/advanced/curved-boundary-conditions.md, uses the + # quadrature-point facet normal mesh.Gamma and a POSITIVE + # coefficient. A negative one is anti-damping and the linear solve + # fails immediately, which is how this was got wrong the first time. + # + # READ THIS COLUMN WITH stress.py BESIDE IT. Imposed facet by facet + # on a curved boundary, this constraint LOCKS: the leak falls + # because the boundary is being frozen, not because the condition is + # being satisfied in the way that was meant. At cell 0.075 the + # velocity field here differs from the rotated one by 20% in l2 at + # a coefficient of 1e4, and stress.py measures the same thing + # against an exact answer. + G = mesh.Gamma + stokes.add_natural_bc(penalty * G.dot(v.sym) * G, boundary) + elif mode == "penalty_node": + # The same penalty against the measure-weighted NODE normal, which + # is one direction per node rather than one per facet, and does not + # lock. The only difference between this and the line above is which + # normal. + G = mesh.boundary_normal(boundary) + stokes.add_natural_bc(penalty * G.dot(v.sym) * G, boundary) + else: + raise ValueError(mode) + return mesh, stokes, v + + +def converged(stokes): + """A diverged solve still leaves numbers in the array, and they look like + measurements. Two runs in the first parameter sweep here had failed the + line search and were about to be tabulated.""" + return stokes.snes.getConvergedReason() > 0 + + +def leaks(mesh, v): + """max |u.n| on the outer boundary, against the facet normal and against + the true radial direction, both normalised by the flow speed.""" + coords = v.coords + rad = np.linalg.norm(coords, axis=1) + on_outer = np.abs(rad - 1.0) < 1.0e-6 + # squeeze: a vector MeshVariable's .array is (N, 1, dim), and the middle + # axis broadcasts SILENTLY against an (N, dim) array of normals, giving a + # projection of ~1e-16 for a velocity of ~1e-2. It does not raise. + allu = np.squeeze(np.asarray(v.array)) + u = allu[on_outer] + xy = coords[on_outer] + speed = np.linalg.norm(allu, axis=1).max() + + # true normal of the circle the mesh approximates + n_true = xy / np.linalg.norm(xy, axis=1)[:, None] + leak_true = np.abs((u * n_true).sum(axis=1)).max() / speed + return leak_true, on_outer.sum(), speed + + +def sweep(cells=(0.15, 0.10, 0.075, 0.05), + modes=("penalty", "penalty_node", "nitsche", "constraint", "rotated")): + """Does the leak fall with the mesh, or is it already at round-off? + + A weakly imposed constraint is satisfied to the accuracy of the + DISCRETISATION, so its leak should fall as the mesh is refined. A strongly + imposed one is satisfied to the accuracy of the ARITHMETIC, so its leak + should sit at round-off and stay there. That difference is the claim, and + the refinement is what tells them apart -- a single resolution cannot. + """ + print("| cell size | %s |" % " | ".join(m for m in modes)) + print("|---|" + "---|" * len(modes)) + for cell in cells: + row = [] + for mode in modes: + mesh, stokes, v = build(mode, cell=cell) + stokes.solve() + row.append(("%.2e" % leaks(mesh, v)[0]) if converged(stokes) + else "diverged") + print("| %.3f | %s |" % (cell, " | ".join(row))) + + +def parameter_sweep(): + """The distinction between a direct penalty and Nitsche, measured. + + Both are bounded above by conditioning: the penalty stops improving past + 1e4 and fails at 1e6, and Nitsche's line search fails from gamma = 1e4. + Neither escapes tuning. What differs is the floor each reaches first -- + 1e-3 for the penalty, 3e-5 for Nitsche -- and that Nitsche is ALSO bounded + below, at gamma = 1, where the form stops being coercive. Its usable window + has a threshold at each end, and gamma = 10 sits in the middle of it on any + mesh because gamma is dimensionless and the term already carries mu/h. + """ + print("\ndirect penalty (FACET normal): leak against penalty magnitude") + print("\n| penalty | leak/|u| |") + print("|---|---|") + for pen in (1.0e2, 1.0e3, 1.0e4, 1.0e5, 1.0e6): + mesh, stokes, v = build("penalty", penalty=pen) + stokes.solve() + cell = ("%.2e" % leaks(mesh, v)[0]) if converged(stokes) else "diverged" + print("| %.0e | %s |" % (pen, cell)) + + print("\ndirect penalty (NODE normal): leak against penalty magnitude") + print("\n| penalty | leak/|u| |") + print("|---|---|") + for pen in (1.0e2, 1.0e3, 1.0e4, 1.0e5, 1.0e6): + mesh, stokes, v = build("penalty_node", penalty=pen) + stokes.solve() + cell = ("%.2e" % leaks(mesh, v)[0]) if converged(stokes) else "diverged" + print("| %.0e | %s |" % (pen, cell)) + + print("\nNitsche: leak against gamma") + print("\n| gamma | leak/|u| |") + print("|---|---|") + for g in (1.0, 10.0, 100.0, 1000.0, 1.0e4, 1.0e5): + mesh, stokes, v = build("nitsche", gamma=g) + stokes.solve() + cell = ("%.2e" % leaks(mesh, v)[0]) if converged(stokes) else "diverged" + print("| %g | %s |" % (g, cell)) + + +if __name__ == "__main__": + if sys.argv[1:2] == ["sweep"]: + sweep() + elif sys.argv[1:2] == ["params"]: + parameter_sweep() + else: + modes = sys.argv[1:] or ["free", "penalty", "penalty_node", "nitsche", + "constraint", "rotated"] + print("%-9s %10s %10s %8s" % ("mode", "leak/|u|", "|u|max", "nodes")) + for mode in modes: + mesh, stokes, v = build(mode) + stokes.solve() + lk, n, speed = leaks(mesh, v) + print("%-9s %10.3e %10.3e %8d" % (mode, lk, speed, n)) diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated-basis-data.json b/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated-basis-data.json new file mode 100644 index 0000000..76c066e --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated-basis-data.json @@ -0,0 +1,1011 @@ +{ + "vertices": [ + [ + -3.5, + -2.9 + ], + [ + -2.625, + -2.9 + ], + [ + -1.75, + -2.9 + ], + [ + -0.875, + -2.9 + ], + [ + 0.0, + -2.9 + ], + [ + 0.875, + -2.9 + ], + [ + 1.75, + -2.9 + ], + [ + 2.625, + -2.9 + ], + [ + 3.5, + -2.9 + ], + [ + -3.0625, + -1.5690676390839076 + ], + [ + -2.1875, + -1.4005781953894207 + ], + [ + -1.3125, + -1.4008544079276464 + ], + [ + -0.4375, + -1.5774930208125677 + ], + [ + 0.4375, + -1.744573187635388 + ], + [ + 1.3125, + -1.9201468750186008 + ], + [ + 2.1875, + -1.9061055410007908 + ], + [ + 3.0625, + -1.7488486663272358 + ], + [ + -3.5, + -0.3632396002330607 + ], + [ + -2.625, + -0.06405866884136113 + ], + [ + -1.75, + 0.1666008481044865 + ], + [ + -0.875, + -0.0676029423016371 + ], + [ + 0.0, + -0.4232124172195797 + ], + [ + 0.875, + -0.7747080243853817 + ], + [ + 1.75, + -0.9985423178094506 + ], + [ + 2.625, + -0.7459268346817169 + ], + [ + 3.5, + -0.5126942107107557 + ], + [ + -3.5, + 0.9051405996504092 + ], + [ + -2.625, + 1.3539119967379585 + ], + [ + -1.75, + 1.6999012721567301 + ], + [ + -0.875, + 1.3485955865475447 + ], + [ + 0.0, + 0.8151813741706304 + ], + [ + 0.875, + 0.2879379634219279 + ], + [ + 1.75, + -0.047813476714175795 + ], + [ + 2.625, + 0.331109747977425 + ], + [ + 3.5, + 0.6809586839338668 + ] + ], + "triangles": [ + [ + 31, + 30, + 21 + ], + [ + 30, + 20, + 21 + ], + [ + 26, + 17, + 18 + ], + [ + 27, + 26, + 18 + ], + [ + 24, + 25, + 33 + ], + [ + 25, + 34, + 33 + ], + [ + 20, + 12, + 21 + ], + [ + 12, + 20, + 11 + ], + [ + 16, + 7, + 8 + ], + [ + 25, + 16, + 8 + ], + [ + 16, + 25, + 24 + ], + [ + 32, + 24, + 33 + ], + [ + 20, + 19, + 11 + ], + [ + 27, + 19, + 28 + ], + [ + 19, + 27, + 18 + ], + [ + 29, + 20, + 30 + ], + [ + 19, + 29, + 28 + ], + [ + 29, + 19, + 20 + ], + [ + 17, + 9, + 18 + ], + [ + 9, + 17, + 0 + ], + [ + 13, + 12, + 4 + ], + [ + 5, + 13, + 4 + ], + [ + 13, + 5, + 14 + ], + [ + 12, + 13, + 21 + ], + [ + 12, + 3, + 4 + ], + [ + 3, + 12, + 11 + ], + [ + 2, + 3, + 11 + ], + [ + 16, + 15, + 7 + ], + [ + 15, + 16, + 24 + ], + [ + 5, + 6, + 14 + ], + [ + 15, + 6, + 7 + ], + [ + 6, + 15, + 14 + ], + [ + 32, + 23, + 24 + ], + [ + 23, + 15, + 24 + ], + [ + 15, + 23, + 14 + ], + [ + 9, + 10, + 18 + ], + [ + 10, + 19, + 18 + ], + [ + 19, + 10, + 11 + ], + [ + 10, + 2, + 11 + ], + [ + 1, + 9, + 0 + ], + [ + 10, + 1, + 2 + ], + [ + 1, + 10, + 9 + ], + [ + 22, + 32, + 31 + ], + [ + 22, + 23, + 32 + ], + [ + 22, + 31, + 21 + ], + [ + 13, + 22, + 21 + ], + [ + 22, + 13, + 14 + ], + [ + 23, + 22, + 14 + ] + ], + "surface": [ + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 34 + ], + "frames": [ + { + "p": [ + -3.5, + 0.9051405996504092 + ], + "n": [ + -0.3058957249037569, + 0.9520650216690062 + ], + "t": [ + 0.9520650216690062, + 0.3058957249037569 + ] + }, + { + "p": [ + -2.625, + 1.3539119967379585 + ], + "n": [ + -0.5302120052673807, + 0.8478651009862024 + ], + "t": [ + 0.8478651009862024, + 0.5302120052673807 + ] + }, + { + "p": [ + -1.75, + 1.6999012721567301 + ], + "n": [ + 0.0005151017428689071, + 0.9999998673350886 + ], + "t": [ + 0.9999998673350886, + -0.0005151017428689071 + ] + }, + { + "p": [ + -0.875, + 1.3485955865475447 + ], + "n": [ + 0.5426587927926065, + 0.8399532335820077 + ], + "t": [ + 0.8399532335820077, + -0.5426587927926065 + ] + }, + { + "p": [ + 0.0, + 0.8151813741706304 + ], + "n": [ + 0.4836463257494078, + 0.8752635212260921 + ], + "t": [ + 0.8752635212260921, + -0.4836463257494078 + ] + }, + { + "p": [ + 0.875, + 0.2879379634219279 + ], + "n": [ + 0.5429867639410579, + 0.8397412543068359 + ], + "t": [ + 0.8397412543068358, + -0.5429867639410579 + ] + }, + { + "p": [ + 1.75, + -0.047813476714175795 + ], + "n": [ + -0.0574834315887112, + 0.9983464604498711 + ], + "t": [ + 0.9983464604498711, + 0.0574834315887112 + ] + }, + { + "p": [ + 2.625, + 0.331109747977425 + ], + "n": [ + -0.5055403695426544, + 0.8628029524536158 + ], + "t": [ + 0.8628029524536158, + 0.5055403695426544 + ] + }, + { + "p": [ + 3.5, + 0.6809586839338668 + ], + "n": [ + -0.1847046296744634, + 0.982794078012693 + ], + "t": [ + 0.982794078012693, + 0.1847046296744634 + ] + } + ], + "curve": [ + [ + -3.5, + 0.9051405996504089 + ], + [ + -3.441666666666667, + 0.924709919910546 + ], + [ + -3.3833333333333333, + 0.9459569983210258 + ], + [ + -3.325, + 0.9689047102584147 + ], + [ + -3.2666666666666666, + 0.9935569625066398 + ], + [ + -3.2083333333333335, + 1.019896408691709 + ], + [ + -3.15, + 1.0478823531124308 + ], + [ + -3.091666666666667, + 1.077448917340561 + ], + [ + -3.033333333333333, + 1.1085035443807112 + ], + [ + -2.975, + 1.1409259131181626 + ], + [ + -2.9166666666666665, + 1.1745673310574358 + ], + [ + -2.8583333333333334, + 1.2092506658746796 + ], + [ + -2.8, + 1.2447708660876644 + ], + [ + -2.7416666666666667, + 1.280896108319149 + ], + [ + -2.6833333333333336, + 1.31736959344396 + ], + [ + -2.625, + 1.353911996737958 + ], + [ + -2.5666666666666664, + 1.3902245584714161 + ], + [ + -2.5083333333333333, + 1.425992781792527 + ], + [ + -2.45, + 1.4608906848912475 + ], + [ + -2.3916666666666666, + 1.4945855350373591 + ], + [ + -2.333333333333333, + 1.526742973893888 + ], + [ + -2.275, + 1.5570324272567517 + ], + [ + -2.216666666666667, + 1.5851326787631834 + ], + [ + -2.158333333333333, + 1.6107374767717508 + ], + [ + -2.1, + 1.633561037067616 + ], + [ + -2.041666666666667, + 1.6533433016773964 + ], + [ + -1.9833333333333334, + 1.6698548161227493 + ], + [ + -1.925, + 1.6829010939636901 + ], + [ + -1.8666666666666667, + 1.6923263483659823 + ], + [ + -1.8083333333333333, + 1.6980164853777744 + ], + [ + -1.75, + 1.6999012721567306 + ], + [ + -1.6916666666666667, + 1.6979556149374022 + ], + [ + -1.6333333333333333, + 1.692199905332323 + ], + [ + -1.575, + 1.682699418789918 + ], + [ + -1.5166666666666666, + 1.6695627748036432 + ], + [ + -1.4583333333333335, + 1.6529394938808053 + ], + [ + -1.4, + 1.6330167104634032 + ], + [ + -1.3416666666666668, + 1.61001512313947 + ], + [ + -1.2833333333333332, + 1.5841842828851125 + ], + [ + -1.225, + 1.5557973361588586 + ], + [ + -1.1666666666666665, + 1.5251453520084297 + ], + [ + -1.1083333333333334, + 1.492531370690367 + ], + [ + -1.0499999999999998, + 1.4582643155616355 + ], + [ + -0.9916666666666667, + 1.4226529102631265 + ], + [ + -0.9333333333333331, + 1.3859997397146477 + ], + [ + -0.875, + 1.348595586547545 + ], + [ + -0.8166666666666664, + 1.3107141647846698 + ], + [ + -0.7583333333333333, + 1.272607360377266 + ], + [ + -0.7000000000000002, + 1.2345010741961506 + ], + [ + -0.6416666666666666, + 1.1965917478187014 + ], + [ + -0.5833333333333335, + 1.1590436364841392 + ], + [ + -0.5249999999999999, + 1.1219868773708934 + ], + [ + -0.4666666666666668, + 1.0855163852548482 + ], + [ + -0.4083333333333332, + 1.0496915919046683 + ], + [ + -0.3500000000000001, + 1.0145370304190584 + ], + [ + -0.2916666666666665, + 0.9800437511642502 + ], + [ + -0.2333333333333334, + 0.9461715419910942 + ], + [ + -0.17499999999999982, + 0.9128519118949112 + ], + [ + -0.1166666666666667, + 0.8799917840856701 + ], + [ + -0.058333333333333126, + 0.8474778314169195 + ], + [ + 0.0, + 0.8151813741706305 + ], + [ + 0.05833333333333357, + 0.7829637472749286 + ], + [ + 0.1166666666666667, + 0.7506820312100971 + ], + [ + 0.17500000000000027, + 0.7181950283314438 + ], + [ + 0.2333333333333334, + 0.6853693544451201 + ], + [ + 0.2916666666666665, + 0.6520855046989421 + ], + [ + 0.3500000000000001, + 0.6182437438110678 + ], + [ + 0.4083333333333332, + 0.5837696640761547 + ], + [ + 0.4666666666666668, + 0.548619251245789 + ], + [ + 0.5250000000000004, + 0.5127832990715608 + ], + [ + 0.583333333333333, + 0.4762910187656371 + ], + [ + 0.6416666666666666, + 0.4392127004940097 + ], + [ + 0.7000000000000002, + 0.4016613007004925 + ], + [ + 0.7583333333333337, + 0.363792851740313 + ], + [ + 0.8166666666666664, + 0.325805618850044 + ], + [ + 0.875, + 0.28793796342192796 + ], + [ + 0.9333333333333336, + 0.2504649100526243 + ], + [ + 0.9916666666666671, + 0.2136934567143398 + ], + [ + 1.0499999999999998, + 0.17795671114661105 + ], + [ + 1.1083333333333334, + 0.14360698042499975 + ], + [ + 1.166666666666667, + 0.11100798268452128 + ], + [ + 1.2249999999999996, + 0.08052638813756408 + ], + [ + 1.2833333333333332, + 0.05252292884049681 + ], + [ + 1.3416666666666668, + 0.02734334129772975 + ], + [ + 1.4000000000000004, + 0.005309421388709379 + ], + [ + 1.458333333333333, + -0.013289523913108559 + ], + [ + 1.5166666666666666, + -0.028204549692437686 + ], + [ + 1.5750000000000002, + -0.03923401019190231 + ], + [ + 1.6333333333333337, + -0.04622932722506001 + ], + [ + 1.6916666666666664, + -0.049099395459559214 + ], + [ + 1.75, + -0.047813476714175684 + ], + [ + 1.8083333333333336, + -0.042402488429535534 + ], + [ + 1.8666666666666671, + -0.032958648592165685 + ], + [ + 1.9249999999999998, + -0.01963349805658543 + ], + [ + 1.9833333333333334, + -0.002634378770072665 + ], + [ + 2.041666666666667, + 0.017780499734667887 + ], + [ + 2.0999999999999996, + 0.04130822544601709 + ], + [ + 2.158333333333333, + 0.06760836083270272 + ], + [ + 2.216666666666667, + 0.09631088206901939 + ], + [ + 2.2750000000000004, + 0.12702461196713855 + ], + [ + 2.333333333333333, + 0.15934587047687743 + ], + [ + 2.3916666666666666, + 0.19286707036237616 + ], + [ + 2.45, + 0.22718500012101506 + ], + [ + 2.5083333333333337, + 0.26190856024390546 + ], + [ + 2.5666666666666664, + 0.2966657509787571 + ], + [ + 2.625, + 0.3311097479774251 + ], + [ + 2.6833333333333336, + 0.36492394452614507 + ], + [ + 2.741666666666667, + 0.39782588332598795 + ], + [ + 2.8, + 0.429570044919817 + ], + [ + 2.8583333333333334, + 0.4599495019114082 + ], + [ + 2.916666666666667, + 0.4887964864025283 + ], + [ + 2.9750000000000005, + 0.5159819512151704 + ], + [ + 3.033333333333333, + 0.5414142324664324 + ], + [ + 3.091666666666667, + 0.5650369413091443 + ], + [ + 3.1500000000000004, + 0.5868262259143289 + ], + [ + 3.208333333333333, + 0.6067875511844959 + ], + [ + 3.2666666666666666, + 0.6249521436990915 + ], + [ + 3.325, + 0.6413732437141221 + ], + [ + 3.3833333333333337, + 0.6561222955694307 + ], + [ + 3.4416666666666664, + 0.6692851936251155 + ], + [ + 3.5, + 0.680958683933867 + ] + ] +} \ No newline at end of file diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated-basis.typ b/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated-basis.typ new file mode 100644 index 0000000..3657d85 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated-basis.typ @@ -0,0 +1,186 @@ +// Rotated boundary conditions: where the basis changes, and where in the +// solver the rotation lives. +// +// Thesis: rotating the degrees of freedom leaves the discrete problem in a +// MIXED basis, but the obligation is CONTAINED -- it lives in the velocity +// block and its multigrid substructure, and the Schur/pressure machinery +// wrapping it never handles a rotated vector. +// +// House style follows the cetz-figures skill in the underworld3 repository +// (cetz 0.3.4, hex colours, helper-function pattern). +// +// Build, from this directory: +// +// python3 generate-rotated-basis.py # writes the JSON geometry +// typst compile --format svg rotated-basis.typ ../figures/rotated-basis.svg +// typst compile --ppi 200 rotated-basis.typ ../figures/rotated-basis.png + +#import "@preview/cetz:0.3.4" + +#set page(width: auto, height: auto, margin: 16pt) +#set text(size: 10pt) + +#cetz.canvas({ + import cetz.draw: * + + // Colours + let cart-fg = rgb("#4a7bf7") + let cart-bg = rgb("#dce8fc") + let rot-fg = rgb("#e57373") + let rot-bg = rgb("#fce4ec") + let plain-bg = rgb("#f2f2f0") + let mg-bg = rgb("#fdf1f4") + let schur-fg = rgb("#49a87c") // the un-rotated half + let schur-bg = rgb("#e8f4e8") + let ink = rgb("#1a1a1a") + let muted = rgb("#7a7a7a") + let hair = rgb("#d0d0d0") + + // ====================================================================== + // LEFT: a deformed surface, meshed. Geometry (nodes, Delaunay triangles, + // surface frames) comes from generate-rotated-basis.py via JSON -- the + // skill's rule, and the reason every node is now in the triangulation. + // An earlier version joined nodes by a distance threshold and left some out. + // ====================================================================== + let data = json("rotated-basis-data.json") + let vtx = data.vertices + let at(i) = (vtx.at(i).at(0), vtx.at(i).at(1)) + let is-surface(i) = data.surface.contains(i) + + // triangles first, then the surface, then the nodes: painter's algorithm + for tri in data.triangles { + line(at(tri.at(0)), at(tri.at(1)), at(tri.at(2)), close: true, + stroke: (paint: hair, thickness: 0.55pt)) + } + + line(..data.curve.map(q => (q.at(0), q.at(1))), + stroke: (paint: ink, thickness: 1.2pt)) + + // A rotated frame at every constrained node: n the outward surface normal, + // t the tangent. Both follow the node, which is the whole point. + for (k, f) in data.frames.enumerate() { + let p = (f.p.at(0), f.p.at(1)) + let n = f.n + let tg = f.t + let ln = 0.72 + let lt = 0.46 + line(p, (p.at(0) + ln * n.at(0), p.at(1) + ln * n.at(1)), + mark: (end: ">", scale: 0.4), stroke: (paint: rot-fg, thickness: 1.1pt)) + line(p, (p.at(0) + lt * tg.at(0), p.at(1) + lt * tg.at(1)), + mark: (end: ">", scale: 0.4), stroke: (paint: rot-fg, thickness: 1.1pt)) + if k == 4 { + content((p.at(0) + 1.20 * ln * n.at(0), p.at(1) + 1.20 * ln * n.at(1)), + text(fill: rot-fg, size: 9pt, $n$)) + content((p.at(0) + 1.55 * lt * tg.at(0) + 0.24 * n.at(0), + p.at(1) + 1.55 * lt * tg.at(1) + 0.24 * n.at(1)), + text(fill: rot-fg, size: 9pt, $t$)) + } + } + + // the Cartesian frame at one interior node -- identical at every other one, + // which is what makes the surface the odd one out + let ip = at(12) // an interior node with room around it + line(ip, (ip.at(0) + 0.66, ip.at(1)), mark: (end: ">", scale: 0.4), + stroke: (paint: cart-fg, thickness: 1.1pt)) + line(ip, (ip.at(0), ip.at(1) + 0.66), mark: (end: ">", scale: 0.4), + stroke: (paint: cart-fg, thickness: 1.1pt)) + content((ip.at(0) + 0.88, ip.at(1)), text(fill: cart-fg, size: 9pt, $x$)) + content((ip.at(0), ip.at(1) + 0.88), text(fill: cart-fg, size: 9pt, $y$)) + + for i in range(vtx.len()) { + if is-surface(i) { + circle(at(i), radius: 0.13, fill: rot-fg, + stroke: (paint: rot-fg, thickness: 1pt)) + } else { + circle(at(i), radius: 0.11, fill: cart-bg, + stroke: (paint: cart-fg, thickness: 1pt)) + } + } + + content((0, 3.05), text(weight: "bold", size: 10pt, fill: ink, + "A deformed surface has no preferred direction")) + + circle((-3.3, -3.55), radius: 0.13, fill: rot-fg, stroke: (paint: rot-fg)) + content((-3.05, -3.55), anchor: "west", + text(size: 9pt, fill: ink, [surface node --- solve for $(v_n, v_t)$, hold $v_n$])) + circle((-3.3, -4.1), radius: 0.11, fill: cart-bg, + stroke: (paint: cart-fg, thickness: 1pt)) + content((-3.05, -4.1), anchor: "west", + text(size: 9pt, fill: ink, [interior node --- solve for $(v_x, v_y)$])) + + line((4.35, -4.9), (4.35, 3.3), stroke: (paint: hair, thickness: 0.8pt)) + + // ====================================================================== + // RIGHT: where the rotation lives. Two blocks ABUTTING, not nested: the + // velocity solve is rotated, the Schur/pressure solve is not, and the + // single un-rotation on the boundary between them feeds both the Schur + // solve and everything outside. + // ====================================================================== + let panel(tl, br, title, subtitle, bg, edge) = { + import cetz.draw: * + rect(tl, br, fill: bg, stroke: (paint: edge, thickness: 1pt), radius: 5pt) + content((tl.at(0) + 0.30, tl.at(1) - 0.36), anchor: "west", + text(weight: "bold", size: 9.5pt, fill: edge, title)) + if subtitle != none { + // cetz content() lays out at natural width and spills over the rect. + // Box it to the panel's own width so the text wraps inside the border. + // The canvas default is 1cm per unit, so the arithmetic is direct. + let tw = (br.at(0) - tl.at(0) - 0.60) * 1cm + content((tl.at(0) + 0.30, tl.at(1) - 0.90), anchor: "north-west", + box(width: tw, text(size: 8.5pt, fill: muted, subtitle))) + } + } + + content((10.3, 3.05), text(weight: "bold", size: 10pt, fill: ink, + "Where the rotation lives")) + + // -- the rotated half --------------------------------------------------- + panel((4.9, 1.65), (11.75, -3.70), [Velocity solve --- rotated], + [$hat(A) = Q^T A Q$, #h(0.7em) $hat(b) = Q^T b$, #h(0.7em) + $v_n$ held strongly at surface nodes], rot-bg, rot-fg) + + panel((5.35, -0.15), (11.50, -3.25), "Multigrid", + [the transfers are the only further obligation], mg-bg, rot-fg) + + let mgrow(y, lhs, rhs) = { + import cetz.draw: * + content((5.65, y), anchor: "west", text(size: 9pt, fill: ink, lhs)) + content((8.05, y), anchor: "west", text(size: 9pt, fill: rot-fg, rhs)) + } + mgrow(-1.85, [prolongation], [$P -> Q^T P$]) + mgrow(-2.40, [coarse operators], [inherit $Q$ via $R A P$]) + mgrow(-2.93, [coarse solve], [SVD (rigid rotations)]) + + // -- the un-rotated half, abutting -------------------------------------- + panel((13.35, 1.65), (17.85, -1.35), [Fieldsplit / Schur solve], + [pressure and constraints. Isotropic, and carrying no boundary + condition of this kind, so it never sees a rotated vector.], + schur-bg, schur-fg) + + // -- one un-rotation on the boundary, feeding both ---------------------- + // The junction sits in the gap BETWEEN the two blocks, low enough to clear + // the velocity block's own subtitle -- one un-rotation, two consumers. + let jx = 12.55 + let jy = -0.40 + line((11.75, jy), (jx, jy), stroke: (paint: rot-fg, thickness: 1pt)) + circle((jx, jy), radius: 0.075, fill: rot-fg, stroke: (paint: rot-fg)) + content((jx, jy + 0.45), text(size: 9pt, fill: rot-fg, $v = Q hat(v)$)) + + // branch 1: into the Schur solve, which needs it un-rotated + line((jx, jy), (13.29, jy), mark: (end: ">", scale: 0.4), + stroke: (paint: rot-fg, thickness: 1pt)) + + // branch 2: out to everything else + line((jx, jy), (jx, -2.70), stroke: (paint: rot-fg, thickness: 1pt)) + line((jx, -2.70), (13.20, -2.70), mark: (end: ">", scale: 0.4), + stroke: (paint: rot-fg, thickness: 1pt)) + content((13.33, -2.70), anchor: "west", text(size: 8.5pt, fill: rot-fg, + [and to everything outside ---\ output, advection, the surface update])) + + // convention, stated -- a reader who assumes the transpose reads the + // whole figure backwards + line((4.75, -5.05), (18.1, -5.05), stroke: (paint: hair, thickness: 0.8pt)) + content((4.75, -5.55), anchor: "west", text(size: 9pt, fill: ink, + [Convention: the columns of $Q$ are the nodal frame, so $hat(v) = Q^T v$, + and $Q = I$ at every unconstrained node.])) +}) diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated_dof_knockout.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated_dof_knockout.py new file mode 100644 index 0000000..bed0900 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/rotated_dof_knockout.py @@ -0,0 +1,129 @@ +"""What the rotated constraint does to the MATRIX, not to the answer. + +On an axis-aligned box the per-node rotation Q is the identity, so the rotated +treatment and a component Dirichlet condition are the same arithmetic: agreeing +answers prove nothing about the machinery. Turning the frame makes Q a real +rotation, but the answers still agree -- correctly, since a rigid rotation is an +orthogonal similarity of the discrete problem. Neither run inspects what was +actually eliminated. + +This does. Two statements about the assembled system: + + * Q is orthogonal, and at 45 degrees it is nowhere near the identity; + * the converged Cartesian reaction r = A.u - b at a constrained node is + PARALLEL TO THE WALL NORMAL. That is the signature of the knockout: the row + that was struck is the wall-normal momentum balance, so the force the + constraint had to supply has no tangential part. If the elimination were + done in the coordinate frame -- the failure a rotated run cannot otherwise + distinguish -- the reaction would lie along a coordinate axis instead. + +The reaction is the solver's own `A.u - b`, stashed at convergence +(`_rotated_freeslip_info["reaction"]`); it is what `boundary_normal_traction` +reads the normal component of. + + python3 rotated_dof_knockout.py + +Run against underworld3 `bugfix/multiplier-traction` (PR #617). +""" +import sys + +import numpy as np + +import underworld3 as uw + +from underworld3.utilities import rotated_bc as RB + +import solcx_rotated as RT + +RES = 16 +CORNER = 1.0e-9 + + +def reaction_vectors(stokes, theta, boundary="Top"): + """The Cartesian nodal reaction at each node of `boundary`, corners dropped. + + Reads the same stashed residual `boundary_normal_traction` uses, but keeps + both components instead of projecting onto the normal. + """ + info = stokes._rotated_freeslip_info + dm = stokes.dm + dim = stokes.mesh.dim + rcl = dm.getLocalVec() + dm.globalToLocal(info["reaction"], rcl) + rca = np.asarray(rcl.getArray()) + lsec = dm.getLocalSection() + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + + coords, vectors = [], [] + for q, _nrm in RB._boundary_velocity_nodes(stokes, boundary, normal=None): + lo = lsec.getFieldOffset(q, RB._VELOCITY_FIELD) + x = RB._point_coord(dm, dim, cvec, csec, v0, v1, q) + coords.append(x) + vectors.append(np.array(rca[lo:lo + dim])) + dm.restoreLocalVec(rcl) + coords, vectors = np.array(coords), np.array(vectors) + + # A corner belongs to two walls and carries both reactions, so its direction + # is a sum of two normals and says nothing about either wall. + material = coords @ RT.rotation(theta) + keep = (material[:, 0] > CORNER) & (material[:, 0] < 1.0 - CORNER) + return coords[keep], vectors[keep] + + +def q_is_orthogonal(stokes): + """||Q^T Q x - x|| / ||x||, and ||Q x - x|| / ||x||. + + The second is reported by the script but is NOT the discriminator it looks + like: even on the aligned box Q reorders the local axes so the wall normal + comes first, so it is never the identity. What separates the two frames is + the DIRECTION the constraint eliminates, which the reaction shows. + """ + info = stokes._rotated_freeslip_info + Q, Qt = info["Q"], info["Qt"] + x = Q.createVecRight() + x.setRandom() + y = Q.createVecLeft() + Q.mult(x, y) + z = x.duplicate() + Qt.mult(y, z) + z.axpy(-1.0, x) + orthogonality = z.norm() / x.norm() + # ||Q x - x||: zero when Q is the identity, O(1) when it turns the frame + w = y.copy() + w.axpy(-1.0, x) + return orthogonality, w.norm() / x.norm() + + +def check(theta_degrees): + theta = np.radians(theta_degrees) + mesh, stokes, v, truth = RT.build("rotated", res=RES, theta=theta) + stokes.solve() + n = RT.rotation(theta) @ np.array([0.0, 1.0]) + t = np.array([-n[1], n[0]]) + + coords, r = reaction_vectors(stokes, theta) + magnitude = np.linalg.norm(r, axis=1) + live = magnitude > 1.0e-30 + tangential = np.abs(r[live] @ t) / magnitude[live] + along_y = np.abs(r[live] @ np.array([0.0, 1.0])) / magnitude[live] + along_x = np.abs(r[live] @ np.array([1.0, 0.0])) / magnitude[live] + orthogonality, _ = q_is_orthogonal(stokes) + return dict(nodes=int(live.sum()), tangential=float(tangential.max()), + along_x=float(along_x.mean()), along_y=float(along_y.mean()), + orthogonality=float(orthogonality)) + + +if __name__ == "__main__": + print("Rotated free slip on every wall, %d x %d. The reaction is the solver's" + % (RES, RES)) + print("own A.u - b at convergence.") + print() + print("| frame | nodes | ||Q^tQ - I|| | reaction . t | reaction . x | reaction . y |") + print("|---|---|---|---|---|---|") + for degrees in (0.0, 45.0): + got = check(degrees) + print("| %.0f° | %d | %.1e | %.1e | %.3f | %.3f |" + % (degrees, got["nodes"], got["orthogonality"], got["tangential"], + got["along_x"], got["along_y"]), flush=True) diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/solcx.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/solcx.py new file mode 100644 index 0000000..756f764 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/solcx.py @@ -0,0 +1,311 @@ +"""The other half of the stress test: a lateral viscosity contrast. + +`stress.py` is the geometry half -- a curved boundary with the rheology +trivial. This is the rheology half, with the geometry trivial. No exact +solution has both, so the two halves are separate tests and the note does not +conflate them. + +On a box, "no flow through this wall" is a single velocity component and every +treatment here reduces to the same constraint. Nothing about this test is about +normals. What it is about is whether a WEAKLY imposed constraint holds the +degrees of freedom it was given when the viscosity beside them jumps by six +orders of magnitude -- and SolCx is the only case of the two halves with an +exact surface stress to check that against. + +SolCx: unit box, viscosity eta_A left of x = 0.5 and eta_B right of it, +forcing (0, cos(pi.x) sin(pi.z)), free slip on all four walls. +`uw.analytic.SolCx` publishes the exact stress, and `topography_top` is +-sigma_zz on the top boundary, which is the quantity here. + +Three walls carry the ordinary component Dirichlet condition. The treatment +under test is on the TOP wall only, so what is measured is attributable to it. + +The metrics +----------- + * the leak, max |u.n| / |u|, reported SEPARATELY over the soft half of the + wall and the stiff half. One coefficient has to hold both, and whether it + can is the whole question; + * the surface stress: relative l2 error of sigma_zz along the top wall + against the exact one, both mean-removed (the box is enclosed, so the + pressure -- and with it the level of sigma_zz -- is fixed only up to a + constant). + + python3 solcx.py sweep # resolution, each treatment + python3 solcx.py contrast # the viscosity ratio, each treatment + python3 solcx.py params # the penalty coefficient, both halves + python3 solcx.py control # the top wall left free + +Run against underworld3 `bugfix/multiplier-traction` (PR #617); the +constrained solver's `traction()` is the fix this note prompted. +""" +import sys + +import numpy as np +import sympy + +import underworld3 as uw + +from underworld3.utilities.boundary_flux import _boundary_field_nodes + +RES = 32 +ETA_B = 1.0e6 +PENALTY = 1.0e4 +GAMMA = 10.0 + + +def build(mode, res=RES, eta_B=ETA_B, penalty=PENALTY, gamma=GAMMA): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res, res), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3) + + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + solver_class = (uw.systems.Stokes_Constrained if mode == "constraint" + else uw.systems.Stokes) + stokes = solver_class(mesh, velocityField=v, pressureField=p) + + exact = uw.analytic.SolCx(mesh, eta_A=1.0, eta_B=eta_B, x_c=0.5, n=1) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = exact.fn_viscosity + if mode != "constraint": + # Stokes_Constrained builds its own Schur preconditioner from the + # operator's blocks and refuses this assignment. + stokes.saddle_preconditioner = 1.0 / exact.fn_viscosity + stokes.bodyforce = exact.fn_bodyforce + stokes.tolerance = 1.0e-9 + # Enclosed: the constant pressure mode is a nullspace and a direct solve on + # the singular saddle returns a quiet, wrong answer without this. With the + # top wall left free (the control) the domain is open and there is no such + # mode -- asserting one there is what makes the control fail to solve + # rather than fail to hold the boundary. + stokes.petsc_use_pressure_nullspace = (mode != "free") + + stokes.add_dirichlet_bc((0.0, None), "Left") + stokes.add_dirichlet_bc((0.0, None), "Right") + stokes.add_dirichlet_bc((None, 0.0), "Bottom") + + if mode == "dirichlet": + # What a box lets you do, and the reference the others are judged by. + stokes.add_dirichlet_bc((None, 0.0), "Top") + elif mode == "constraint": + stokes.add_constraint_bc(0.0, "Top") + elif mode == "rotated": + stokes.add_rotated_freeslip_bc(0.0, "Top") + elif mode == "nitsche": + stokes.add_nitsche_bc(0.0, "Top", gamma=gamma, theta=1) + elif mode == "penalty": + # A CONSTANT coefficient, deliberately. Scaling it by the local + # viscosity is the fair thing to want here -- one number cannot be large + # against 1e6 and moderate against 1 -- and it does not solve: a + # Piecewise viscosity inside the boundary term fails the line search at + # every magnitude tried, from 1 to 1e3 times mu, against both normals. + # A constant 1e4 solves, and that is what this column is. + stokes.add_natural_bc(penalty * mesh.Gamma.dot(v.sym) * mesh.Gamma, "Top") + elif mode == "free": + pass + else: + raise ValueError(mode) + + return mesh, stokes, v, exact + + +def converged(stokes): + return stokes.snes.getConvergedReason() > 0 + + +def trace(solver, field_id, var, boundary="Top"): + """(coords, values) for `var` at the nodes it carries on `boundary`, + selected by the mesh boundary LABEL.""" + nodes, *_ = _boundary_field_nodes(solver, boundary, field_id) + coords = np.array([node[2] for node in nodes]) + tree = uw.kdtree.KDTree(np.ascontiguousarray(var.coords)) + index = np.asarray(tree.query(np.ascontiguousarray(coords), 1)[1]).flatten() + return coords, np.squeeze(np.asarray(var.array))[index] + + +def recovered_traction(mesh, stokes): + """sigma_zz on the top wall, projected out of the solved fields.""" + field = uw.discretisation.MeshVariable("Szz", mesh, 1, degree=2) + projection = uw.systems.Projection(mesh, field) + projection.uw_function = stokes.stress[1, 1] + projection.solve() + return trace(projection, 0, field) + + +def penalty_traction(stokes, v, kappa): + """The traction a direct penalty holds the wall with. + + The penalty term IS that traction, so sigma_zz is -kappa (u.n) on the top + wall, where n = (0, 1). Returned with the sign of a reaction, as the other + two routes are, so `measure` negates all three alike. + """ + coords = np.asarray(v.coords) + top = np.abs(coords[:, 1] - 1.0) < 1.0e-9 + u = np.squeeze(np.asarray(v.array)) + return coords[top], kappa * u[top, 1] + + +def reaction_traction(stokes, mode, v=None, penalty=PENALTY): + """The traction the solve returned, for the three methods that return one. + + Nitsche is the exception: its term is written in sigma(u), so reading it + back means differentiating the velocity, which is what + `recovered_traction` does. + """ + if mode == "penalty": + return penalty_traction(stokes, v, penalty) + if mode == "rotated": + return stokes.boundary_normal_traction("Top") + if mode == "constraint": + # `traction()` is h + r(u.n - g), the WHOLE boundary term. Reading the + # multiplier alone here was wrong by an order of magnitude and a sign at + # a 1e6 viscosity contrast, because the default r is viscosity-weighted + # (underworld3#607, fixed in #617). The box is flat, so evaluating the + # expression at the trace nodes is safe -- on a convex curved boundary it + # would extrapolate (#605), which is why stress.py reads arrays instead. + coords, _h = trace(stokes, 2, stokes.multiplier("Top")) + return coords, np.asarray( + uw.function.evaluate(stokes.traction("Top"), coords)).reshape(-1) + return None + + +def stress_error(coords, values, exact, trim=0.0): + """Relative l2 error along the wall, both mean-removed. + + The mean is the gauge: the box is enclosed, so the pressure and with it the + level of sigma_zz is fixed only up to a constant. The deviation is what + topography is built from and the only part that is determined. + + `trim` drops nodes within that distance of the two ends of the wall. The + corners are where the treatment under test meets the side walls' component + conditions, so a node there is constrained twice and by two different + mechanisms; trimming separates what the wall does from what the corner does. + """ + coords = np.asarray(coords) + values = np.asarray(values) + if trim > 0.0: + keep = (coords[:, 0] > trim) & (coords[:, 0] < 1.0 - trim) + coords, values = coords[keep], values[keep] + truth = -exact.topography_top(coords) # topography_top is -sigma_zz + order = np.argsort(coords[:, 0]) + got = values[order] - values.mean() + truth = truth[order] - truth.mean() + return float(np.linalg.norm(got - truth) / np.linalg.norm(truth)) + + +def leaks(v): + """max |u_z| / |u| on the top wall, over the soft half and the stiff half. + + Reported separately because one penalty coefficient has to hold both, and + the two sides do not ask the same thing of it. + """ + coords = v.coords + top = np.abs(coords[:, 1] - 1.0) < 1.0e-9 + u = np.squeeze(np.asarray(v.array)) + speed = np.linalg.norm(u, axis=1).max() + out = [] + for side in (coords[:, 0] < 0.5, coords[:, 0] > 0.5): + mask = top & side + out.append(np.abs(u[mask, 1]).max() / speed) + return out + + +def measure(mode, **kwargs): + mesh, stokes, v, exact = build(mode, **kwargs) + stokes.solve() + if not converged(stokes): + return None + soft, stiff = leaks(v) + out = {"soft": soft, "stiff": stiff, "velocity": exact.velocity_error(v)} + coords, values = recovered_traction(mesh, stokes) + out["recovered"] = stress_error(coords, values, exact) + out["trimmed"] = stress_error(coords, values, exact, trim=2.0 / kwargs.get("res", RES)) + read = reaction_traction(stokes, mode, v=v, penalty=kwargs.get("penalty", PENALTY)) + if read is not None: + coords, values = read + # The reaction is the traction holding the wall, opposite in sign to + # sigma_zz -- the same convention as in stress.py. + out["reaction"] = stress_error(coords, -np.asarray(values), exact) + return out + + +MODES = ("dirichlet", "penalty", "nitsche", "constraint", "rotated") + + +def sweep(resolutions=(16, 32, 64), modes=MODES): + print("surface stress error, eta_B/eta_A = %.0e" % ETA_B) + print() + print("| elements | " + " | ".join(modes) + " | rotated (reaction) | multiplier |") + print("|---" * (len(modes) + 3) + "|") + for res in resolutions: + row, extra = [], {} + for mode in modes: + got = measure(mode, res=res) + row.append("diverged" if got is None else "%.2e" % got["recovered"]) + if got and "reaction" in got: + extra[mode] = "%.2e" % got["reaction"] + print("| %d | %s | %s | %s |" + % (res, " | ".join(row), extra.get("rotated", "-"), + extra.get("constraint", "-")), flush=True) + + +def contrast(ratios=(1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e6), modes=MODES): + """Does the treatment survive the viscosity jump getting bigger? + + `dirichlet` is the control, and it is the one to read first. It is the + ordinary component condition a box allows, it holds u.n exactly, and its + velocity error is 1e-5 at every contrast here -- so whatever it reads in + the stress column is the RECOVERY's error and not a boundary condition's. + A treatment can only be said to be worse than the reference where it is + worse than that. + """ + print("whole wall / trimmed by two elements at each end") + print() + print("| eta_B/eta_A | " + " | ".join(modes) + " |") + print("|---" * (len(modes) + 1) + "|") + for ratio in ratios: + row = [] + for mode in modes: + got = measure(mode, eta_B=ratio) + row.append("diverged" if got is None + else "%.2e / %.2e" % (got["recovered"], got["trimmed"])) + print("| %.0e | %s |" % (ratio, " | ".join(row)), flush=True) + + +def parameters(): + """The leak on each half of the wall, against the parameter. + + A dimensionless gamma carries mu/h with it and so asks the same of both + halves. A bare penalty coefficient does not, and the stiff half is where + that shows. + """ + for mode, values, label in ( + ("penalty", (1e2, 1e4, 1e6, 1e8), "penalty coefficient"), + ("nitsche", (1.0, 10.0, 100.0, 1000.0), "Nitsche gamma")): + print("\n%s, eta_B/eta_A = %.0e" % (label, ETA_B)) + print("\n| %s | leak, soft half | leak, stiff half | stress error |" % label) + print("|---|---|---|---|") + for value in values: + kwargs = {"gamma": value} if mode == "nitsche" else {"penalty": value} + got = measure(mode, **kwargs) + if got is None: + print("| %g | diverged | | |" % value) + continue + print("| %g | %.2e | %.2e | %.2e |" + % (value, got["soft"], got["stiff"], got["recovered"]), flush=True) + + +def control(): + for mode in ("free", "dirichlet"): + got = measure(mode) + if got is None: + print("%-10s diverged" % mode, flush=True) + continue + print("%-10s leak soft %.2e stiff %.2e velocity %.2e stress %.2e" + % (mode, got["soft"], got["stiff"], got["velocity"], + got["recovered"]), flush=True) + + +if __name__ == "__main__": + command = sys.argv[1:2] or ["sweep"] + {"sweep": sweep, "contrast": contrast, + "params": parameters, "control": control}[command[0]]() diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/solcx_rotated.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/solcx_rotated.py new file mode 100644 index 0000000..69ffe11 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/solcx_rotated.py @@ -0,0 +1,277 @@ +"""SolCx in a rotated frame: an exact solution on walls that are not axis-aligned. + +`solcx.py` runs SolCx on the unit box, where "no flow through the top wall" is a +single velocity component and every treatment reduces to the same constraint. +That is deliberate -- it isolates the viscosity contrast from the geometry -- but +it also means the box says nothing about normals. + +This rotates the whole problem instead. The mesh is turned through THETA, and the +viscosity, the forcing and the exact solution are turned with it, so the physics +is identical and no wall is aligned with a coordinate direction. The exact +solution is still SolCx's, read at the pre-image of each point. + +With x the rotated coordinate and xi = R^T x the material coordinate, + + eta(x) = eta_solcx(xi) + f(x) = R f_solcx(xi) + u(x) = R u_solcx(xi) + sigma(x) = R sigma_solcx(xi) R^T + +and because sigma_nn is frame-invariant, the surface stress on the rotated top +wall is exactly sigma_zz of the unrotated solution at the pre-image. Nothing has +to be un-rotated to make the comparison. + +Every wall carries the treatment under test: a component Dirichlet condition +cannot express u.n = 0 on a wall at 45 degrees, which is the point. The corners +are where two walls meet, and there both constraints apply -- u.n = 0 against two +orthogonal normals, which fixes the node at zero. That is what the exact solution +says a free-slip corner does, so it is a constraint and not a choice. + + python3 solcx_rotated.py # every treatment, at THETA + python3 solcx_rotated.py angles # the rotated constraint against the angle + +Run against underworld3 `bugfix/multiplier-traction` (PR #617). +""" +import sys + +import numpy as np +import sympy + +import underworld3 as uw + +from underworld3.utilities.boundary_flux import _boundary_field_nodes + +RES = 32 +ETA_B = 1.0e6 +PENALTY = 1.0e4 +GAMMA = 10.0 +THETA = np.pi / 4.0 +WALLS = ("Left", "Right", "Bottom", "Top") + + +def rotation(theta): + c, s = np.cos(theta), np.sin(theta) + return np.array([[c, -s], [s, c]]) + + +def build(mode, res=RES, eta_B=ETA_B, penalty=PENALTY, gamma=GAMMA, theta=THETA): + """The rotated problem, and one boundary treatment on every wall.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res, res), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3) + + R = rotation(theta) + mesh.deform(np.asarray(mesh.X.coords) @ R.T) + + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + solver_class = (uw.systems.Stokes_Constrained if mode == "constraint" + else uw.systems.Stokes) + stokes = solver_class(mesh, velocityField=v, pressureField=p) + + # Built on the rotated mesh, but only for its SYMBOLS: every expression below + # is composed with xi = R^T x before it is used, which is what turns the + # solution with the mesh. + exact = uw.analytic.SolCx(mesh, eta_A=1.0, eta_B=eta_B, x_c=0.5, n=1) + x, z = mesh.X + c, s = sympy.cos(theta), sympy.sin(theta) + material = {x: c * x + s * z, z: -s * x + c * z} # xi = R^T x + Rs = sympy.Matrix([[c, -s], [s, c]]) + + def turn(expr): + """An expression of the material coordinate, written in the rotated one.""" + return expr.subs(material, simultaneous=True) + + fn_viscosity = turn(exact.fn_viscosity) + fn_bodyforce = (Rs * turn(exact.fn_bodyforce).T).T + fn_velocity = (Rs * turn(sympy.Matrix(exact.fn_velocity)).reshape(2, 1)).T + # sigma_nn is invariant, so the exact surface stress on the rotated top wall + # is the unrotated sigma_zz at the pre-image -- no rotation of the stress. + fn_sigma_nn = turn(exact.fn_stress[1, 1]) + + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = fn_viscosity + if mode != "constraint": + stokes.saddle_preconditioner = 1.0 / fn_viscosity + stokes.bodyforce = fn_bodyforce + stokes.tolerance = 1.0e-9 + stokes.petsc_use_pressure_nullspace = True + + if mode == "aligned": + # The NEGATIVE CONTROL, and the reason the rotated run means anything: + # the condition the unrotated box would use, imposed unchanged on walls + # that have been turned. It pins the wrong component of the velocity at + # every boundary node, so if the rotation matrix were quietly the + # identity somewhere in the machinery, this is the answer we would get. + stokes.add_dirichlet_bc((0.0, None), "Left") + stokes.add_dirichlet_bc((0.0, None), "Right") + stokes.add_dirichlet_bc((None, 0.0), "Bottom") + stokes.add_dirichlet_bc((None, 0.0), "Top") + return mesh, stokes, v, dict(velocity=fn_velocity, sigma_nn=fn_sigma_nn, + viscosity=fn_viscosity, R=R, theta=theta) + + for wall in WALLS: + if mode == "rotated": + stokes.add_rotated_freeslip_bc(0.0, wall) + elif mode == "constraint": + stokes.add_constraint_bc(0.0, wall) + elif mode == "nitsche": + stokes.add_nitsche_bc(0.0, wall, gamma=gamma, theta=1) + elif mode == "penalty": + n = mesh.boundary_normal(wall) + stokes.add_natural_bc(penalty * n.dot(v.sym) * n, wall) + else: + raise ValueError(mode) + + return mesh, stokes, v, dict(velocity=fn_velocity, sigma_nn=fn_sigma_nn, + viscosity=fn_viscosity, R=R, theta=theta) + + +def converged(stokes): + return stokes.snes.getConvergedReason() > 0 + + +def velocity_error(v, truth): + got = np.squeeze(np.asarray(v.array)) + want = np.asarray(uw.function.evaluate(truth["velocity"], v.coords)).reshape(got.shape) + return float(np.linalg.norm(got - want) / np.linalg.norm(want)) + + +def leak(v, mesh, theta=THETA): + """max |u.n| / |u| on the rotated top wall, against the TRUE wall normal.""" + n = rotation(theta) @ np.array([0.0, 1.0]) + coords = np.asarray(v.coords) + material = coords @ rotation(theta) # xi = R^T x + top = np.abs(material[:, 1] - 1.0) < 1.0e-9 + u = np.squeeze(np.asarray(v.array)) + return float(np.abs(u[top] @ n).max() / np.linalg.norm(u, axis=1).max()) + + +def corner_speed(v, theta=THETA): + """|u| at the four corners, which two free-slip walls fix at zero.""" + coords = np.asarray(v.coords) + material = coords @ rotation(theta) + u = np.squeeze(np.asarray(v.array)) + corner = ((np.abs(material[:, 0] - 0.0) < 1.0e-9) | (np.abs(material[:, 0] - 1.0) < 1.0e-9)) \ + & ((np.abs(material[:, 1] - 0.0) < 1.0e-9) | (np.abs(material[:, 1] - 1.0) < 1.0e-9)) + if not corner.any(): + return float("nan") + return float(np.linalg.norm(u[corner], axis=1).max() / np.linalg.norm(u, axis=1).max()) + + +def trace_top(solver, mesh, theta=THETA): + """Velocity-degree nodes on the rotated top wall.""" + nodes, *_ = _boundary_field_nodes(solver, "Top", 0) + return np.array([node[2] for node in nodes]) + + +def stress_error(coords, values, truth, trim=0.0, theta=THETA): + """Relative l2 error in sigma_nn along the rotated top wall, mean removed.""" + coords = np.asarray(coords) + values = np.asarray(values).reshape(-1) + material = coords @ rotation(theta) + if trim > 0.0: + keep = (material[:, 0] > trim) & (material[:, 0] < 1.0 - trim) + coords, values, material = coords[keep], values[keep], material[keep] + want = np.asarray(uw.function.evaluate(truth["sigma_nn"], coords)).reshape(-1) + order = np.argsort(material[:, 0]) + got = values[order] - values.mean() + want = want[order] - want.mean() + return float(np.linalg.norm(got - want) / np.linalg.norm(want)) + + +def recovered_sigma_nn(mesh, stokes, theta=THETA): + """sigma_nn on the rotated top wall, projected out of the solved fields.""" + n = rotation(theta) @ np.array([0.0, 1.0]) + normal = sympy.Matrix([[n[0], n[1]]]) + field = uw.discretisation.MeshVariable("Snn", mesh, 1, degree=2) + projection = uw.systems.Projection(mesh, field) + projection.uw_function = (normal * stokes.stress * normal.T)[0, 0] + projection.solve() + coords = trace_top(projection, mesh, theta) + tree = uw.kdtree.KDTree(np.ascontiguousarray(field.coords)) + index = np.asarray(tree.query(np.ascontiguousarray(coords), 1)[1]).flatten() + return coords, np.squeeze(np.asarray(field.array))[index] + + +def reaction_sigma_nn(mesh, stokes, v, mode, penalty=PENALTY, theta=THETA): + """The traction the solve already has, for the three treatments that have one.""" + n = rotation(theta) @ np.array([0.0, 1.0]) + if mode == "rotated": + coords, values = stokes.boundary_normal_traction("Top") + return coords, -np.asarray(values) + if mode == "constraint": + coords = trace_top(stokes, mesh, theta) + return coords, -np.asarray( + uw.function.evaluate(stokes.traction("Top"), coords)).reshape(-1) + if mode == "penalty": + coords = trace_top(stokes, mesh, theta) + tree = uw.kdtree.KDTree(np.ascontiguousarray(v.coords)) + index = np.asarray(tree.query(np.ascontiguousarray(coords), 1)[1]).flatten() + u = np.squeeze(np.asarray(v.array))[index] + return coords, -penalty * (u @ n) + return None + + +MODES = ("rotated", "constraint", "penalty", "nitsche", "aligned") + + +def measure(mode, **kwargs): + mesh, stokes, v, truth = build(mode, **kwargs) + stokes.solve() + if not converged(stokes): + return None + theta = kwargs.get("theta", THETA) + out = { + "leak": leak(v, mesh, theta), + "corner": corner_speed(v, theta), + "velocity": velocity_error(v, truth), + } + coords, values = recovered_sigma_nn(mesh, stokes, theta) + out["recovered"] = stress_error(coords, values, truth, theta=theta) + read = reaction_sigma_nn(mesh, stokes, v, mode, + penalty=kwargs.get("penalty", PENALTY), theta=theta) + if read is not None: + coords, values = read + out["reaction"] = stress_error(coords, values, truth, theta=theta) + # Two elements off each end. A corner node belongs to two walls, and the + # measure-weighted normal there is the bisector of two orthogonal ones, + # so u.n at a corner is not the quantity either wall asked for. + out["reaction_trimmed"] = stress_error( + coords, values, truth, trim=2.0 / kwargs.get("res", RES), theta=theta) + return out + + +def table(modes=MODES, **kwargs): + print("SolCx rotated through %.1f degrees, %d x %d." + % (np.degrees(kwargs.get("theta", THETA)), RES, RES)) + print() + print("| treatment | leak | corner |u| | velocity | sigma_nn, recovered | sigma_nn, from the solve |") + print("|---|---|---|---|---|---|") + for mode in modes: + got = measure(mode, **kwargs) + if got is None: + print("| %s | did not converge | | | | |" % mode, flush=True) + continue + print("| %s | %.1e | %.1e | %.2e | %.3f | %s |" + % (mode, got["leak"], got["corner"], got["velocity"], got["recovered"], + "%.3f" % got["reaction"] if "reaction" in got else "—"), flush=True) + + +def angles(degrees=(0.0, 15.0, 30.0, 45.0), mode="rotated"): + print("%s, against the angle of the frame" % mode) + print() + print("| degrees | leak | velocity | sigma_nn, from the solve |") + print("|---|---|---|---|") + for d in degrees: + got = measure(mode, theta=np.radians(d)) + if got is None: + print("| %.0f | did not converge | | |" % d, flush=True) + continue + print("| %.0f | %.1e | %.2e | %s |" + % (d, got["leak"], got["velocity"], + "%.3f" % got["reaction"] if "reaction" in got else "—"), flush=True) + + +if __name__ == "__main__": + command = sys.argv[1] if len(sys.argv) > 1 else "table" + {"table": table, "angles": angles}[command]() diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/stress.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/stress.py new file mode 100644 index 0000000..8476cf6 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/stress.py @@ -0,0 +1,448 @@ +"""Does the surface stress come out right? + +`leak.py` measures the CONSTRAINT -- how much flow each boundary treatment lets +through. This measures the thing the constraint is wanted for: the wall-normal +stress on the boundary, which is dynamic topography once divided by the +buoyancy scale. Against an exact answer, so it is a measurement rather than a +comparison of two methods that might both be wrong. + +The oracle +---------- +Kramer, Davies & Wilson (2021) give exact Stokes solutions in a cylindrical +annulus, and their `assess` package publishes the RADIAL STRESS directly -- +`radial_stress = tau_rr - p` -- which is the quantity here. Underworld wraps +the package as `uw.analytic.CylindricalStokes`. + +The case is the smooth one: density (r/r_o)^k cos(n.theta) driving the flow, +free slip on both radii, viscosity 1. On the outer boundary the exact +sigma_rr is a pure cos(n.theta) (checked, residual 1e-16), so the measurement +is the amplitude of that harmonic and the metric is its relative error. + +The inner boundary carries the EXACT velocity as a Dirichlet condition rather +than a free-slip treatment of its own. The exact solution satisfies both, so +the problem is unchanged; what it buys is that the treatment under test is the +only free-slip condition in the model, and that the rigid-rotation nullspace of +two free-slip circles is gone. + +What is compared +---------------- +Two routes to the surface stress, and the difference between them is the +note's argument: + + * RECOVERED -- project r.sigma.r out of the solved velocity and pressure. + Available for every treatment, and the only route the weak ones have. + * REACTION -- the constraint reaction of the rotated method + (`boundary_normal_traction`), and the multiplier field of the constraint + method. Not recovered from the solution: an unknown the solve returned. + +Both are compared against the same exact amplitude, using the TRUE radial +direction, so no treatment is being scored against its own normal. + +Sign: `boundary_normal_traction` and the multiplier both return the constraint +reaction, which is opposite in sign to sigma_rr as `assess` publishes it (it is +the traction that holds the boundary, not the traction the fluid exerts). The +dynamic-topography formula h = -(sigma_nn - mean)/(rho.g) carries the sign +back. Amplitudes are compared unsigned and the sign is printed so the +convention stays visible. + + python3 stress.py sweep # the refinement table + python3 stress.py params # penalty coefficient and Nitsche gamma + python3 stress.py locking # the facet normal does not converge + python3 stress.py control # the metric fires when the BC is removed + +Run against underworld3 `bugfix/multiplier-traction` (PR #617); the +constrained solver's `traction()` is the fix this note prompted. +""" +import sys + +import numpy as np +import sympy + +import underworld3 as uw + +# The boundary trace of a field, selected by the mesh boundary LABEL. An earlier +# version of this comparison selected it by a radius band that narrowed with the +# mesh, which quietly admitted a different node set at each resolution. This is +# the same selector the solver's own reaction recovery uses. +from underworld3.utilities.boundary_flux import _boundary_field_nodes + +N = 2 # azimuthal wavenumber of the density anomaly +K = 3 # its radial power +R_I, R_O = 0.5, 1.0 +CELL = 0.075 +PENALTY = 1.0e4 +GAMMA = 10.0 + + +def build(mode, cell=CELL, penalty=PENALTY, gamma=GAMMA): + """The annulus, the forcing, and one boundary treatment on the outer arc.""" + mesh = uw.meshing.Annulus(radiusInner=R_I, radiusOuter=R_O, cellSize=cell) + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + theta = sympy.atan2(y, x) + + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + solver_class = (uw.systems.Stokes_Constrained if mode == "constraint" + else uw.systems.Stokes) + stokes = solver_class(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + + # The forcing `assess` solves for: rho = (r/r_o)^k cos(n.theta), gravity + # inward. This is the convention the Kramer benchmark in underworld3's + # docs/examples uses, and the velocity error below is what checks it. + rho = ((r / R_O) ** K) * sympy.cos(N * theta) + stokes.bodyforce = -rho * sympy.Matrix([[x / r, y / r]]) + + exact = uw.analytic.CylindricalStokes( + mesh, n=N, k=K, r_inner=R_I, r_outer=R_O, + density="smooth", boundary="free") + + v_exact = uw.discretisation.MeshVariable("Uex", mesh, mesh.dim, degree=2) + v_exact.data[:] = exact.evaluate("velocity", v_exact.coords) + stokes.add_dirichlet_bc(v_exact.sym, "Lower") + + if mode == "constraint": + stokes.add_constraint_bc(0.0, "Upper") + elif mode == "rotated": + stokes.add_rotated_freeslip_bc(0.0, "Upper") + elif mode == "nitsche": + # Nitsche's default normal is mesh.boundary_normal(boundary) -- the + # measure-weighted node normal, not the per-facet one. + stokes.add_nitsche_bc(0.0, "Upper", gamma=gamma, theta=1) + elif mode == "penalty": + # The documented direct penalty, with the quadrature-point FACET normal. + stokes.add_natural_bc(penalty * mesh.Gamma.dot(v.sym) * mesh.Gamma, "Upper") + elif mode == "penalty_node": + # The same penalty against the measure-weighted NODE normal. The only + # difference between this and the line above is which normal, and it is + # the difference between a method that converges and one that locks. + node_n = mesh.boundary_normal("Upper") + stokes.add_natural_bc(penalty * node_n.dot(v.sym) * node_n, "Upper") + elif mode == "free": + pass # the control: no condition on the outer arc + else: + raise ValueError(mode) + + return mesh, stokes, v, exact + + +def converged(stokes): + """A diverged solve leaves numbers in the array that look like measurements.""" + return stokes.snes.getConvergedReason() > 0 + + +def exact_amplitude(exact): + """Amplitude of cos(n.theta) in the exact sigma_rr on the outer boundary. + + Sampled densely on the true circle and least-squares fitted. The residual + is returned too: it is 1e-16, which is what says the exact surface stress + really is a single harmonic and the amplitude is the whole of it. + """ + phi = np.linspace(0.0, 2.0 * np.pi, 720, endpoint=False) + points = np.c_[R_O * np.cos(phi), R_O * np.sin(phi)] + sigma = np.array([exact._above.radial_stress_cartesian(pt) for pt in points]) + (mean, c, s), residual = _fit(phi, sigma) + return c, residual + + +def _fit(angles, values): + """Least squares fit of a constant plus the degree-N harmonic. + + A least-squares fit rather than a quadrature: boundary nodes are not evenly + spaced in theta, and on a P2 trace vertices and edge midpoints alternate + with unequal gaps. Returns ((mean, cos, sin), max residual). + """ + basis = np.c_[np.ones_like(angles), np.cos(N * angles), np.sin(N * angles)] + coefficients, *_ = np.linalg.lstsq(basis, values, rcond=None) + residual = np.abs(values - basis @ coefficients).max() + return coefficients, residual + + +def trace(solver, field_id, var, boundary="Upper"): + """(coords, values) for `var` at the nodes it carries on `boundary`.""" + nodes, *_ = _boundary_field_nodes(solver, boundary, field_id) + coords = np.array([node[2] for node in nodes]) + tree = uw.kdtree.KDTree(np.ascontiguousarray(var.coords)) + index = np.asarray(tree.query(np.ascontiguousarray(coords), 1)[1]).flatten() + # squeeze: a MeshVariable's .array carries a middle axis that broadcasts + # silently against an (N,) index -- see leak.py. + return coords, np.squeeze(np.asarray(var.array))[index] + + +def recovered_traction(mesh, stokes): + """sigma_rr on the outer boundary, projected out of the solved fields. + + The recovery every treatment can do, and the only one the weak forms have. + Against the TRUE radial direction, which is also the direction the oracle + publishes, so the comparison does not depend on the solver's normal. + """ + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + radial = sympy.Matrix([[x / r, y / r]]) + sigma_rr = (radial * stokes.stress * radial.T)[0, 0] + + field = uw.discretisation.MeshVariable("Srr", mesh, 1, degree=2) + projection = uw.systems.Projection(mesh, field) + projection.uw_function = sigma_rr + projection.solve() + return trace(projection, 0, field) + + +# The augmented-Lagrangian parameter this problem gets by default: the base +# (1e4) times the local viscosity, which is 1 here. +AUGMENTATION = 1.0e4 + + +def multiplier_traction(stokes, v, boundary="Upper"): + """The WHOLE traction a multiplier constraint holds the boundary with. + + The momentum row carries `h + r(u.n - g)`, so `h` alone is short by `r` times + the discrete constraint residual. `stokes.traction(boundary)` is that sum as + an expression; this reads it off the node arrays instead, because evaluating + an expression at points sitting exactly on a convex curved boundary + extrapolates from the containing cell (underworld3#605) and that error would + land on top of the measurement. + + `u.n` is taken against the true radial direction rather than the constraint's + node normal. The two differ at O(h^2) and they multiply a term that is itself + a correction. + """ + nodes, *_ = _boundary_field_nodes(stokes, boundary, 2) + coords = np.array([node[2] for node in nodes]) + # The multiplier is carried at the VELOCITY degree, so h and u share a node + # set and one lookup serves both. Building it is the only part of this that + # is not arithmetic, and it depends on the mesh rather than on the solution: + # a time-stepping consumer builds it once and reuses it every step. + index = _trace_index(v, coords) + h = np.squeeze(np.asarray(stokes.multiplier(boundary).array))[index] + u = np.squeeze(np.asarray(v.array))[index] + normal = coords / np.linalg.norm(coords, axis=1)[:, None] + return coords, h + AUGMENTATION * (u * normal).sum(axis=1) + + +_TRACE_INDEX = {} + + +def _trace_index(var, coords): + """Row of `var`'s arrays for each trace coordinate, cached per (mesh, trace).""" + key = (id(var), coords.shape[0], float(coords.sum())) + if key not in _TRACE_INDEX: + tree = uw.kdtree.KDTree(np.ascontiguousarray(var.coords)) + _TRACE_INDEX[key] = np.asarray( + tree.query(np.ascontiguousarray(coords), 1)[1]).flatten() + return _TRACE_INDEX[key] + + +def penalty_traction(stokes, v, kappa=PENALTY, boundary="Upper"): + """The traction a direct penalty holds the boundary with. + + The term the weak form carries IS that traction, so sigma_nn is + -kappa (u.n) on the boundary and there is nothing to recover. Read at the + boundary nodes against the TRUE radial direction, as `recovered_traction` + is, so the comparison does not depend on the solver's normal. + + It is only as good as the constraint it holds: the leak and this number are + the same quantity scaled by kappa. + """ + nodes, *_ = _boundary_field_nodes(stokes, boundary, 0) + coords = np.array([node[2] for node in nodes]) + index = _trace_index(v, coords) + u = np.squeeze(np.asarray(v.array))[index] + normal = coords / np.linalg.norm(coords, axis=1)[:, None] + return coords, -kappa * (u * normal).sum(axis=1) + + +def reaction_traction(stokes, mode, v=None): + """The traction each treatment returns without recovering one. + + Three of the four do. Nitsche is the exception: its term is written in + sigma(u), so reading it back means differentiating the velocity, which is + what `recovered_traction` does. + """ + if mode == "rotated": + return stokes.boundary_normal_traction("Upper") + if mode == "constraint": + return multiplier_traction(stokes, v) + if mode in ("penalty", "penalty_node"): + return penalty_traction(stokes, v) + return None + + +def amplitude_error(coords, values, reference): + """Relative error in the harmonic amplitude, unsigned, and the sign.""" + angles = np.arctan2(coords[:, 1], coords[:, 0]) + (_mean, c, _s), _residual = _fit(angles, values) + return abs(abs(c) - abs(reference)) / abs(reference), np.sign(c) + + +def split(coords): + """Vertices from edge midpoints. A vertex of the annulus mesh sits exactly + on the circle; a P2 edge midpoint sits on the chord, inside it by the + sagitta. They are not interchangeable: vertex values of sigma_nn carry the + O(h) facet-geometry error and midpoint values are superconvergent + (underworld3#414).""" + on_circle = np.abs(np.linalg.norm(coords, axis=1) - R_O) < 1.0e-9 + return {"vertex": on_circle, "midpoint": ~on_circle} + + +def leak(v): + """max |u.n| on the outer boundary against the TRUE radial direction, + normalised by the flow speed -- the measurement leak.py tabulates.""" + coords = v.coords + radius = np.linalg.norm(coords, axis=1) + outer = np.abs(radius - R_O) < 1.0e-6 + u = np.squeeze(np.asarray(v.array)) + normal = coords[outer] / radius[outer][:, None] + return np.abs((u[outer] * normal).sum(axis=1)).max() / np.linalg.norm(u, axis=1).max() + + +def measure(mode, cell=CELL, **kwargs): + """One solve, and everything read off it.""" + mesh, stokes, v, exact = build(mode, cell=cell, **kwargs) + stokes.solve() + if not converged(stokes): + return None + reference, _residual = exact_amplitude(exact) + out = { + "leak": leak(v), + "velocity": exact.error("velocity", v), + "exact": reference, + } + coords, values = recovered_traction(mesh, stokes) + out["recovered"], out["sign"] = amplitude_error(coords, values, reference) + read = reaction_traction(stokes, mode, v=v) + if read is not None: + coords, values = read + out["reaction"], out["reaction_sign"] = amplitude_error(coords, values, reference) + for name, mask in split(coords).items(): + out["reaction_" + name] = amplitude_error( + coords[mask], values[mask], reference)[0] + return out + + +MODES = ("penalty_node", "nitsche", "constraint", "rotated") + + +def both(cells=(0.15, 0.10, 0.075, 0.05), modes=MODES): + """The constraint and the stress, from the SAME solves. + + The two questions -- does the boundary hold, and is the answer right -- are + one experiment. Reporting them from separate runs invited the reader to + compare numbers that came from different problems. + """ + exact, _residual = exact_amplitude(build("free", cell=0.2)[3]) + print("exact sigma_rr on r = %.2f: %.10f cos(%d.theta)" % (R_O, exact, N)) + print() + print("leak / stress error, from one solve each") + print() + print("| cell size | " + " | ".join(modes) + " |") + print("|---" * (len(modes) + 1) + "|") + for cell in cells: + row = [] + for mode in modes: + got = measure(mode, cell=cell) + row.append("diverged" if got is None + else "%.1e / %.1e" % (got["leak"], got["recovered"])) + print("| %.3f | %s |" % (cell, " | ".join(row)), flush=True) + print() + print("the reaction routes, same solves") + print() + print("| cell size | rotated reaction | multiplier traction | penalty term |") + print("|---|---|---|---|") + for cell in cells: + entries = [] + for mode in ("rotated", "constraint", "penalty_node"): + got = measure(mode, cell=cell) + entries.append("-" if not got or "reaction" not in got + else "%.1e" % got["reaction"]) + print("| %.3f | %s |" % (cell, " | ".join(entries)), flush=True) + + +def sweep(cells=(0.15, 0.10, 0.075, 0.05), modes=MODES): + """Does the surface stress converge to the exact one, and how fast?""" + exact, residual = exact_amplitude(build("free", cell=0.2)[3]) + print("exact sigma_rr on r = %.2f: %.10f cos(%d.theta), residual %.1e" + % (R_O, exact, N, residual)) + print() + print("relative error in the surface stress amplitude") + print() + print("| cell size | " + " | ".join(modes) + " | rotated (reaction) | multiplier |") + print("|---" * (len(modes) + 3) + "|") + for cell in cells: + row, extra = [], {} + for mode in modes: + got = measure(mode, cell=cell) + row.append("diverged" if got is None else "%.2e" % got["recovered"]) + if got and "reaction" in got: + extra[mode] = "%.2e" % got["reaction"] + print("| %.3f | %s | %s | %s |" + % (cell, " | ".join(row), + extra.get("rotated", "-"), extra.get("constraint", "-")), flush=True) + + +def parameters(cell=CELL): + """The two weak methods against their own parameter, with the leak beside + the stress -- which is where they part company.""" + for mode, values, label in ( + ("penalty", (1e3, 1e4, 1e5, 1e6, 1e8), "penalty (facet normal)"), + ("penalty_node", (1e3, 1e4, 1e5, 1e6), "penalty (node normal)"), + ("nitsche", (1.0, 10.0, 100.0, 1000.0, 1e4), "Nitsche gamma")): + print("\n%s, cell %.3f" % (label, cell)) + print("\n| %s | leak | velocity error | stress error |" + % ("gamma" if mode == "nitsche" else "coefficient")) + print("|---|---|---|---|") + for value in values: + kwargs = {"gamma": value} if mode == "nitsche" else {"penalty": value} + got = measure(mode, cell=cell, **kwargs) + if got is None: + print("| %g | diverged | | |" % value) + continue + print("| %g | %.2e | %.2e | %.2e |" + % (value, got["leak"], got["velocity"], got["recovered"]), flush=True) + + +def locking(cells=(0.15, 0.10, 0.075, 0.05, 0.035), penalty=1.0e6): + """The facet normal, pushed hard, does not converge to free slip. + + Imposing u.n = 0 facet by facet on a curved boundary constrains a corner + node in two directions at once, and the discrete limit is not the smooth + problem. Refine it and it stays wrong -- while the leak, which is the + metric that would normally be trusted, reads 1e-5. + """ + print("direct penalty at %.0e, against the FACET normal" % penalty) + print("\n| cell size | leak | velocity error | stress error |") + print("|---|---|---|---|") + for cell in cells: + got = measure("penalty", cell=cell, penalty=penalty) + if got is None: + print("| %.3f | diverged | | |" % cell) + continue + print("| %.3f | %.2e | %.2e | %.2e |" + % (cell, got["leak"], got["velocity"], got["recovered"]), flush=True) + print("\nthe same coefficient against the measure-weighted NODE normal") + print("\n| cell size | leak | velocity error | stress error |") + print("|---|---|---|---|") + for cell in cells: + got = measure("penalty_node", cell=cell, penalty=penalty) + if got is None: + print("| %.3f | diverged | | |" % cell) + continue + print("| %.3f | %.2e | %.2e | %.2e |" + % (cell, got["leak"], got["velocity"], got["recovered"]), flush=True) + + +def control(cell=CELL): + """Take the boundary condition away. A metric that cannot see that is not + measuring anything.""" + for mode in ("free", "rotated"): + got = measure(mode, cell=cell) + print("%-9s leak %.2e velocity error %.2e stress error %.2e" + % (mode, got["leak"], got["velocity"], got["recovered"]), flush=True) + + +if __name__ == "__main__": + command = sys.argv[1:2] or ["sweep"] + {"sweep": sweep, "both": both, "params": parameters, + "locking": locking, "control": control}[command[0]]() diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/timing.py b/articles/boundary-conditions-on-non-planar-boundaries/examples/timing.py new file mode 100644 index 0000000..e068184 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/timing.py @@ -0,0 +1,101 @@ +"""What each treatment costs. + +The note compares four ways of imposing free slip on accuracy. This measures the +other half of the choice: what each one costs to solve, and what its surface +traction costs to recover once solved. + +Two numbers per treatment, because they are paid at different times: + + * the SOLVE. The rotated constraint changes the operator; the multiplier adds a + field and enlarges the saddle point; the weak forms add a boundary term to a + system that is otherwise the plain Stokes one. + * the RECOVERY of the surface traction. The weak forms have to project + `n.sigma.n` out of the solution, which is a second (scalar, symmetric) solve. + The rotated constraint and the multiplier read theirs off the state the solve + already returned, which is arithmetic on the boundary trace and no solve at + all. + +Method: each configuration is built and solved once UNTIMED (JIT compilation, +PETSc setup and the first-touch allocations are not what is being measured), then +timed `repeats` times. Runs are sequential by construction -- concurrent PETSc +solves contend for memory bandwidth and inflate each other by tens of per cent. + + python3 timing.py # the table + python3 timing.py 0.05 # at one cell size + +Run against underworld3 `bugfix/multiplier-traction`. +""" +import sys +import time + +import numpy as np +import sympy + +import underworld3 as uw + +import stress as S + +REPEATS = 3 +# Big enough that a solve is seconds rather than hundredths: at 10k nodes the +# four treatments were separated by less than the run-to-run spread. +CELLS = (0.02, 0.0125) +MODES = ("penalty_node", "nitsche", "constraint", "rotated") + + +def _time(call, repeats=REPEATS): + """Median of `repeats` timings, and the spread, in seconds. + + One UNTIMED call first. The first recovery in a process compiles its + projection: timed cold it read 1.145 s where the same call reads 0.038 s + warm, which is a measurement of the JIT and not of the method. + """ + call() + got = [] + for _ in range(repeats): + start = time.perf_counter() + call() + got.append(time.perf_counter() - start) + return float(np.median(got)), float(np.max(got) - np.min(got)) + + +def solve_cost(mode, cell): + """Seconds to solve, and the nonlinear/linear iteration counts.""" + mesh, stokes, v, exact = S.build(mode, cell=cell) + stokes.solve() # untimed: JIT, setup, first touch + assert S.converged(stokes) + iterations = stokes.snes.getIterationNumber() + median, spread = _time(stokes.solve) + return mesh, stokes, v, exact, median, iterations + + +def recovery_cost(mesh, stokes, mode): + """Seconds to get sigma_nn on the boundary, by the route that treatment has.""" + if mode == "rotated": + return _time(lambda: stokes.boundary_normal_traction("Upper")) + if mode == "constraint": + # The whole traction, h + r(u.n - g). Reading it is an expression build + # plus a boundary-trace evaluation -- no solve. + return _time(lambda: S.multiplier_traction(stokes, stokes.u)) + # The weak forms: project n.sigma.n and read its trace. A scalar solve. + return _time(lambda: S.recovered_traction(mesh, stokes)) + + +def table(cells=CELLS, modes=MODES): + print("seconds, median of %d timed repeats after one untimed warm-up" % REPEATS) + print() + print("| cell size | velocity nodes | " + " | ".join( + "%s: solve / recover" % m for m in modes) + " |") + print("|---" * (len(modes) + 2) + "|") + for cell in cells: + row, nodes = [], None + for mode in modes: + mesh, stokes, v, exact, solve, its = solve_cost(mode, cell) + nodes = len(v.coords) if nodes is None else nodes + recover, _s = recovery_cost(mesh, stokes, mode) + row.append("%.2f / %.3f" % (solve, recover)) + print("| %.3f | %d | %s |" % (cell, nodes, " | ".join(row)), flush=True) + + +if __name__ == "__main__": + cells = (float(sys.argv[1]),) if sys.argv[1:] else CELLS + table(cells=cells) diff --git a/articles/boundary-conditions-on-non-planar-boundaries/examples/topography-data.json b/articles/boundary-conditions-on-non-planar-boundaries/examples/topography-data.json new file mode 100644 index 0000000..299a819 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/examples/topography-data.json @@ -0,0 +1 @@ +{"contrasts": [100.0, 1000000.0], "truth": {"100.0": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2794438305496565, 0.2792336766663841, 0.27860509329600136, 0.27756371664641966, 0.2761189449799955, 0.27428394554731295, 0.2720756644146153, 0.2695148393052133, 0.26662601560986043, 0.2634375657559502, 0.25998171216052346, 0.25629455402751367, 0.25241609828546413, 0.24839029499817417, 0.2442650776174339, 0.24009240848423946, 0.23592833002271932, 0.23183302210949724, 0.22787086614045374, 0.22411051635688298, 0.22062497903397357, 0.21749170017642513, 0.21479266240896497, 0.2126144917936135, 0.21104857535088614, 0.21019119010879417, 0.21014364455164997, 0.21101243339038472, 0.21290940662749624, 0.2159519539429761, 0.22026320548276043, 0.22597225018855668, 0.03791009878192678, 0.013354106362604068, -0.010480844032193426, -0.03360399506264387, -0.05602180902156315, -0.07773822478757993, -0.09875490409895744, -0.11907146745092187, -0.138685719891038, -0.15759386696074398, -0.17579072100665732, -0.1932698980625954, -0.21002400548230854, -0.22604482048363547, -0.2413234597470689, -0.2558505401954737, -0.2696163310668727, -0.282610897378695, -0.29482423486964227, -0.30624639649422947, -0.3168676105350941, -0.32667839038921354, -0.3356696360761876, -0.34383272750965166, -0.3511596095665777, -0.35764286898384134, -0.36327580310611063, -0.3680524805059488, -0.3719677934919559, -0.3750175025187017, -0.377198272508923, -0.37850770109619103, -0.37894433879410594]], "1000000.0": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.27879056265826546, 0.27858570337229305, 0.2779730314553374, 0.2769582662024702, 0.27555094469691077, 0.27376442925837235, 0.2716159179932234, 0.26912645856914497, 0.26632096537230676, 0.2632282402406444, 0.2598809970026598, 0.2563158900873196, 0.2525735475071686, 0.24869860855375314, 0.24473976658191288, 0.24074981729753853, 0.2367857130020395, 0.23290862328611167, 0.22918400270550734, 0.22568166601244988, 0.22247587155821263, 0.21964541352524677, 0.21727372369121403, 0.21544898347243532, 0.21426424704071523, 0.21381757635534415, 0.21421218900143502, 0.2155566197767385, 0.21796489702181682, 0.22155673474309964, 0.22645774163500287, 0.23279964816615384, 0.03680939726302895, 0.012188222315832912, -0.011705092371091049, -0.034880357082689994, -0.05734456400948344, -0.07910214645220835, -0.10015522720934285, -0.12050385645768474, -0.14014623940740678, -0.15907895398616861, -0.17729715878197405, -0.19479479145141543, -0.21156475777865164, -0.22759911155083346, -0.24288922539765032, -0.2574259527261132, -0.2711997808665712, -0.28420097553215773, -0.29641971668134703, -0.3078462258619563, -0.3184708851047058, -0.32828434742526713, -0.3372776389855125, -0.3454422529573666, -0.3527702351261416, -0.3592542612646509, -0.3648877063038704, -0.3696647053225625, -0.3735802063730181, -0.3766300151578111, -0.37881083156900464, -0.38012027809879434, -0.38055692012825537]]}, "data": {"100.0": {"dirichlet": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2782418279108243, 0.27795090847582005, 0.2774011632642354, 0.27628016945796835, 0.2749096468589217, 0.27299898546768503, 0.2708573623719712, 0.2682280887166753, 0.2653954547431387, 0.26214836101975664, 0.2587356150452108, 0.25500264306868486, 0.2511511347800553, 0.24709565503897396, 0.24297848312505035, 0.23879529911563188, 0.2346181571206335, 0.2305334155643952, 0.2265348800525844, 0.22280950764873517, 0.21925945709397013, 0.2161923282455665, 0.21339339993378814, 0.2113188002657265, 0.20957556074635947, 0.20893919527764163, 0.20845420730063483, 0.20991187692890423, 0.20988190560080627, 0.21602237921120912, 0.2104292535820529, 0.23167127253185973, 0.18066180199740106, 0.03642722690571339, -0.034080033793114924, -0.049773364839333034, -0.0525673479605493, -0.08007006850432118, -0.10507379400592183, -0.12408430104956213, -0.14240250492325265, -0.16131478457439166, -0.17984906936951062, -0.19689712203809437, -0.21362162042299337, -0.22932787700518756, -0.24471782303496337, -0.2589115976681528, -0.272785584594719, -0.2854656331625424, -0.29781844254478496, -0.3089326508225187, -0.31971372925302616, -0.32922487059443944, -0.33839662665456205, -0.3462664870648695, -0.35379208667626927, -0.359988876889968, -0.36583653193352944, -0.37033425753757104, -0.3744789565980029, -0.3772569301478814, -0.3796799521989641, -0.3807266465265126, -0.3814158520790764]], "penalty": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2734434607089738, 0.27315597627556965, 0.2726038442577027, 0.27148703738510366, 0.2701149994470591, 0.2682091816534774, 0.26606702424066286, 0.2634426222743534, 0.2606103655050831, 0.2573681099208353, 0.25395625931606136, 0.25022746546257624, 0.24637736004361915, 0.24232433432861228, 0.23820862042263533, 0.23402512285081123, 0.229848256457529, 0.22576037392496515, 0.22175932772350343, 0.21802579094864466, 0.21446893310131174, 0.2113865408140821, 0.2085699684063549, 0.20647553424539836, 0.204670110879541, 0.20404778701302284, 0.20327761403502864, 0.2050285359612285, 0.20336878393219462, 0.21181301227580251, 0.19654386554297246, 0.23198695784655315, 0.1198370861771435, 0.01778744460117032, -0.03363799805624428, -0.057102447999933796, -0.06974987031571535, -0.08961535955844206, -0.1103728073242995, -0.12818860374591623, -0.1450130704149113, -0.16188256603773057, -0.178385536907645, -0.1940369578885589, -0.2091185764480828, -0.22356862562770735, -0.2374692891837668, -0.25063006763403217, -0.26319913141129864, -0.275016235669273, -0.28621978600999415, -0.29663438588479835, -0.3064108477333798, -0.3153723223859676, -0.3236752228497963, -0.3311383039371429, -0.33792488038100404, -0.3438513139234628, -0.34908661958144444, -0.35344545190816357, -0.3571019270694604, -0.3598701030211166, -0.36192832677838466, -0.3630910620089356, -0.3635400102829753]], "constraint": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2785268916874456, 0.2782362072813296, 0.2776869780008344, 0.2765668146312993, 0.27519730099262374, 0.2732881720138809, 0.2711481412282027, 0.2685207428606257, 0.26569026483136077, 0.26244567675318464, 0.2590353886996808, 0.25530537671028253, 0.2514568541420485, 0.247404308395994, 0.24329057485016536, 0.23911003822493224, 0.23493615229162712, 0.23085474121392618, 0.2268586626864397, 0.22313658339294942, 0.2195889229988021, 0.21652394204303937, 0.2137274697673687, 0.21165448540580842, 0.20991417750006322, 0.20927536538669214, 0.2087959418075856, 0.21024626081018533, 0.2102304777196122, 0.21634526556566497, 0.21083195418937337, 0.2319343235246382, 0.18140933737821205, 0.03710021137842936, -0.033767913006413165, -0.049437718903158484, -0.05225749017332368, -0.07984288850937456, -0.10493596233475519, -0.12400773437050396, -0.1423879035076937, -0.16136552321687428, -0.17996047779102095, -0.19706771481013513, -0.21384716452295602, -0.22960824730856488, -0.24504901704304766, -0.25929194388038135, -0.2732125512354152, -0.28593517421703807, -0.2983306462766995, -0.3094817777056305, -0.3202999145530307, -0.3298423374087925, -0.339045655091293, -0.34694179522966867, -0.35449298560787923, -0.36071072452809877, -0.3665783491096681, -0.3710910656541535, -0.37524951088462055, -0.3780373337494884, -0.38046797907811386, -0.3815183855556869, -0.38221012110082725]], "constraint+r": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2780250586999775, 0.27773438682315654, 0.2771851800188121, 0.2760650636888428, 0.27469560674976967, 0.2727865576611125, 0.27064661571218995, 0.2680193271236989, 0.26518896671483305, 0.26194451412250713, 0.25853436786333756, 0.25480451157574835, 0.2509561490686256, 0.24690377245854922, 0.242790210027737, 0.23860984786903577, 0.23443613555923398, 0.2303548948167916, 0.2263589825052833, 0.22263705857531693, 0.2190895457148127, 0.21602469274269523, 0.21322833670477423, 0.21115543899248954, 0.20941520331826666, 0.2087764181058338, 0.20829701423008262, 0.2097472730488394, 0.20973149037932884, 0.21584602355777155, 0.21033294181381051, 0.2314344319069696, 0.18016536669378158, 0.036457528857207225, -0.034115311497949524, -0.04971982646962084, -0.052527848749710386, -0.079998307970861, -0.10498682736757695, -0.12397913371335657, -0.14228271879829957, -0.16118126515147083, -0.17969874071626651, -0.19673469765546833, -0.2134442329529046, -0.2291396446364526, -0.2445160777858872, -0.2586996591664178, -0.2725622639532868, -0.2852318760667844, -0.297575700296849, -0.3086803687308042, -0.31945342998369725, -0.3289560927880209, -0.3381210632905599, -0.345984302885642, -0.3535040299590128, -0.35969586164334605, -0.36553903779089797, -0.3700329514100264, -0.374174069745771, -0.3769502766834633, -0.3793707943045827, -0.3804168241192964, -0.38110567741693985]], "rotated": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2780836286127151, 0.27771539632297365, 0.2772428305102198, 0.27604699263763605, 0.2747504930631134, 0.2727703314117768, 0.27069687518109836, 0.2680058830064543, 0.26523270629810397, 0.26193480783388373, 0.25856965089573136, 0.2547995240068448, 0.25098098127662455, 0.2469045168334859, 0.2428025264997572, 0.23861737559553833, 0.23443380033075087, 0.23037029000125492, 0.2263398706608663, 0.2226613613295593, 0.21905211615813608, 0.21605840296978365, 0.2131750807386179, 0.21119556380202675, 0.20937285640589257, 0.20879921200153523, 0.20843453486593808, 0.2096075178676945, 0.21104770070592893, 0.21469010147701914, 0.21866498925289288, 0.22427876033188035, 0.18039142951778595, 0.043777232514659314, -0.04270028156733962, -0.048433405359016087, -0.05405013382132685, -0.07974361695642906, -0.10529693172802608, -0.12390339680101578, -0.14238283266565382, -0.16113615447979718, -0.17976235748776243, -0.19669488095931784, -0.21350152330931566, -0.2291004440016308, -0.24457269294554113, -0.25866014444221747, -0.27261945765994294, -0.2851917842944213, -0.2976338672581672, -0.3086396132955738, -0.31951271238821494, -0.32891466057816465, -0.3381814738183109, -0.34594223273495883, -0.3535654865484508, -0.35965323304082075, -0.36560138395742436, -0.36998987602458144, -0.37423709186296483, -0.3769068901616813, -0.37943423816495914, -0.38037327805168714, -0.38116926455346495]]}, "1000000.0": {"dirichlet": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.27752274369103486, 0.27724027321740596, 0.27670345943583224, 0.2756124375114777, 0.2742765303892372, 0.27241795635766675, 0.27033338912412286, 0.26777936416287207, 0.265027448559996, 0.2618802838820674, 0.2585736139450435, 0.2549672768107007, 0.2512494206790332, 0.24734978614149308, 0.24339659542086076, 0.23940161875775093, 0.235422103173992, 0.23156153892293302, 0.22779814123708836, 0.22433768100133936, 0.22106524838005803, 0.21830751667913043, 0.21583244039896896, 0.2141232426341522, 0.21275992257475634, 0.2125349352506103, 0.21246462671910016, 0.21447426697403696, 0.21492221328688862, 0.2215784218665878, 0.21596015497183804, 0.2391933134896972, 0.18689817251492444, 0.03654286149267027, -0.03677064670728274, -0.052000462632792066, -0.05376432467841929, -0.08161405497668446, -0.10686726794956597, -0.12583712299055602, -0.14409872060269907, -0.16303945310787604, -0.18160045058577629, -0.19864824208265927, -0.21537333894205118, -0.23108210021829773, -0.24647600643344675, -0.26066837976758317, -0.2745426518951873, -0.28722023127896706, -0.29957239363151517, -0.31068308804462297, -0.32146262111321877, -0.3309698476907814, -0.3401397662558658, -0.34800567050601244, -0.3555294883396407, -0.36172261850795273, -0.36756884655580196, -0.3720634566732739, -0.3762073113521223, -0.37898288075599185, -0.3814058016747456, -0.38245091115743074, -0.3831408425465471]], "penalty": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.036850510598519674, 0.03657177549854867, 0.03603314327254125, 0.03494749202791725, 0.033611528911999526, 0.03175991680018898, 0.029677224692320876, 0.027131223240968894, 0.024383144740652344, 0.021245155304743375, 0.01794405386879113, 0.014347632788891562, 0.010637289697611955, 0.006747599697386004, 0.0028037445321102044, -0.0011818489296434931, -0.00515085421151638, -0.009002017465669082, -0.012754311851438402, -0.01620678502147323, -0.019465940435647328, -0.022220346759982405, -0.024674013417018838, -0.02639446334952045, -0.027692813749661833, -0.02800520543345714, -0.027758297081010935, -0.02624937202987564, -0.024024084591645806, -0.020242420741043418, -0.015542505027065397, -0.009034359490511169, -0.0012647965330273674, -0.0012804468890246557, -0.001296159332320159, -0.0013119716899741962, -0.001327935384908286, -0.0013439785832496574, -0.001360091775117622, -0.0013762037231243732, -0.001392274487265316, -0.0014082275214679645, -0.0014240132422142894, -0.0014395622179975143, -0.001454820378448915, -0.0014697255638460595, -0.001484222833585092, -0.0014982579663364013, -0.0015117770343765174, -0.001524733387699722, -0.0015370758552288918, -0.0015487656497557732, -0.0015597549884327384, -0.0015700125750107478, -0.0015794947517246793, -0.0015881775625823501, -0.0015960220286219936, -0.001603011410432691, -0.001609111802216845, -0.0016143134903136365, -0.0016185880230317794, -0.0016219325584775602, -0.0016243244122855186, -0.0016257674389771929, -0.0016262449959342446]], "constraint": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.023484259415375942, 0.0232020202711758, 0.02266571018525724, 0.02157549803804424, 0.02024057302121921, 0.018383491150676036, 0.016300474461657145, 0.013748266085288709, 0.010998446344867498, 0.007853709897319767, 0.004549423301431272, 0.0009459401198293258, -0.0027690294801543326, -0.0066658401469344925, -0.010615725615785464, -0.014608205604802131, -0.018584558628493973, -0.022442049465185125, -0.026202827658879733, -0.029660757409433962, -0.03292898935429569, -0.035688075444630256, -0.03814912599179709, -0.03987610828192619, -0.04116771160847842, -0.041507184303919996, -0.04116958345804227, -0.03981866615117834, -0.03701046147997857, -0.03417396172279246, -0.02603671437742066, -0.025055396908712883, 0.002842712802135537, 0.004564911787143319, 0.006884826327662872, 0.006940148907527011, 0.007078792429550493, 0.007702410701266821, 0.008335927445931493, 0.008792982524708611, 0.009244054863143598, 0.009707128272565209, 0.010162760464173751, 0.010581189525329281, 0.010991618637396445, 0.011377344184664313, 0.011754999875756839, 0.012103601444532975, 0.012443992894976816, 0.012755394406333755, 0.013058447303511023, 0.01333140423948434, 0.01359583675196847, 0.013829442671310137, 0.01405434683024299, 0.014247684523135613, 0.014432149934686913, 0.014584436157581693, 0.014727709183049822, 0.014838293503167059, 0.01493976415433261, 0.015008169622097958, 0.015067384526976983, 0.01509330391127834, 0.015109989993076806]], "constraint+r": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.27827295232270866, 0.27799072510310396, 0.27745443720660745, 0.27636427064045477, 0.2750294011185695, 0.2731723967515867, 0.2710894667243674, 0.2685373648332267, 0.26578765955658706, 0.26264305426825607, 0.2593389052436307, 0.2557355723214629, 0.2520207574277194, 0.24812410922470843, 0.24417438823560103, 0.24018207473569753, 0.23620588723992864, 0.23234855731648824, 0.22858793567393898, 0.22513015012911847, 0.22186205430977357, 0.21910308321731584, 0.21664213512909292, 0.21491522497103693, 0.2136236750906849, 0.21328421703098704, 0.21362180320186464, 0.21497266484680805, 0.21778075195790858, 0.22061713396262916, 0.2287540420553322, 0.22973531876440098, 0.11736021141805625, 0.04732052773252468, -0.04702147239602905, -0.049272870985694495, -0.05490970755874752, -0.08027137695252728, -0.10603327572315013, -0.12462136379886418, -0.14296382090562548, -0.16179686958535233, -0.18032452820810507, -0.19734223078574897, -0.21403167270841106, -0.22971917816819126, -0.24507620003831437, -0.2592537139202632, -0.2730954345838461, -0.28576010971311017, -0.29808349511917365, -0.3091844056706605, -0.319937573371295, -0.32943784544782484, -0.33858382981943347, -0.3464461797440681, -0.3539480236370066, -0.36014080522941005, -0.3659676175949522, -0.3704646839739383, -0.3745913860235228, -0.3773730433634762, -0.3797812751972055, -0.3808350600318977, -0.38151385950633276]], "rotated": [[0.0, 0.01562500000006506, 0.03125000000013012, 0.04687500000019512, 0.06250000000026013, 0.07812500000032518, 0.09375000000039024, 0.10937500000045527, 0.1250000000005203, 0.14062500000058537, 0.1562500000006504, 0.1718750000007154, 0.1875000000007804, 0.20312500000084543, 0.2187500000009105, 0.2343750000009755, 0.2500000000010405, 0.26562500000110556, 0.2812500000011706, 0.2968750000012357, 0.3125000000013007, 0.3281250000013657, 0.3437500000014307, 0.3593750000014958, 0.3750000000015609, 0.3906250000016259, 0.4062500000016909, 0.42187500000175593, 0.437500000001821, 0.453125000001886, 0.468750000001951, 0.48437500000200523, 0.5000000000020595, 0.5156250000019619, 0.5312500000018643, 0.5468750000017559, 0.5625000000016476, 0.5781250000015392, 0.5937500000014309, 0.6093750000013225, 0.625000000001214, 0.6406250000011057, 0.6562500000009973, 0.6718750000008888, 0.6875000000007805, 0.7031250000006721, 0.7187500000005638, 0.7343750000004554, 0.7500000000003471, 0.7656250000002386, 0.7812500000001301, 0.7968750000000216, 0.8124999999999132, 0.8281249999998048, 0.8437499999996964, 0.8593749999995988, 0.8749999999995012, 0.8906249999995555, 0.9062499999996096, 0.9218749999996747, 0.9374999999997398, 0.9531249999998048, 0.9687499999998699, 0.9843749999999349, 1.0], [0.2773645142220426, 0.2770060688807109, 0.2765450633890458, 0.2753805502872182, 0.27411721394250915, 0.27219055431830586, 0.270172569837212, 0.2675583549103214, 0.264864124809081, 0.26166785067861936, 0.2584067587750641, 0.25476517535677, 0.25107796436136676, 0.2471595536987249, 0.24321883493795293, 0.23922443966955664, 0.23523535645694285, 0.23139898110803986, 0.22760004918816168, 0.22418987272589588, 0.2208540370406222, 0.21817363982107119, 0.21560955151308286, 0.21399951597823821, 0.21255319001604114, 0.21239267713569004, 0.2124491323385774, 0.21415914374546688, 0.21614573306724275, 0.22018794453786592, 0.22457317893160572, 0.2314676491421537, 0.1866197430100396, 0.04422423108761631, -0.0457764628855414, -0.050601873930320895, -0.055312151321664345, -0.08127631924355917, -0.10710005268126711, -0.12565292652483745, -0.14407914993973828, -0.16285900112620205, -0.18151211254900068, -0.19844439117336296, -0.21525131124255012, -0.23085311952189194, -0.24632886003262447, -0.26041539688483256, -0.27437447721509206, -0.2869448623421712, -0.29938574686970465, -0.3103885533909812, -0.3212595067422845, -0.33065815280881056, -0.339922498263013, -0.3476799486158023, -0.35530076019385726, -0.36138551651158, -0.3673315679605391, -0.37171762380416096, -0.37596331172230896, -0.3786313961517514, -0.38115794999324015, -0.3820960984472562, -0.38289211882890845]]}}} \ No newline at end of file diff --git a/articles/boundary-conditions-on-non-planar-boundaries/figures/banner.jpg b/articles/boundary-conditions-on-non-planar-boundaries/figures/banner.jpg new file mode 100644 index 0000000..25179a9 Binary files /dev/null and b/articles/boundary-conditions-on-non-planar-boundaries/figures/banner.jpg differ diff --git a/articles/boundary-conditions-on-non-planar-boundaries/figures/locking.png b/articles/boundary-conditions-on-non-planar-boundaries/figures/locking.png new file mode 100644 index 0000000..e31caca Binary files /dev/null and b/articles/boundary-conditions-on-non-planar-boundaries/figures/locking.png differ diff --git a/articles/boundary-conditions-on-non-planar-boundaries/figures/rotated-basis.png b/articles/boundary-conditions-on-non-planar-boundaries/figures/rotated-basis.png new file mode 100644 index 0000000..9d9ca00 Binary files /dev/null and b/articles/boundary-conditions-on-non-planar-boundaries/figures/rotated-basis.png differ diff --git a/articles/boundary-conditions-on-non-planar-boundaries/figures/rotated-basis.svg b/articles/boundary-conditions-on-non-planar-boundaries/figures/rotated-basis.svg new file mode 100644 index 0000000..e72529d --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/figures/rotated-basis.svg @@ -0,0 +1,2118 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/articles/boundary-conditions-on-non-planar-boundaries/figures/topography.png b/articles/boundary-conditions-on-non-planar-boundaries/figures/topography.png new file mode 100644 index 0000000..83e3d4b Binary files /dev/null and b/articles/boundary-conditions-on-non-planar-boundaries/figures/topography.png differ diff --git a/articles/boundary-conditions-on-non-planar-boundaries/metadata.yml b/articles/boundary-conditions-on-non-planar-boundaries/metadata.yml new file mode 100644 index 0000000..a745314 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/metadata.yml @@ -0,0 +1,39 @@ +# Validated in CI against schemas/article-metadata.schema.json. +# `pixi run validate` checks this and the cross-file invariants a schema cannot +# express -- that the article file is named .md, that canonical_path +# matches the slug, and that no legacy DOI is ever paired with a new registrant. +id: UWTN 2026-016 +slug: boundary-conditions-on-non-planar-boundaries +title: Boundary conditions on non-planar boundaries +article_type: technical-note +status: published +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliation: Australian National University +publication_date: 2026-08-23 +version: 1.0.0 +# The deposit writes archive_doi and repository_record_id when the note is +# published; leave them out until then. `doi` and `doi_registrant` were here +# once and are not fields the schema knows -- every note made from this +# template failed `pixi run validate` on all three of them. +license: CC-BY-4.0 +canonical_path: /boundary-conditions-on-non-planar-boundaries/ +legacy_paths: [] +# Facets, from vocabulary.yml. Both keys must be present even when empty: a +# note with no subject is normal -- many are purely about method. +subjects: + - surface-processes +methods: + - solvers + - finite-elements + - meshing +ghost_tags: + - Underworld Code +figures: 3 +# `figures` counts the figures in the body; the banner is not one. Unsplash's +# licence asks for the credit wherever the photograph appears, so it lives here +# and is rebuilt into the page with the banner. +banner: figures/banner.jpg +banner_credit: "Photo by Alexander Londoño / Unsplash" +source: native diff --git a/articles/boundary-conditions-on-non-planar-boundaries/references.bib b/articles/boundary-conditions-on-non-planar-boundaries/references.bib new file mode 100644 index 0000000..d8c2b29 --- /dev/null +++ b/articles/boundary-conditions-on-non-planar-boundaries/references.bib @@ -0,0 +1,39 @@ +% Pinned here rather than fetched at build time, following the convention in the +% MMPDE note: a build that depends on doi.org answering can publish a note with a +% broken citation, and a deposited PDF cannot be repaired afterwards. +% +% Every citation in this note is pinned, not only the awkward one. The preview +% build failed on CI with "Citation data from doi.org was not available or +% malformed" for two of them while resolving fine locally: a build that depends +% on doi.org answering will eventually publish a note with a broken citation, and +% a deposited PDF cannot be repaired afterwards. +% +% The Zhong entry has a second reason. The DOI contains parentheses, +% 10.1016/0031-9201(93)90078-N, and MyST's inline DOI citation form stops at the +% first "(" -- the build reported: Could not find DOI "10.1016/0031-9201". A +% bibtex key sidesteps that. +% +% Do not write an at-sign in these comments. The bibtex reader scans for one +% wherever it appears and then expects an entry: an at-sign inside a comment +% here failed the whole file with "expected lbrace, got label". + +@article{Zhong_1993, + title = {Accurate determination of surface normal stress in viscous flow from a consistent boundary flux method}, + volume = {78}, + ISSN = {0031-9201}, + url = {http://dx.doi.org/10.1016/0031-9201(93)90078-N}, + DOI = {10.1016/0031-9201(93)90078-n}, + number = {1-2}, + journal = {Physics of the Earth and Planetary Interiors}, + publisher = {Elsevier BV}, + author = {Zhong, Shijie and Gurnis, Michael and Hulbert, Gregory}, + year = {1993}, + month = jun, + pages = {1--8} +} + +@article{Nitsche_1971, title={Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind}, volume={36}, ISSN={1865-8784}, url={http://dx.doi.org/10.1007/BF02995904}, DOI={10.1007/bf02995904}, number={1}, journal={Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg}, publisher={Springer Science and Business Media LLC}, author={Nitsche, J.}, year={1971}, month=July, pages={9–15} } + +@article{Engelman_1982, title={The implementation of normal and/or tangential boundary conditions in finite element codes for incompressible fluid flow}, volume={2}, ISSN={1097-0363}, url={http://dx.doi.org/10.1002/fld.1650020302}, DOI={10.1002/fld.1650020302}, number={3}, journal={International Journal for Numerical Methods in Fluids}, publisher={Wiley}, author={Engelman, M. S. and Sani, R. L. and Gresho, P. M.}, year={1982}, month=July, pages={225–238} } + +@article{Kramer_2021, title={Analytical solutions for mantle flow in cylindrical and spherical shells}, volume={14}, ISSN={1991-9603}, url={http://dx.doi.org/10.5194/gmd-14-1899-2021}, DOI={10.5194/gmd-14-1899-2021}, number={4}, journal={Geoscientific Model Development}, publisher={Copernicus GmbH}, author={Kramer, Stephan C. and Davies, D. Rhodri and Wilson, Cian R.}, year={2021}, month=Apr, pages={1899–1919} } diff --git a/classification.yml b/classification.yml index 65201be..4925dd5 100644 --- a/classification.yml +++ b/classification.yml @@ -322,3 +322,7 @@ setting-up-full-multigrid: article_type: technical-note subjects: [mantle-convection] methods: [solvers, meshing, parallel-hpc] +boundary-conditions-on-non-planar-boundaries: + article_type: technical-note + subjects: [surface-processes] + methods: [solvers, finite-elements, meshing] diff --git a/scripts/build_reader_pages.py b/scripts/build_reader_pages.py index bd3d771..b6dfa75 100644 --- a/scripts/build_reader_pages.py +++ b/scripts/build_reader_pages.py @@ -49,9 +49,13 @@ } } * { box-sizing: border-box; } +/* height, not just min-height: a percentage height inside a flex column needs + a DEFINITE container, and without it the viewer below fell back to its + min-height (512px) however tall the window was -- a letterbox showing a + few lines of a page scaled to the full window width. */ body { margin: 0; background: var(--uwtn-paper); color: var(--uwtn-ink); font-family: var(--uwtn-sans); font-size: 16px; line-height: 1.5; - display: flex; flex-direction: column; min-height: 100vh; } + display: flex; flex-direction: column; height: 100vh; min-height: 100vh; } header { border-bottom: 1px solid var(--uwtn-rule); padding: 1.1rem 1.4rem; } .wrap { max-width: 68rem; margin: 0 auto; width: 100%; } .kicker { font-size: .78rem; letter-spacing: .09em; text-transform: uppercase; @@ -66,9 +70,11 @@ .actions a:hover { border-color: var(--uwtn-accent); color: var(--uwtn-accent); } .actions a.primary { background: var(--uwtn-accent); border-color: var(--uwtn-accent); color: #fff; } -.reader { flex: 1 1 auto; min-height: 32rem; } -.reader object, .reader iframe { display: block; width: 100%; height: 100%; - min-height: 32rem; border: 0; } +/* min-height: 0 lets a flex item shrink below its content; without it the + viewer cannot give the window's height back and the page scrolls instead. */ +.reader { flex: 1 1 auto; min-height: 0; display: flex; } +.reader object, .reader iframe { flex: 1 1 auto; display: block; width: 100%; + height: 100%; min-height: 20rem; border: 0; } .fallback { padding: 2rem 1.4rem; color: var(--uwtn-muted); } """ @@ -94,7 +100,7 @@
- +

This browser will not display a PDF here — most phones will not.

Open the PDF or diff --git a/scripts/new_article.py b/scripts/new_article.py index a010c81..c60324f 100644 --- a/scripts/new_article.py +++ b/scripts/new_article.py @@ -17,6 +17,7 @@ import datetime import pathlib import re +import subprocess import sys ROOT = pathlib.Path(__file__).resolve().parent.parent @@ -44,6 +45,37 @@ def load_authors(): return registry +def ids_on_branches(year): + """Article numbers claimed on a branch, merged or not. + + The rest of the allocator reads the working tree, which cannot see a number + claimed by a note still in review on its own branch. Two notes drafted in + parallel were therefore both offered the same number, and UWTN 2026-012 was + claimed twice before anyone noticed. + + Every local and remote-tracking ref is searched, so this covers open pull + requests as far as they have been fetched. A note on a branch that has never + been pushed to a remote this checkout tracks is still invisible, which is why + `pixi run validate` checks for duplicates as well. + """ + def git(*args): + return subprocess.run(("git",) + args, cwd=str(ROOT), + capture_output=True, text=True) + + refs = git("for-each-ref", "--format=%(refname)", "refs/heads", "refs/remotes") + if refs.returncode != 0: + return set() # no git, or not a checkout: the tree is all we have + + used = set() + for ref in refs.stdout.split(): + found = git("grep", "-h", "-E", r"^id: +UWTN +[0-9]{4}-[0-9]{3}", + ref, "--", "articles/*/metadata.yml") + for match in re.finditer(r"UWTN\s+(\d{4})-(\d{3})", found.stdout): + if match.group(1) == str(year): + used.add(int(match.group(2))) + return used + + def next_article_id(year): """Allocate an ID that no existing article uses. @@ -52,7 +84,7 @@ def next_article_id(year): number already present in that year keeps a new note clear of anything the backfill will produce, and an ID that has been published never moves. """ - used = set() + used = ids_on_branches(year) for meta in ARTICLES.glob("*/metadata.yml"): match = re.search(r"^id:\s*(UWTN\s+(\d{4})-(\d{3}))\s*$", meta.read_text(encoding="utf-8"), re.M) diff --git a/static/uwtn.css b/static/uwtn.css index 63106de..710aae6 100644 --- a/static/uwtn.css +++ b/static/uwtn.css @@ -273,6 +273,7 @@ article.content figcaption p { height: 100%; max-height: 260px; object-fit: cover; + margin: 0; /* Decorative, so it should not compete with the words. */ filter: saturate(0.82); transition: filter 0.25s ease, transform 0.4s ease; @@ -297,6 +298,11 @@ article.content figcaption p { width: 100%; max-height: 300px; object-fit: cover; + /* The theme's prose rule gives every image that is not inside a

+ `margin: 2em 0`. Inside the banner's frame that reads as a white band + above and below the picture -- invisible on a banner with a white + background, obvious on a photograph. */ + margin: 0; } @media (max-width: 720px) { diff --git a/templates/article-template/ARTICLE.md b/templates/article-template/ARTICLE.md index 000571b..3865815 100644 --- a/templates/article-template/ARTICLE.md +++ b/templates/article-template/ARTICLE.md @@ -26,12 +26,18 @@ but not this corner of Underworld. ## A section -Maths is plain LaTeX, inline as $\nabla \cdot \boldsymbol{\sigma} = 0$ or -displayed: +Maths is plain LaTeX, inline as $\nabla \cdot \boldsymbol{\sigma} = 0$. Display +maths goes in a `math` directive rather than a `$$` block: both build, but the +fence survives an editor that is not a MyST editor, and it takes a label +without the two traps `$$` has (see CONTRIBUTING.md). -$$ +```{math} +:label: eq-momentum -\nabla \cdot \boldsymbol{\sigma}(u, \nabla u) - \mathbf{f}(u, \nabla u) = 0 -$$ +``` + +Refer to it as {eq}`eq-momentum`. Drop the `:label:` line if nothing refers to +the equation. Code blocks are fenced and language-tagged. Keep lines under ~84 characters: that is the 95th percentile of this corpus and it is what fits the PDF measure diff --git a/tests/test_migration.py b/tests/test_migration.py index 35e7cce..561edb4 100644 --- a/tests/test_migration.py +++ b/tests/test_migration.py @@ -2214,3 +2214,66 @@ def test_reader_link_is_visible_and_click_safe(): style = (pathlib.Path(__file__).resolve().parent.parent / "static" / "uwtn.css").read_text(encoding="utf-8") assert ".uwtn-pdf-link" in style, "the injected link has no style" + + +def test_every_eq_reference_has_a_labelled_equation(): + """`{eq}`name`` needs `$$ (name)` on the CLOSING delimiter line. + + A label written on a line of its own renders as literal text, the reference + resolves to nothing, and `myst build --typst` fails the article outright + with "label does not exist" -- an HTML-only check would pass it. That is + how it reached a pushed branch once already. + """ + for path in sorted((ROOT / "articles").glob("*/*.md")): + text = path.read_text(encoding="utf-8") + # [ \t]* and not \s*: \s crosses newlines, so the very form this test + # exists to catch -- the label on its own line -- would match. + labels = set(re.findall(r"^\$\$[ \t]*\(([^)]+)\)[ \t]*$", text, re.M)) + labels |= set(re.findall(r"^\s*:label:\s*(\S+)\s*$", text, re.M)) + labels |= set(re.findall(r"\\label\{([^}]+)\}", text)) + for name in set(re.findall(r"\{eq\}`([^`]+)`", text)): + assert name in labels, ( + "%s references {eq}`%s` but no equation carries that label. " + "The label belongs on the closing $$ line: `$$ (%s)`." + % (path.name, name, name)) + + +# A line inside display math that markdown would claim for itself if the block +# is not a block of its own: a list bullet, a heading, a quote, a fence. +_MARKDOWN_LINE = re.compile(r"^\s*([-+*>#]|\d+[.)])(\s|$)") + + +def test_display_math_that_needs_a_blank_line_has_one(): + """`$$` opening on the line after a paragraph is a lazy continuation. + + For one-line arithmetic that renders anyway, which is why the corpus has + examples of it. It breaks when a line INSIDE the block is something + markdown wants: `- \\int ...` opens a list, the paragraph ends, the math + never closes, and the page shows raw LaTeX with half the equation in a + bullet. `myst build` reports nothing at all -- no warning, exit 0 -- so + this test is the only guard. It cost the boundary-conditions note a + published-looking draft once. + """ + for path in sorted((ROOT / "articles").glob("*/*.md")): + lines = path.read_text(encoding="utf-8").split("\n") + in_code = False + opened = None + for i, line in enumerate(lines): + if line.startswith("```"): + in_code = not in_code + continue + if in_code: + continue + if opened is None: + if line.rstrip() == "$$": + opened = i + elif line.startswith("$$"): + body = lines[opened + 1:i] + risky = [b for b in body if _MARKDOWN_LINE.match(b)] + if risky and opened > 0 and lines[opened - 1].strip(): + raise AssertionError( + "%s line %d: display math opens straight after text and " + "contains %r, which markdown will take for itself. Put a " + "blank line before the `$$`." + % (path.name, opened + 1, risky[0].strip()[:40])) + opened = None