solid-rocket-grain-reverse-design

Reverse-design a 2D solid-propellant rocket grain from a target chamber-pressure curve, by evolving a neural network that is the grain shape.

An independent Python reimplementation of Li, Li, He & Liang, "Reverse Design of Solid Propellant Grain for a Performance-Matching Goal: Shape Optimization via Evolutionary Neural Network", Aerospace 2022, 9(10), 552 (doi:10.3390/aerospace9100552).

---

The idea

Classical grain design is forward and parametric: pick a star/wagon-wheel/dogbone template, tune its handful of dimensions, simulate the burn-back, compare the thrust curve, repeat. The template you start from bounds the shapes you can ever reach.

This method inverts that. The grain shape is not parameterised by geometry at all — it is the zero level set of a phase field φ(x,y) emitted by a tiny neural network. A genetic algorithm evolves the network's 81 weights; whatever iso-contour those weights happen to produce is the candidate grain. There is no shape template, so there is no template to be trapped by.

Four coupled subsystems, one direction of data flow:

target p_c(t) ──[M1 reverse ballistics]──→ target (l̂_b·C_T) ~ ŵ
                                                        │
Δp ─[M3 FNN]→ φ(x,y) → r,s ─[M2 φ-PEF FEM]→ W ─[M4 contour]→ design l̂_b ~ ŵ
 ↑                                                             │
 └──────────────[M6 genetic algorithm]←──[M5 COD fitness]←─────┘

	Module	Role
M1	ballistics.py	Reverse internal ballistics. Inverts equilibrium chamber pressure into the dimensionless target burn-perimeter curve. Pure quadrature.
M2	solver.py, mesh.py	φ-PEF burn-back. A nonlinear diffusion–eikonal PDE whose solution W(x,y) is the burned web — i.e. when each point ignites. Solved by hand-coded Newton on a fixed mesh. The bottleneck.
M3	fnn.py, pretrain.py	The 2→20→1 tanh network (81 params) that emits φ. Forward pass only — no backprop in the loop.
M4	objective.py	Burn perimeter = arc length of the W = w iso-contour, by marching triangles.
M5	objective.py	Fitness. Least-squares slope C_D, coefficient of determination COD, and FV = 1 − COD + P.
M6	ga.py	pymoo GA over the 81 weights. ~10⁴ φ-PEF solves per run.

Why a phase field. Burn-back normally requires re-meshing the propellant as it recedes. φ-PEF instead marks solid (φ>0) and gas (φ<0) on one fixed mesh and gives the gas a huge virtual burn rate, so the flame front crosses it instantly. Mesh and the diffusion skeleton are assembled once and reused across every GA evaluation — the single biggest reason the search is tractable.

---

Status

Built stage by stage; every stage had a hard quantitative gate that had to pass before the next was written. All gates are checked-in scripts you can re-run.

Stage	Gate	Result
0 Environment	scikit-fem Poisson < 1e-6	4.4e-16 ✅
1 Reverse ballistics	dual-thrust two-plateau ≈30 → ≈16	27.81 / 17.12 ✅
2a Fan-sector mesh	1000 elements, boundary tags partition	1000; 20/20/50/50 facets ✅
2b Diffusion skeleton	symmetric, sparse, φ-independent	‖S·1‖ = 1e-16 ✅
2c Newton + K-continuation	‖F‖ < 1e-8 in < 15 steps	2 steps, 1e-10 ✅
2d Square ∩ circle benchmark	φ-PEF vs body-fitted PEF < 1% nodal	0.63% max, 0.37% L2 ✅
3 Phase-field FNN	81 params, φ ∈ [−1,1]	✅
4 Pre-training	MSE ~1e-7…1e-8 (paper: 3.89e-9)	5.998e-9 ✅
5 Objective	identical→COD=1; mismatch→COD<0.9	✅
6 GA loop	monotone best-fitness; tube→star morphology	✅

Headline benchmark — dual-thrust motor

Quantity	This work	Paper
COD	0.967	0.952
C_D (regression slope)	32.58	31.77

Reconstructed from the two-plateau (10 MPa / 5 MPa) pressure curve. The method reproduces the paper's result and slightly exceeds its reported fit quality.

On that COD. The GA's own scoring at n_nodes=20 reported ~0.89 and appeared to plateau below the paper. That turned out to be a metric artifact, not a worse grain: COD is sensitive to the interpolation node count, because sparse nodes over-weight a small late-web residual at the curve overlap edge. Re-scoring the same saved winner across node counts (diag_nodes.py) gives 0.891 (n=20) → 0.947 (n=80) → 0.960 (n=200) → 0.967 (converged, n≥1600). The converged, density-independent number is the honest one. An earlier "0.918 @ n=200" note in the working log was wrong and has been discarded.

Star benchmark — deferred, deliberately

