Skip to content

Commit 9c8d950

Browse files
lmoresilmoresi
committed
Add publications scaffold: blog posts, GMD paper outlines, and figures
Blog posts (4 published on underworldcode.org, 1 draft): - Our Journey from Underworld2 to Underworld3 - AI and Scientific Software: What We Learned Rebuilding Underworld3 - How Underworld3 Turns SymPy into C (with companion notebook) - Mesh Variables and PETSc Vectors: Keeping Arrays in Sync (draft) - Blog post index (README.md) with 25 planned topics GMD paper scaffolds: - sympy-machinery: myst.yml, outline, references (Paper 1) - particles-surfaces: myst.yml, outline, references (Paper 2) Also includes Typst/cetz diagram source for the arrays sync flow, and the solver unification design document. Underworld development team with AI support from Claude Code
1 parent 01351a3 commit 9c8d950

22 files changed

Lines changed: 3112 additions & 0 deletions

‎docs/developer/design/SOLVER_UNIFICATION_DESIGN.md‎

Lines changed: 146 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
# Solver Unification Design
2+
3+
> Status: **Proposed** — for implementation after VEP validation is complete
4+
5+
## Goal
6+
7+
Eliminate the need for separate `VE_Stokes` and (future) `VE_NavierStokes` solver
8+
classes. The constitutive model declares what infrastructure it needs; the solver
9+
creates it lazily.
10+
11+
## Current Architecture
12+
13+
| Solver | DuDt (velocity) | DFDt (flux/stress) | Constitutive models |
14+
|--------|-----------------|-------------------|-------------------|
15+
| `Stokes` | — | — | Viscous, VP |
16+
| `VE_Stokes` | — | SemiLagrangian (stress history) | VEP |
17+
| `NavierStokes` | SemiLagrangian (velocity) | SemiLagrangian (AM flux) | Viscous, VP |
18+
| `VE_NavierStokes` | does not exist | — | — |
19+
20+
Problem: user must choose the correct solver class based on the constitutive model.
21+
Using VEP on plain Stokes silently drops stress history (now caught by barrier in PR #95).
22+
23+
## Proposed Architecture
24+
25+
Two solver classes: `Stokes` and `NavierStokes`. Each detects whether the
26+
constitutive model requires stress history and creates DFDt infrastructure lazily.
27+
No separate VE variants needed.
28+
29+
### Constitutive model contract
30+
31+
```python
32+
class Constitutive_Model:
33+
@property
34+
def requires_stress_history(self):
35+
return False # Viscous, VP
36+
37+
class ViscoElasticPlasticFlowModel(ViscousFlowModel):
38+
@property
39+
def requires_stress_history(self):
40+
return True # VEP
41+
```
42+
43+
### Solver behaviour
44+
45+
```python
46+
@constitutive_model.setter
47+
def constitutive_model(self, model):
48+
# ... existing setup ...
49+
if model.requires_stress_history and self.Unknowns.DFDt is None:
50+
self._create_stress_history_ddt(order=model.order)
51+
```
52+
53+
### Constitutive model is assigned once
54+
55+
Changing parameters (viscosity, yield stress, modulus) is fine — they flow through
56+
UWexpressions and PetscDS constants[]. Swapping the constitutive model class after
57+
the first solve is not supported (DFDt allocation, JIT structure changes).
58+
59+
### VE_Stokes becomes a backward-compat alias
60+
61+
```python
62+
class VE_Stokes(Stokes):
63+
"""Deprecated: use Stokes directly with VEP constitutive model."""
64+
def __init__(self, mesh, order=2, **kwargs):
65+
super().__init__(mesh, **kwargs)
66+
self._create_stress_history_ddt(order=order)
67+
```
68+
69+
### solve() hooks
70+
71+
```python
72+
def solve(self, timestep=None, ...):
73+
if self.Unknowns.DFDt is not None:
74+
if timestep is None:
75+
raise ValueError("timestep required for viscoelastic solve")
76+
self.constitutive_model._update_bdf_coefficients()
77+
self.DFDt.update_pre_solve(timestep, store_result=False)
78+
79+
# PETSc solve
80+
self._snes_solve(...)
81+
82+
if self.Unknowns.DFDt is not None:
83+
self._post_solve_stress_history(timestep)
84+
```
85+
86+
### tau property
87+
88+
```python
89+
@property
90+
def tau(self):
91+
if self.Unknowns.DFDt is not None:
92+
return self.DFDt.psi_star[0] # stored actual stress
93+
else:
94+
return self._lazy_tau_projection() # on-demand projection
95+
```
96+
97+
## NavierStokes Considerations
98+
99+
NS already uses DFDt for Adams-Moulton (Crank-Nicolson) stabilisation of the
100+
viscous flux. The AM scheme stores `η·∇u` at previous timesteps for flux averaging:
101+
102+
$$F^{n+1/2} = \theta \cdot F^{n+1} + (1-\theta) \cdot F^{n*}$$
103+
104+
For VE-NS, the DFDt must serve both purposes:
105+
- AM flux averaging for time integration stability
106+
- VE stress history for the Maxwell constitutive law
107+
108+
### Possible unification
109+
110+
If DFDt stores the **actual deviatoric stress** (as PR #89 implements for VE_Stokes),
111+
the AM scheme can be reformulated to read from it:
112+
113+
$$F^{n+1/2} = \theta \cdot \sigma^{n+1} + (1-\theta) \cdot \sigma^{n*}$$
114+
115+
where σ^{n*} is the advected stress from `psi_star[0]`. This unifies the two uses:
116+
the DFDt stores actual stress, and both AM stabilisation and VE constitutive law
117+
read from the same history chain.
118+
119+
### Open questions for VE-NS
120+
121+
- Is order-1 AM sufficient alongside VE stress history? The VE history already
122+
provides temporal accuracy for the elastic part; AM only needs to stabilise the
123+
viscous/advection part.
124+
- Can the AM flux and VE stress contributions simply be added symbolically in
125+
the F0/F1 expressions? SymPy handles the algebra; the JIT compiles the combined
126+
expression. This might avoid needing a separate DFDt entirely — the NS solver's
127+
existing DFDt carries the AM flux, and the VE stress history is a separate DFDt
128+
created by the constitutive model.
129+
- This needs to be driven by physics requirements (a problem that demands VE-NS),
130+
not implemented speculatively.
131+
132+
### Recommendation
133+
134+
Implement VE-NS only after:
135+
1. VEP on Stokes is fully validated (current priority)
136+
2. Solver unification for Stokes is complete and tested
137+
3. NS benchmarks are passing as a baseline
138+
4. A physics problem demands VE-NS
139+
140+
## Implementation Order
141+
142+
1. ~~PR #95: Barrier + is_viscoplastic fix~~ (done)
143+
2. Move DFDt creation from VE_Stokes.__init__ to Stokes.constitutive_model setter
144+
3. Move pre/post solve hooks from VE_Stokes.solve() to Stokes.solve()
145+
4. VE_Stokes becomes backward-compat alias
146+
5. Same pattern for NavierStokes (when physics demands it)

‎publications/blog-posts/README.md‎

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,202 @@
1+
# Blog Posts for underworldcode.org (Ghost)
2+
3+
Building blocks for the two GMD papers and the 3.0.0 release announcement.
4+
Each post is self-contained but feeds into larger pieces.
5+
6+
## Status Key
7+
- **draft**: outline or raw notes exist
8+
- **ready**: content written, needs review
9+
- **published**: live on underworldcode.org
10+
11+
---
12+
13+
## The Release Anchor
14+
15+
1. **"Underworld3 Reaches 3.0.0"**
16+
Status: draft (uw3-release-announcement.md — to be written last)
17+
Feeds into: index for all other posts
18+
Content: What changed since 0.99/JOSS. New features, fixed bugs, API changes.
19+
Links out to the explanatory posts below.
20+
21+
## Origin Story
22+
23+
2. **"Our Journey from Underworld2 to Underworld3"**
24+
Status: **published** (2026-03-23)
25+
Feeds into: Paper 1 introduction, release post
26+
Content: UW1 and UW2 as the same engine in different clothes. Why we rewrote:
27+
PETSc maturity, StGermain liability, hex-only meshes, SWIG agony, SNES unlock.
28+
What SymPy gave us. What UW3 can do that UW2 could not (curved BCs, coordinates,
29+
constitutive models without C, swappable time derivatives).
30+
31+
## Symbolic Machinery (→ Paper 1)
32+
33+
3. **"How Underworld3 Turns SymPy into C"**
34+
Status: **published** (2026-04-01)
35+
Feeds into: Paper 1 core (sections 4–6)
36+
Content: The six-stage JIT pipeline: strong form templates, automatic Jacobians,
37+
symbolic wrappers, unwrapping/non-dimensionalisation, C code generation,
38+
per-function caching, PETSc callback registration.
39+
URL: https://www.underworldcode.org/how-underworld3-turns-sympy-into-c/
40+
41+
4. **"Automatic Jacobians for Free"**
42+
Status: not started
43+
Feeds into: Paper 1 (solver integration)
44+
Content: How symbolic differentiation eliminates hand-coded Newton derivatives.
45+
Why this was impossible in UW2 (opaque C Functions). SNES integration.
46+
47+
5. **"Constitutive Models as Symbolic Objects"**
48+
Status: not started
49+
Feeds into: Paper 1 (constitutive models)
50+
Content: Composable rheology: viscous, plastic, elastic, transverse isotropic.
51+
How SymPy lets you mix and simplify material laws interactively.
52+
Contrast with UW2's fn.rheology (one C class per model).
53+
54+
6. **"Constants That Aren't Constant"**
55+
Status: not started
56+
Feeds into: Paper 1 (expression system)
57+
Content: The PetscDS constants mechanism for routing UWexpressions to C variables.
58+
Why this matters for time-dependent coefficients (BDF/AM ramp, continuation).
59+
60+
7. **"Time Derivatives You Can See"**
61+
Status: not started
62+
Feeds into: Paper 1 (time discretisation)
63+
Content: The DDt hierarchy: Lagrangian, Semi-Lagrangian, Eulerian, Symbolic.
64+
How BDF/AM schemes appear as symbolic rewrites. Contrast with UW2's fixed RK4.
65+
66+
8. **"Natural Mathematical Syntax in Scientific Python"**
67+
Status: not started
68+
Feeds into: Paper 1 (mathematical objects)
69+
Content: The MathematicalMixin: why `density * velocity` works.
70+
How operator overloading connects to SymPy's Matrix API.
71+
72+
## Units & Scaling (→ Paper 1)
73+
74+
9. **"Physical Units in Computational Geodynamics"**
75+
Status: not started
76+
Feeds into: Paper 1 (units system)
77+
Content: The Pint-based units system. String input, object storage,
78+
transparent container principle. Why units not dimensionality.
79+
80+
10. **"Non-dimensionalisation Without Tears"**
81+
Status: not started
82+
Feeds into: Paper 1 (units system)
83+
Content: Reference quantities, nondimensional unwrapping mode,
84+
solver works in scaled space while user thinks in physical units.
85+
86+
## Particles & Surfaces (→ Paper 2)
87+
88+
11. **"Finding Particles in a Parallel Mesh"**
89+
Status: not started
90+
Feeds into: Paper 2 core
91+
Content: How particle-to-processor assignment works: DMSwarm migration,
92+
spatial indexing, the handshake between PETSc's mesh decomposition and
93+
particle ownership. What happens when particles cross processor boundaries.
94+
95+
12. **"Particles That Know Calculus"**
96+
Status: not started
97+
Feeds into: Paper 2 (proxy variables)
98+
Content: Swarm variables as first-class symbolic objects. Proxy mesh
99+
variables via RBF projection. How particle data participates in weak forms.
100+
101+
13. **"Stress Has a History"**
102+
Status: not started
103+
Feeds into: Paper 2 (material history)
104+
Content: Viscoelastic stress storage and advection. Explicit stress
105+
history architecture. Order-1 vs order-2 validation.
106+
107+
14. **"Ghost Boundaries Done Right"**
108+
Status: not started
109+
Feeds into: Paper 2 (boundary integrals)
110+
Content: Internal and external boundary integrals in parallel.
111+
The ownership problem, PETSc patches, MPI-correct integration.
112+
113+
15. **"Tracking Interfaces in Large Deformation"**
114+
Status: not started
115+
Feeds into: Paper 2 (surfaces)
116+
Content: Material interfaces, level-set weighted composites,
117+
population control.
118+
119+
## Geometry & Meshing (→ Paper 2)
120+
121+
16. **"Meshing for Planetary Scale"**
122+
Status: not started
123+
Feeds into: Paper 2 (meshing)
124+
Content: Cubed-sphere construction, RegionalSphericalBox,
125+
geodetic projections. Why lat-lon grids have singularities.
126+
127+
17. **"Geographic Meshes for Regional Geodynamics"**
128+
Status: not started
129+
Feeds into: Paper 2 (meshing)
130+
Content: Ellipsoidal Earth geometry with topography,
131+
RegionalGeographicBox, boundary labelling, coordinate transforms.
132+
Real-world regional modelling with proper geodesy.
133+
134+
18. **"Symbolic Geometry: Differential Operators in Curvilinear Coordinates"**
135+
Status: not started
136+
Feeds into: Paper 1 or 2 (coordinates)
137+
Content: How gradient/divergence/curl auto-adjust for
138+
spherical/cylindrical via the coordinate system factory.
139+
140+
19. **"Adaptive Meshes That Follow the Physics"**
141+
Status: not started
142+
Feeds into: Paper 2 (adaptivity)
143+
Content: Metric-tensor-based adaptation, swarm-mediated variable
144+
transfer during remeshing.
145+
146+
## Infrastructure & Craft
147+
148+
20. **"Mesh Variables and PETSc Vectors: Keeping Arrays in Sync"**
149+
Status: not started
150+
Feeds into: Paper 1 & 2
151+
Content: The self-validating data cache, NDArray_With_Callback,
152+
why `with mesh.access()` is gone and direct `.data[...]` works now.
153+
Contrast with UW2's context managers.
154+
155+
21. **"Safe Parallelism Without MPI Expertise"**
156+
Status: not started
157+
Feeds into: Paper 2
158+
Content: uw.pprint(), selective_ranks(), why PETSc handles the
159+
hard parts. The MPICH-on-macOS scaling bug story.
160+
161+
22. **"From Notebook to Supercomputer in Zero Edits"**
162+
Status: not started
163+
Feeds into: Paper 1 (introspection)
164+
Content: How the same Jupyter notebook runs on a laptop and
165+
12,000 cores. Mathematical display, literate computing philosophy.
166+
167+
23. **"Interactive 3D Geodynamics in a Browser"**
168+
Status: not started
169+
Feeds into: neither paper directly
170+
Content: PyVista/trame integration, P2 field visualisation,
171+
server proxy detection. Contrast with UW2's LavaVu.
172+
173+
## Development Process
174+
175+
24. **"AI and Scientific Software: What We Learned Rebuilding Underworld3"**
176+
Status: **published** (2026-03-23)
177+
Feeds into: standalone (high interest)
178+
Content: Co-evolution of code and AI tools. Four phases from first contact
179+
to productivity jumps. What works, what doesn't, the units system slog.
180+
URL: https://www.underworldcode.org/ai-and-scientific-software-what-we-learned-rebuilding-underworld3/
181+
182+
25. **"Notation Gymnastics in Continuum Mechanics"**
183+
Status: not started
184+
Feeds into: Paper 1 (tensors)
185+
Content: Voigt, Mandel, full tensor. Why Mandel preserves inner
186+
products. Dimension-independent indexing.
187+
188+
---
189+
190+
## Suggested Writing Order
191+
192+
Priority based on: feeds into papers, standalone interest, dependency chain.
193+
194+
| # | Post | Why first |
195+
|---|------|-----------|
196+
| 1 | Journey from UW2 to UW3 (#2) | Sets up everything else |
197+
| 2 | SymPy into C (#3) | Paper 1 centrepiece |
198+
| 3 | Finding particles (#11) | Paper 2 centrepiece |
199+
| 4 | Units (#9) | Foundation, standalone interest |
200+
| 5 | AI strategy (#24) | Standalone, high external interest |
201+
| 6 | Arrays in sync (#20) | Practical, bridges both papers |
202+
| 7 | Release announcement (#1) | Written last, links to all others |

0 commit comments

Comments
 (0)