diff --git a/docs/advanced/figures/material_index.png b/docs/advanced/figures/material_index.png new file mode 100644 index 000000000..3fcb3826f Binary files /dev/null and b/docs/advanced/figures/material_index.png differ diff --git a/docs/advanced/figures/material_share.png b/docs/advanced/figures/material_share.png new file mode 100644 index 000000000..0a1b28a40 Binary files /dev/null and b/docs/advanced/figures/material_share.png differ diff --git a/docs/advanced/figures/repopulation_counts.png b/docs/advanced/figures/repopulation_counts.png new file mode 100644 index 000000000..a1a6545b8 Binary files /dev/null and b/docs/advanced/figures/repopulation_counts.png differ diff --git a/docs/advanced/figures/repopulation_particles.png b/docs/advanced/figures/repopulation_particles.png new file mode 100644 index 000000000..c028586dd Binary files /dev/null and b/docs/advanced/figures/repopulation_particles.png differ diff --git a/docs/advanced/index.md b/docs/advanced/index.md index 3a2605939..51e5c081e 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -75,6 +75,13 @@ Dynamic remeshing and adaptive refinement strategies. **[→ Mesh Adaptation](mesh-adaptation.md)** +### Particles: Population Control and Materials +Keep every cell sampled as the flow deforms the swarm, and declare materials +on a `MaterialSwarm` so an interface stays where the particles put it — read +at the integration points, with no level sets in the model script. + +**[→ Particle Population and Materials](particle-population-and-materials.md)** + ### Semi-Lagrangian Time Integration (SLCN / SL-BDF2) How `AdvDiffusionSLCN` discretizes advection–diffusion in time: the BDF time-derivative and Adams-Moulton/θ flux knobs, and how to pair them @@ -138,6 +145,7 @@ stress-visualisation custom-meshes curved-boundary-conditions mesh-adaptation +particle-population-and-materials semi-lagrangian-time-integration eulerian-advection-diffusion eulerian-navier-stokes diff --git a/docs/advanced/particle-population-and-materials.md b/docs/advanced/particle-population-and-materials.md new file mode 100644 index 000000000..2cdf5d2a2 --- /dev/null +++ b/docs/advanced/particle-population-and-materials.md @@ -0,0 +1,508 @@ +--- +title: "Particles: population control and materials" +--- + +# Particles: population control and materials + +Two things a particle method has to get right, and how Underworld3 does them. + +**Population control** keeps every cell holding enough particles to be worth +integrating, however hard the flow deforms the swarm. + +**A material index read where the assembler actually looks** keeps a material +interface where the particles put it, instead of smearing it over a cell. + +Both are demonstrated below with runnable scripts. Neither changes how you +write anything: a particle field reaches the mathematics through a proxy, and +what the proxy changes is where the sampling happens, not what you can say. + +## Population control + +A swarm is not conserved cell by cell. Particles leave through outflow +boundaries and nothing arrives through inflow boundaries unless you put it +there; strong deformation sweeps them out of some cells and piles them into +others. A cell that ends up with too few particles cannot support the +reconstruction that gives the mesh its material properties, and one with none +carries no information at all. + +`Swarm.repopulate()` takes a census of the particles per cell and refills the +starved ones from the cell's own lattice, choosing the lattice points furthest +from the particles already there. A new particle takes each swarm variable +from its neighbours: a bounded reconstruction for a continuous field, and its +nearest neighbour's value whole for an integer field, because the average of +two material labels is not a material label. + +Set `swarm.population_control` and every `advection()` ends with a +repopulation, before anything reads the swarm again: + +```python +swarm.population_control = dict() # back to the populate() density +swarm.population_control = dict(min_per_cell=8) # or a floor you choose +swarm.population_control = dict(min_per_cell=8, values={T: 0.0}) # an inflow datum +``` + +Or call it yourself, which is what you want if the swarm moves by some route +of your own: + +```python +added, removed = swarm.repopulate(min_per_cell=8) +``` + +`values` overrides what a new particle receives, per variable, with a constant +or a callable of the coordinates. That is how an inflow boundary gets the +right material or temperature rather than a copy of whatever drifted past. +`max_per_cell` will also thin over-full cells, by dropping the particles +closest to a neighbour; it is off by default because discarding particles +costs accuracy. + +### Demonstration + +`docs/examples/utilities/intermediate/Ex_Swarm_Population_Control.py`. Pure-shear extension +$\mathbf{v} = (x, -y)$ on a fixed mesh, so the side walls are outflow and the +top and bottom are inflow, with a marker layer through the middle. Particles +that leave are deleted (`mesh.return_coords_to_bounds = None`), which is what +an open boundary means. + +```{figure} figures/repopulation_particles.png +:alt: Particles in an extending box, with and without population control + +Without population control (top) the box drains: after two time units 1020 of +the original 7640 particles remain and 626 of the 764 cells are empty, and the +marker layer has shredded into streaks. With it (bottom) no cell is ever +empty, and the layer thins as the extension requires rather than falling apart. +``` + +```{figure} figures/repopulation_counts.png +:alt: Particle count and starved-cell count against time + +The count in the box and the number of starved cells. Population control is +not fighting the outflow, which is physical and correct; it is refilling the +inflow side, where the flow brings nothing. +``` + +### What it is worth, in a number + +The figures above are particle scatters. The field-level cost depends on how +the material is read, and it is worth knowing which case needs the refilling. +Pure shear thins a marker layer by $e^{-t}$, so after $t = 2$ its area should +be 13.5% of where it started; anything else is the representation failing: + +| `proxy_location` | population control | layer area vs exact | +|---|---|---| +| `"cells"` | on | 1.03x | +| `"cells"` | off | **5.03x** | +| `"integration_points"` | on | 1.06x | +| `"integration_points"` | off | 1.07x | + +Starvation wrecks the per-cell fit, which needs particles *in that cell*. The +nearest-particle mapping degrades gracefully — an empty cell still finds a +plausible particle nearby — so with the default mapping the field-level answer +barely moves, and losing particles shows up only in the swarm's own +bookkeeping. Population control is cheap and always safe; `"cells"` is where +it is *necessary*. + +## Materials + +Name the materials, say where they are, and stop. + +```python +materials = uw.swarm.MaterialSwarm(mesh, fill_param=3) + +materials.add("mantle", shear_viscosity_0=1.0, density=3300) +materials.add("slab", shear_viscosity_0=1.0e3, density=3400) + +materials["slab"] = mesh.X[1] > 0.53 + +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.materials = materials +stokes.bodyforce = -materials.density * mesh.CoordinateSystem.unit_e_1 +``` + +That is the whole interface. `stokes.materials = materials` sets every +constitutive-model parameter the materials declare *and* the model recognises, +by name — here `shear_viscosity_0`. A property the model does not own, +`density`, is not pushed anywhere; it is available as a blended symbol for the +model script to use where it belongs. Properties can be changed and regions +repainted afterwards: the blend is symbolic and the push repeats. + +A `MaterialSwarm` **is** a `Swarm`, so everything on the previous page still +applies to it. Declare any extra per-particle state **before** the materials +are first read, though — the first read allocates the particles, and a swarm +cannot gain variables once it holds them: + +```python +materials = uw.swarm.MaterialSwarm(mesh, fill_param=3) +materials.add("mantle", shear_viscosity_0=1.0) +materials.add("slab", shear_viscosity_0=1.0e3) + +strain = uw.swarm.SwarmVariable("eps", materials, 1, # per-particle state + proxy_location="integration_points", + proxy_sampling="share") + +materials["slab"] = mesh.X[1] > 0.53 # reads nothing yet +stokes.materials = materials # this allocates + +materials.population_control = dict(min_per_cell=8) +materials.advection(v.sym, dt) +``` + +```{warning} +Allocating the level sets is **collective**: every rank must reach it +together. It happens on the first *read* of any material property — +`materials.density`, `materials["x"].mask`, `materials.index`, or +`stokes.materials = ...` — so a read inside a rank-local branch +(`if rank == 0: ...`) deadlocks. Call `materials.build()` at a point where +every rank is together if the ordering is ever in doubt. +``` + +Regions can be a symbolic condition on the mesh coordinates (and `&`, `|`, `~` +combinations of them), a boolean array over the particles, or a callable of the +coordinate array. `materials.add(..., where=...)` is the same thing said in one +line. + +Assignment **replaces**: `materials["slab"] = A` followed by +`materials["slab"] = B` leaves the material at B, and anything it held outside +B reverts to the first material declared. Across *different* materials the +later assignment wins where they overlap. + +A few things are refused rather than half-done, because each of them used to +produce a plausible-looking wrong answer: + +| | | +|---|---| +| adding or deleting a material after the first read | the index of a material *is* which level set is its | +| two distributions sharing an explicit `name` on one mesh | they would share level sets, silently | +| `mixing(...)` naming a property no material declares | it would have no effect | +| harmonic mixing of a property some material sets to zero | `1/Σ(φᵢ/vᵢ)` divides by it | +| a mesh label value that does not exist | asking PETSc for it aborts the run | + +Dimensional values are non-dimensionalised when they are **read**, not when +they are declared, so a registry written before the model's reference +quantities are set means the same thing as one written after. + +### What a material is, and where it is + +Those are two questions, and Underworld3 keeps them apart. + +**What** is a `uw.MaterialRegistry` entry: a name and a table of properties, +with optional description and reference. A registry knows nothing about +geometry, so it can be written before there is a mesh, shared between models, +and exported and read back as plain data. + +```python +rocks = uw.MaterialRegistry() +rocks.add("mantle", shear_viscosity_0=1.0, density=3300, + description="upper mantle", reference="Turcotte & Schubert (2014)") +rocks.add("slab", shear_viscosity_0=1.0e3, density=3400) +``` + +**Where** is a *distribution*, and there are two. `MaterialSwarm` carries the +materials on particles, so they advect with the flow. `MaterialRegions` ties +them to the mesh: + +```python +materials = uw.MaterialRegions(mesh, registry=rocks) +materials["slab"] = "Slab" # a gmsh physical group +materials["slab"] = mesh.X[1] > 0.53 # or a geometric condition + +stokes.materials = materials +``` + +Both take the same registry, declare materials the same way, and hand the same +thing to the solver. `materials.add(...)` on either is shorthand that writes +into its registry, so a two-material script never has to mention one. + +| | `MaterialSwarm` | `MaterialRegions` | +|---|---|---| +| the material moves | yes | no | +| level sets are | sampled from the particles | exact, 0 or 1 | +| needs population control | yes | no | +| interface position | sub-cell, to the particle density | sub-cell, to the rule | +| cost per step | a proxy fill | nothing; built once | + +Use regions when the geometry is fixed — a layered model, an inclusion, a +basin. It is exact and free. Use the swarm when the material is carried by +the flow, or when it has to carry per-particle history as well. + +### A property can be a law, not just a number + +This is the whole reason for the machinery underneath: + +```python +materials.add("crust", shear_viscosity_0=eta_0 * sympy.exp(-T.sym[0])) +``` + +There is no number to store at an integration point for that material, so it +cannot be handled by sampling a viscosity field. It is carried symbolically and +combined with the other materials' laws by the level sets described below. + +### There is no route that puts the property on the particles + +Carrying the viscosity itself on the particles and sampling it looks simpler +and is the wrong shape. PETSc calls the compiled pointwise functions with +values tabulated at the rule points, so a solver handed a sampled *answer* +cannot reach the parts of the constitutive law it needs — the tangent, the +yield surface, the history update. Keeping the law symbolic and letting the +materials weight it is what makes a material model composable. (The mechanical +symptom is easy to see too: the linear-exact reconstruction a smooth field +wants overshoots a factor-1000 viscosity jump to −219, and a negative viscosity +is not a viscosity. Asking for nearest-particle sampling on a plain +`SwarmVariable` raises, and the error says so.) + +### What is underneath + +A `MaterialSwarm` carries one integer label per particle and presents the mesh +with one *level set* per material — a partition of unity — so that a property +is `Σ φᵢ · valueᵢ`. Where the masks are 0 or 1 that sum is a select, and if +every property were a number a single stored coefficient field would do the +same job; the sum earns its place the moment a property is a law, because there +is then no other way to combine N expressions into one symbol the assembler can +compile. + +That machinery is `uw.swarm.IndexSwarmVariable` and its `createMask`. Models +written before `MaterialSwarm` use it directly and still work; new models +should not need to see it. The two things worth knowing about it are the two +arguments `MaterialSwarm` passes through, below. + +#### `proxy_location` — where the level sets live + +Where the level sets are stored is what the assembler actually reads: + +```python +materials = uw.swarm.MaterialSwarm(mesh, proxy_location="cells") +``` + +| `proxy_location` | the level sets are | at an interface | +|---|---|---| +| `"integration_points"` **(default)** | a value at each point of the quadrature rule | exactly 0 or 1, and the interface keeps its sub-cell position | +| `"cells"` | a polynomial material fraction per cell | a sharp step at cell edges, with a gradient inside the cell | +| `"nodes"` | a continuous field per material | a node on the interface averages both, so the cells either side see a viscosity that is neither | + +The default is the integration points: it is the classic particle-in-cell +material mapping of Ellipsis and Underworld, and it is measurably the best of +the three. The nodal option is kept for continuity with existing models, but +its smear is about one cell wide *however many particles you add* — that is a +property of the basis, not of the swarm — and the current distance-weighted +fill actually *widens* the band as the swarm is refined. + +#### `proxy_sampling` — what each point reads + +How the particles in a cell become the value at each of its integration +points. Applies only to `proxy_location="integration_points"`. + +| `proxy_sampling` | each point takes | masks | +|---|---|---| +| `"nearest"` **(default)** | the material of its nearest particle | exactly 0 or 1 | +| `"share"` | the material fractions of the particles it speaks for — those whose nearest integration point *within their own cell* is this one | fractional where a cell is crossed | + +`"nearest"` is sharp and assumption-free: nothing is mixed, so no mixing rule +is implied. It sub-samples, though — at ten particles per cell and six rule +points, most particles never reach the assembly. + +`"share"` is the cell-restricted Voronoi share, and it is the closest thing to +the true Voronoi integration that a fixed quadrature rule allows: the rule +points partition their own cell between them, every particle lands in exactly +one part, and nothing crosses a cell boundary. For an *identity* that buys +fractional masks in the crossed cells, and then the answer depends on how you +blend — see the trap below. + +#### What the choice is worth + +Interface at $y = 0.53$ on an irregular mesh, so no scheme can be exact (the +P2 velocity cannot hold a kink inside a cell), viscosity contrast 1000, +against the exact layered Couette profile: + +| particles per cell | `"nearest"` | `"share"` + `createMask` | `"share"`, blended harmonically | +|---|---|---|---| +| 3 (2 420) | 3.48e-2 | 3.53e-2 | 1.83e-2 | +| 8 (10 890) | **1.85e-2** | 3.53e-2 | 1.27e-2 | +| 15 (32 912) | 1.85e-2 | 3.61e-2 | **1.24e-2** | + +```{figure} figures/material_share.png +:alt: Identity error against particle density, and the history overshoot + +Left: only `"nearest"` and the harmonically-blended share improve as particles +are added; the arithmetic (Voigt) blend of fractional masks is flat. Right: the +share is bounded by the particle values because it is an average of them, while +the reconstruction's overshoot grows with the swarm. +``` + +Three things to read off it. + +**`"nearest"` converges to a floor and then stops.** By about eight particles +per cell the error stops moving: what limits it is the quadrature rule and the +velocity space, not the swarm. Past that point, refine the mesh — more +particles buy nothing. + +**Fractional masks need a mixing rule you have chosen.** `createMask` blends +*arithmetically*, which for a flux at a common strain rate is a Voigt average, +and on a sharp contrast a Voigt average does not converge with particle +density at all — the middle column is flat. The same masks blended +harmonically (Reuss, `1.0 / material.createMask([1/η₀, 1/η₁])`) converge, and +past the `"nearest"` floor. So fractions are worth having when the material +genuinely *is* a sub-cell mixture and you have picked the rule on physical +grounds; for an interface, `"nearest"` says the true thing (nothing is mixed) +without you having to. + +**Cost is not the discriminator.** Filling the level sets took 1.0 / 1.7 / +3.7 ms for `"nearest"` and 0.7 / 1.3 / 3.9 ms for `"share"` at the three +densities. + +With the interface on mesh edges instead, the exact velocity lies in the P2 +space and the only error left is the material representation: + +| `proxy_location` | assembled $\int \eta$ (exact 500.5) | velocity $L_2$ error | +|---|---|---| +| `"nodes"` | 500.5000 | 8.0e-2 | +| `"integration_points"` | 500.5000 | **1.8e-7** | +| `"cells"` | 500.5000 | **1.8e-7** | + +All three integrate the viscosity correctly *in the mean*, which is exactly +why a bulk diagnostic cannot see the difference. Only the placement differs, +and the placement is what the solve feels. + +```{figure} figures/material_index.png +:alt: The material mask across the interface and the resulting velocity error + +Left: the upper-material mask along a line crossing the interface, as the weak +form sees it. The nodal level set ramps linearly across a whole cell; the other +two are a step in the right place and lie on top of each other. Right: the +resulting error in the velocity against the exact layered flow. +``` + +`docs/examples/utilities/intermediate/Ex_Swarm_Material_Index.py` runs it. + +### Mixing, when the masks are fractional + +`materials.mixing(shear_viscosity_0="harmonic")` chooses how a property is +blended. With the default sampling exactly one mask is 1 at every integration +point, so every rule gives the same answer and there is nothing to choose. With +`proxy_sampling="share"` there is, and the table above is the reason to choose +it on physical grounds rather than by trying both. + +## Material state and history + +Identity is not the only thing particles carry. A viscoelastic stress, an +accumulated strain, a damage variable — these are *state*, earned by being +advected, and every particle's is different. They ride on ordinary +`SwarmVariable`s, and they get the same choice of where the assembler reads +them and how: + +```python +stress = uw.swarm.SwarmVariable( + "tau", materials, (2, 2), # the MaterialSwarm is the swarm + proxy_location="integration_points", + proxy_sampling="share") # every particle's history contributes +``` + +`Lagrangian_Swarm` takes the same argument, so a viscoelastic stress history +carried on the particles is read the same way: + +```python +DFDt = uw.systems.ddt.Lagrangian_Swarm( + swarm=materials, psi_fn=sympy.Matrix.zeros(2, 2), + vtype=uw.VarType.SYM_TENSOR, degree=1, continuous=False, + order=2, step_averaging=1, + proxy_location="integration_points", proxy_sampling="share") +``` + +Build it before anything reads the materials: it adds swarm variables, and a +swarm cannot gain variables once it holds particles. + +For state, `"share"` is the one to reach for, and for a different reason than +it was rejected for identity: + +- **It uses the whole swarm.** Each point averages the particles it represents + rather than adopting one of them, so a history that varies within a cell is + represented by all of it. This is the "more PIC than not" part: the + quadrature rule is fixed, but what it reads can still be a weighted account + of every particle. +- **It cannot cross a material boundary.** The default `"reconstruct"` gathers + from the nearest particles by distance, and that stencil ignores cell walls; + across a jump it both smears and overshoots. Measured on a discontinuous + particle field read at the integration points (interface at $y = 0.53$, + $h = 0.1$): + + | particles per cell | `"reconstruct"` range | `"share"` range | + |---|---|---| + | 3 | −0.051 … +1.137 | 0.000 … 1.000 | + | 8 | −0.068 … +1.108 | 0.000 … 1.000 | + | 15 | −0.109 … +1.097 | 0.000 … 1.000 | + + The share is bounded by the particle values by construction — it is an + average of them — so it cannot invent a stress the swarm never held. The + reconstruction's overshoot *grows* as particles are added. + +Keep `"reconstruct"` (the default) for a field that really is smooth: it is +exact for linear fields, where the share carries a small averaging error. + +**Cost.** The share needs to know each particle's cell, and locating particles +is the expensive half: 18.8 ms for 32 912 particles against 3.2 ms for the +share itself, on a mesh of 242 cells. That location is cached on the swarm +until the particles move and is reused by every share variable and by +`repopulate`, so a model with several history fields and population control +pays it once per step. For comparison, the `"reconstruct"` path costs 8.1 ms +per step on the same swarm (its cached operator is geometry-only, so it is +rebuilt whenever the particles move). + +## What can be said about a particle field + +A particle field is a first-class citizen of the symbolic algebra: it carries a +symbol, and that symbol goes wherever a mesh variable's symbol goes. The one +exception is a derivative of the integration-point form. + +| | `"nodes"` | `"integration_points"` | `"cells"` | +|---|---|---|---| +| arithmetic with mesh variables and `sympy` | yes | yes | yes | +| `uw.maths.Integral` | yes | yes | yes | +| `uw.function.evaluate` anywhere | yes | yes | yes | +| projection onto a mesh variable | yes | yes | yes | +| viscosity, body force, any solver term | yes | yes | yes | +| a **gradient** of the expression | yes | refused | yes | + +The element that holds a value at each integration point has no gradient to +give, so the compiler refuses one rather than returning the silent zero its +tabulation would produce. The gradient is still available from a projection of +the same data, and `"cells"` is that projection: its level sets are a +least-squares polynomial per cell, so they differentiate directly and with no +global solve. Recovering $\partial_x$ of a quadratic particle field: + +| where the proxy lives | gradient error | +|---|---| +| `"nodes"`, degree 1 | 2.1e-3 | +| `"nodes"`, degree 2 | 1.1e-4 | +| `"cells"`, degree 1 | 1.3e-3 | +| `"cells"`, degree 2 | **2.4e-7** | +| a global L2 projection onto P2, then differentiate | 1.1e-4 | +| `"integration_points"` | refused | + +The per-cell fit is the most accurate of them because it is local and exact for +polynomials up to its degree, and no further projection follows it. + +**Two paths, and the difference is deliberate.** A weak form and a query are +not the same thing: + +- **In a weak form**, a derivative of an integration-point field is *refused*. + Answering would mean either the silent zero its own tabulation gives, or a + reconstruction chosen behind your back and paid for at every assembly. + Which discretisation the gradient comes from is a modelling decision, so it + stays yours: build the variable with `proxy_location="cells"` and the level + sets are already polynomials. +- **`uw.function.evaluate`** is a query, and it *answers*. It fits the + integration-point values cell by cell and differentiates that, once, for + this call. The result converges (2.4e-3, 6.1e-4, 2.6e-4 for a quadratic + field as the cell size halves from 1/5 to 1/20) but it is a *recovered* + gradient, so treat it as a diagnostic rather than as the field's own + derivative. The `"cells"` route reaches 2.4e-7 on the same field. + +So the rule is: sample at the integration points when you want a value placed +exactly, fit per cell when you want to differentiate, and expect `evaluate` to +help you look at a gradient either way. + +## What this rests on + +The reconstruction behind `"cells"`, the delta element behind +`"integration_points"`, the share, and the guards that keep them well posed +are described in {doc}`../developer/subsystems/integration-point-variables`. +The same machinery carries semi-Lagrangian and fully Lagrangian histories, +including the viscoelastic stress history. diff --git a/docs/api/materials.md b/docs/api/materials.md index ad8476be5..5ca84b32d 100644 --- a/docs/api/materials.md +++ b/docs/api/materials.md @@ -1,27 +1,72 @@ # Materials +Underworld3 keeps two questions apart: **what** a material is, and **where** it +is. + +*What* is a {py:class}`~underworld3.MaterialRegistry` entry — a name and a +table of properties, where a property may be a number, a quantity, or a *law* +such as `eta_0 * sympy.exp(-T.sym[0])`. A registry knows nothing about +geometry, so it can be built before there is a mesh and shared between models. + +*Where* is a **distribution**, and there are two: + +- {py:class}`~underworld3.swarm.MaterialSwarm` — carried by particles, so the + materials advect with the flow. +- {py:class}`~underworld3.MaterialRegions` — tied to the mesh, from gmsh + physical groups or from a geometric condition. Exact, and needs no particles. + +Both present the same face to a solver — `stokes.materials = materials`, which +sets every constitutive-model parameter the materials declare and the model +recognises — and both build it from the same partition of unity, one level set +per material. + +See {doc}`../advanced/particle-population-and-materials` for the user guide. + ```{eval-rst} .. automodule:: underworld3.materials :no-members: ``` -## Material Property +## Defining materials ```{eval-rst} +.. autoclass:: underworld3.MaterialRegistry + :members: + :show-inheritance: + +.. autoclass:: underworld3.MaterialDefinition + :members: + :show-inheritance: + .. autoclass:: underworld3.MaterialProperty :members: :show-inheritance: ``` -## Material Registry +## Distributing materials ```{eval-rst} -.. autoclass:: underworld3.MaterialRegistry +.. autoclass:: underworld3.swarm.MaterialSwarm + :members: + :show-inheritance: + +.. autoclass:: underworld3.MaterialRegions :members: :show-inheritance: + +.. autoclass:: underworld3.materials.MaterialDistribution + :members: + :show-inheritance: + +.. autoclass:: underworld3.materials.BoundMaterial + :members: ``` -## Multi-Material Models +## Different laws per material + +A distribution blends per-material *parameter values* into one constitutive +model. When the materials need genuinely different constitutive **laws** — +one viscous, one viscoelastic — compose the models instead: ```{eval-rst} .. autoclass:: underworld3.MultiMaterialConstitutiveModel diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index d313da6ce..a8e106e8e 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -69,6 +69,32 @@ symbol onto a `MeshVariable` explicitly with `SNES_Projection`; the projection of the field is an ordinary weak form and is exact for data that the target space can represent. + +### The derivative: refused in a weak form, recovered by `evaluate` + +An integration-point variable's tabulated gradient is identically zero, so a +derivative of its symbol would be a silent zero. The two paths are handled +differently on purpose: + +- **Code generation for a weak form** (`utilities/_jitextension.py`, + `_no_derivative`) raises. A hidden reconstruction inside a residual would be + a per-assembly cost and would decide a discretisation on the user's behalf. + The message names the remedy: `proxy_location="cells"`, whose level sets are + per-cell polynomials and differentiate directly. +- **`uw.function.evaluate`** (`function/_function.pyx`, + `_integration_point_sources_to_cell_fit`) substitutes any integration-point + source appearing under a derivative by a per-cell least-squares fit of its + own values, then lets the ordinary derivative machinery run. The fit is + allowed to be exactly determined (`nmin = Nb`) because the rule is unisolvent + for that degree; the default `Nb + 2` would send every cell to the linear + patch and leave the recovered gradient first order. + +Measured on `x^2 + 2y` carried at the integration points, the recovered +gradient converges: 2.4e-3, 6.1e-4, 2.6e-4 at cell sizes 1/5, 1/10, 1/20. The +direct `"cells"` route (degree 2, fitted from particles) gives 2.4e-7 on the +same field, because it is exact for a quadratic and nothing is projected +afterwards. + ## Guards The field has no gradient (its tabulated derivative is identically zero), so @@ -335,12 +361,20 @@ nodes and back. ```python swarm = uw.swarm.Swarm(mesh) -M = uw.swarm.SwarmVariable("M", swarm, 1, proxy_location="integration_points") +tau = uw.swarm.SwarmVariable("tau", swarm, (2, 2), + proxy_location="integration_points") swarm.populate(fill_param=3) -M.data[:, 0] = ... # per particle -stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_0 * M.sym[0] + eta_1 * (1 - M.sym[0]) +tau.data[...] = ... # per particle ``` +For a **material**, this is not the entry point: use `uw.swarm.MaterialSwarm`, +which owns an `IndexSwarmVariable` whose level sets live at the integration +points by default and blends the declared properties by the resulting +partition of unity. See +{doc}`../../advanced/particle-population-and-materials` for why the direct +route is not offered — sampling a property field hands the solver an answer +where it needs a constitutive law. + The reconstruction itself is unchanged (a linear-exact RBF over the nearest particles, `rbf_interpolate`); only its target moved. A particle-carried material step is reproduced at the integration points with less than half @@ -350,6 +384,41 @@ the interface (`tests/test_0067_integration_point_proxy.py`). has no gradient, so a derivative of the swarm variable's symbol is refused. Vector and tensor swarm variables get a multi-component proxy. +### `proxy_sampling`: what each point reads + +`proxy_location` is where; `proxy_sampling` is what. + +| | `"reconstruct"` (default) | `"share"` | +|---|---|---| +| gathers from | the `nnn` nearest particles, by distance | the particles whose nearest integration point *in their own cell* is this one | +| respects cell walls | no | yes | +| linear fields | exact | small averaging error | +| bounded by the particle values | no (overshoots a jump) | yes, it is a mean of them | +| particles used | the stencil's | all of them, each exactly once | + +`"share"` is the cell-restricted Voronoi share +(`underworld3/utilities/particle_share.py`): `share_assignment` maps each +particle to one flat index in the cell-major `(ncells, Nq)` layout, +`share_average` reduces by `np.bincount`. A rule point whose share is empty +falls back to the nearest particle anywhere on the rank, and the count of +those is left on `var._share_empty` — persistently non-zero means the swarm is +too thin for the rule, and `Swarm.repopulate` is the fix. + +The assignment needs each particle's owning cell. UW3 swarms are +`DMSWARM_BASIC`, so PETSc holds no cell id and the locator has to run; +`Swarm._owning_cells()` caches the result and drops it wherever `_kdtree` is +dropped, so the share proxies, the population census and anything else +cell-local pay for one location per step between them. On 32 912 particles / +242 cells the location is 18.8 ms and the share itself 3.2 ms, against 8.1 ms +for the `"reconstruct"` path (whose cached operator is geometry-only, so it is +rebuilt every time the particles move). + +There is no `"nearest"` here. Sampling one particle's value whole is the +material mapping, and materials go through `MaterialSwarm` / its +`IndexSwarmVariable` (`proxy_sampling="nearest"` there, or `"share"` for +fractional masks); asking for it on a plain `SwarmVariable` raises and names +the alternative. + `Lagrangian_Swarm(..., proxy_location="integration_points")` applies the same to the fully Lagrangian history: the slots carried on the particles are reconstructed at the integration points and the weak form reads them diff --git a/docs/examples/utilities/intermediate/Ex_Swarm_Material_Index.py b/docs/examples/utilities/intermediate/Ex_Swarm_Material_Index.py new file mode 100644 index 000000000..49dd2fb3c --- /dev/null +++ b/docs/examples/utilities/intermediate/Ex_Swarm_Material_Index.py @@ -0,0 +1,149 @@ +# --- +# jupyter: +# jupytext: +# formats: py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% [markdown] +""" +# Materials on a swarm + +**PHYSICS:** fluid_mechanics +**DIFFICULTY:** intermediate +**PURPOSE:** demonstration + +## Description + +Two viscosity layers, 1 and 1000, carried by particles and driven from the +top. The model names its materials and says where they are; it never writes +a level set, a mask, or a blend. + +With the interface on mesh edges the exact velocity is piecewise linear and +lies in the P2 velocity space, so the only error in the solve is how the +material is represented. Read at the integration points (the default), where +each point takes the material of its nearest particle, the problem solves to +2e-7. Run with `-uw_proxy_location nodes` to watch the nodal level set smear +the interface across a cell and leave an L2 error of 8e-2. +""" + +# %% +import numpy as np +import sympy + +import underworld3 as uw + +params = uw.Params( + uw_proxy_location="integration_points", # or "nodes", or "cells" + uw_proxy_sampling="nearest", # or "share" + uw_cell_size=0.1, + uw_eta_top=1000.0, + uw_interface=0.5, + uw_fill_param=3, +) + +mesh = uw.meshing.UnstructuredSimplexBox( + cellSize=params.uw_cell_size, qdegree=2, regular=True +) +v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1) + +# %% [markdown] +""" +## The materials + +`MaterialSwarm` is a `Swarm` that carries materials. Declare them with the +property names the constitutive model knows, then say where each one is: a +symbolic condition on the mesh coordinates, a boolean array over the +particles, or a callable. Every `add` must come before the first read, which +is what allocates the particles and the level sets. +""" + +# %% +materials = uw.swarm.MaterialSwarm( + mesh, + fill_param=params.uw_fill_param, + proxy_location=params.uw_proxy_location, + proxy_sampling=( + params.uw_proxy_sampling + if params.uw_proxy_location == "integration_points" + else None + ), +) +materials.add("lower", shear_viscosity_0=1.0, density=3300) +materials.add("upper", shear_viscosity_0=params.uw_eta_top, density=3400) + +materials["upper"] = mesh.X[1] > params.uw_interface + +# %% [markdown] +""" +## The solve + +`stokes.materials = materials` sets every parameter the constitutive model +recognises — here `shear_viscosity_0`. `density` is not a viscous-model +parameter, so it stays a blended symbol for the model script to use; this +problem is driven by the boundary, so it goes unused. +""" + +# %% +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.materials = materials + +stokes.add_dirichlet_bc((1.0, 0.0), "Top") +stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") +stokes.tolerance = 1e-8 +stokes.solve() + +# %% [markdown] +""" +## Against the exact layered Couette profile + +The error is a global integral rather than a nodal norm, so the number is the +same however the mesh is partitioned. +""" + +# %% +h, eta_top = params.uw_interface, params.uw_eta_top +gradient = 1.0 / (h + (1.0 - h) / eta_top) +y = mesh.X[1] +exact = sympy.Piecewise( + (gradient * y, y < h), + (gradient * h + gradient / eta_top * (y - h), True), +) + +error = uw.maths.Integral(mesh, (v.sym[0] - exact) ** 2).evaluate() ** 0.5 +viscosity = uw.maths.Integral(mesh, materials.shear_viscosity_0).evaluate() +exact_viscosity = 1.0 * h + eta_top * (1.0 - h) + +uw.pprint( + f"proxy_location={params.uw_proxy_location} " + f"sampling={params.uw_proxy_sampling} fill={params.uw_fill_param}: " + f"assembled int(eta) {viscosity:.4f} (exact {exact_viscosity:.4f}) | " + f"velocity L2 {error:.3e}" +) + +# %% [markdown] +""" +## What the weak form sees + +The upper-material mask along a line crossing the interface. At the +integration points it is a step in the right place; the nodal level set ramps +across a whole cell. +""" + +# %% +line = np.column_stack( + [np.full(201, 0.5), np.linspace(max(0.0, h - 0.25), min(1.0, h + 0.25), 201)] +) +upper = np.asarray(uw.function.evaluate(materials["upper"].mask, line)).reshape(-1) +uw.pprint(f"mask along the line: {upper.min():+.4f} .. {upper.max():+.4f}") diff --git a/docs/examples/utilities/intermediate/Ex_Swarm_Population_Control.py b/docs/examples/utilities/intermediate/Ex_Swarm_Population_Control.py new file mode 100644 index 000000000..2cdbeba41 --- /dev/null +++ b/docs/examples/utilities/intermediate/Ex_Swarm_Population_Control.py @@ -0,0 +1,141 @@ +# --- +# jupyter: +# jupytext: +# formats: py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% [markdown] +""" +# Population control in an extending box + +**PHYSICS:** fluid_mechanics +**DIFFICULTY:** intermediate +**PURPOSE:** demonstration + +## Description + +Pure-shear extension $\\mathbf{v} = (x, -y)$ on a fixed mesh. The side walls +are outflow and the top and bottom are inflow, and nothing arrives through an +inflow boundary unless you put it there — so the cells along the top and +bottom starve. + +What that costs depends on how the material is read. This example measures it +on the area of a marker layer, a global integral of the layer's own material +mask, so the number is the same however the mesh is partitioned. + +Incompressible pure shear thins the layer by $e^{-t}$, so after $t = 2$ the +area should be $e^{-2} = 13.5\\%$ of where it started. Anything else is the +material representation failing, not physics. + +Run with `-uw_population_control 0` to switch the refilling off, and with +`-uw_proxy_location integration_points` to see how much the choice of +mapping matters. +""" + +# %% +import math + +import sympy + +import underworld3 as uw + +params = uw.Params( + uw_population_control=1, + uw_proxy_location="cells", # or "integration_points" + uw_cell_size=0.08, + uw_steps=40, + uw_dt=0.05, + uw_fill_param=3, + uw_min_per_cell=6, +) + +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5), + cellSize=params.uw_cell_size, qdegree=2, +) +x, y = mesh.X +velocity = sympy.Matrix([[x, -y]]) # incompressible pure shear + +# A true outflow: particles that leave the box are deleted rather than +# clamped onto the wall. +mesh.return_coords_to_bounds = None + +# %% [markdown] +""" +## A marker layer on a material swarm + +`proxy_location="cells"` fits a polynomial to the particles **of each cell**, +so a cell that runs out of particles has nothing to fit. That is the mapping +population control exists for. +""" + +# %% +materials = uw.swarm.MaterialSwarm( + mesh, fill_param=params.uw_fill_param, + proxy_location=params.uw_proxy_location, +) +materials.add("matrix", density=1.0) +materials.add("layer", density=1.0) +materials["layer"] = sympy.Abs(y) < 0.2 + +if params.uw_population_control: + # Refill starved cells at the end of every advection. A new particle takes + # the material of its nearest neighbour: repopulate() does that for any + # INTEGER variable without being asked, because the average of two + # material labels is not a label. + materials.population_control = dict(min_per_cell=params.uw_min_per_cell) + +layer_area = uw.maths.Integral(mesh, materials["layer"].mask) + +# %% [markdown] +""" +## Extend the box +""" + +# %% +initial_area = layer_area.evaluate() + +for _ in range(params.uw_steps): + materials.advection(velocity, params.uw_dt, order=2) + +elapsed = params.uw_steps * params.uw_dt +final_area = layer_area.evaluate() +expected = initial_area * math.exp(-elapsed) + +uw.pprint( + f"proxy_location={params.uw_proxy_location} " + f"population_control={bool(params.uw_population_control)}: " + f"layer area {initial_area:.4f} -> {final_area:.4f} at t={elapsed:g} " + f"(exact thinning gives {expected:.4f}, " + f"so this is {final_area / expected:.2f}x the right answer)" +) + +# %% [markdown] +""" +## What the numbers say + +Measured on this rig, `cellSize=0.08`, 40 steps: + +| `proxy_location` | population control | layer area | vs exact | +|---|---|---|---| +| `"cells"` | on | 0.1125 | 1.03x | +| `"cells"` | off | 0.5505 | **5.03x** | +| `"integration_points"` | on | 0.1169 | 1.06x | +| `"integration_points"` | off | 0.1175 | 1.07x | + +Starvation wrecks the per-cell fit, which needs particles *in that cell*. The +nearest-particle mapping degrades gracefully by comparison — an empty cell +still finds a plausible particle nearby — so with `"integration_points"` the +field-level answer barely moves and the cost of losing particles shows up +only in the swarm's own bookkeeping. + +Population control is cheap and always safe; this is where it is *necessary*. +""" diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py index 4ed69e815..1469569c0 100644 --- a/src/underworld3/__init__.py +++ b/src/underworld3/__init__.py @@ -222,7 +222,12 @@ def view(): create_thermal_convection_model, ) from .parameters import ParameterRegistry, ParameterType -from .materials import MaterialRegistry, MaterialProperty +from .materials import ( + MaterialRegistry, + MaterialProperty, + MaterialDefinition, + MaterialRegions, +) from .constitutive_models import MultiMaterialConstitutiveModel # uw.quantity is THE quantity factory (returns UWQuantity, exposed alongside # for isinstance checks); uw.create_quantity is deprecated (see units.py). diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f7b1b1951..c7d53b68c 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -160,6 +160,7 @@ class SolverBaseClass(uw_object): self._order = 0 self._constitutive_model = None + self._materials = None self._rebuild_after_mesh_update = self._build self.name = "Solver_{}_".format(self.instance_number) @@ -2076,6 +2077,8 @@ class SolverBaseClass(uw_object): ): self._check_expression_meshes() + if self._materials is not None: + self._materials.check() if self.is_setup: return @@ -2727,6 +2730,50 @@ class SolverBaseClass(uw_object): return + @property + def materials(self): + """The materials this solver's coefficients come from. + + Assigning a :class:`~underworld3.swarm.MaterialSwarm` sets every + constitutive-model parameter the materials declare *and* the model + recognises, by name — so a model script names its materials and their + properties, and never writes a level set or a mask:: + + materials = uw.swarm.MaterialSwarm(mesh, fill_param=3) + materials.add("mantle", shear_viscosity_0=1.0) + materials.add("slab", shear_viscosity_0=1.0e3, density=3400) + materials["slab"] = mesh.X[1] > 0.53 + + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.materials = materials # sets shear_viscosity_0 + + A declared property the model does not recognise (``density`` here) is + not pushed anywhere; it is available as a blended symbol, + ``materials.density``, for the model script to use where it belongs. + One that is neither recognised nor read is reported at solve time, + because a misspelled viscosity is silently the default one. + + Properties may be changed, and materials repainted, after assignment: + the blend is symbolic and the push repeats on every change. + """ + return self._materials + + @materials.setter + def materials(self, material_swarm): + if material_swarm is None: + previous = self._materials + self._materials = None + if previous is not None: + previous._detach(self) # or it keeps pushing to this solver + return + if not hasattr(material_swarm, "_attach"): + raise TypeError( + "solver.materials expects a MaterialSwarm (uw.swarm.MaterialSwarm), " + f"not {type(material_swarm).__name__}" + ) + self._materials = material_swarm + material_swarm._attach(self) + @property def constitutive_model(self): """ @@ -2785,6 +2832,11 @@ class SolverBaseClass(uw_object): if self._constitutive_model.requires_stress_history and self.Unknowns.DFDt is None: self._create_stress_history_ddt(order=self._constitutive_model.order) + # Materials assigned before the constitutive model still have to + # reach it: the push is a no-op while there is no model to push to. + if getattr(self, "_materials", None) is not None: + self._materials._push_to(self) + # May not work due to flux being incomplete if self.Unknowns.DFDt is not None: self.Unknowns.DFDt.psi_fn = self._constitutive_model.flux.T diff --git a/src/underworld3/discretisation/enhanced_variables.py b/src/underworld3/discretisation/enhanced_variables.py index 7c3b8e3cc..9a1388fa3 100644 --- a/src/underworld3/discretisation/enhanced_variables.py +++ b/src/underworld3/discretisation/enhanced_variables.py @@ -957,6 +957,14 @@ class IntegrationPointVariable(EnhancedMeshVariable): >>> eta_q.cell_data[...] = 1.0 # (ncells, Nq, 1) >>> eta_q.coords # the physical integration points >>> stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_q.sym + + Gradients: this variable has none of its own (it stores a value at each + integration point and nothing between them, so its tabulated derivative is + identically zero). A derivative of its symbol in a WEAK FORM is refused + rather than answered with that zero; use a discontinuous mesh variable, or + ``SwarmVariable(proxy_location="cells")``, when a solve needs the gradient. + ``uw.function.evaluate`` of the same derivative does answer, by fitting the + values per cell first: a recovered gradient, for inspection. """ _base_variable_class = _BaseIntegrationPointVariable diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 3fb4e9165..96540e25f 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -787,6 +787,88 @@ def _project_to_work_variable(expr, mesh, smoothing=1e-6): return work_var +def _integration_point_sources_to_cell_fit(expr, mesh, derivfns): + """Replace integration-point variables that appear under a derivative by a + per-cell least-squares fit of their own values, which does have a gradient. + + An integration-point variable holds a value at each integration point and + nothing between them: its tabulated gradient is identically zero, so both + the residual and an L2 projection of the derivative would return a silent + zero. The values themselves are good data, though, and fitting them cell by + cell (the reconstruction of + :class:`~underworld3.utilities.cell_polynomial_projection.CellPolynomialProjector`, + the same one behind ``SwarmVariable(proxy_location="cells")``) gives a + polynomial per cell that differentiates directly, with no projection solve. + + **Two paths, deliberately different.** ``uw.function.evaluate`` is a query: + it answers, by way of this recovery, and the answer converges (measured on + a quadratic particle field: 2.4e-3, 6.1e-4, 2.6e-4 as the cell size halves + from 1/5 to 1/20). Assembly of a weak form still REFUSES, because there a + hidden reconstruction would be a per-assembly cost and would quietly decide + a discretisation the user should choose. If a solve needs the gradient, + build the variable with ``proxy_location="cells"``, whose level sets are + already polynomials: that is both cheaper and sharper (2.4e-7 on the same + field, exact for a quadratic, because no further projection follows). + + Returns ``(expr, derivfns)`` with the substitutions applied. + """ + from underworld3.utilities.cell_polynomial_projection import CellPolynomialProjector + + ip_sources = [v for v in derivfns + if getattr(v, "is_integration_point", False)] + if not ip_sources: + return expr, derivfns + + subs = {} + for source_var in ip_sources: + # Degree: enough to carry a gradient, and no more than the rule can + # support (the fit falls back per cell when a cell is short of points). + degree = max(1, min(2, mesh.qdegree)) + # Cached on the mesh, like the other evaluate work variables: adding a + # field rebuilds the mesh DM, so one per (source, degree) and no more. + cache_key = f"_eval_ipgrad_{source_var.clean_name}_{degree}" + work = getattr(mesh, cache_key, None) + if work is None: + work = uw.discretisation.MeshVariable( + cache_key, mesh, source_var.shape, source_var.vtype, + degree=degree, continuous=False, + varsymbol=rf"{{ \widehat{{{source_var.symbol}}} }}", + ) + setattr(mesh, cache_key, work) + setattr(mesh, cache_key + "_projector", CellPolynomialProjector(work)) + projector = getattr(mesh, cache_key + "_projector") + if projector.mesh_version != mesh._mesh_version: + projector = CellPolynomialProjector(work) + setattr(mesh, cache_key + "_projector", projector) + # The rule has exactly as many points per cell as a fit of this degree + # has coefficients (that is what "the rule is unisolvent for P_k" + # means), so allow an exactly-determined fit: nmin = Nb rather than the + # default Nb + 2, which would send every cell to the linear patch and + # leave the gradient first order. The projector's own conditioning + # guard still catches a degenerate cell. + fitted = projector.fit( + np.asarray(source_var.coords_nd), + np.asarray(source_var.data).reshape(-1, source_var.num_components), + nmin=projector.Nb, + ) + work.data[...] = fitted.reshape(work.data.shape) + + src_flat, work_flat = source_var.sym_1d, work.sym_1d + for k in range(source_var.num_components): + subs[src_flat[k]] = work_flat[k] + for deriv_expr, diffindex in derivfns[source_var]: + subs[deriv_expr] = work_flat[deriv_expr.component].diff(mesh.X[diffindex]) + + expr = expr.subs(subs) if hasattr(expr, "subs") else expr + derivfns = {v: d for v, d in derivfns.items() if v not in ip_sources} + # the substituted derivatives are now derivatives of ordinary mesh + # variables; let the caller re-extract them + _, _, new_derivs = uw.function.fn_mesh_vars_in_expression(expr) + for v, d in new_derivs.items(): + derivfns.setdefault(v, d) + return expr, derivfns + + def _clement_to_work_variable(expr, mesh, derivfns): """ Evaluate expression at nodes using Clement gradient recovery (no solve). @@ -979,6 +1061,12 @@ def evaluate_nd( expr, # Two modes: # - Quick (rbf=True, force_l2=False): Clement gradient at nodes, no solve # - Accurate (force_l2=True or rbf=False): L2 projection, requires solve + if derivfns and mesh is not None: + # An integration-point source has no gradient of its own: fit its + # values per cell first, then the ordinary derivative machinery below + # differentiates a polynomial (see the helper). + expr, derivfns = _integration_point_sources_to_cell_fit(expr, mesh, derivfns) + if derivfns and mesh is not None: if evalf: raise RuntimeError( diff --git a/src/underworld3/materials.py b/src/underworld3/materials.py index c2f0ca817..2ecb4a974 100644 --- a/src/underworld3/materials.py +++ b/src/underworld3/materials.py @@ -1,394 +1,1091 @@ -""" -Material Management System for Underworld3 Models +r"""Materials: what they are, and where they are. + +Two questions, kept apart on purpose. + +**What a material is** is a :class:`MaterialRegistry` entry — a name and a table +of properties, where a property may be a number, a dimensional quantity, or a +*law* (``eta_0 * sympy.exp(-T.sym[0])``). Definitions carry description and +reference metadata, export and import as plain dictionaries, and say nothing +about geometry, so one registry can serve several models. + +**Where a material is** is a *distribution*. Underworld3 has two: + +- :class:`~underworld3.swarm.MaterialSwarm` — carried by particles, so it + advects with the flow. The one to use when the material moves. +- :class:`MaterialRegions` — tied to the mesh, from gmsh physical groups or + from a geometric condition. Exact, needs no particles and no population + control. The one to use when the material does not move. + +Both present the same face to a solver:: + + stokes.materials = materials + +which sets every constitutive-model parameter the materials declare and the +model recognises, by name. Both build that from the same partition of unity — +one level set per material — because a property that is a law cannot be stored +as a value and has to be combined symbolically. -This module provides structured material management with property definitions, -region assignments, and automatic propagation to constitutive models and solvers. +Example +------- +>>> rocks = uw.MaterialRegistry() +>>> rocks.add("mantle", shear_viscosity_0=1.0, density=3300) +>>> rocks.add("slab", shear_viscosity_0=1.0e3, density=3400) +>>> +>>> materials = uw.swarm.MaterialSwarm(mesh, registry=rocks, fill_param=3) +>>> materials["slab"] = mesh.X[1] > 0.53 +>>> stokes.materials = materials + +or, for a mesh that carries the geometry itself, + +>>> materials = uw.MaterialRegions(mesh, registry=rocks) +>>> materials["slab"] = "Slab" # a gmsh physical group +>>> stokes.materials = materials """ -import weakref -from typing import Any, Dict, List, Optional, Union, Callable -from dataclasses import dataclass, field +import warnings from enum import Enum +from typing import Any, Dict, List, Optional + import numpy as np +import sympy + +import underworld3 as uw + +__all__ = [ + "MaterialProperty", + "MaterialDefinition", + "MaterialRegistry", + "MaterialDistribution", + "BoundMaterial", + "MaterialRegions", + "create_standard_mantle_material", + "create_standard_crust_material", + "create_high_viscosity_material", +] + +_MIXING_RULES = ("arithmetic", "harmonic") class MaterialProperty(Enum): - """Standard material properties for geodynamic models""" + """Names for the properties a geodynamic material usually carries. + + The *value* of each member is the property name, and it is the name a + constitutive model knows it by — ``MaterialProperty.VISCOSITY`` is + ``"viscosity"``, which is the established alias for + ``shear_viscosity_0``. Using the enum is optional; a plain keyword is + equivalent and is what most models do. + """ - # Mechanical properties + # Mechanical VISCOSITY = "viscosity" DENSITY = "density" YIELD_STRESS = "yield_stress" COHESION = "cohesion" FRICTION_ANGLE = "friction_angle" - # Thermal properties + # Thermal THERMAL_CONDUCTIVITY = "thermal_conductivity" THERMAL_DIFFUSIVITY = "thermal_diffusivity" HEAT_CAPACITY = "heat_capacity" THERMAL_EXPANSION = "thermal_expansion" - # Elastic properties + # Elastic YOUNGS_MODULUS = "youngs_modulus" POISSONS_RATIO = "poissons_ratio" SHEAR_MODULUS = "shear_modulus" - # Flow properties + # Flow PERMEABILITY = "permeability" POROSITY = "porosity" - # Custom properties - CUSTOM = "custom" +def _property_name(prop): + return prop.value if isinstance(prop, MaterialProperty) else str(prop) -@dataclass -class MaterialDefinition: + +def _as_symbolic(value, name, material): + """Coerce a declared property value into something sympy can carry. + + A dimensional quantity is non-dimensionalised HERE, so this is called when + the value is *read*, not when it is declared: the scaling is a property of + the model, and a registry is meant to be writable before the model exists. + Resolving at declaration time made the same declaration mean different + numbers depending on whether the reference quantities had been set yet. """ - Definition of a material with its properties and metadata. + if isinstance(value, sympy.Basic): + return value + try: + return sympy.sympify(value) + except (TypeError, ValueError, AttributeError): + # Not a plain number or expression. Fall through to the quantity path + # below, which raises with the offending value if that fails too. + pass + try: # a dimensional quantity + return sympy.sympify(uw.scaling.non_dimensionalise(value)) + except Exception as exc: # noqa: BLE001 - name the value + raise TypeError( + f"material {material!r}: property {name!r} = {value!r} is neither a " + "number, a sympy expression, nor a quantity that can be " + "non-dimensionalised" + ) from exc + + +class MaterialDefinition: + """What a material is: a name, a property table, and provenance. + + Created by :meth:`MaterialRegistry.add`. Says nothing about where the + material is — that is a distribution's job. - Attributes: - ----------- + Attributes + ---------- name : str - Material name (e.g., 'mantle', 'crust', 'plume') + index : int + Position in the registry, and therefore which level set is this + material's in any distribution built from it. properties : dict - Dictionary of property_name -> value mappings - description : str - Human-readable description - reference : str - Literature reference or source - temperature_dependent : dict - Temperature-dependent property functions - pressure_dependent : dict - Pressure-dependent property functions - constitutive_model : object - Associated constitutive model instance + Property name -> value. Values are sympy objects: a number, or a law. + description, reference : str + Free text, carried through export/import. """ - name: str - properties: Dict[str, Any] = field(default_factory=dict) - description: str = "" - reference: str = "" - temperature_dependent: Dict[str, Callable] = field(default_factory=dict) - pressure_dependent: Dict[str, Callable] = field(default_factory=dict) - constitutive_model: Optional[Any] = None - - def set_property(self, prop: Union[MaterialProperty, str], value: Any): - """Set a material property value""" - prop_name = prop.value if isinstance(prop, MaterialProperty) else prop - self.properties[prop_name] = value - - def get_property(self, prop: Union[MaterialProperty, str], default=None): - """Get a material property value""" - prop_name = prop.value if isinstance(prop, MaterialProperty) else prop - return self.properties.get(prop_name, default) - - def has_property(self, prop: Union[MaterialProperty, str]) -> bool: - """Check if material has a specific property""" - prop_name = prop.value if isinstance(prop, MaterialProperty) else prop - return prop_name in self.properties - - def evaluate_property( - self, prop: Union[MaterialProperty, str], temperature=None, pressure=None - ): + def __init__(self, registry, name, index, properties=None, + description="", reference=""): + self._registry = registry + self.name = name + self.index = index + self.description = description + self.reference = reference + self.properties = {} + if properties: + self.set(**properties) + + def set(self, **properties): + """Set or replace properties, and re-push them to attached solvers. + + The value is stored AS DECLARED. A dimensional quantity is + non-dimensionalised when it is read (see :meth:`resolved`), because + the scaling belongs to the model and a registry may be written before + there is one. """ - Evaluate a material property, accounting for temperature/pressure dependence. - - Parameters: - ----------- - prop : MaterialProperty or str - Property to evaluate - temperature : float or array, optional - Temperature for evaluation - pressure : float or array, optional - Pressure for evaluation - - Returns: - -------- - Property value (scalar or array) - """ - prop_name = prop.value if isinstance(prop, MaterialProperty) else prop - - # Get base value - base_value = self.get_property(prop_name) - if base_value is None: - raise ValueError(f"Material '{self.name}' does not have property '{prop_name}'") - - # Apply temperature dependence - if temperature is not None and prop_name in self.temperature_dependent: - temp_func = self.temperature_dependent[prop_name] - base_value = temp_func(base_value, temperature) - - # Apply pressure dependence - if pressure is not None and prop_name in self.pressure_dependent: - pressure_func = self.pressure_dependent[prop_name] - base_value = pressure_func(base_value, pressure) + for key, value in properties.items(): + key = _property_name(key) + _as_symbolic(value, key, self.name) # validate now, resolve later + self.properties[key] = value + self._registry._changed() + return self + + def resolved(self, prop): + """The property as sympy, non-dimensionalised against the model's + current scaling.""" + name = _property_name(prop) + return _as_symbolic(self.properties[name], name, self.name) + + # -- the pre-2026 spelling, kept working ------------------------------ + def set_property(self, prop, value): + """Set one property. ``set(**{name: value})`` is the shorter form.""" + return self.set(**{_property_name(prop): value}) + + def get_property(self, prop, default=None): + return self.properties.get(_property_name(prop), default) + + def has_property(self, prop) -> bool: + return _property_name(prop) in self.properties - return base_value + def __repr__(self): + props = ", ".join(f"{k}={v}" for k, v in self.properties.items()) + return f"" class MaterialRegistry: - """ - Central registry for material definitions and region assignments. + """A set of material definitions, independent of any mesh or swarm. - Features: - --------- - - Material property database with validation - - Region-based material assignments - - Temperature/pressure dependent properties - - Integration with constitutive models - - Automatic property propagation to solvers + A registry is the *what*: it can be built before there is a mesh, shared + between models, exported to a dictionary and read back. Attach it to a + distribution — a :class:`~underworld3.swarm.MaterialSwarm` or a + :class:`MaterialRegions` — to say where each material is. - Example: + Examples -------- - >>> registry = MaterialRegistry() - >>> mantle = registry.create_material('mantle') - >>> mantle.set_property('viscosity', 1e21) - >>> mantle.set_property('density', 3300) - >>> registry.assign_to_region('mantle', region_id=1) + >>> rocks = uw.MaterialRegistry() + >>> rocks.add("mantle", shear_viscosity_0=1.0, density=3300) + >>> rocks.add("slab", shear_viscosity_0=1.0e3, density=3400) + >>> rocks.list_materials() + ['mantle', 'slab'] """ - def __init__(self): + def __init__(self, materials=None): self._materials: Dict[str, MaterialDefinition] = {} - self._region_assignments: Dict[int, str] = {} # region_id -> material_name - self._callbacks: List[Callable] = [] # Material change callbacks + self._order: List[str] = [] + self._callbacks = [] + self._built = [] # distributions that have built their level sets self._version = 0 + if materials: + for name, properties in dict(materials).items(): + self.add(name, **properties) - def create_material( - self, name: str, description: str = "", reference: str = "" - ) -> MaterialDefinition: - """ - Create a new material definition. + # -- declaring --------------------------------------------------------- - Parameters: - ----------- + def add(self, name, description="", reference="", **properties): + """Declare a material. + + Parameters + ---------- name : str - Material name - description : str - Human-readable description - reference : str - Literature reference - - Returns: - -------- + description, reference : str, optional + Free text carried through export/import. + **properties + Property values, by the name a constitutive model knows them by + (``shear_viscosity_0``, ``density``, ...). A value may be a number, + a quantity, or a symbolic law. + + Returns + ------- MaterialDefinition - New material instance """ + self._refuse_if_built(f"declare material {name!r}") if name in self._materials: - raise ValueError(f"Material '{name}' already exists") - - material = MaterialDefinition(name=name, description=description, reference=reference) - + raise ValueError(f"material {name!r} has already been declared") + material = MaterialDefinition( + self, name, len(self._order), properties, description, reference + ) self._materials[name] = material - self._version += 1 - + self._order.append(name) + self._changed() return material - def get_material(self, name: str) -> Optional[MaterialDefinition]: - """Get a material by name""" + def create_material(self, name, description="", reference=""): + """Declare a material with no properties yet (the pre-2026 spelling).""" + return self.add(name, description=description, reference=reference) + + def get_material(self, name) -> Optional[MaterialDefinition]: return self._materials.get(name) def list_materials(self) -> List[str]: - """List all material names""" - return list(self._materials.keys()) + """Material names, in declaration order — which is level-set order.""" + return list(self._order) - def delete_material(self, name: str): - """Delete a material definition""" - if name in self._materials: - del self._materials[name] - # Remove any region assignments - self._region_assignments = { - region_id: mat_name - for region_id, mat_name in self._region_assignments.items() - if mat_name != name - } - self._version += 1 - - def assign_to_region(self, material_name: str, region_id: int): + def delete_material(self, name): + """Remove a material, and reindex the rest. + + Only safe before a distribution has been built from the registry: the + index of a material IS which level set is its, so removing one + renumbers the others. """ - Assign a material to a mesh region. - - Parameters: - ----------- - material_name : str - Name of material to assign - region_id : int - Mesh region identifier + self._refuse_if_built(f"delete material {name!r}") + if name not in self._materials: + raise KeyError( + f"no material {name!r}; declared: {self.list_materials()}" + ) + del self._materials[name] + self._order.remove(name) + for i, key in enumerate(self._order): + self._materials[key].index = i + self._changed() + + # -- reading ----------------------------------------------------------- + + @property + def materials(self): + """The definitions, in declaration order.""" + return tuple(self._materials[name] for name in self._order) + + def declared_properties(self): + """Every property name declared by any material.""" + names = set() + for material in self.materials: + names.update(material.properties) + return names + + def __len__(self): + return len(self._order) + + def __iter__(self): + return iter(self.materials) + + def __contains__(self, name): + return name in self._materials + + def __getitem__(self, name): + try: + return self._materials[name] + except KeyError: + raise KeyError( + f"no material {name!r}; declared: {self.list_materials()}" + ) from None + + def __repr__(self): + return f"MaterialRegistry({self.list_materials()})" + + # -- change notification ---------------------------------------------- + + def _refuse_if_built(self, what): + """A material's index IS which level set is its, so the set of + materials cannot change once a distribution has allocated them. + + Without this, ``registry.add`` after a build produced a material with + no level set (an ``IndexError`` from ``blend``, swallowed into a + warning by the change callback, leaving the solver holding a stale + blend), and ``delete_material`` silently re-pointed every later + material at its neighbour's level set. """ - if material_name not in self._materials: - raise ValueError(f"Material '{material_name}' does not exist") - - self._region_assignments[region_id] = material_name - self._notify_callbacks("region_assignment", region_id, material_name) - - def get_region_material(self, region_id: int) -> Optional[str]: - """Get the material assigned to a region""" - return self._region_assignments.get(region_id) - - def get_material_regions(self, material_name: str) -> List[int]: - """Get all regions assigned to a material""" - return [ - region_id - for region_id, mat_name in self._region_assignments.items() - if mat_name == material_name - ] + if self._built: + names = ", ".join(sorted(type(d).__name__ for d in self._built)) + raise RuntimeError( + f"cannot {what}: this registry is already in use by {names}, " + "whose level sets are allocated one per material. Declare " + "every material before the first read." + ) + + def add_callback(self, callback): + """Register ``callback()`` to run whenever a definition changes.""" + self._callbacks.append(callback) + + def _changed(self): + self._version += 1 + for callback in list(self._callbacks): + try: + callback() + except Exception as exc: # noqa: BLE001 + warnings.warn(f"material registry callback failed: {exc}") + + # -- serialisation ----------------------------------------------------- + + def export_config(self) -> Dict[str, Any]: + """Export as plain data. Property values are stringified sympy.""" + return { + "materials": { + m.name: { + "properties": {k: str(v) for k, v in m.properties.items()}, + "description": m.description, + "reference": m.reference, + } + for m in self.materials + }, + "version": self._version, + } + + def import_config(self, config: Dict[str, Any]): + """Read back an :meth:`export_config` dictionary.""" + for name, entry in config.get("materials", {}).items(): + material = self.add( + name, + description=entry.get("description", ""), + reference=entry.get("reference", ""), + ) + material.set(**{ + key: sympy.sympify(value) + for key, value in entry.get("properties", {}).items() + }) + return self + + +class BoundMaterial: + """A material seen through one distribution. + + What :meth:`MaterialDistribution.__getitem__` returns: the definition plus + the things that only mean something once you know *where* the material is — + its level set, and how to paint it. + """ + + def __init__(self, distribution, definition): + self._distribution = distribution + self._definition = definition + + @property + def name(self): + return self._definition.name - def evaluate_property_field( - self, - prop: Union[MaterialProperty, str], - region_field: np.ndarray, - temperature: Optional[np.ndarray] = None, - pressure: Optional[np.ndarray] = None, - ) -> np.ndarray: + @property + def index(self): + return self._definition.index + + @property + def properties(self): + return self._definition.properties + + @property + def mask(self): + r"""This material's level set, :math:`\phi_i`, as a symbol. + + 1 where the material is, 0 where it is not. Rarely needed — a property + blend is what a model usually wants — but it is the right thing for a + per-material diagnostic, e.g. ``uw.maths.Integral(mesh, slab.mask)`` + for the area the material occupies. """ - Evaluate a material property over a field of region IDs. - - Parameters: - ----------- - prop : MaterialProperty or str - Property to evaluate - region_field : array - Array of region IDs - temperature : array, optional - Temperature field for evaluation - pressure : array, optional - Pressure field for evaluation - - Returns: - -------- - array - Property values corresponding to each region + return self._distribution._mask_of(self._definition) + + def occupies(self, region): + """Put this material wherever ``region`` is true. + + What ``region`` may be depends on the distribution — see + :meth:`MaterialDistribution.__setitem__`. """ - prop_name = prop.value if isinstance(prop, MaterialProperty) else prop + self._distribution[self.name] = region + return self - # Initialize output array - result = np.zeros_like(region_field, dtype=float) + def set(self, **properties): + """Change properties, and re-push them to attached solvers.""" + self._definition.set(**properties) + return self - # Evaluate property for each unique region - unique_regions = np.unique(region_field) + def __repr__(self): + return f"<{type(self._distribution).__name__} {self._definition!r}>" - for region_id in unique_regions: - material_name = self.get_region_material(region_id) - if material_name is None: - raise ValueError(f"No material assigned to region {region_id}") - material = self.get_material(material_name) - if material is None: - raise ValueError(f"Material '{material_name}' not found") +class MaterialDistribution: + """Shared behaviour of everything that says *where* materials are. - # Get mask for this region - mask = region_field == region_id + A distribution owns a :class:`MaterialRegistry` and supplies one level set + per material. Everything else — blending a property into a symbol, pushing + it to a solver, the mixing rule, the unused-property check — is the same + whether the materials ride on particles or on mesh regions, and lives here. - # Extract temperature/pressure for this region if provided - region_temp = temperature[mask] if temperature is not None else None - region_pressure = pressure[mask] if pressure is not None else None + Subclasses provide ``_ensure_built()`` and ``_level_sets()``. + """ - # Evaluate property for this region - prop_value = material.evaluate_property(prop_name, region_temp, region_pressure) + #: How many distributions have taken the default name, so that two on one + #: mesh do not collide over their level-set variable names. + _default_name_count = 0 - # Assign to result - result[mask] = prop_value + def _check_level_sets_are_new(self, mesh, before, count): + """Refuse a name whose level sets already existed on this mesh. - return result + Two distributions sharing a name silently shared their level sets: + creating the second one's variables printed "Variable ... already + exists - Skipping" to stdout and handed back the FIRST distribution's + variables, so painting the second changed the first one's answers with + no error anywhere. The default names are counter-suffixed to avoid + this; an explicit ``name=`` was unguarded. - def add_callback(self, callback: Callable): + Checked by counting what the mesh actually gained rather than by + predicting the sanitised variable names, which are not the strings + passed in — ``"M^{[0]}"`` becomes ``"M0"``. + """ + gained = len(mesh.vars) - before + if gained != count: + raise ValueError( + f"name={self._distribution_name!r} is already in use on this " + f"mesh: {count} level sets were requested and {gained} were " + "created, so this distribution would share another's storage. " + "Give it a different name." + ) + + def _init_distribution(self, registry=None, name=None): + self._registry = registry if registry is not None else MaterialRegistry() + if name is None: + n = MaterialDistribution._default_name_count + MaterialDistribution._default_name_count += 1 + name = "material" if n == 0 else f"material_{n}" + self._distribution_name = name + self._solvers = [] + self._mixing_rules = {} + self._read_properties = set() + self._reported = set() # check() speaks once per finding + self._registry.add_callback(self._push_all) + + # -- declaring (shorthand onto the registry) -------------------------- + + @property + def registry(self): + """The :class:`MaterialRegistry` these materials are defined in.""" + return self._registry + + @property + def materials(self): + """The declared materials, in order, bound to this distribution.""" + return tuple(BoundMaterial(self, m) for m in self._registry.materials) + + def _warn_on_shadowed_properties(self, properties): + """A property whose name is also an attribute of this object is + unreachable as ``materials.`` — normal lookup wins and returns + the attribute, silently, with no blend in sight. The two + distributions have different attributes, so a name that works on one + can be shadowed on the other.""" + # NB not hasattr(self, n): that goes through __getattr__, which + # resolves a property by BUILDING the level sets, mid-add(). + own = set(vars(self)) + for klass in type(self).__mro__: + own |= set(vars(klass)) + shadowed = sorted(n for n in properties if n in own) + if shadowed: + warnings.warn( + f"material properties {shadowed} share a name with an " + f"attribute of {type(self).__name__}, so materials. will " + "return the attribute, not the blend. Read them with " + "materials.blend(name) instead, or rename them.", + stacklevel=3, + ) + + def add(self, name, where=None, description="", reference="", **properties): + """Declare a material here — shorthand for ``registry.add(...)``. + + Returns a :class:`BoundMaterial`, so the definition and its place in + this distribution are reachable from one handle. ``where`` is the same + as assigning a region afterwards. """ - Add a callback function for material changes. + self._check_can_declare(name) + self._warn_on_shadowed_properties(properties) + self._registry.add(name, description=description, reference=reference, + **properties) + bound = self[name] + if where is not None: + self[name] = where + return bound + + def _check_can_declare(self, name): + """Subclass hook: raise if it is too late to add a material.""" + + def __getitem__(self, name): + return BoundMaterial(self, self._registry[name]) + + def __len__(self): + return len(self._registry) + + # -- level sets (subclass) -------------------------------------------- + + def build(self): + """Allocate the level sets now, rather than on first read. + + **Collective**: it creates mesh variables, so every rank must call it + together. That is also true of the first *read* of any material + property — ``materials.density``, ``materials["x"].mask``, + ``materials.index``, or ``solver.materials = ...`` — because the read + builds them. Putting such a read inside a rank-local branch + (``if rank == 0: ...``) deadlocks. Call this at a point where every + rank is together if the ordering is ever in doubt. + """ + self._ensure_built() + return self - Parameters - ---------- - callback : callable - Function called as ``callback(event_type, *args)`` + def _ensure_built(self): + raise NotImplementedError + + def _level_sets(self): + """The per-material level-set variables, in registry order.""" + raise NotImplementedError + + def _mask_of(self, definition): + self._ensure_built() + return self._level_sets()[definition.index].sym[0] + + # -- properties -------------------------------------------------------- + + def mixing(self, **rules): + """Choose how a property is blended where masks are fractional. + + ``materials.mixing(shear_viscosity_0="harmonic")``. Only meaningful + where a mask can take a value strictly between 0 and 1; where exactly + one mask is 1 every rule gives the same answer. + + ``"arithmetic"`` (the default) is a Voigt average: for a flux at a + common gradient it is the flux blend. On a sharp contrast it does not + converge with sampling density. ``"harmonic"`` is the Reuss average, + which does. Pick on physical grounds. """ - self._callbacks.append(callback) + declared = self._registry.declared_properties() + for name, rule in rules.items(): + name = _property_name(name) + if rule not in _MIXING_RULES: + raise ValueError( + f"mixing rule for {name!r} must be one of {_MIXING_RULES}, " + f"not {rule!r}" + ) + if name not in declared: + raise KeyError( + f"no material declares {name!r}, so a mixing rule for it " + f"would have no effect; declared: {sorted(declared)}" + ) + self._mixing_rules[name] = rule + self._push_all() + return self - def _notify_callbacks(self, event_type: str, *args): - """Notify all callbacks of a material change""" - for callback in self._callbacks: - try: - callback(event_type, *args) - except Exception as e: - print(f"Warning: Material callback failed: {e}") + def blend(self, name, mixing=None): + """The material-weighted symbol for property ``name``. - def export_config(self) -> Dict[str, Any]: - """Export material configuration""" - return { - "materials": { - name: { - "properties": mat.properties, - "description": mat.description, - "reference": mat.reference, - } - for name, mat in self._materials.items() - }, - "region_assignments": dict(self._region_assignments), - "version": self._version, - } + Usually reached as ``materials.``; call it directly to override + the mixing rule for one use. + """ + name = _property_name(name) + missing = [m.name for m in self._registry.materials + if name not in m.properties] + if missing: + raise KeyError( + f"property {name!r} is not declared by {missing}. Every " + "material must declare a property that is blended — a material " + "with no viscosity is not a material with zero viscosity." + ) + + self._ensure_built() + self._read_properties.add(name) + values = [m.resolved(name) for m in self._registry.materials] + masks = self._level_sets() + rule = mixing or self._mixing_rules.get(name, "arithmetic") + + if rule == "arithmetic": + return sum(masks[i].sym[0] * value for i, value in enumerate(values)) + + # A harmonic blend is 1 / sum(phi_i / v_i), so a material whose value + # is zero poisons the symbol at declaration time: sympy folds 1/0 to + # ComplexInfinity and the JIT then dies with a bare C-printer + # traceback naming neither the material nor the property. A value that + # can VANISH (a law in T, say) is the same defect deferred to run + # time, where it surfaces as Integral = nan behind a RuntimeWarning. + zeros = [ + m.name for m, value in zip(self._registry.materials, values) + if value.is_zero + ] + if zeros: + raise ValueError( + f"harmonic mixing of {name!r} divides by its value in " + f"{zeros}, which is zero. A material with no {name} cannot be " + "blended harmonically; use arithmetic mixing, or give it a " + "small non-zero value." + ) + return 1 / sum(masks[i].sym[0] / value for i, value in enumerate(values)) + + def __getattr__(self, name): + # Only reached when normal lookup fails, so this cannot shadow a real + # attribute. Private names never resolve to a material property. + if name.startswith("_"): + raise AttributeError(name) + try: + registry = object.__getattribute__(self, "_registry") + except AttributeError: + raise AttributeError(name) from None + if not any(name in m.properties for m in registry.materials): + raise AttributeError( + f"{type(self).__name__!r} object has no attribute {name!r}, and " + "no material declares a property of that name" + ) + try: + return self.blend(name) + except KeyError as exc: + # hasattr() swallows AttributeError and nothing else, so a KeyError + # escaping here breaks hasattr / getattr(o, n, default) for every + # caller. Keep the diagnosis, honour the contract. + raise AttributeError(str(exc.args[0] if exc.args else exc)) from None + + # -- the handoff to a solver ------------------------------------------ + + def _attach(self, solver): + """Called by ``solver.materials = ...``.""" + if not any(s is solver for s in self._solvers): + self._solvers.append(solver) + self._push_to(solver) + + def _detach(self, solver): + """Stop pushing to a solver that no longer wants these materials.""" + self._solvers = [s for s in self._solvers if s is not solver] + + def _push_all(self): + for solver in self._solvers: + self._push_to(solver) + + def _push_to(self, solver): + """Set every property the solver's constitutive model recognises.""" + constitutive_model = getattr(solver, "_constitutive_model", None) + if constitutive_model is None: + return # pushed again when the constitutive model is set + if len(self._registry) == 0: + return + parameters = constitutive_model.Parameters + recognised = set(type(parameters)._list_valid_parameters(type(parameters))) + for name in sorted(self._registry.declared_properties() & recognised): + setattr(parameters, name, self.blend(name)) + + def _recognised(self): + names = set() + for solver in self._solvers: + constitutive_model = getattr(solver, "_constitutive_model", None) + if constitutive_model is None: + continue + parameters = constitutive_model.Parameters + names |= set( + type(parameters)._list_valid_parameters(type(parameters)) + ) + return names + + def unclaimed(self): + """Declared properties no attached model recognises and nobody read. + + Not in itself a problem: ``density`` is unclaimed until the model + script reads ``materials.density`` for its body force. It is the list + :meth:`check` looks through for misspellings. + """ + return sorted( + self._registry.declared_properties() + - self._recognised() + - self._read_properties + ) + + def check(self): + """Report a declared property that looks like a misspelled parameter. + + Called from the solver build. A property no constitutive model + recognises is perfectly normal — that is how ``density`` reaches a body + force. What is not normal is one that *nearly* matches a parameter the + model does have and that nothing has read: ``viscocity`` is not + ``viscosity``, it is silently the default viscosity, and nothing else + will ever say so. + """ + import difflib + + recognised = sorted(self._recognised()) + suspicious = [] + for name in self.unclaimed(): + close = difflib.get_close_matches(name, recognised, n=1, cutoff=0.8) + if close: + suspicious.append((name, close[0])) + for name, suggestion in suspicious: + if (name, suggestion) in self._reported: + continue # a nonlinear solve rebuilds repeatedly + self._reported.add((name, suggestion)) + declared = self._registry.declared_properties() + fate = ( + f"{suggestion!r} is set from the materials that DO declare it" + if suggestion in declared + else f"{suggestion!r} keeps its default value" + ) + warnings.warn( + f"material property {name!r} is declared but nothing reads it, " + f"and the constitutive model has a similarly-named parameter " + f"{suggestion!r} — did you mean that? As it stands {fate}, and " + f"nothing reads {name!r}.", + stacklevel=3, + ) + return suspicious + + # -- display ----------------------------------------------------------- + + def _object_viewer(self): + uw.pprint(f"{type(self).__name__}: {len(self._registry)} materials") + for material in self._registry.materials: + uw.pprint(f" [{material.index}] {material.name}") + for key, value in material.properties.items(): + uw.pprint(f" {key} = {value}") + + +class MaterialRegions(MaterialDistribution): + """Materials tied to the mesh rather than to particles. + + For a model whose materials do not move: layers, inclusions, a basin, a + gmsh model with physical groups. The level sets are exact — a material + either owns an integration point or it does not — and there are no + particles to populate, advect or repopulate. + + A region may be given as a mesh label (a gmsh physical group, reaching + Underworld3 as ``mesh.regions``), or as a geometric condition, which is + resolved at the integration points and so keeps sub-cell position. + + Parameters + ---------- + mesh : Mesh + registry : MaterialRegistry, optional + Where the materials are defined. A fresh one is made if omitted, and + :meth:`~MaterialDistribution.add` writes into it. + name : str, optional + Base name for the level-set variables. Defaults to ``material``, then + ``material_1`` and so on, so two distributions on one mesh do not + collide. + + Examples + -------- + >>> materials = uw.MaterialRegions(mesh) + >>> materials.add("mantle", shear_viscosity_0=1.0) + >>> materials.add("crust", shear_viscosity_0=1.0e3) + >>> materials["crust"] = "Crust" # a gmsh physical group + >>> materials["crust"] = mesh.X[1] > 0.8 # or a condition + >>> stokes.materials = materials + + See Also + -------- + underworld3.swarm.MaterialSwarm : the same interface, carried by particles. + """ - def import_config(self, config: Dict[str, Any]): - """Import material configuration from exported dict""" - if "materials" in config: - for name, mat_config in config["materials"].items(): - material = self.create_material( - name, mat_config.get("description", ""), mat_config.get("reference", "") + def __init__(self, mesh, registry=None, name=None): + self.mesh = mesh + self._level_set_vars = None + self._pending = [] + self._init_distribution(registry=registry, name=name) + + # -- building ---------------------------------------------------------- + + def _check_can_declare(self, name): + if self._level_set_vars is not None: + raise RuntimeError( + f"cannot add material {name!r}: the level sets are already " + "built. Declare every material before the first read." + ) + + def _ensure_built(self): + if self._level_set_vars is not None: + return + if len(self._registry) == 0: + raise RuntimeError( + "no materials have been declared — call add() before reading " + "a material property" + ) + _vars_before = len(self.mesh.vars) + self._level_set_vars = [ + uw.discretisation.IntegrationPointVariable( + f"{self._distribution_name}^{{[{i}]}}", self.mesh, + ) + for i in range(len(self._registry)) + ] + self._check_level_sets_are_new( + self.mesh, _vars_before, len(self._registry)) + self._registry._built.append(self) + # material 0 owns everything not claimed by anyone else + self._level_set_vars[0].array[...] = 1.0 + for i in range(1, len(self._level_set_vars)): + self._level_set_vars[i].array[...] = 0.0 + + pending, self._pending = self._pending, [] + for definition, region in pending: + self._paint(definition, region) + + def _level_sets(self): + return self._level_set_vars + + # -- painting ---------------------------------------------------------- + + def __setitem__(self, name, region): + """Assign a region to a material. + + ``region`` may be + + - a mesh label name (``"Crust"``), optionally ``(name, value)`` — the + cells in that label, and every integration point in them; + - a symbolic condition on the mesh coordinates (``mesh.X[1] > 0.8``, + with ``&``, ``|``, ``~``), resolved at each integration point; + - a boolean array over the integration points, or a callable of their + coordinates. + + Assignment is ordered and cumulative: a later region overwrites an + earlier one where they overlap. + """ + definition = self._registry[name] + if self._level_set_vars is None: + self._pending.append((definition, region)) + else: + self._paint(definition, region) + + def _paint(self, definition, region): + """Give ``definition`` exactly ``region``, and nothing else. + + Assignment replaces: ``materials["x"] = A`` then ``materials["x"] = B`` + leaves the material at B, not at A union B. Anything this material + held outside the new region reverts to material 0, which is what + ``x[k] = v`` means in Python and what a user correcting a region in a + notebook cell expects. + """ + selected = self._region_mask(region) + held = np.asarray(self._level_set_vars[definition.index].array).reshape(-1) > 0.5 + released = held & ~selected + for i, var in enumerate(self._level_set_vars): + values = np.asarray(var.array).reshape(-1).copy() + values[selected] = 1.0 if i == definition.index else 0.0 + if released.any(): # back to the default material + values[released] = 1.0 if i == 0 else 0.0 + var.array[:, 0, 0] = values + + def _integration_points(self): + self._ensure_built() + return np.asarray(self._level_set_vars[0].integration_points) + + def _region_mask(self, region): + """A boolean array over the flattened integration points.""" + points = self._integration_points() + ncells, nq, _ = points.shape + flat = points.reshape(-1, points.shape[-1]) + + if isinstance(region, str): + return self._label_mask(region, None, ncells, nq) + if (isinstance(region, tuple) and len(region) == 2 + and isinstance(region[0], str)): + return self._label_mask(region[0], region[1], ncells, nq) + if isinstance(region, np.ndarray): + selected = np.asarray(region).reshape(-1) + if selected.dtype != bool or selected.shape[0] != flat.shape[0]: + raise ValueError( + "a region given as an array must be a boolean array of " + f"length {flat.shape[0]} (the integration points), not " + f"{selected.dtype} of length {selected.shape[0]}" + ) + return selected + if isinstance(region, sympy.Basic): + return _condition_mask(region, flat) + if callable(region): + return np.asarray(region(flat), dtype=bool).reshape(-1) + + raise TypeError( + "a region must be a mesh label name, a symbolic condition on the " + "mesh coordinates, a boolean array over the integration points, or " + f"a callable of their coordinates — not {type(region).__name__}" + ) + + def _label_mask(self, label_name, label_value, ncells, nq): + """Integration points of the cells carried by a mesh label. + + Two PETSc hazards govern the shape of this method. + + ``getStratumIS(v)`` for a value that is not in the label's live value + set **hard-aborts** the process — no exception, no traceback, every + rank gone (the repo warns about this twice; cf. the "Centre" + pseudo-label). And on a rank where the stratum is empty, petsc4py + returns a live ``IS`` wrapping a NULL handle rather than ``None``, so + ``getIndices()`` on it segfaults too. Both are probed for below. + + The value set is also **rank-local**: a label live on one rank can be + absent on another, and a label with two values can look single-valued + to each rank separately. Resolving the value from the local set would + raise on some ranks and not others (a hang), or silently paint the + union of two regions. The set is therefore reduced across ranks + before anything branches on it. + """ + dm = self.mesh.dm + if not dm.hasLabel(label_name): + available = [dm.getLabelName(i) for i in range(dm.getNumLabels())] + regions = getattr(self.mesh, "regions", None) + if regions is not None: + available += [r.name for r in regions] + raise KeyError( + f"the mesh has no label {label_name!r}. Available: " + f"{sorted(set(available))}" + ) + + label = dm.getLabel(label_name) + values_is = label.getValueIS() + local_values = ( + {int(v) for v in values_is.getIndices()} if values_is is not None else set() + ) + + # Reduce first, then branch: every rank must resolve the same value + # and raise the same errors. + global_values = set() + for rank_values in uw.mpi.comm.allgather(local_values): + global_values |= rank_values + + if label_value is None: + regions = getattr(self.mesh, "regions", None) + if regions is not None and label_name in regions.__members__: + label_value = int(regions[label_name].value) + elif len(global_values) == 1: + label_value = int(next(iter(global_values))) + else: + raise ValueError( + f"label {label_name!r} carries values " + f"{sorted(global_values)}; say which with " + f"materials[name] = ({label_name!r}, value)" ) - material.properties = mat_config.get("properties", {}) - if "region_assignments" in config: - self._region_assignments = config["region_assignments"] + label_value = int(label_value) + if label_value not in global_values: + raise ValueError( + f"label {label_name!r} has no stratum with value {label_value} " + f"anywhere on the mesh (live values: {sorted(global_values)}). " + "Asking PETSc for it would abort the run." + ) + + selected = np.zeros((ncells, nq), dtype=bool) + if label_value not in local_values: + return selected.reshape(-1) # live elsewhere, no cells here + + stratum = label.getStratumIS(label_value) + try: + if stratum is None or stratum.getSize() == 0: + return selected.reshape(-1) + cells = np.asarray(stratum.getIndices()) + finally: + if stratum is not None: + stratum.destroy() + + c_start, c_end = dm.getHeightStratum(0) + cells = cells[(cells >= c_start) & (cells < c_end)] - c_start + selected[cells] = True + return selected.reshape(-1) + + +def _condition_mask(condition, coords): + """Evaluate a symbolic condition at ``coords`` -> boolean array.""" + import sympy.logic.boolalg as boolalg + from sympy.core.relational import Relational + + n = coords.shape[0] + + if isinstance(condition, boolalg.BooleanTrue): + return np.ones(n, dtype=bool) + if isinstance(condition, boolalg.BooleanFalse): + return np.zeros(n, dtype=bool) + if isinstance(condition, sympy.And): + return np.logical_and.reduce( + [_condition_mask(a, coords) for a in condition.args]) + if isinstance(condition, sympy.Or): + return np.logical_or.reduce( + [_condition_mask(a, coords) for a in condition.args]) + if isinstance(condition, sympy.Not): + return ~_condition_mask(condition.args[0], coords) + + if isinstance(condition, Relational): + values = np.asarray( + uw.function.evaluate(condition.lhs - condition.rhs, coords) + ).reshape(-1) + if isinstance(condition, sympy.StrictGreaterThan): + return values > 0.0 + if isinstance(condition, sympy.GreaterThan): + return values >= 0.0 + if isinstance(condition, sympy.StrictLessThan): + return values < 0.0 + if isinstance(condition, sympy.LessThan): + return values <= 0.0 + if isinstance(condition, sympy.Eq): + return np.isclose(values, 0.0) + if isinstance(condition, sympy.Ne): + return ~np.isclose(values, 0.0) + + raise TypeError( + f"cannot read {condition!r} as a region: expected a comparison of mesh " + "coordinates, or an &/|/~ combination of them" + ) - def __repr__(self): - return f"MaterialRegistry({len(self._materials)} materials, {len(self._region_assignments)} assignments)" +# --------------------------------------------------------------------------- +# Convenience definitions. The values are SI; non-dimensionalise them, or set +# reference quantities on the model, before using them in a solve. +# --------------------------------------------------------------------------- -# Common material definitions for geodynamics -def create_standard_mantle_material(registry: MaterialRegistry) -> MaterialDefinition: - """Create a standard mantle material with typical properties""" - material = registry.create_material( +def create_standard_mantle_material(registry: MaterialRegistry): + """A standard upper-mantle material (Turcotte & Schubert, 2014).""" + return registry.add( "mantle", description="Standard upper mantle material", reference="Turcotte & Schubert (2014)", + viscosity=1e21, density=3300, + thermal_conductivity=3.0, thermal_diffusivity=1e-6, + thermal_expansion=3e-5, ) - material.set_property(MaterialProperty.VISCOSITY, 1e21) # Pa·s - material.set_property(MaterialProperty.DENSITY, 3300) # kg/m³ - material.set_property(MaterialProperty.THERMAL_CONDUCTIVITY, 3.0) # W/m/K - material.set_property(MaterialProperty.THERMAL_DIFFUSIVITY, 1e-6) # m²/s - material.set_property(MaterialProperty.THERMAL_EXPANSION, 3e-5) # K⁻¹ - - return material - -def create_standard_crust_material(registry: MaterialRegistry) -> MaterialDefinition: - """Create a standard crustal material with typical properties""" - material = registry.create_material( +def create_standard_crust_material(registry: MaterialRegistry): + """A standard continental-crust material (Turcotte & Schubert, 2014).""" + return registry.add( "crust", description="Standard continental crust material", reference="Turcotte & Schubert (2014)", + viscosity=1e22, density=2700, + thermal_conductivity=2.5, thermal_diffusivity=1e-6, + thermal_expansion=3e-5, ) - material.set_property(MaterialProperty.VISCOSITY, 1e22) # Pa·s - material.set_property(MaterialProperty.DENSITY, 2700) # kg/m³ - material.set_property(MaterialProperty.THERMAL_CONDUCTIVITY, 2.5) # W/m/K - material.set_property(MaterialProperty.THERMAL_DIFFUSIVITY, 1e-6) # m²/s - material.set_property(MaterialProperty.THERMAL_EXPANSION, 3e-5) # K⁻¹ - - return material - -def create_high_viscosity_material( - registry: MaterialRegistry, name: str = "high_visc", viscosity_contrast: float = 1000 -) -> MaterialDefinition: - """Create a high viscosity material for inclusion studies""" - material = registry.create_material( +def create_high_viscosity_material(registry: MaterialRegistry, + name: str = "high_visc", + viscosity_contrast: float = 1000): + """A stiff inclusion, for viscosity-contrast studies.""" + return registry.add( name, description=f"High viscosity material (contrast {viscosity_contrast}x)", reference="User defined", + viscosity=1e21 * viscosity_contrast, density=3300, + thermal_conductivity=3.0, thermal_diffusivity=1e-6, ) - - # Base properties similar to mantle - material.set_property(MaterialProperty.VISCOSITY, 1e21 * viscosity_contrast) - material.set_property(MaterialProperty.DENSITY, 3300) - material.set_property(MaterialProperty.THERMAL_CONDUCTIVITY, 3.0) - material.set_property(MaterialProperty.THERMAL_DIFFUSIVITY, 1e-6) - - return material diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 76fefc2e1..7fc16b643 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -121,6 +121,112 @@ class SwarmType(Enum): from underworld3.utilities.dimensionality_mixin import DimensionalityMixin +#: Which sampling rules each proxy target can offer. ``proxy_location`` chooses +#: WHERE the assembler reads the particle data; ``proxy_sampling`` chooses WHAT +#: it reads there. The two axes are not fully independent: a per-cell +#: least-squares fit is a reconstruction by construction, and the Voronoi share +#: needs the cell-major integration-point layout to restrict itself to a cell. +_PROXY_SAMPLING_BY_LOCATION = { + "nodes": ("reconstruct",), + "integration_points": ("reconstruct", "share"), + "cells": ("reconstruct",), +} + +#: ``"nearest"`` is deliberately absent above. Reading one particle's value +#: whole is the MATERIAL mapping, and a material is an +#: :class:`IndexSwarmVariable` — it carries the label, presents one level set +#: per material, and lets the properties be blended by that partition of unity. +#: Carrying a property field directly on the particles and sampling it is the +#: route this does not offer: it puts the constitutive relationship on the +#: particles, where the solver cannot reach the pieces it needs. +_NO_HAND_ROLLED_MATERIALS = ( + "a nearest-particle read is the material mapping, and a material is an " + "IndexSwarmVariable: it carries the label, presents one level set per " + "material, and blends the properties through that partition of unity " + "(see createMask). Building the property field itself on the particles is " + "deliberately not offered — the solver needs the constitutive law, not a " + "sampled answer to it." +) + +_PROXY_SAMPLING_WHY = { + ("nodes", "nearest"): _NO_HAND_ROLLED_MATERIALS, + ("integration_points", "nearest"): _NO_HAND_ROLLED_MATERIALS, + ("cells", "nearest"): _NO_HAND_ROLLED_MATERIALS, + ("cells", "share"): ( + "a 'cells' proxy already uses every particle in the cell, by fitting " + "them. The share is the integration-point equivalent: use " + "proxy_location='integration_points'." + ), + ("nodes", "share"): ( + "the share partitions a cell between its integration points, and a " + "NODE is shared between cells, so there is no cell to restrict to. Use " + "proxy_location='integration_points'." + ), +} + + +def _validate_proxy_sampling(proxy_location, proxy_sampling): + """Check a (location, sampling) pair and return the sampling rule. + + Raises rather than silently ignoring an unavailable combination: a + material mapping that quietly did something other than what was asked is + exactly the failure this machinery exists to prevent. + """ + allowed = _PROXY_SAMPLING_BY_LOCATION.get(proxy_location) + if allowed is None: + raise ValueError( + "proxy_location must be 'nodes', 'integration_points' or 'cells', " + f"not {proxy_location!r}" + ) + if proxy_sampling in allowed: + return proxy_sampling + if proxy_sampling not in ("reconstruct", "nearest", "share"): + raise ValueError( + "proxy_sampling must be 'reconstruct', 'nearest' or 'share', " + f"not {proxy_sampling!r}" + ) + why = _PROXY_SAMPLING_WHY.get((proxy_location, proxy_sampling), "") + raise ValueError( + f"proxy_sampling={proxy_sampling!r} is not available for " + f"proxy_location={proxy_location!r} (offered: {', '.join(allowed)})" + + (" — " + why if why else "") + ) + + +def _validate_index_proxy_sampling(proxy_location, proxy_sampling): + """Resolve ``proxy_sampling`` for a level-set (index) variable. + + The level sets of an :class:`IndexSwarmVariable` are built differently at + each target — a distance-weighted nodal fill, a per-cell least-squares + fraction, or a direct read at the integration points — so the sampling + axis only opens up for the last of those. ``None`` means "the sensible + one for this location". + """ + if proxy_location not in _PROXY_SAMPLING_BY_LOCATION: + raise ValueError( + "proxy_location must be 'nodes', 'integration_points' or 'cells', " + f"not {proxy_location!r}" + ) + if proxy_location != "integration_points": + if proxy_sampling not in (None, "reconstruct"): + raise ValueError( + f"proxy_sampling={proxy_sampling!r} is only available for " + "proxy_location='integration_points'. A 'nodes' level set is a " + "distance-weighted fill and a 'cells' level set is a " + "least-squares fraction; neither takes a sampling rule." + ) + return "reconstruct" + if proxy_sampling is None: + return "nearest" + if proxy_sampling not in ("nearest", "share"): + raise ValueError( + "proxy_sampling for an IndexSwarmVariable at the integration points " + f"must be 'nearest' (a material label, sharp) or 'share' (material " + f"fractions from every particle), not {proxy_sampling!r}" + ) + return proxy_sampling + + class SwarmVariable(DimensionalityMixin, MathematicalMixin, Stateful, uw_object): r""" Variable supported by a particle swarm (point cloud). @@ -152,13 +258,19 @@ class SwarmVariable(DimensionalityMixin, MathematicalMixin, Stateful, uw_object) reconstructed from the nearest particles at every integration point and read there directly, with no second interpolation (the Ellipsis / Underworld PIC-LIP mapping); a material interface keeps - its sub-cell position, and the proxy has no gradient (a derivative - of its symbol is refused). ``proxy_degree`` / ``proxy_continuous`` - are ignored in that case. ``"cells"``: a discontinuous mesh variable + its sub-cell position, and the proxy has no gradient of its own: a + derivative of its symbol in a WEAK FORM is refused (see ``"cells"`` + below for the remedy), while ``uw.function.evaluate`` of the same + derivative answers by fitting the values per cell first. + ``proxy_degree`` / ``proxy_continuous`` are ignored in that case. ``"cells"``: a discontinuous mesh variable of ``proxy_degree`` holding, in every cell, the least-squares polynomial through the particles that cell holds (a thin cell takes a linear fit to the particles nearest its centroid, an empty cell - keeps its previous value). Exact for + keeps its previous value). THIS is the target to choose when a solve + needs a gradient: the level sets are polynomials, so they + differentiate directly in a weak form, with no projection solve + (degree 2 recovers the gradient of a quadratic particle field to + 2e-7). Exact for polynomial particle fields up to ``proxy_degree``, integrated exactly by the default rule, sharp at cell edges, with a gradient, and no neighbour search across ranks; see @@ -169,6 +281,31 @@ class SwarmVariable(DimensionalityMixin, MathematicalMixin, Stateful, uw_object) proxy_continuous : bool, default=True Whether the proxy uses continuous (True) or discontinuous (False) interpolation. + proxy_sampling : {"reconstruct", "share"}, default="reconstruct" + How the particle values become a value at each proxy point. + ``proxy_location`` says *where* the assembler reads; this says *what* + it reads there. + + - ``"reconstruct"``: a weighted fit over the ``nnn`` nearest + particles, exact for linear fields. What a smooth field wants. + - ``"share"``: the mean over the particles this point speaks for — + those whose nearest integration point *within their own cell* is + this one (the cell-restricted Voronoi share). Every particle reaches + the assembly, weighted by the sub-region it represents, which is + what a HISTORY wants: stress, damage, accumulated strain, anything + whose value was earned by being advected. Requires + ``proxy_location="integration_points"``. + + There is no nearest-particle option here. Sampling one particle's + value whole is the *material* mapping, and materials are built from + :class:`IndexSwarmVariable`, which carries the label and presents one + level set per material for the properties to be blended by. Putting a + property field on the particles and sampling it hands the solver an + answer where it needs a constitutive law, so it is not offered. + + ``"cells"`` fits a polynomial per cell and takes ``"reconstruct"`` + only; ``"share"`` needs the cell-major integration-point layout. An + unavailable combination raises rather than being quietly ignored. varsymbol : str, optional LaTeX symbol for display. Defaults to ``name``. rebuild_on_cycle : bool, default=True @@ -214,6 +351,7 @@ def __init__( proxy_degree=1, proxy_continuous=True, proxy_location="nodes", + proxy_sampling="reconstruct", _register=True, _proxy=True, varsymbol=None, @@ -405,6 +543,7 @@ def __init__( ) self._cell_projector = None self._proxy_location = proxy_location + self._proxy_sampling = _validate_proxy_sampling(proxy_location, proxy_sampling) self._create_proxy_variable() # Inert: kept for backward compatibility with the removed @@ -1351,6 +1490,60 @@ def _cells_to_meshVar(self, meshVar): return # Maybe rbf_interpolate for this one and meshVar is a special case + def _share_to_integration_points(self, meshVar, values): + r"""Values at the integration points by the cell-restricted Voronoi share. + + Each particle is assigned to the nearest integration point *of its own + cell*, and each point reads the mean over the particles assigned to + it. Every particle reaches the assembly exactly once, weighted by the + sub-region of the cell it represents — the nearest-particle rule, by + contrast, discards whichever particles are not closest to a rule point, + which at ten per cell is most of them. + + A rule point whose share is empty (an under-filled or empty cell) falls + back to the nearest particle anywhere on this rank, which is the same + answer ``proxy_sampling="nearest"`` would have given it. The count of + such points is kept on ``self._share_empty`` — a persistently non-zero + value means the swarm is too thin for the rule, and + :meth:`Swarm.repopulate` is the fix. + + Parameters + ---------- + meshVar : IntegrationPointVariable + The proxy being filled. + values : ndarray, shape (Np, ncomp) + Particle values, non-dimensional. + + Returns + ------- + ndarray, shape (ncells * Nq, ncomp) + """ + from underworld3.utilities.particle_share import ( + share_assignment, + share_average, + ) + + ipc = np.asarray(meshVar.integration_points) # (ncells, Nq, cdim) + npoints = ipc.shape[0] * ipc.shape[1] + values = np.asarray(values, dtype=float) + values = values.reshape(values.shape[0], -1) + + flat = share_assignment( + ipc, np.asarray(self.swarm._particle_coordinates.data), + self.swarm._owning_cells(), + ) + means, counts = share_average(flat, values, npoints) + + empty = counts == 0 + self._share_empty = int(empty.sum()) + if self._share_empty: + _, nearest = self.swarm._get_kdtree().query( + ipc.reshape(-1, ipc.shape[-1])[empty], k=1, sqr_dists=False + ) + means[empty] = values[np.asarray(nearest).reshape(-1)] + + return means + def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False, order=1, monotone=False): """ @@ -1407,6 +1600,10 @@ def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False, order=1, stacklevel=2, ) Values = current_values + elif getattr(self, "_proxy_sampling", "reconstruct") == "share": + Values = self._share_to_integration_points( + meshVar, self.unpack_raw_data_from_petsc(squeeze=False) + ) elif monotone: # The limiter is data-dependent, so it cannot ride on a cached # geometry-only operator; take the direct path. @@ -2442,6 +2639,49 @@ class IndexSwarmVariable(SwarmVariable): Polynomial degree for mesh projection (default 1). proxy_continuous : bool Whether mesh proxy is continuous (default True). + proxy_location : {"integration_points", "nodes", "cells"}, default="integration_points" + Where the level sets live. + + ``"integration_points"`` (the default): the level sets are stored at + the points where the assembler evaluates the weak form, so a material + interface keeps its sub-cell position and each mask is exactly 0 or 1 + (the Ellipsis / Underworld particle-in-cell material mapping). These + level sets have no gradient of their own, so a derivative of a mask in + a weak form is refused rather than silently answered with zero; + ``uw.function.evaluate`` still answers, by fitting per cell. + + ``"nodes"``: a continuous field per material. A node on an interface + averages both materials, so the cells either side see a property that + is neither, and the smear is about one cell wide however many + particles you add — it is a property of the basis, not of the swarm. + Kept for continuity with existing models; measurably the worst of the + three at an interface (layered Couette: 8.0e-2 against 1.8e-7). + + ``"cells"``: a polynomial material fraction per cell, clamped to + [0, 1] and renormalised. Sharp at cell edges and DIFFERENTIABLE — the + one to choose when a solve needs the gradient of a material property. + + ``proxy_continuous`` applies only to ``"nodes"``. + proxy_sampling : {"nearest", "share"}, optional + Only for ``proxy_location="integration_points"``; defaults to + ``"nearest"``. + + ``"nearest"``: each integration point takes the material of its + nearest particle, so the masks are 0 or 1 and sum to 1 by + construction. Nothing is mixed, so no mixing assumption is made. The + error falls as particles are added until the quadrature rule limits + it — about eight particles per cell on P2 velocity, after which the + mesh, not the swarm, is what to refine. + + ``"share"``: each integration point takes the material FRACTIONS of + the particles it speaks for (the cell-restricted Voronoi share), so a + cell the interface crosses carries fractional masks. Use it when the + material genuinely is a sub-cell mixture rather than an interface, and + note the trap: :meth:`createMask` blends properties arithmetically, + which is a Voigt (equal-strain-rate) average, and on a sharp contrast + that does not converge with particle density — 3.5e-2 flat, against + 1.2e-2 for the same masks blended harmonically. Fractions are only + worth having when you have chosen the mixing rule deliberately. Examples -------- @@ -2463,6 +2703,8 @@ def __init__( indices=1, proxy_degree=1, proxy_continuous=True, + proxy_location="integration_points", + proxy_sampling=None, update_type=0, npoints=5, radius=0.5, @@ -2471,6 +2713,18 @@ def __init__( varsymbol=None, ): self.indices = indices + proxy_sampling = _validate_index_proxy_sampling(proxy_location, proxy_sampling) + if update_type != 0 and proxy_location != "nodes": + import warnings + + warnings.warn( + f"update_type={update_type} selects between two NODAL fill " + f"algorithms and has no effect at proxy_location=" + f"{proxy_location!r}; pass proxy_location='nodes' if that is " + "what you meant.", + stacklevel=2, + ) + self._cell_projector = None self.nnn = npoints self.radius_s = radius # **2 # changed to radius self.update_type = update_type @@ -2490,16 +2744,12 @@ def __init__( _proxy=False, varsymbol=varsymbol, ) - """ - vtype = (None,) - dtype = (float,) - proxy_degree = (1,) - proxy_continuous = (True,) - _register = (True,) - _proxy = (True,) - varsymbol = (None,) - rebuild_on_cycle = (True,) - """ + # AFTER super().__init__, which sets _proxy_location / _proxy_sampling + # from its own defaults (this class does not forward the arguments, + # since the base single-proxy _meshVar is not built here at all). + self._proxy_location = proxy_location + self._proxy_sampling = proxy_sampling + # The indices variable defines how many "level set" maps we create as components in the proxy variable import sympy @@ -2508,13 +2758,21 @@ def __init__( self._meshLevelSetVars = [None] * self.indices for i in range(indices): - self._meshLevelSetVars[i] = uw.discretisation.MeshVariable( - name + R"^{[" + str(i) + R"]}", - self.swarm.mesh, - num_components=1, - degree=proxy_degree, - continuous=proxy_continuous, - ) + lname = name + R"^{[" + str(i) + R"]}" + if proxy_location == "integration_points": + self._meshLevelSetVars[i] = uw.discretisation.IntegrationPointVariable( + lname, self.swarm.mesh, + ) + elif proxy_location == "cells": + self._meshLevelSetVars[i] = uw.discretisation.MeshVariable( + lname, self.swarm.mesh, num_components=1, + degree=proxy_degree, continuous=False, + ) + else: + self._meshLevelSetVars[i] = uw.discretisation.MeshVariable( + lname, self.swarm.mesh, num_components=1, + degree=proxy_degree, continuous=proxy_continuous, + ) self._MaskArray[0, i] = self._meshLevelSetVars[i].sym[0, 0] # Initialize lazy evaluation state @@ -2726,6 +2984,95 @@ def view(self): uw.pprint(f"IndexSwarmVariable {self}") uw.pprint(f"Numer of indices {self.indices}") + def _update_index_proxies_from_particles(self): + r"""Fill the level sets directly from the particles, with no nodal step. + + ``proxy_location="integration_points"`` with ``proxy_sampling="nearest"`` + (the default): every integration point takes the material of its + NEAREST PARTICLE, so each level set is exactly 0 or 1 there and the + masks sum to 1 by construction. That is the Ellipsis / Underworld + particle-in-cell material mapping: the interface keeps its sub-cell + position, no node averages two materials, and no reconstruction can + overshoot into a negative viscosity. + + ``proxy_sampling="share"``: every integration point takes the material + FRACTIONS of the particles it speaks for — those whose nearest + integration point within their own cell is this one. Every particle + contributes; a cell the interface crosses carries fractional masks. + See the class docstring for when that is what you want, and for the + mixing-rule trap that comes with it. + + ``proxy_location="cells"``: each level set is the least-squares + polynomial through the cell's own particle indicators + (:class:`~underworld3.utilities.cell_polynomial_projection.CellPolynomialProjector`), + clamped to :math:`[0, 1]` and renormalised so the masks still sum to 1. + A material fraction per cell, with a gradient, sharp at cell edges. + """ + from underworld3.utilities.cell_polynomial_projection import CellPolynomialProjector + + # Collective read/write sequence: every rank walks the same variables + # (only the values differ), as in the nodal path's starved-rank guard. + # One particle is enough here, unlike the nodal path's weighted + # average: the nearest-particle answer is well defined from a single + # particle, and the cell fit falls back to its patch. Only a rank with + # NO particles has nothing to say. + starved = self.swarm.local_size < 1 + if starved: + if self.swarm._population_generation > 0: + import warnings + + warnings.warn( + f"IndexSwarmVariable proxy update: rank {uw.mpi.rank} holds " + f"{max(self.swarm.local_size, 0)} particles; level-set " + f"variables for '{self.clean_name}' left unchanged on this rank.", + stacklevel=2, + ) + for var in self._meshLevelSetVars: + var.data[:, 0] = var.data[:, 0] # keep, but write collectively + return + + Xp = np.asarray(self.swarm._particle_coordinates.data) + idx = np.asarray(self.data).reshape(-1).astype(int) + indicator = (idx[:, None] == np.arange(self.indices)[None, :]).astype(float) + + if self._proxy_location == "integration_points": + if self._proxy_sampling == "share": + # A mean of rows that each sum to 1 still sums to 1, and the + # empty-share fallback is one particle's 0/1 row, so the masks + # remain a partition of unity whatever the sampling density. + U = self._share_to_integration_points( + self._meshLevelSetVars[0], indicator + ) + else: + # Every level set is stored at the same points, so the + # nearest-particle lookup is done once for all of them. + tree = uw.kdtree.KDTree(Xp) + _, nearest = tree.query( + np.asarray(self._meshLevelSetVars[0].coords_nd), + k=1, sqr_dists=False, + ) + U = indicator[np.asarray(nearest).reshape(-1)] + for ii, var in enumerate(self._meshLevelSetVars): + var.data[:, 0] = U[:, ii] + return + + # "cells": one fit for every index at once (the level sets share a basis) + projector = self._cell_projector + if (projector is None or projector.var is not self._meshLevelSetVars[0] + or projector.mesh_version != self.swarm.mesh._mesh_version): + projector = CellPolynomialProjector(self._meshLevelSetVars[0]) + self._cell_projector = projector + old = np.column_stack([np.asarray(v.data[:, 0]) for v in self._meshLevelSetVars]) + U = projector.fit(Xp, indicator, old=old) + # A fitted indicator can leave [0, 1]; clamp, then renormalise so the + # masks remain a partition of unity (createMask stays a weighted mean). + U = np.clip(U, 0.0, 1.0) + total = U.sum(axis=1) + good = total > 1.0e-12 + U[good] /= total[good, None] + for ii, var in enumerate(self._meshLevelSetVars): + var.data[:, 0] = U[:, ii] + def _update_proxy_variables(self): """ This method updates the proxy mesh (vector) variable for the index variable on the current swarm locations @@ -2744,7 +3091,12 @@ def _update_proxy_variables(self): update_type 0: assign the particles to the nearest mesh_levelset nodes, and calculate the value on nodes from them. update_type 1: calculate the material property value on mesh_levelset nodes from the nearest N particles directly. + ``proxy_location`` other than ``"nodes"`` takes neither route: see + :meth:`_update_index_proxies_from_particles`. """ + if self._proxy_location != "nodes": + self._update_index_proxies_from_particles() + return # Starved-rank guard (SWARM-07): with <= 1 local particles the # nearest-neighbour machinery cannot run — KDTree construction on an # empty coordinate array raises IndexError, aborting/hanging the @@ -3178,6 +3530,7 @@ def _invalidate_canonical_data(self): # Invalidate cached spatial index self._kdtree = None + self._owning_cells_cache = None def _flush_pending_petsc_sync(self): """Pack canonical arrays written while migration was suppressed. @@ -3330,6 +3683,31 @@ def _proxy_interpolation_operator(self, meshVar, nnn, p, order): self._proxy_interpolation_cache[key] = (kdtree, operator) return operator + def _owning_cells(self): + """Local owning cell of every particle, cached until they move. + + A UW3 swarm is ``DMSWARM_BASIC``: PETSc keeps no cell id for us, so + the cell has to be located. That is the expensive half of any + cell-local particle operation, and several of them (the Voronoi share, + a population census, a per-cell fit) want the same answer within one + update, so it is cached here and dropped wherever ``_kdtree`` is. + + Returns + ------- + ndarray, shape (local_size,), int64 + Local cell index, or ``-1`` for a particle the local mesh does not + contain (one in flight between ranks, or just outside an open + boundary). + """ + if getattr(self, "_owning_cells_cache", None) is None: + X = np.asarray(self._particle_coordinates.data) + self._owning_cells_cache = ( + np.asarray(self.mesh._robust_owning_cells(X), dtype=np.int64) + if X.shape[0] + else np.zeros(0, dtype=np.int64) + ) + return self._owning_cells_cache + def _get_kdtree(self): """ Return a cached KDTree for the swarm particle coordinates. @@ -5040,6 +5418,7 @@ def repopulate( values=None, nnn=None, order=0, + nearest=None, verbose=False, ): """Add particles to cells that hold too few, remove from cells that hold @@ -5077,6 +5456,12 @@ def repopulate( order : {0, 1}, optional RBF reconstruction order for new particles: 0 bounded (default), 1 linear-exact. + nearest : variable, name, or list of them, optional + Variables (or names) a new particle takes whole from its nearest + existing neighbour instead of by reconstruction. INTEGER-valued + variables are always in this set: a material index is a label, and + the average of two labels is not a label. Use it for any other + field that must stay one of its own values. Returns ------- @@ -5100,7 +5485,13 @@ def repopulate( self._flush_pending_petsc_sync() X = np.array(self._particle_coordinates.data, copy=True) if self.local_size > 0 \ else np.zeros((0, dim)) - cells = np.asarray(mesh._robust_owning_cells(X), dtype=np.int64) if X.shape[0] else np.zeros(0, np.int64) + # The census and the Voronoi-share proxy fill want the same answer, and + # locating particles is the expensive half of both, so take the swarm's + # cached assignment (dropped whenever the particles move). + cells = self._owning_cells() + if cells.shape[0] != X.shape[0]: # cache raced the copy above + cells = (np.asarray(mesh._robust_owning_cells(X), dtype=np.int64) + if X.shape[0] else np.zeros(0, np.int64)) npc = np.bincount(cells[cells >= 0], minlength=ncells) lat_cells = np.asarray(mesh._robust_owning_cells(lattice), dtype=np.int64) owned = np.zeros(ncells, dtype=bool) @@ -5157,9 +5548,19 @@ def repopulate( nnn = min(nnn, max(n_old, 1)) rbf_order = order if nnn >= dim + 2 else 0 operator = None + nearest_row = None if n_old > 0: - operator = uw.kdtree.KDTree(X).interpolation_matrix( + tree = uw.kdtree.KDTree(X) + operator = tree.interpolation_matrix( np.asarray(new_coords), nnn=nnn, p=2, order=rbf_order) + _, nearest_row = tree.query(np.asarray(new_coords), k=1) + nearest_row = np.asarray(nearest_row).reshape(-1) + nearest_spec = nearest or () + if isinstance(nearest_spec, str) or not hasattr(nearest_spec, "__iter__"): + nearest_spec = (nearest_spec,) # a bare name or variable + nearest_set = set() + for item in nearest_spec: + nearest_set.add(item if isinstance(item, str) else getattr(item, "clean_name", item)) # raw values of every variable at the old particles, BEFORE the add raw_old = {} for name, var in self._vars.items(): @@ -5192,15 +5593,24 @@ def repopulate( continue spec = values.get(var, values.get(name, values.get(var.clean_name))) ncomp = raw_old[name].shape[1] if n_old > 0 else var.num_components + # A label cannot be averaged: integer variables (a material + # index) take their nearest neighbour's value whole. + take_nearest = ( + name in nearest_set + or var.clean_name in nearest_set + or np.issubdtype(np.dtype(getattr(var, "_petsc_dtype", float)), np.integer) + ) if spec is not None: vals = spec(np.asarray(new_coords)) if callable(spec) else spec vals = np.broadcast_to(np.asarray(vals, dtype=float).reshape(n_new, -1) if np.ndim(vals) > 0 else vals, (n_new, ncomp)) + elif take_nearest and nearest_row is not None: + vals = raw_old[name][nearest_row] elif operator is not None: vals = operator @ raw_old[name] else: vals = np.zeros((n_new, ncomp)) f = self.dm.getField(var.clean_name).reshape((-1, ncomp)) - f[n_old:, :] = np.asarray(vals).reshape(n_new, ncomp) + f[n_old:, :] = np.asarray(vals).reshape(n_new, ncomp).astype(f.dtype, copy=False) self.dm.restoreField(var.clean_name) added = n_new self._invalidate_canonical_data() @@ -5713,3 +6123,9 @@ def advection( ## - PIC layouts of particles are not directly available / must be done by hand ## - No automatic migration - must compute ranks for the particle swarms ## - No automatic definition of coordinate fields (need to add by hand) + + +# Materials live in their own module (this one is long enough) but belong to +# the swarm namespace: a MaterialSwarm IS a Swarm. Imported at the end so the +# submodule can import Swarm and IndexSwarmVariable from here. +from underworld3.swarm_materials import MaterialSwarm # noqa: E402,F401 diff --git a/src/underworld3/swarm_materials.py b/src/underworld3/swarm_materials.py new file mode 100644 index 000000000..0538c4f57 --- /dev/null +++ b/src/underworld3/swarm_materials.py @@ -0,0 +1,252 @@ +r"""Materials carried by particles. + +The particle half of the material system: a :class:`MaterialSwarm` is a +:class:`~underworld3.swarm.Swarm` that carries a +:class:`~underworld3.materials.MaterialRegistry`, so the materials advect with +the flow. Use it when the material moves; use +:class:`~underworld3.materials.MaterialRegions` when it does not. + + materials = uw.swarm.MaterialSwarm(mesh, fill_param=3) + + materials.add("mantle", shear_viscosity_0=1.0, density=3300) + materials.add("slab", shear_viscosity_0=1.0e3, density=3400) + + materials["slab"] = mesh.X[1] > 0.53 + + stokes.materials = materials + stokes.bodyforce = -materials.density * mesh.CoordinateSystem.unit_e_1 + +Everything about *what* a material is — declaring it, its properties, blending +them into the symbol a solver reads — lives in +:mod:`underworld3.materials` and is shared with the mesh-region distribution. +What is here is *where*: the particles, the integer label they carry, and the +level sets built from it. + +Why a partition of unity, when the masks are 0 or 1 +--------------------------------------------------- +If every material's property were a *number*, the level sets would be +redundant: exactly one mask is 1 at each integration point, so +:math:`\sum_i \phi_i \eta_i` is a select, and one stored coefficient field +would do the same job. The masks earn their place the moment a property is a +*law* — + + materials.add("crust", shear_viscosity_0=eta_0 * sympy.exp(-T.sym[0])) + +— because there is then no number to store, and the only way to combine N +expressions into one symbol the assembler can compile is the weighted sum. +That is also why a material property must not be evaluated on the particles +and sampled: the solver needs the law, not an answer to it. +""" + +import numpy as np +import sympy + +import underworld3 as uw +from underworld3.materials import MaterialDistribution, _condition_mask +from underworld3.swarm import IndexSwarmVariable, Swarm + +__all__ = ["MaterialSwarm"] + + +class MaterialSwarm(MaterialDistribution, Swarm): + """A swarm that carries materials. + + Parameters + ---------- + mesh : Mesh + The mesh the particles live on. + registry : MaterialRegistry, optional + Where the materials are defined. A fresh one is made if omitted, and + :meth:`add` writes into it; pass one to share definitions with another + model or with a :class:`~underworld3.materials.MaterialRegions`. + fill_param : int, default 3 + Particles per cell at population. The swarm is populated lazily, the + first time anything needs particles, so that every ``add`` can happen + first; call :meth:`populate` to force it. + proxy_sampling : {"nearest", "share"}, default "nearest" + How each integration point reads the particles. ``"nearest"`` gives + each point the material of its nearest particle — masks are exactly 0 + or 1, nothing is mixed and no mixing rule is implied. ``"share"`` + gives fractional masks in a cell the interface crosses; use it only + when the material genuinely *is* a sub-cell mixture, and set the + mixing rule deliberately (see + :meth:`~underworld3.materials.MaterialDistribution.mixing`). + proxy_location : {"integration_points", "cells", "nodes"} + Where the level sets live. The default reads the material where the + assembler evaluates the weak form. ``"cells"`` is the one to choose if + a solve needs the *gradient* of a material property. + name : str, optional + Base name for the underlying index variable and its level sets. + Defaults to ``material``, then ``material_1`` and so on, so two + distributions on one mesh do not collide. + + Examples + -------- + >>> materials = uw.swarm.MaterialSwarm(mesh, fill_param=3) + >>> materials.add("mantle", shear_viscosity_0=1.0, density=3300) + >>> materials.add("slab", shear_viscosity_0=1.0e3, density=3400) + >>> materials["slab"] = mesh.X[1] > 0.53 + >>> stokes.materials = materials + + It is a swarm, so it also advects and repopulates: + + >>> materials.population_control = dict(min_per_cell=8) + >>> materials.advection(v.sym, dt) + + See Also + -------- + underworld3.materials.MaterialRegions : the same interface, tied to the mesh. + underworld3.materials.MaterialRegistry : where the definitions live. + """ + + def __init__( + self, + mesh, + registry=None, + fill_param=3, + proxy_sampling="nearest", + proxy_location="integration_points", + name=None, + recycle_rate=0, + verbose=False, + clip_to_mesh=True, + ): + # Material bookkeeping BEFORE the Swarm constructor: Swarm touches + # attributes that __getattr__ would otherwise try to resolve as a + # material property. + self._material_pending = [] + self._material_fill_param = fill_param + self._material_proxy_sampling = proxy_sampling + self._material_proxy_location = proxy_location + self._index_var = None + self._init_distribution(registry=registry, name=name) + + Swarm.__init__( + self, mesh, recycle_rate=recycle_rate, verbose=verbose, + clip_to_mesh=clip_to_mesh, + ) + + # -- the index variable, and the particles ---------------------------- + + @property + def index(self): + """The underlying :class:`~underworld3.swarm.IndexSwarmVariable`. + + The machinery, exposed for the cases this interface does not cover. + A model should not normally need it. + """ + self._ensure_built() + return self._index_var + + def _check_can_declare(self, name): + if self._index_var is not None: + raise RuntimeError( + f"cannot add material {name!r}: the swarm is already populated " + f"with {len(self._registry)} materials and their level sets. " + "Declare every material before the first read (or before " + "calling populate)." + ) + + def populate(self, fill_param=None): + """Create the particles and the level sets, and apply any painting. + + Called automatically the first time anything reads the materials; call + it explicitly to fix the moment, or to override ``fill_param``. + """ + if fill_param is not None: + self._material_fill_param = fill_param + self._ensure_built() + return self + + def _ensure_built(self): + if self._index_var is not None: + return + if len(self._registry) == 0: + raise RuntimeError( + "no materials have been declared — call add() before reading " + "a material property" + ) + + _vars_before = len(self.mesh.vars) + self._index_var = IndexSwarmVariable( + self._distribution_name, + self, + indices=len(self._registry), + proxy_location=self._material_proxy_location, + proxy_sampling=( + self._material_proxy_sampling + if self._material_proxy_location == "integration_points" + else None + ), + ) + # local_size is -1, not 0, on a swarm that has never been populated. + if self.local_size <= 0: + Swarm.populate(self, fill_param=self._material_fill_param) + + self._check_level_sets_are_new( + self.mesh, _vars_before, len(self._registry)) + self._registry._built.append(self) + + pending, self._material_pending = self._material_pending, [] + for definition, region in pending: + self._paint(definition, region) + + def _level_sets(self): + return self._index_var._meshLevelSetVars + + # -- painting --------------------------------------------------------- + + def __setitem__(self, name, region): + """Put a material wherever ``region`` is true. + + ``region`` may be a symbolic condition on the mesh coordinates + (``mesh.X[1] > 0.53``, with ``&``, ``|``, ``~``), a boolean array over + the particles, or a callable of the particle coordinate array. + Painting is ordered and cumulative: a later region overwrites an + earlier one where they overlap. + """ + definition = self._registry[name] + if self._index_var is None: + # Deferred so that every add() can still happen; applied in order + # at population time. + self._material_pending.append((definition, region)) + else: + self._paint(definition, region) + + def _paint(self, definition, region): + """Give ``definition`` exactly ``region``, and nothing else — see + ``MaterialRegions._paint``; assignment replaces rather than unions.""" + selected = self._region_mask(region) + labels = np.asarray(self._index_var.array).reshape(-1) + released = (labels == definition.index) & ~selected + with uw.synchronised_array_update(): + self._index_var.array[selected, 0, 0] = definition.index + if released.any(): + self._index_var.array[released, 0, 0] = 0 + + def _region_mask(self, region): + """A boolean array over the particles.""" + coords = np.asarray(self._particle_coordinates.data) + n = coords.shape[0] + + if isinstance(region, np.ndarray): + selected = np.asarray(region).reshape(-1) + if selected.dtype != bool or selected.shape[0] != n: + raise ValueError( + f"a region given as an array must be a boolean array of " + f"length {n} (the local particle count), not " + f"{selected.dtype} of length {selected.shape[0]}" + ) + return selected + + if isinstance(region, sympy.Basic): + return _condition_mask(region, coords) + + if callable(region): + return np.asarray(region(coords), dtype=bool).reshape(-1) + + raise TypeError( + "a region must be a symbolic condition on the mesh coordinates " + "(e.g. mesh.X[1] > 0.53), a boolean array over the particles, or " + f"a callable of the coordinate array — not {type(region).__name__}" + ) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 9f29b09b4..46309ccdb 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -3915,6 +3915,15 @@ class Lagrangian_Swarm(_DDtBase): ``"cells"`` fits a least-squares polynomial of degree ``degree`` per cell, exact for polynomial histories and integrated exactly by the default rule. + proxy_sampling : {"reconstruct", "share"}, optional + How each slot's proxy reads the particles (default ``"reconstruct"``). + ``"share"`` requires ``proxy_location="integration_points"`` and gives + every integration point the mean of the particles it speaks for — the + cell-restricted Voronoi share. For a HISTORY that is usually what is + wanted: every particle's state contributes, the average is bounded by + the particle values so it cannot invent a stress the swarm never held, + and the stencil cannot reach across a cell wall into another material. + See :doc:`../../advanced/particle-population-and-materials`. step_averaging : int, optional Number of steps for history averaging (default ``2``). @@ -3971,6 +3980,7 @@ def __init__( smoothing=0.0, step_averaging=2, proxy_location="nodes", + proxy_sampling="reconstruct", particle_update="pic", residual_retention=1.0, ): @@ -4000,6 +4010,11 @@ def __init__( # (discontinuous, degree `degree`), exact for polynomial histories # and integrated exactly by the default rule. self.proxy_location = proxy_location + # "share": each integration point averages the particles whose nearest + # point within their own cell it is, so every particle's history + # reaches the assembly and the average stays inside the range the + # particles hold. + self.proxy_sampling = proxy_sampling self._init_history_tracking(order) @@ -4033,6 +4048,7 @@ def _initialise_before_first_move(): proxy_degree=degree, proxy_continuous=continuous, proxy_location=proxy_location, + proxy_sampling=proxy_sampling, varsymbol=rf"{varsymbol}^{{ {'*'*(i+1)} }}", ) ) diff --git a/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index 2d1ebbe11..baa3ccf6f 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -957,8 +957,13 @@ def _no_derivative(self, printer): raise RuntimeError( f"{self.__class__.__name__}: derivative of an integration-point " "variable has no meaning (the field is defined only at the " - "quadrature points). Remove the derivative or project the " - "variable onto a nodal MeshVariable first." + "quadrature points), so the gradient here would be a silent " + "zero. This is refused in a WEAK FORM only, where the " + "discretisation is yours to choose: build the variable with " + "proxy_location='cells' instead, whose level sets are a " + "least-squares polynomial per cell and differentiate directly. " + "uw.function.evaluate() of the same derivative does answer: as " + "a query it recovers the gradient from a per-cell fit for you." ) for var in varlist: diff --git a/src/underworld3/utilities/particle_share.py b/src/underworld3/utilities/particle_share.py new file mode 100644 index 000000000..f843d1822 --- /dev/null +++ b/src/underworld3/utilities/particle_share.py @@ -0,0 +1,118 @@ +r"""The cell-restricted Voronoi share: which particles an integration point speaks for. + +A particle method has to answer one question before it can assemble anything: +what value does the quadrature rule read at each of its points? Taking the +nearest particle answers it sharply but throws most of the swarm away — at ten +particles per cell and six rule points, better than half of them never reach +the assembly at all. That is fine for a *label*, where sub-sampling costs +nothing but sharpness, and wrong for a *history*, where every particle carries +state that was earned by being advected. + +The share is the middle course, and it is what the classic Voronoi +particle-in-cell integration was reaching for. Each particle is assigned to the +nearest integration point **of its own cell**, so the rule points partition the +cell between them and every particle lands in exactly one part. An integration +point then reads the mean over the sub-region it represents. No particle is +discarded, nothing is smeared across a cell boundary, and the whole thing is a +gather with no tree: the owning cell is already known, and the choice within a +cell is over :math:`N_q` candidates. + +The restriction to the particle's own cell is not a detail. Without it a +particle just across a cell wall is often nearer to a rule point on the far +side, and a material interface leaks into the neighbouring cell — which +defeats the reason for sampling at the integration points at all. + +See Also +-------- +underworld3.utilities.cell_polynomial_projection : the fitted alternative, + which gives a differentiable field per cell rather than a value per point. +""" + +import numpy as np + +__all__ = ["share_assignment", "share_average"] + + +def share_assignment(integration_points, coords, cells, chunk_bytes=64 << 20): + """Assign each particle to the nearest integration point of its own cell. + + Parameters + ---------- + integration_points : ndarray, shape (ncells, Nq, cdim) + The physical rule points, in local cell order — the layout + :attr:`IntegrationPointVariable.integration_points` returns. + coords : ndarray, shape (Np, cdim) + Particle coordinates, non-dimensional (the same frame as + ``integration_points``). + cells : ndarray, shape (Np,) + Local owning cell per particle; a negative entry means "not on this + rank's mesh" and is left unassigned. + chunk_bytes : int + Working-set bound for the distance evaluation. The gather is + ``(chunk, Nq, cdim)``, so this caps memory rather than particle count. + + Returns + ------- + flat : ndarray, shape (Np,), int64 + Index into the flattened ``(ncells * Nq,)`` point ordering, or ``-1`` + for a particle whose cell is negative. + """ + integration_points = np.asarray(integration_points, dtype=float) + coords = np.asarray(coords, dtype=float) + cells = np.asarray(cells).reshape(-1) + + ncells, Nq, cdim = integration_points.shape + flat = np.full(coords.shape[0], -1, dtype=np.int64) + live = np.flatnonzero((cells >= 0) & (cells < ncells)) + if live.size == 0 or ncells == 0: + return flat + + # (chunk, Nq, cdim) doubles per pass + step = max(1, int(chunk_bytes // max(Nq * cdim * 8, 1))) + for start in range(0, live.size, step): + sel = live[start : start + step] + cc = cells[sel] + d2 = ((integration_points[cc] - coords[sel][:, None, :]) ** 2).sum(axis=-1) + flat[sel] = cc * Nq + d2.argmin(axis=1) + + return flat + + +def share_average(flat, values, npoints): + """Mean of each integration point's share, and how many particles it holds. + + Parameters + ---------- + flat : ndarray, shape (Np,) + The assignment from :func:`share_assignment`; ``-1`` entries are + ignored. + values : ndarray, shape (Np, ncomp) + Particle values. + npoints : int + Total number of integration points on this rank (``ncells * Nq``). + + Returns + ------- + means : ndarray, shape (npoints, ncomp) + The share mean; zero where a point holds no particles. + counts : ndarray, shape (npoints,) + Particles per integration point. A zero here is the caller's problem + to fill — a rule point in an empty or badly-sampled cell. + """ + flat = np.asarray(flat).reshape(-1) + values = np.asarray(values, dtype=float).reshape(flat.shape[0], -1) + ncomp = values.shape[1] + + ok = flat >= 0 + idx = flat[ok] + counts = np.bincount(idx, minlength=npoints)[:npoints] + + means = np.zeros((npoints, ncomp), dtype=float) + if idx.size: + vals = values[ok] + for k in range(ncomp): + means[:, k] = np.bincount(idx, weights=vals[:, k], minlength=npoints)[:npoints] + nz = counts > 0 + means[nz] /= counts[nz, None] + + return means, counts diff --git a/tests/test_0071_index_swarm_proxy_location.py b/tests/test_0071_index_swarm_proxy_location.py new file mode 100644 index 000000000..ff22d35c9 --- /dev/null +++ b/tests/test_0071_index_swarm_proxy_location.py @@ -0,0 +1,397 @@ +"""A material index read at the integration points, or as a fraction per cell. + +``IndexSwarmVariable(..., proxy_location=...)`` chooses where the level sets +live. At the integration points each one is exactly 0 or 1 (the material of +the nearest particle), which is what keeps a layered viscosity exact; per +cell they are a polynomial fraction, clamped and renormalised so the masks +still sum to one. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _layered(tag, location, cell_size=0.1, fill=3, degree=1, h=0.5): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=cell_size, qdegree=2, regular=True) + swarm = uw.swarm.Swarm(mesh) + mat = uw.swarm.IndexSwarmVariable(tag, swarm, indices=2, proxy_degree=degree, + proxy_location=location) + swarm.populate(fill_param=fill) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + mat.data[:, 0] = (X[:, 1] > h).astype(int) + return mesh, swarm, mat + + +@pytest.mark.parametrize("location", ["nodes", "integration_points", "cells"]) +def test_masks_are_bounded_and_sum_to_one(location): + mesh, swarm, mat = _layered(f"P{location[:4]}", location) + for var in mat._meshLevelSetVars: + d = np.asarray(var.data[:, 0]) + assert d.min() >= -1e-12 and d.max() <= 1 + 1e-12, (d.min(), d.max()) + total = sum(np.asarray(v.data[:, 0]) for v in mat._meshLevelSetVars) + assert np.allclose(total, 1.0, atol=1e-12) + # and through the weak form + assert abs(uw.maths.Integral(mesh, mat.sym[0] + mat.sym[1]).evaluate() - 1.0) < 1e-10 + + +def test_integration_point_masks_are_exactly_zero_or_one(): + """Each integration point takes the material of its nearest particle: a + label, not an average. That is what a nodal level set cannot do.""" + mesh, swarm, mat = _layered("Q", "integration_points") + for var in mat._meshLevelSetVars: + d = np.asarray(var.data[:, 0]) + assert set(np.unique(d).tolist()) <= {0.0, 1.0}, np.unique(d) + assert mat._meshLevelSetVars[0].is_integration_point + # the nodal proxy, by contrast, averages across the interface + _, _, nodal = _layered("N", "nodes") + d = np.asarray(nodal._meshLevelSetVars[0].data[:, 0]) + assert ((d > 1e-6) & (d < 1 - 1e-6)).any() + + +def test_layered_couette_is_exact_at_the_integration_points(): + """Interface on mesh edges: the exact velocity is in the P2 space, so the + only error is the material representation.""" + eta_top, h = 1.0e3, 0.5 + results = {} + for location in ("nodes", "integration_points"): + mesh, swarm, mat = _layered(f"C{location[:4]}", location, h=h) + v = uw.discretisation.MeshVariable(f"v{location[:4]}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{location[:4]}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = mat.createMask([1.0, eta_top]) + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1e-8 + stokes.solve() + A = 1.0 / (h + (1.0 - h) / eta_top) + Xv = np.asarray(v.coords) + exact = np.where(Xv[:, 1] < h, A * Xv[:, 1], A * h + A / eta_top * (Xv[:, 1] - h)) + results[location] = np.abs(np.asarray(v.data[:, 0]) - exact).max() + assert results["integration_points"] < 1e-5, results + assert results["nodes"] > 1e-2, results # the control + + +def test_repopulation_keeps_a_material_index_a_label(): + """A new particle takes its nearest neighbour's index whole: averaging two + labels does not give a label.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.15, qdegree=2) + swarm = uw.swarm.Swarm(mesh) + mat = uw.swarm.IndexSwarmVariable("R", swarm, indices=3, proxy_location="integration_points") + scalar = uw.swarm.SwarmVariable("Rs", swarm, 1) + swarm.populate(fill_param=3) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + mat.data[:, 0] = np.digitize(X[:, 1], [0.33, 0.66]) + scalar.data[:, 0] = X[:, 0] + for i in np.sort(np.nonzero(X[:, 0] < 0.4)[0])[::-1]: + swarm.dm.removePointAtIndex(int(i)) + swarm._invalidate_canonical_data() + added, _ = swarm.repopulate(order=1) + from mpi4py import MPI + assert uw.mpi.comm.allreduce(added, op=MPI.SUM) > 0 # added is rank-local + Xn = np.asarray(swarm._particle_coordinates.data) + idx = np.asarray(mat.data[:, 0]) + assert set(np.unique(idx).tolist()) <= {0, 1, 2}, np.unique(idx) + # right material for all but the particles nearest an interface + if idx.shape[0] > 0: + assert np.mean(idx == np.digitize(Xn[:, 1], [0.33, 0.66])) > 0.9 + # a float field is still reconstructed, not copied from one neighbour + assert np.abs(np.asarray(scalar.data[:, 0]) - Xn[:, 0]).max() < 1e-10 + + +def test_population_control_keeps_an_extending_box_sampled(): + """Pure shear with outflow sides and inflow top and bottom: without + population control the inflow cells empty.""" + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5), + cellSize=0.15, qdegree=2) + mesh.return_coords_to_bounds = None + x, y = mesh.X + V = sympy.Matrix([[x, -y]]) + c0, c1 = mesh.dm.getHeightStratum(0) + + def empty_cells(swarm): + P = np.asarray(swarm._particle_coordinates.data) + cells = np.asarray(mesh._robust_owning_cells(P)) + return int((np.bincount(cells[cells >= 0], minlength=c1 - c0) == 0).sum()) + + counts = {} + for control in (False, True): + swarm = uw.swarm.Swarm(mesh) + swarm.populate(fill_param=3) + if control: + swarm.population_control = dict(min_per_cell=6) + for _ in range(10): + swarm.advection(V, 0.1, order=2) + counts[control] = empty_cells(swarm) + from mpi4py import MPI + assert uw.mpi.comm.allreduce(counts[True], op=MPI.SUM) == 0, counts + assert uw.mpi.comm.allreduce(counts[False], op=MPI.SUM) > 0, counts + + +def test_nearest_accepts_a_bare_name_and_a_single_particle_rank_is_not_starved(): + """`nearest="F"` names one variable, not five characters; and one particle + is enough for the nearest-particle mapping (only a rank with none is).""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2) + swarm = uw.swarm.Swarm(mesh) + flag = uw.swarm.SwarmVariable("Flag", swarm, 1) + swarm.populate(fill_param=3) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + flag.data[:, 0] = (X[:, 0] > 0.5).astype(float) + for i in np.sort(np.nonzero(X[:, 0] < 0.4)[0])[::-1]: + swarm.dm.removePointAtIndex(int(i)) + swarm._invalidate_canonical_data() + swarm.repopulate(nearest="Flag") + vals = np.asarray(flag.data[:, 0]) + if vals.shape[0] > 0: # a float field, kept as a label + assert set(np.unique(vals).tolist()) <= {0.0, 1.0}, np.unique(vals) + + # one particle: the nearest-particle mapping is still well defined + mesh2 = uw.meshing.UnstructuredSimplexBox(cellSize=0.35, qdegree=2) + swarm2 = uw.swarm.Swarm(mesh2) + mat = uw.swarm.IndexSwarmVariable("One", swarm2, indices=2, + proxy_location="integration_points") + swarm2.populate(fill_param=2) + keep = 1 if swarm2.local_size > 0 else 0 + for i in range(swarm2.local_size - 1, keep - 1, -1): + swarm2.dm.removePointAtIndex(int(i)) + swarm2._invalidate_canonical_data() + if swarm2.local_size == 1: + with uw.synchronised_array_update(): + mat.data[:, 0] = 1 + mat._proxy_stale = True + mat._update_proxy_if_stale() + d = np.asarray(mat._meshLevelSetVars[1].data[:, 0]) + if swarm2.local_size == 1: + assert np.allclose(d, 1.0), d[:5] # every point takes that particle + + +@pytest.mark.parametrize("location", ["nodes", "integration_points", "cells"]) +def test_the_symbol_participates_in_expressions(location): + """First class means the symbol goes where a mesh variable's symbol goes: + arithmetic, integrals, evaluation, projection, a solver term. Only a + gradient of the integration-point form is refused, because the element + that holds a value at each integration point has no gradient.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2) + x, y = mesh.X + tag = location[:4] + T = uw.discretisation.MeshVariable(f"T{tag}", mesh, 1, degree=2) + T.data[:, 0] = np.asarray(T.coords)[:, 0] + swarm = uw.swarm.Swarm(mesh) + mat = uw.swarm.IndexSwarmVariable(f"E{tag}", swarm, indices=2, proxy_location=location) + fld = uw.swarm.SwarmVariable(f"F{tag}", swarm, 1, proxy_location=location) + swarm.populate(fill_param=3) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + mat.data[:, 0] = (X[:, 1] > 0.5).astype(int) + fld.data[:, 0] = X[:, 0] ** 2 + + # arithmetic with a mesh variable and with sympy, under an integral + assert np.isfinite(uw.maths.Integral( + mesh, mat.sym[1] * T.sym[0] + sympy.sin(x) * fld.sym[0]).evaluate()) + # evaluation at arbitrary points + got = uw.function.evaluate(mat.sym[1] + fld.sym[0], np.array([[0.31, 0.42], [0.7, 0.8]])) + assert np.isfinite(np.asarray(got)).all() + # projection of a material-weighted property + P = uw.discretisation.MeshVariable(f"P{tag}", mesh, 1, degree=1) + proj = uw.systems.solvers.SNES_Projection(mesh, P) + proj.uw_function = mat.createMask([1.0, 5.0]) + proj.solve() + # both material values are recovered; a continuous projection of a sharp + # mask overshoots at the interface, which is the projection's business + Pv = np.asarray(P.data) + assert Pv.max() > 4.5 and Pv.min() < 1.5, (Pv.min(), Pv.max()) + # a solver term: viscosity and body force at once + v = uw.discretisation.MeshVariable(f"v{tag}", mesh, mesh.dim, degree=2) + q = uw.discretisation.MeshVariable(f"q{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=q) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = mat.createMask([1.0, 10.0]) + stokes.bodyforce = sympy.Matrix([[0, -mat.sym[1] * (1 + T.sym[0])]]) + for b in ("Top", "Bottom"): + stokes.add_dirichlet_bc((0.0, 0.0), b) + for b in ("Left", "Right"): + stokes.add_dirichlet_bc((sympy.oo, 0.0), b) + stokes.solve() + assert np.abs(np.asarray(v.data)).max() > 0 + + # the gradient: available except at the integration points + grad = uw.maths.Integral(mesh, fld.sym[0].diff(x) ** 2) + if location == "integration_points": + with pytest.raises(RuntimeError, match="integration-point"): + grad.evaluate() + else: + assert grad.evaluate() > 0 + + +def test_the_gradient_is_available_from_the_cell_fit(): + """The integration-point form has no gradient, but the same particle data + fitted per cell does, and it is the most accurate of the routes: the fit is + local and exact for polynomials up to its degree.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2) + x, y = mesh.X + all_pts = np.array([[0.31, 0.42], [0.62, 0.25], [0.5, 0.5]]) + # uw.function.evaluate is rank-local: only ask for points this rank owns. + pts = all_pts[np.asarray(mesh._robust_owning_cells(all_pts)) >= 0] + if pts.shape[0] == 0: + pytest.skip("no sample point on this rank") + exact = 2 * pts[:, 0] # d/dx of x^2 + 2y + + def gradient_error(tag, location, degree): + swarm = uw.swarm.Swarm(mesh) + var = uw.swarm.SwarmVariable(tag, swarm, 1, proxy_location=location, + proxy_degree=degree) + swarm.populate(fill_param=3) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + var.data[:, 0] = X[:, 0] ** 2 + 2 * X[:, 1] + got = np.asarray(uw.function.evaluate(var.sym[0].diff(x), pts)).reshape(-1) + return np.abs(got - exact).max() + + cells2 = gradient_error("Gc2", "cells", 2) + nodes2 = gradient_error("Gn2", "nodes", 2) + assert cells2 < 1e-5, cells2 # exact for a quadratic + assert cells2 < nodes2 / 10, (cells2, nodes2) + + # The integration-point form: refused in a WEAK FORM, recovered by + # evaluate (which fits the values per cell for the one call). + swarm = uw.swarm.Swarm(mesh) + ip = uw.swarm.SwarmVariable("Gq", swarm, 1, proxy_location="integration_points") + swarm.populate(fill_param=3) + Xp = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + ip.data[:, 0] = Xp[:, 0] ** 2 + 2 * Xp[:, 1] + with pytest.raises(RuntimeError, match="proxy_location='cells'"): + uw.maths.Integral(mesh, ip.sym[0].diff(x)).evaluate() # the weak form + recovered = np.asarray(uw.function.evaluate(ip.sym[0].diff(x), pts)).reshape(-1) + assert np.abs(recovered - exact).max() < 5e-3, recovered # the query answers + # and the direct route is the sharper of the two + assert cells2 < np.abs(recovered - exact).max() + + +def _layered_solve(tag, viscosity, swarm_setup, eta_top=1.0e3, h=0.5): + """The layered Couette solve, given a viscosity expression.""" + mesh, swarm, viscosity = swarm_setup + v = uw.discretisation.MeshVariable(f"v{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = viscosity + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1e-8 + stokes.solve() + A = 1.0 / (h + (1.0 - h) / eta_top) + Xv = np.asarray(v.coords) + exact = np.where(Xv[:, 1] < h, A * Xv[:, 1], A * h + A / eta_top * (Xv[:, 1] - h)) + return np.sqrt(np.mean((np.asarray(v.data[:, 0]) - exact) ** 2)) + + +def test_a_property_on_particles_is_not_a_material(): + """The hand-rolled material: carry the viscosity itself on the particles + and sample it. It is refused, and the refusal names the thing to use + instead — because the only reconstruction on offer overshoots a jump into + a NEGATIVE viscosity, and because a solver needs the constitutive law, not + a sampled answer to it.""" + eta_top, h = 1.0e3, 0.5 + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2, regular=True) + swarm = uw.swarm.Swarm(mesh) + + with pytest.raises(ValueError, match="IndexSwarmVariable"): + uw.swarm.SwarmVariable("etaN", swarm, 1, + proxy_location="integration_points", + proxy_sampling="nearest") + + # why: the reconstruction that IS on offer cannot hold a jump + eta = uw.swarm.SwarmVariable("etaR", swarm, 1, proxy_location="integration_points") + swarm.populate(fill_param=3) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + eta.data[:, 0] = np.where(X[:, 1] > h, eta_top, 1.0) + eta._update_proxy_if_stale() + assert np.asarray(eta._meshVar.data[:, 0]).min() < 0.0 # a negative viscosity + + # the supported route is exact on the same problem + mesh, swarm, mat = _layered("Mok", "integration_points", h=h) + err = _layered_solve("ok", None, (mesh, swarm, mat.createMask([1.0, eta_top])), + eta_top, h) + assert err < 1e-5, err + + +def test_the_share_uses_every_particle_and_is_linear_exact(): + """proxy_sampling="share" assigns every particle to one integration point + of its OWN cell, and each point reads the mean of its share. Nothing is + discarded (unlike a nearest-particle read) and nothing crosses a cell + boundary.""" + from underworld3.utilities.particle_share import share_assignment + + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2, regular=True) + swarm = uw.swarm.Swarm(mesh) + f = uw.swarm.SwarmVariable("fsh", swarm, 1, proxy_location="integration_points", + proxy_sampling="share") + swarm.populate(fill_param=4) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + f.data[:, 0] = 2.0 * X[:, 0] - 3.0 * X[:, 1] + 1.0 + f._update_proxy_if_stale() + + ipc = np.asarray(f._meshVar.integration_points) + flat = share_assignment(ipc, X, swarm._owning_cells()) + assert (flat >= 0).sum() == X.shape[0] # every particle used + cells = swarm._owning_cells() + assert np.array_equal(flat // ipc.shape[1], cells) # and only in its own cell + + # the share mean of a linear field is close to its value at the point + pts = ipc.reshape(-1, mesh.dim) + exact = 2.0 * pts[:, 0] - 3.0 * pts[:, 1] + 1.0 + got = np.asarray(f._meshVar.data[:, 0]) + assert np.sqrt(np.mean((got - exact) ** 2)) < 0.1 * mesh.get_min_radius() + + +def test_share_level_sets_are_fractional_but_still_a_partition(): + """The index share gives material FRACTIONS where a cell is crossed, which + nearest sampling never does — and the masks still sum to one, so + createMask stays a weighted mean.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.15, qdegree=2, regular=False) + fractional = {} + for tag, sampling in (("Snear", "nearest"), ("Sshar", "share")): + swarm = uw.swarm.Swarm(mesh) + mat = uw.swarm.IndexSwarmVariable(tag, swarm, indices=2, + proxy_location="integration_points", + proxy_sampling=sampling) + swarm.populate(fill_param=4) + X = np.asarray(swarm._particle_coordinates.data) + with uw.synchronised_array_update(): + mat.data[:, 0] = (X[:, 1] > 0.53).astype(int) # cuts through cells + U = np.column_stack([np.asarray(v.data[:, 0]) for v in mat._meshLevelSetVars]) + assert np.allclose(U.sum(axis=1), 1.0, atol=1e-12), U.sum(axis=1) + assert U.min() >= -1e-12 and U.max() <= 1 + 1e-12 + fractional[sampling] = int(((U > 1e-9) & (U < 1 - 1e-9)).any(axis=1).sum()) + + assert fractional["nearest"] == 0, fractional + assert fractional["share"] > 0, fractional + + +def test_an_unavailable_sampling_combination_raises(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=2, regular=True) + swarm = uw.swarm.Swarm(mesh) + with pytest.raises(ValueError, match="integration_points"): + uw.swarm.SwarmVariable("bad1", swarm, 1, proxy_location="nodes", + proxy_sampling="share") + with pytest.raises(ValueError, match="integration_points"): + uw.swarm.SwarmVariable("bad2", swarm, 1, proxy_location="cells", + proxy_sampling="share") + with pytest.raises(ValueError, match="integration_points"): + uw.swarm.IndexSwarmVariable("bad3", swarm, indices=2, + proxy_location="nodes", proxy_sampling="share") diff --git a/tests/test_0072_material_swarm.py b/tests/test_0072_material_swarm.py new file mode 100644 index 000000000..8c3255991 --- /dev/null +++ b/tests/test_0072_material_swarm.py @@ -0,0 +1,254 @@ +"""Materials named on a swarm, with no level sets in sight. + +``MaterialSwarm`` is the user-facing material interface: declare the +materials and their properties, say where each one is, hand the swarm to the +solver. The index variable, the level sets and the blend are machinery. +""" + +import warnings + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +ETA_TOP, H = 1.0e3, 0.5 + + +def _box(cell_size=0.1, regular=True): + return uw.meshing.UnstructuredSimplexBox( + cellSize=cell_size, qdegree=2, regular=regular + ) + + +def _layered(mesh, tag, fill=3, **kwargs): + materials = uw.swarm.MaterialSwarm(mesh, fill_param=fill, name=tag, **kwargs) + materials.add("lower", shear_viscosity_0=1.0, density=3300) + materials.add("upper", shear_viscosity_0=ETA_TOP, density=3400) + materials["upper"] = mesh.X[1] > H + return materials + + +def _couette(mesh, tag, materials): + v = uw.discretisation.MeshVariable(f"v{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.materials = materials + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1e-8 + stokes.solve() + A = 1.0 / (H + (1.0 - H) / ETA_TOP) + X = np.asarray(v.coords) + exact = np.where(X[:, 1] < H, A * X[:, 1], A * H + A / ETA_TOP * (X[:, 1] - H)) + return np.sqrt(np.mean((np.asarray(v.data[:, 0]) - exact) ** 2)) + + +def test_a_model_names_materials_and_never_writes_a_mask(): + """The whole interface: two named materials, a region, one assignment.""" + mesh = _box() + materials = _layered(mesh, "Ma") + + assert _couette(mesh, "a", materials) < 1e-5 + # and the blend is exactly what the assembler integrated + assert abs( + float(uw.maths.Integral(mesh, materials.shear_viscosity_0).evaluate()) + - (1.0 * H + ETA_TOP * (1.0 - H)) + ) < 1e-8 + assert abs( + float(uw.maths.Integral(mesh, materials["upper"].mask).evaluate()) - (1 - H) + ) < 1e-8 + + +def test_a_property_may_be_a_law_not_a_number(): + """The reason the partition of unity exists: a material whose viscosity is + an expression cannot be stored as a value at the integration points, but it + blends symbolically.""" + mesh = _box(cell_size=0.2) + T = uw.discretisation.MeshVariable("Tlaw", mesh, 1, degree=1) + with uw.synchronised_array_update(): + T.data[:, 0] = np.asarray(T.coords)[:, 1] + + materials = uw.swarm.MaterialSwarm(mesh, fill_param=3, name="Mlaw") + materials.add("a", shear_viscosity_0=1.0) + materials.add("b", shear_viscosity_0=10 * sympy.exp(-T.sym[0])) + materials["b"] = mesh.X[1] > H + + blended = materials.shear_viscosity_0 + assert T.sym[0] in blended.atoms(type(T.sym[0])) # the law survived + # exact: 0.5 * 1 + integral over the top half of 10 exp(-y) + exact = 0.5 + 10.0 * (np.exp(-0.5) - np.exp(-1.0)) + assert abs(float(uw.maths.Integral(mesh, blended).evaluate()) - exact) < 2e-3 + + +@pytest.mark.parametrize("kind", ["symbolic", "combined", "array", "callable"]) +def test_regions_can_be_written_four_ways(kind): + mesh = _box(cell_size=0.2) + materials = uw.swarm.MaterialSwarm(mesh, fill_param=3, name=f"Mr{kind[:4]}") + materials.add("bg", shear_viscosity_0=1.0) + materials.add("box", shear_viscosity_0=2.0) + materials.populate() + + x, y = mesh.X + if kind == "symbolic": + region, area = y > H, 0.5 + elif kind == "combined": + region, area = (y > H) & (x < 0.5), 0.25 + elif kind == "array": + coords = np.asarray(materials._particle_coordinates.data) + region, area = coords[:, 1] > H, 0.5 + else: + region, area = (lambda c: c[:, 1] > H), 0.5 + + materials["box"] = region + assert abs( + float(uw.maths.Integral(mesh, materials["box"].mask).evaluate()) - area + ) < 0.02 + + +def test_only_recognised_properties_reach_the_constitutive_model(): + """shear_viscosity_0 is pushed; density is not a viscous-model parameter + and stays a symbol for the model script to use in the body force.""" + mesh = _box(cell_size=0.2) + materials = _layered(mesh, "Mp", fill=2) + v = uw.discretisation.MeshVariable("vp", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("pp", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.materials = materials + + pushed = stokes.constitutive_model.Parameters.shear_viscosity_0 + assert abs( + float(uw.maths.Integral(mesh, pushed.sym if hasattr(pushed, "sym") else pushed + ).evaluate()) - (1.0 * H + ETA_TOP * (1.0 - H)) + ) < 1e-8 + + assert not hasattr(stokes.constitutive_model.Parameters, "density") + assert abs( + float(uw.maths.Integral(mesh, materials.density).evaluate()) + - (3300 * H + 3400 * (1 - H)) + ) < 1e-8 + + +def test_a_property_nothing_uses_is_reported(): + """A misspelled viscosity is silently the default one, so say so.""" + mesh = _box(cell_size=0.25) + materials = uw.swarm.MaterialSwarm(mesh, fill_param=2, name="Mw") + materials.add("a", shear_viscosity_0=1.0, viscocity=99.0) + materials.add("b", shear_viscosity_0=2.0, viscocity=99.0) + materials["b"] = mesh.X[1] > H + + v = uw.discretisation.MeshVariable("vw", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("pw", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.materials = materials + stokes.add_dirichlet_bc((0.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + stokes._build() + assert any("viscocity" in str(w.message) for w in caught), [ + str(w.message) for w in caught + ] + + +def test_a_property_missing_from_one_material_is_an_error(): + """Two surfaces, two exception types on purpose. ``blend`` raises KeyError + naming the materials that are missing the property; attribute access + converts that to AttributeError, because hasattr() swallows AttributeError + and nothing else -- a KeyError escaping __getattr__ breaks hasattr() and + getattr(o, n, default) for every caller.""" + mesh = _box(cell_size=0.25) + materials = uw.swarm.MaterialSwarm(mesh, fill_param=2, name="Mm") + materials.add("a", shear_viscosity_0=1.0, density=3300) + materials.add("b", shear_viscosity_0=2.0) # no density + + with pytest.raises(KeyError, match="density"): + materials.blend("density") + + with pytest.raises(AttributeError, match="density"): + materials.density + + assert hasattr(materials, "density") is False + + +def test_materials_must_be_declared_before_the_swarm_is_populated(): + mesh = _box(cell_size=0.25) + materials = uw.swarm.MaterialSwarm(mesh, fill_param=2, name="Md") + materials.add("a", shear_viscosity_0=1.0) + materials.populate() + with pytest.raises(RuntimeError, match="already populated"): + materials.add("b", shear_viscosity_0=2.0) + + +def test_a_property_changed_after_attachment_reaches_the_solver(): + mesh = _box(cell_size=0.2) + materials = _layered(mesh, "Mc", fill=2) + v = uw.discretisation.MeshVariable("vc", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("pc", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.materials = materials + + materials["upper"].set(shear_viscosity_0=7.0) + pushed = stokes.constitutive_model.Parameters.shear_viscosity_0 + expr = pushed.sym if hasattr(pushed, "sym") else pushed + assert abs( + float(uw.maths.Integral(mesh, expr).evaluate()) - (1.0 * H + 7.0 * (1 - H)) + ) < 1e-8 + + +def test_the_mixing_rule_only_matters_for_fractional_masks(): + """With the default sampling exactly one mask is 1, so arithmetic and + harmonic blending give the same integral; with the share they do not.""" + mesh = _box(cell_size=0.15, regular=False) + sharp = _layered(mesh, "Msh", fill=4) + both = [ + float(uw.maths.Integral(mesh, sharp.blend("shear_viscosity_0", m)).evaluate()) + for m in ("arithmetic", "harmonic") + ] + assert abs(both[0] - both[1]) < 1e-6 * abs(both[0]) + + shared = _layered(mesh, "Msr", fill=4, proxy_sampling="share") + both = [ + float(uw.maths.Integral(mesh, shared.blend("shear_viscosity_0", m)).evaluate()) + for m in ("arithmetic", "harmonic") + ] + assert abs(both[0] - both[1]) > 1.0 + + +def test_it_is_a_swarm(): + """Advection, population control and extra state variables all still work: + a MaterialSwarm is a Swarm.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5), cellSize=0.15, qdegree=2 + ) + materials = uw.swarm.MaterialSwarm(mesh, fill_param=3, name="Mad") + materials.add("bg", shear_viscosity_0=1.0) + materials.add("layer", shear_viscosity_0=10.0) + strain = uw.swarm.SwarmVariable( + "eps", materials, 1, proxy_location="integration_points", + proxy_sampling="share", + ) + materials["layer"] = sympy.Abs(mesh.X[1]) < 0.2 + materials.population_control = dict(min_per_cell=6) + + area0 = float(uw.maths.Integral(mesh, materials["layer"].mask).evaluate()) + x, y = mesh.X + for _ in range(6): + materials.advection(sympy.Matrix([[x, -y]]), 0.05, order=2) + + labels = np.unique(np.asarray(materials.index.data).reshape(-1)) + assert set(labels.tolist()) <= {0, 1} # still labels + area1 = float(uw.maths.Integral(mesh, materials["layer"].mask).evaluate()) + assert area1 < area0 # the layer thinned + assert area1 > 0.3 * area0 # but did not fall apart + assert np.asarray(strain.data).shape[0] == materials.local_size diff --git a/tests/test_0073_material_regions.py b/tests/test_0073_material_regions.py new file mode 100644 index 000000000..6881c91a6 --- /dev/null +++ b/tests/test_0073_material_regions.py @@ -0,0 +1,312 @@ +"""Materials tied to the mesh rather than to particles. + +``MaterialRegions`` is the other distribution: same declaration, same handoff +to a solver, but the level sets come from mesh labels or from geometry, so +there are no particles to populate or advect. A registry can be shared between +the two. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +ETA_TOP, H = 1.0e3, 0.5 + + +def _box(cell_size=0.1, regular=True): + return uw.meshing.UnstructuredSimplexBox( + cellSize=cell_size, qdegree=2, regular=regular + ) + + +def _couette(mesh, tag, materials, h=H): + v = uw.discretisation.MeshVariable(f"v{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.materials = materials + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1e-8 + stokes.solve() + A = 1.0 / (h + (1.0 - h) / ETA_TOP) + X = np.asarray(v.coords) + exact = np.where(X[:, 1] < h, A * X[:, 1], A * h + A / ETA_TOP * (X[:, 1] - h)) + return np.sqrt(np.mean((np.asarray(v.data[:, 0]) - exact) ** 2)) + + +def test_regions_need_no_particles(): + """The same model as the swarm case, with the geometry on the mesh.""" + mesh = _box() + materials = uw.MaterialRegions(mesh, name="Rc") + materials.add("lower", shear_viscosity_0=1.0, density=3300) + materials.add("upper", shear_viscosity_0=ETA_TOP, density=3400) + materials["upper"] = mesh.X[1] > H + + assert _couette(mesh, "rc", materials) < 1e-5 + assert abs( + float(uw.maths.Integral(mesh, materials.shear_viscosity_0).evaluate()) + - (1.0 * H + ETA_TOP * (1.0 - H)) + ) < 1e-8 + + +def test_region_masks_are_a_sharp_partition_of_unity(): + mesh = _box(cell_size=0.15, regular=False) + materials = uw.MaterialRegions(mesh, name="Rp") + materials.add("a", shear_viscosity_0=1.0) + materials.add("b", shear_viscosity_0=2.0) + materials.add("c", shear_viscosity_0=3.0) + materials["b"] = mesh.X[1] > 0.33 + materials["c"] = mesh.X[1] > 0.66 + + materials._ensure_built() + U = np.column_stack( + [np.asarray(v.data[:, 0]) for v in materials._level_sets()] + ) + assert np.allclose(U.sum(axis=1), 1.0, atol=1e-12), U.sum(axis=1) + assert set(np.unique(U).tolist()) <= {0.0, 1.0} # no fractions, ever + assert abs(float(uw.maths.Integral(mesh, materials["c"].mask).evaluate()) + - (1 - 0.66)) < 0.02 + + +def test_a_region_can_come_from_a_mesh_label(): + """The gmsh-physical-group route, exercised with a label built by hand.""" + mesh = _box(cell_size=0.2) + dm = mesh.dm + c0, c1 = dm.getHeightStratum(0) + centroids = np.array([dm.computeCellGeometryFVM(c)[1] for c in range(c0, c1)]) + + dm.createLabel("TopHalf") + upper_cells = [c for c in range(c0, c1) if centroids[c - c0][1] > H] + for c in upper_cells: + dm.setLabelValue("TopHalf", c, 1) + + materials = uw.MaterialRegions(mesh, name="Rl") + materials.add("lower", shear_viscosity_0=1.0) + materials.add("upper", shear_viscosity_0=ETA_TOP) + materials["upper"] = "TopHalf" + + area = float(uw.maths.Integral(mesh, materials["upper"].mask).evaluate()) + assert abs(area - (1 - H)) < 0.05, area + + with pytest.raises(KeyError, match="no label"): + materials["upper"] = "NotALabel" + + +def test_a_registry_is_shared_between_distributions(): + """Definitions live in the registry, so two distributions can use one set + of materials — and a property change reaches both.""" + mesh = _box(cell_size=0.2) + rocks = uw.MaterialRegistry() + rocks.add("lower", shear_viscosity_0=1.0) + rocks.add("upper", shear_viscosity_0=ETA_TOP) + + regions = uw.MaterialRegions(mesh, registry=rocks, name="Rsh") + regions["upper"] = mesh.X[1] > H + swarm = uw.swarm.MaterialSwarm(mesh, registry=rocks, fill_param=3, name="Ssh") + swarm["upper"] = mesh.X[1] > H + + exact = 1.0 * H + ETA_TOP * (1.0 - H) + for distribution in (regions, swarm): + got = float( + uw.maths.Integral(mesh, distribution.shear_viscosity_0).evaluate() + ) + assert abs(got - exact) < 1e-8, (distribution, got) + + rocks["upper"].set(shear_viscosity_0=7.0) + exact = 1.0 * H + 7.0 * (1.0 - H) + for distribution in (regions, swarm): + got = float( + uw.maths.Integral(mesh, distribution.shear_viscosity_0).evaluate() + ) + assert abs(got - exact) < 1e-8, (distribution, got) + + +def test_a_registry_stands_alone_and_round_trips(): + """No mesh, no swarm: a material library is just definitions.""" + rocks = uw.MaterialRegistry() + rocks.add("mantle", shear_viscosity_0=1.0, density=3300, + description="upper mantle", reference="T&S (2014)") + rocks.add("slab", shear_viscosity_0=1.0e3, density=3400) + + assert rocks.list_materials() == ["mantle", "slab"] + assert rocks["mantle"].index == 0 and rocks["slab"].index == 1 + assert rocks.declared_properties() == {"shear_viscosity_0", "density"} + + config = rocks.export_config() + clone = uw.MaterialRegistry().import_config(config) + assert clone.list_materials() == rocks.list_materials() + assert clone["mantle"].description == "upper mantle" + assert float(clone["slab"].get_property("shear_viscosity_0")) == 1.0e3 + + # the enum spelling reaches the same place as the keyword + rocks["slab"].set_property(uw.MaterialProperty.DENSITY, 3500) + assert float(rocks["slab"].get_property("density")) == 3500 + + +def test_materials_must_be_declared_before_the_level_sets_are_built(): + mesh = _box(cell_size=0.25) + materials = uw.MaterialRegions(mesh, name="Rd") + materials.add("a", shear_viscosity_0=1.0) + materials["a"] = mesh.X[1] > H + materials._ensure_built() + with pytest.raises(RuntimeError, match="already built"): + materials.add("b", shear_viscosity_0=2.0) + + +def test_a_property_can_be_a_law_for_regions_too(): + mesh = _box(cell_size=0.2) + T = uw.discretisation.MeshVariable("Trl", mesh, 1, degree=1) + with uw.synchronised_array_update(): + T.data[:, 0] = np.asarray(T.coords)[:, 1] + + materials = uw.MaterialRegions(mesh, name="Rlaw") + materials.add("a", shear_viscosity_0=1.0) + materials.add("b", shear_viscosity_0=10 * sympy.exp(-T.sym[0])) + materials["b"] = mesh.X[1] > H + + exact = 0.5 + 10.0 * (np.exp(-0.5) - np.exp(-1.0)) + got = float(uw.maths.Integral(mesh, materials.shear_viscosity_0).evaluate()) + assert abs(got - exact) < 2e-3, got + + +# --------------------------------------------------------------------------- +# Adversarial review, 2026-09-10. Every test below reproduces a defect that +# was confirmed against the previous revision. +# --------------------------------------------------------------------------- + + +def test_a_label_value_that_does_not_exist_is_refused_not_dereferenced(): + """PETSc HARD-ABORTS on getStratumIS() for a value outside a label's live + set -- no exception, no traceback, every rank gone. The API's own error + message steered users straight into it ("say which with (label, value)"), + and on a UW3 box the only valid value is 99999, so every guess crashed.""" + mesh = _box(cell_size=0.25) + materials = uw.MaterialRegions(mesh, name="Lv") + materials.add("a", shear_viscosity_0=1.0) + materials.add("b", shear_viscosity_0=2.0) + materials["b"] = ("Elements", 1) # not the live value + + with pytest.raises(ValueError, match="no stratum with value"): + materials.build() + + # the real value still selects the cells + good = uw.MaterialRegions(mesh, name="Lv2") + good.add("a", shear_viscosity_0=1.0) + good.add("b", shear_viscosity_0=2.0) + live = int(mesh.dm.getLabelIdIS("Elements").getIndices()[0]) + good["b"] = ("Elements", live) + assert abs(float(uw.maths.Integral(mesh, good.shear_viscosity_0).evaluate()) + - 2.0) < 1e-9 + + +def test_assigning_a_region_replaces_the_previous_one(): + """``m["x"] = A`` then ``m["x"] = B`` used to leave the material at + A union B -- not what ``=`` means, and a trap for anyone correcting a + region in a notebook cell.""" + mesh = _box(cell_size=0.2) + materials = uw.MaterialRegions(mesh, name="Rp2") + materials.add("bg", rho=1.0) + materials.add("x", rho=2.0) + + materials["x"] = mesh.X[1] > 0.3 + materials["x"] = mesh.X[1] > 0.6 + got = float(uw.maths.Integral(mesh, materials.rho).evaluate()) + assert abs(got - 1.4) < 0.05, f"{got} looks like the union (1.7)" + + +def test_hasattr_does_not_raise_and_does_not_build(): + """``blend`` raises KeyError for a partly-declared property, and + ``__getattr__`` used to let it out -- breaking hasattr() for every caller. + hasattr must also not trigger the (collective) level-set build.""" + mesh = _box(cell_size=0.25) + materials = uw.MaterialRegions(mesh, name="Ha") + materials.add("a", rho=1.0) + materials.add("b", rho=2.0, sigma=3.0) # only b declares sigma + + assert hasattr(materials, "sigma") is False + assert hasattr(materials, "not_a_property") is False + assert materials._level_set_vars is None, "hasattr built the level sets" + + +def test_two_distributions_cannot_share_a_name(): + """They shared their level sets silently: the second one's variables were + skipped, it was handed the first one's, and painting it changed the FIRST + one's answers. The only diagnostic was a print to stdout.""" + mesh = _box(cell_size=0.25) + first = uw.MaterialRegions(mesh, name="Same") + first.add("a", rho=1.0) + first.add("b", rho=2.0) + first["b"] = mesh.X[1] > 0.5 + first.build() + + second = uw.MaterialRegions(mesh, name="Same") + second.add("a", rho=1.0) + second.add("b", rho=2.0) + with pytest.raises(ValueError, match="already in use"): + second.build() + + +def test_the_registry_is_frozen_once_a_distribution_is_built(): + """A material's index IS which level set is its, so adding or deleting one + after allocation re-points the blend. ``registry.add`` bypassed the + distribution's own guard entirely, and the resulting IndexError was + swallowed into a warning while the solver kept a stale blend.""" + mesh = _box(cell_size=0.25) + materials = uw.MaterialRegions(mesh, name="Fz") + materials.add("a", rho=1.0) + materials.add("b", rho=2.0) + materials["b"] = mesh.X[1] > 0.5 + materials.build() + + with pytest.raises(RuntimeError, match="already in use"): + materials.registry.add("late", rho=9.0) + with pytest.raises(RuntimeError, match="already in use"): + materials.registry.delete_material("a") + + +def test_harmonic_mixing_will_not_divide_by_a_material_with_no_value(): + """1 / sum(phi_i / v_i) with v_i == 0 folds to ComplexInfinity at + declaration time and the JIT then dies with a bare C-printer traceback + naming neither the material nor the property. An air layer with + density=0 is the obvious way in.""" + mesh = _box(cell_size=0.25) + materials = uw.MaterialRegions(mesh, name="Hz") + materials.add("air", eta=0.0) + materials.add("rock", eta=1.0) + materials["rock"] = mesh.X[1] < 0.5 + materials.mixing(eta="harmonic") + + with pytest.raises(ValueError, match="divides by its value"): + materials.eta + + # a mixing rule for a property nobody declares is a typo, not a no-op + with pytest.raises(KeyError, match="no material declares"): + materials.mixing(viscocity="harmonic") + + +def test_a_quantity_means_the_same_whenever_it_was_declared(): + """Units were non-dimensionalised at add() time, so the same declaration + gave different numbers depending on whether the model's reference + quantities had been set yet -- 28 orders of magnitude apart, no warning. + A registry is meant to be writable before the model exists.""" + rocks = uw.MaterialRegistry() + rocks.add("early", shear_viscosity_0=uw.quantity(1.0e21, "Pa*s")) + + model = uw.get_default_model() + model.set_reference_quantities( + length=uw.quantity(1.0e6, "m"), + viscosity=uw.quantity(1.0e21, "Pa*s"), + time=uw.quantity(1.0e15, "s"), + ) + rocks.add("late", shear_viscosity_0=uw.quantity(1.0e21, "Pa*s")) + + early = float(rocks["early"].resolved("shear_viscosity_0")) + late = float(rocks["late"].resolved("shear_viscosity_0")) + assert early == late, (early, late) diff --git a/tests/test_0115_index_swarm_vectorized.py b/tests/test_0115_index_swarm_vectorized.py index bf88c29a3..198dee4f5 100644 --- a/tests/test_0115_index_swarm_vectorized.py +++ b/tests/test_0115_index_swarm_vectorized.py @@ -118,6 +118,7 @@ def _build_mesh_and_swarm(fill_param=5, nnn=5, radius=0.5, update_type=0, swarm = uw.swarm.Swarm(mesh) material = uw.swarm.IndexSwarmVariable( "M_test", swarm, indices=2, proxy_degree=1, proxy_continuous=True, + proxy_location="nodes", # update_type only means anything here update_type=update_type, npoints=nnn, radius=radius, npoints_bc=nnn_bc if nnn_bc is not None else 2, ind_bc=ind_bc,