The paper's second benchmark reconstructs a given star grain (its Fig. 15). That target is not recoverable from the published figure: the readable dimensions underdetermine the grain. phipef/star.py builds the star's phase field analytically, and the inversion scripts (invert_star.py, invert_pen.py, explore*_star.py) demonstrate the contradiction — pinning the burnout web to the figure's 0.595 forces the initial perimeter 50% too high (13.97 vs 9.3) and overshoots peak-web (0.467 vs 0.52). No parameter triple in the figure-readable family satisfies all three features at once.

This is a target-acquisition limit, not a method failure — the same machinery reproduces the peak amplitude to ~3% throughout, and reproduces the dual-thrust case outright. Guessing a textbook star would have silently substituted a different grain. The clean reopen is the exact geometry from the paper's cited source (Jiang's thesis, ref [1]).

---

Quickstart

Requires Python 3.11 (pymoo 0.6.1 ships compiled cp311 wheels).

uv venv --python 3.11 .venv
uv pip install -r requirements.txt

.venv/Scripts/python verify_stage0.py      # Windows
# .venv/bin/python verify_stage0.py        # Linux/macOS

Each gate script prints its pass/fail criterion and writes a stage*.png:

python verify_stage1.py             # reverse-ballistics target curve
python verify_stage2_benchmark.py   # the keystone φ-PEF correctness check
python verify_stage5.py             # burn perimeter + COD objective
python verify_stage6_ga.py          # GA wiring + parallel speedup

Run the full evolution (hours; parallel over cores):

python run_stage6_ga.py

---

Layout

phipef/                 the library
  ballistics.py         M1  reverse internal ballistics
  mesh.py               M2  fixed 1/N fan-sector mesh
  solver.py             M2  φ-PEF Newton solve + α/K-continuation
  benchmark.py          M2  square∩circle validation case
  fnn.py                M3  81-parameter phase-field network
  pretrain.py           M3  L-BFGS fit of the initial weights p0
  objective.py          M4+M5  burn perimeter + COD fitness
  ga.py                 M6  pymoo problem (the full forward chain)
  star.py               analytic star grain (Stage 7 target constructor)

verify_stage*.py        stage gates — the executable specification
run_stage6_ga.py        full GA driver
multi_seed.py           seed-scatter diagnostic (search-limited vs representation-limited)
restart_search.py       restart-based exploration search
nodes_only.py           node-count experiment
diag_nodes.py           the COD metric-artifact proof
explore*_star.py        star geometry sweeps
invert_star.py          burnout inversion, bluntness DOF
invert_pen.py           burnout inversion, penetration DOF
experiment_2d.py        terminal-K choice for the benchmark
render_pdf.py           high-DPI figure crops from the source paper

day1progress.md         Stages 0–4 log
day2progress.md         Stages 5–7 log, incl. the star deferral proof
day2review.md           mid-project review
day2todo.md             planning notes
CLAUDE.md               the binding specification this was built against
p0.npy                  frozen pre-trained network weights

---

Deviations from the paper

Each of these was needed to make the published method actually run.

1. α-continuation with a backtracking line search (not in the paper). A cold Newton solve at the target α = 0.15·δl overshoots catastrophically (δW ~ 1/α). Ramping a global α scale 100 → 1 at K=1 first, then running K-continuation, makes the solve robust across ~10⁴ GA evaluations.
2. Terminal K = 1e5, where the paper illustrates K ~ 1e3. The implicit- interface offset scales as ~1/√K; the <1% benchmark gate needs 1e5, at which point the mesh-error floor (0.63%) dominates. The mechanism is unchanged.
3. Divergence-form discretisation of the governing equation rather than the literal non-divergence form. The two differ by a term that vanishes in solid and is annihilated in gas. Divergence form is symmetric and well-posed — and this choice is gated by the square∩circle benchmark against an independent body-fitted reference.
4. Structured polar MeshQuad instead of gmsh for the fan sector. The φ-PEF domain is a plain annular wedge, so a structured grid gives exact element-count control and clean boundary tags in ~20 lines. gmsh is still used for the benchmark's body-fitted lens — the case it is actually for.
5. L-BFGS + weight decay instead of MATLAB's Bayesian-regularisation trainbr for pre-training. Non-critical: those weights are only the GA's initial guess.
6. tanh radial tube profile for the pre-training target, not a linear ramp. A linear ramp saturates to exactly ±1, which a tanh-output network cannot represent (MSE floors ~3e-6); the tanh profile reaches 6e-9.

---

Source paper

The paper is open access (CC BY 4.0) but is not redistributed here. Download it from doi:10.3390/aerospace9100552 and place it at NN_PAPER.pdf if you want to regenerate figure crops with render_pdf.py.

License

Apache-2.0. See NOTICE for attribution of the method's original authors.