Project Context: φ-PEF Evolutionary-Neural-Network Grain Reverse Design

Paste this file as CLAUDE.md (or the initial context message) into Claude Code. It is the binding specification for this project. Treat it as source of truth alongside the paper NN_PAPER.pdf (Li et al., Aerospace 2022, 9, 552).

0. Mission

Programmatically reconstruct, in Python, the reverse-design system that maps a target chamber-pressure curve p_c(t) to a 2D solid-propellant grain shape. The system replaces semi-empirical parameter optimization with shape optimization driven by an evolutionary neural network. Four coupled subsystems:

φ-PEF burn-back FEM — nonlinear stationary solve on a fixed unstructured mesh.
Phase-field FNN — a 2→20→1 network (81 params) that emits the phase field φ(x,y).
Genetic algorithm — evolves the 81 FNN parameters to match the target curve.
COD fitness — linear-regression coefficient of determination as the objective.

The deliverable is the minimum code that reproduces the paper's quantitative benchmarks (Section 9), not a general-purpose framework.

1. Operating Principles (binding)

Think before coding. State assumptions; if multiple interpretations exist, surface them — do not pick silently. If something is unclear, stop and ask.
Simplicity first. Minimum code that passes the gate. No speculative features, no abstractions for single-use code, no configurability that wasn't requested. If 200 lines could be 50, rewrite to 50.
Surgical changes. Touch only what the current stage requires. Match existing style. Remove only orphans your own change created; flag pre-existing dead code, don't delete it.
Goal-driven execution. Every stage below has a → verify: gate with a quantitative pass condition. Do not promote a stage until its gate passes. Do not build ahead — produce the minimum to clear the current gate, then stop for confirmation.

2. Non-Dimensionalization Contract (never violate)

All optimization happens in 2D non-dimensional "cyber space":

Outer radius normalized to a standard value R0 (paper uses 1 m). Physical grains map in by scaling; the solver never reasons about physical R or length L directly.
Grain length enters only through the scalar C_T = L·R / A_t. Geometry and C_T are decoupled.
Dimensionless burn perimeter l̂_b = l_b / R; dimensionless web ŵ = w / R.
The optimizer never sees L; it is recovered post-hoc from the designed C_D = L_D/A_{t,D}.

If any module starts handling absolute lengths inside the GA loop, that is a contract violation — stop and reconsider.

3. Pre-Resolved Technical Decisions (do not re-litigate without flagging)

#	Decision	Rationale
D1	Discretize Eq. 11 in divergence form `∫ αr ∇W·∇v` rather than the literal non-divergence `αr∇²W`.	The two differ by `α∇r·∇W`, which vanishes in solid (`r=1`) and is immaterial in gas (`s=0`). Divergence form is symmetric and well-posed. Gate it on the Eq. 15 benchmark. If that fails, revisit.
D2	`r` and `s` are constant during each W-solve.	They depend only on `φ` (fixed per GA individual), not on `W`. The sole nonlinearity is the eikonal `r²
D3	Hand-coded Newton (scikit-fem `ex10` pattern), no JAX.	Honors minimalism. Escape hatch: `NonlinearForm`+`skfem.autodiff` (adds JAX) only if the hand-derived Jacobian proves error-prone — flag before switching.
D4	`K`-continuation is mandatory. Ramp `K: 1 → 10³` (e.g. `[1, 10, 100, 1000]`), Newton-converging at each step, warm-starting `W`.	Eq. 13's `K ≫ 1` makes the system stiff; a cold solve at `K=1000` will not converge. The paper omits this; it is required in practice.
D5	Constraints via penalty-in-objective `FV = 1 − COD + P` (Eq. 29).	Matches the paper exactly. pymoo-native `n_ieq_constr`/`out["G"]` is the cleaner alternative — note it, but default to the penalty for fidelity.

4. Technical Assumptions (A1–A10)

A1 All work is 2D, non-dimensional (see §2).
A2 Physics holds only inside the paper's envelope: uniform parallel-layer burn rate (∇W·∇W=1 RHS), no erosive combustion, equilibrium pressure (Eq. 5), pressure-independent ρ_p, c*. Not re-validated at runtime.
A3 N (geometry units per half-angle) is a fixed integer per run. Optimal N is found by an outer enumeration loop, not by the GA. Each N → a distinct fixed mesh (a 1/N fan sector, inner radius 0.2R, outer radius R).
A4 Eq. 11 is a nonlinear stationary diffusion–eikonal problem → Newton + K-continuation (D3, D4).
A5 Fixed-mesh invariance: mesh, connectivity, and the geometry-independent diffusion skeleton are assembled once and reused across all ~10⁴ GA evaluations. Only nodal r(x,y), s(x,y) change per individual. This is the primary efficiency lever.
A6 BCs: Dirichlet W=0 on the inner-arc initial burning surface; zero-flux ∂W/∂n=0 on flame-retardant and symmetry edges.
A7 φ ∈ [−1,1] enforced by the tanh output layer and by Δp ∈ [−1,1]; assumed sufficient to suppress interior holes (Requirement 2).
A8 The objective is non-differentiable (isocontour extraction on W) → gradient-free GA; pre-training (not backprop) sets the operating point.
A9 Pre-training substitutes MATLAB Bayesian regularization (trainbr) with L-BFGS + weight decay (or a PyTorch L2 prior). p_0 is only the GA's initial guess, so the substitution is non-critical.
A10 l̂_{b,D}(ŵ) = ℓ_contour(W=w) / (N·R); tail-off accuracy is bounded by the linear-interpolation node density (Eq. 24).

5. System Decomposition (M1–M6 + outer N loop)

Data flow within a single fitness evaluation is strictly one-directional:

Target p_c~t ─[M1]→ (l̂_b,T·C_T)~ŵ_T   (target curve f2)
                                              │
Δp ─[M3 FNN]→ φ(nodes) → r,s ─[M2 φ-PEF]→ W ──┤
                                  │           │
                          [M4 contour→ℓ]      │
                                  ▼           ▼
                       l̂_b,D~ŵ_D (f1) ─[M5 COD]→ FV ─[M6 GA]→ updates Δp
                                                              │
                                            outer loop over N ┘

M1 Reverse internal ballistics (Eq. 4–7): p_c(t) → target (l̂_b,T·C_T)~ŵ_T. Pure quadrature; runs once per request.
M2 φ-PEF FEM solver (Eq. 11–14): nonlinear W-solve on the fixed mesh. Bottleneck.
M3 Phase-field FNN (Eq. 16–18): 2→20→1, tanh, 81 params. Forward-pass only.
M4 Burn-perimeter extraction: isocontours of W at equispaced web levels → l̂_b,D(ŵ).
M5 Objective (Eq. 23–29): regression slope C_D, COD, FV = 1 − COD + P.
M6 GA: real-coded over Δp ∈ [−1,1]^81; population 200, 50 generations.

6. Verified Technology Stack (pinned)

Core (required):

Concern	Library	Confirmed idiom
Language	Python 3.11+	—
Linear algebra / quadrature / least-squares	NumPy, SciPy	`scipy.sparse`, `scipy.optimize`
FEM assembly	scikit-fem ≥ 12	`Basis(m, ElementQuad2())`, `@BilinearForm def f(u,v,w):`, `@LinearForm def g(v,w):`, prev solution via `w['prev']` from `basis.interpolate(W)`, `condense(K, b, D=dirichlet_dofs)`, `solve(...)`
Mesh	gmsh + meshio	one-time fan-sector mesh per `N` → meshio → scikit-fem
Contour length	contourpy (or skimage `find_contours`)	deterministic per-level perimeter
FNN forward pass	NumPy (vectorized)	no DL framework in the inner loop
Pre-training (one-time)	SciPy `L-BFGS-B` + weight decay (or PyTorch)	A9
Optimizer	pymoo 0.6.1	`from pymoo.core.problem import ElementwiseProblem`; `from pymoo.algorithms.soo.nonconvex.ga import GA`; `from pymoo.optimize import minimize`; `from pymoo.parallelization.starmap import StarmapParallelization`

Deployment / ops (defer until core gates pass):

Concern	Library
Parallel fitness	joblib / `multiprocessing` (1 node); Ray (multi-node)
Service	FastAPI + Uvicorn
Async jobs	Redis + RQ (or Ray Jobs)
Container / config / tracking	Docker; Hydra; MLflow

Rejected on purpose (do not introduce): FEniCS/dolfinx (native-binary weight), a deep-learning framework inside the GA loop, any abstraction layer over scikit-fem or pymoo.

7. Confirmed API Idiom Anchors

These are correct against the verified versions. They are anchors, not implementations — write the real modules yourself, but do not deviate from these signatures.

M2 — φ-PEF nonlinear solve (scikit-fem v12, divergence form per D1, Newton per D3, D4):

# r, s are CONSTANT fields per individual (D2). Pass φ as a DiscreteField and threshold
# at quadrature points (cleaner than interpolating a discontinuous nodal s).
@BilinearForm                      # Newton tangent (Jacobian)
def tangent(u, v, w):
    return (w['alpha_r'] * dot(grad(u), grad(v))
            + 2.0 * w['s'] * w['r']**2 * dot(grad(w['Wk']), grad(u)) * v)

@LinearForm                        # residual  F(W_k)
def residual(v, w):
    gW = grad(w['Wk'])
    return (w['alpha_r'] * dot(gW, grad(v))
            + w['s'] * (w['r']**2 * dot(gW, gW) - 1.0) * v)
# Newton step:  J = tangent.assemble(basis, Wk=basis.interpolate(W), ...)
#               F = residual.assemble(basis, Wk=basis.interpolate(W), ...)
#               dW = solve(*condense(J, -F, D=inner_arc_dofs));  W += dW
# Wrap in a K-continuation loop (D4), rebuilding r each K, warm-starting W.

M3 — FNN forward pass (NumPy, 81 params):

def phase_field(p, XY):            # p:(81,)  XY:(n_nodes,2)
    W1 = p[:40].reshape(20, 2);  b1 = p[40:60]
    W2 = p[60:80].reshape(1, 20); b2 = p[80:81]
    a1 = np.tanh(XY @ W1.T + b1)
    return np.tanh(a1 @ W2.T + b2).ravel()    # φ ∈ [-1, 1]
# param count check: (2+1)*20 + (20+1)*1 == 81

M5 — COD objective (NumPy least-squares, Eq. 26–28):

C_D = (y1 @ y2) / (y1 @ y1)                                    # Eq. 26
COD = 1 - ((y2 - C_D*y1) @ (y2 - C_D*y1)) / ((y2 - y2.mean()) @ (y2 - y2.mean()))  # Eq. 27
FV  = 1.0 - COD + P                                            # Eq. 28/29 ; GA minimizes

M6 — GA wrapper (pymoo 0.6.1):

class GrainProblem(ElementwiseProblem):
    def __init__(self, ctx, **kw):
        super().__init__(n_var=81, n_obj=1, xl=-1.0, xu=1.0, **kw)
        self.ctx = ctx                          # holds p0, fixed mesh/basis, target f2
    def _evaluate(self, dp, out, *a, **k):
        out["F"] = fitness(self.ctx.p0 + dp, self.ctx)   # p = p0 + Δp  (Eq. 20)

runner  = StarmapParallelization(pool.starmap)
problem = GrainProblem(ctx, elementwise_runner=runner)
res = minimize(problem, GA(pop_size=200), termination=("n_gen", 50), seed=1)

8. Staged Roadmap (each gate is a hard promotion barrier)

Stage 0 — Environment

Pin core stack; build base image. → verify: import skfem, pymoo, contourpy succeeds; a trivial scikit-fem Poisson solve matches its analytic solution to < 1e-6.

Stage 1 — Reverse ballistics (M1)

Implement Eq. 6 quadrature for l̂_b,T·C_T and ŵ_T; resample to the target curve. → verify: dual-thrust p_c(t) (Eq. 34) + Table 5 params reproduce the two-plateau shape of Fig. 25 (≈30 → ≈16 step) within plotting tolerance.

Stage 2 — Fixed-mesh φ-PEF solver (M2) (critical path — budget the most effort here)

Generate the 1/N fan-sector mesh (inner 0.2R, outer R); tag inner-arc / radial-symmetry / outer-arc boundaries. → verify: element count hits target (star case: 1000 ElementQuad2); boundary tag sets are mutually exclusive and exhaustive.
Assemble the geometry-invariant diffusion skeleton once; set α = 0.15·δl (Eq. 10). → verify: skeleton is symmetric, sparse, independent of any φ.
Newton + K-continuation (D4). → verify: residual ‖F‖₂ < 1e-8 in < 15 Newton steps at terminal K; no NaNs at the gas–solid interface.
Validate against pure-PEF mode on the square∩circle benchmark, analytic φ (Eq. 15). → verify: PEF and φ-PEF W fields agree to < 1% nodal error (also validates D1).

Stage 3 — FNN (M3)

Vectorized forward pass over all nodes. → verify: param count == 81; output ∈ [−1,1] for arbitrary input.

Stage 4 — Pre-training (one-time, M3 init)

Sample analytic tube-grain φ_0 at mesh nodes (Eq. 21). → verify: sampled field shows the concentric tube isobands of Fig. 12.
Fit p_0 by minimizing MSE (Eq. 22), L-BFGS + weight decay (A9). → verify: training MSE reaches the paper's order (≈1e-7 to 1e-8). Freeze p_0; reuse as the GA's affine offset (Eq. 20).

Stage 5 — Objective (M4 + M5)

Extract W isocontours at n equispaced web levels → l̂_b,D(ŵ). → verify: a known star geometry reproduces the rising-then-collapsing C_T·l̂_b/R profile of Fig. 16.
Compute C_D, COD, FV with separation/hole and loading-fraction penalties. → verify: identical curves give COD=1, FV=0; a deliberately mismatched curve gives COD < 0.9.

Stage 6 — GA loop (M6)

Wrap M3→M2→M4→M5 as a pymoo ElementwiseProblem; pop 200, 50 gens, parallel. → verify: one fitness call returns finite FV in bounded wall-time; population evaluates in parallel with ~linear speedup to core count.
Full evolution; log best-FV per generation. → verify: monotone-non-increasing best-fitness trajectory; morphology collapses tube → irregular → smooth star by ~G10 (Fig. 18).

Stage 7 — Benchmark validation

Star, N=12, C_T=30. → verify: C_D ≈ 29.65, COD ≈ 0.995, ε ≈ −1.16% (Eq. 31–32).
Dual-thrust, Table 5. → verify: C_D ≈ 31.77, COD ≈ 0.952 (Eq. 35); reconstructed p_c~t shows stable 10 MPa / 5 MPa stages (Fig. 29).

Stage 8 — Deployment (only after Stage 7 passes)

FastAPI endpoint: input (p_c~t, ρ_p, c*, a, n, R, N-range, C_D-range) → enqueue async job. → verify: endpoint returns a job ID immediately; status polling reflects GA generation.
N-enumeration outer loop as parallel jobs; return per-N Pareto set of (loading fraction, FV). → verify: dual-thrust sweep reproduces the N∈{8,10,12,16,20} morphology transition (dog bone → combined dendrite → wagon wheel, Fig. 30).
Containerize; persist phase fields + metrics to MLflow. → verify: a cold-start container reproduces a stored run's C_D under a fixed RNG seed.

9. Acceptance Gates (summary table)

Case	Source	Required output
Square∩circle, analytic φ	Eq. 15 / Fig. 8	PEF and φ-PEF `W` identical (< 1% nodal)
Star, `N=12`	Eq. 30–32	`C_D≈29.65`, `COD≈0.995`, `ε≈−1.16%`
Star evolution	Fig. 18	tube→star by ~G10; monotone best-`FV`
Dual-thrust, `N=12`	Eq. 35 / Fig. 29	`C_D≈31.77`, `COD≈0.952`; 10/5 MPa plateaus
`N`-sweep	Fig. 30	dog-bone / dendrite / wagon-wheel transition

Failure at any gate blocks promotion of the dependent stage. Highest-risk gates: Stage 2 (Newton convergence under K≫1) and Stage 5 (tail-off perimeter accuracy).

10. What NOT To Do

Do not introduce FEniCS/dolfinx, a DL framework in the inner loop, or any wrapper over scikit-fem/pymoo.
Do not re-assemble the mesh or diffusion skeleton inside the GA loop (violates A5).
Do not let L or absolute lengths enter the optimizer (violates §2).
Do not attempt a cold K=1000 solve (violates D4).
Do not build later stages before the current gate passes. Stop at the gate.
Do not silently choose between the decisions in §3 — if evidence contradicts a default, flag it and present the alternative before switching.

FilesExpand file tree

CLAUDE.md

Latest commit

History