Atomizer/docs/plans/atomizer_fast_solver_technologies.md

# Atomizer Future Development: Fast Solver Technologies

## Document Information
- **Created:** December 2024
- **Updated:** December 28, 2025 (L-BFGS implementation)
- **Purpose:** Technical reference for integrating fast solver technologies into Atomizer
- **Status:** Partially Implemented (L-BFGS complete, ROM/meshfree planned)
- **Priority:** Phase 2+ (after core parametric optimization is stable)

---

## Executive Summary

This document captures research on alternative solver technologies that could dramatically accelerate Atomizer's optimization capabilities. Two main approaches are explored:

1. **Meshfree/Immersed FEA** - Eliminates meshing bottleneck entirely (Intact Solutions, Gridap.jl)
2. **Reduced Order Models (ROM/POD)** - Accelerates existing Nastran by 1000x for parametric problems

**Recommendation:**
- **Short-term:** Implement ROM/POD for parametric optimization (2-4 weeks, 1000x speedup)
- **Medium-term:** Prototype with Gridap.jl for topology optimization capability (2-3 months)
- **Long-term:** Evaluate Intact licensing or productionize custom meshfree solver

---

## Part 1: Meshfree / Immersed FEA Methods

### 1.1 The Problem with Traditional FEA

Traditional FEA requires **body-fitted meshes** that conform to geometry:

```
Traditional FEA Pipeline:
                                                
CAD Geometry → Clean/Heal → Defeature → Mesh Generation → Solve → Post-process
                                              ↑
                                         BOTTLENECK
                                         - 80% of time
                                         - Breaks on complex geometry
                                         - Must remesh for topology opt
```

**Pain points:**
- Meshing complex geometry (lattices, AM parts) often fails
- Small features require fine meshes → huge DOF count
- CAD errors (gaps, overlaps) break meshing
- Topology optimization requires remeshing every iteration
- Automation breaks when meshing fails

### 1.2 The Meshfree Solution: Immersed Methods

**Core idea:** Don't mesh the geometry. Mesh the SPACE.

```
Immersed FEA Pipeline:

CAD Geometry → Convert to Implicit → Immerse in Grid → Solve → Post-process
              (any format works)     (regular grid)    (always works)
```

**How it works:**

```
Step 1: Background Grid              Step 2: Classify Cells
                                    
╔══╤══╤══╤══╤══╗                    ╔══╤══╤══╤══╤══╗
║  │  │  │  │  ║                    ║○ │○ │○ │○ │○ ║   ○ = Outside
╠══╪══╪══╪══╪══╣                    ╠══╪══╪══╪══╪══╣   ● = Inside
║  │  │  │  │  ║    +  Geometry  →  ║○ │◐ │● │◐ │○ ║   ◐ = Cut by boundary
╠══╪══╪══╪══╪══╣       (implicit)   ╠══╪══╪══╪══╪══╣
║  │  │  │  │  ║                    ║○ │◐ │● │◐ │○ ║
╚══╧══╧══╧══╧══╝                    ╚══╧══╧══╧══╧══╝

Step 3: Special Integration for Cut Cells

    Full cell: Standard Gauss quadrature
    
    Cut cell:  Moment-fitting quadrature
               (custom weights that integrate correctly
                over the irregular "inside" portion)
```

### 1.3 Implicit Geometry Representation

**Signed Distance Function (SDF):**
- φ(x) < 0 → inside geometry
- φ(x) > 0 → outside geometry
- φ(x) = 0 → on boundary (the surface)

```
Example: Plate with hole

        φ > 0 (outside)
           ↓
 ─────────────────────
│    ·    │    ·      │
│  -0.3   │  +0.2     │     Numbers show signed distance
│─────────┼───────────│  ←  φ = 0 (surface)
│  -0.5   │  -0.1     │
│    ·    │    ·      │
 ─────────────────────
           ↑
        φ < 0 (inside)
```

**Why implicit is powerful:**
- Boolean operations are trivial: union = min(φ_A, φ_B), intersection = max(φ_A, φ_B)
- Works with ANY input: CAD, STL, CT scans, procedural geometry
- Resolution-independent
- VDB format stores efficiently (sparse, only near surface)

### 1.4 Key Technical Challenges

#### Challenge 1: Cut Cell Integration

When boundary cuts through a cell, standard Gauss quadrature doesn't work.

**Solution: Moment Fitting**

Find quadrature weights w_i such that polynomial integrals are exact:

```
∫ xᵃ yᵇ zᶜ dV = Σ wᵢ · xᵢᵃ yᵢᵇ zᵢᶜ   for all a+b+c ≤ p
Ω_cut           i
```

#### Challenge 2: Boundary Condition Enforcement

Grid doesn't have nodes ON the boundary.

**Solution: Nitsche's Method (Weak Enforcement)**

Add penalty term to variational form:
```
a(u,v) = ∫ ε(u):σ(v) dΩ 
       + ∫ β/h (u-g)·(v) dΓ      ← penalty term
       - ∫ σ(u)·n·v dΓ           ← consistency
       - ∫ σ(v)·n·u dΓ           ← symmetry
```

#### Challenge 3: Small Cut Cells (Ill-Conditioning)

Tiny slivers cause numerical instability.

**Solution: Cell Aggregation (AgFEM)**

Merge small cuts with neighboring full cells:
```
Before:                    After (aggregated):
┌───┬───┐                  ┌───┬───┐
│   │░░░│ ← tiny sliver    │   │▓▓▓│ ← merged with neighbor
└───┴───┘                  └───┴───┘
```

### 1.5 Intact Solutions: Commercial Implementation

**Company:** Intact Solutions, Inc. (Madison, WI)
**Technology:** Immersed Method of Moments (IMM)
**Funding:** DARPA, NSF, NIST
**History:** 10+ years (started with Scan&Solve for Rhino)

**Products:**
| Product | Platform | Use Case |
|---------|----------|----------|
| Intact.Simulation for Grasshopper | Rhino | Design exploration |
| Intact.Simulation for nTop | nTopology | Lattice/AM parts |
| Intact.Simulation for Synera | Synera | Automotive |
| Intact.Simulation for Automation | Headless/API | Batch processing |
| PyIntact | Python API | Custom workflows |
| LevelOpt | All platforms | Topology optimization |

**Capabilities:**
- Physics: Linear static, thermal, modal
- Geometry: STL, VDB, BREP, CT scans, G-code, hybrid
- Materials: Isotropic, orthotropic
- Multi-load case support
- Assembly handling (bonded)

**Performance (from their case study):**
```
Steering knuckle topology optimization:
- 5,600 simulations
- 1,400 design iterations
- 4 load cases
- ~200,000 DOF each
- Total time: 45 hours
- Per simulation: ~30 seconds
- No meshing, no manual intervention
```

**Limitations (inferred):**
- No nonlinear mechanics (yet?)
- No dynamic response (frequency response, random vib)
- Not certified for aerospace (Nastran still required for certification)
- Proprietary (black box)

### 1.6 Open-Source Alternative: Gridap.jl

**Language:** Julia
**Repository:** github.com/gridap/Gridap.jl, github.com/gridap/GridapEmbedded.jl
**License:** MIT
**Maturity:** Research+ (active development, used in publications)

**Why Gridap:**
- Native AgFEM support (handles cut cells properly)
- Level set geometry built-in
- Clean mathematical API
- Good performance (Julia compiles to native code)
- Active community

**Installation:**
```julia
using Pkg
Pkg.add("Gridap")
Pkg.add("GridapEmbedded")
```

**Basic Example:**
```julia
using Gridap
using GridapEmbedded

# 1. Background mesh
domain = (-1, 1, -1, 1, -1, 1)
partition = (20, 20, 20)
bgmodel = CartesianDiscreteModel(domain, partition)

# 2. Implicit geometry (level set)
R = 0.7
geo = sphere(R)  # Built-in primitive

# 3. Cut the mesh
cutgeo = cut(bgmodel, geo)
Ω = Triangulation(cutgeo, PHYSICAL)  # Inside domain
Γ = EmbeddedBoundary(cutgeo)          # Boundary

# 4. FE space with aggregation
reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, 2)
V = TestFESpace(Ω, reffe, conformity=:H1)
U = TrialFESpace(V)

# 5. Weak form (linear elasticity)
const E, ν = 210e9, 0.3
const λ = E*ν / ((1+ν)*(1-2ν))
const μ = E / (2*(1+ν))

ε(u) = 0.5 * (∇(u) + transpose(∇(u)))
σ(u) = λ * tr(ε(u)) * one(ε(u)) + 2*μ * ε(u)

a(u,v) = ∫(ε(v) ⊙ σ(u)) * dΩ
l(v) = ∫(v ⋅ VectorValue(0,0,-1e6)) * dΓ

# 6. Solve
op = AffineFEOperator(a, l, U, V)
uh = solve(op)

# 7. Export
writevtk(Ω, "results", cellfields=["u"=>uh])
```

**Parameterized Solver for Optimization:**
```julia
function solve_bracket(hole_radius, thickness, load)
    # Background mesh
    domain = (0, 1, 0, 0.5, 0, thickness)
    partition = (40, 20, 4)
    bgmodel = CartesianDiscreteModel(domain, partition)
    
    # Parameterized geometry
    function φ(x)
        d_box = max(max(-x[1], x[1]-1), max(-x[2], x[2]-0.5), max(-x[3], x[3]-thickness))
        d_hole = sqrt((x[1]-0.5)^2 + (x[2]-0.25)^2) - hole_radius
        return max(d_box, -d_hole)
    end
    
    geo = AnalyticalGeometry(φ)
    cutgeo = cut(bgmodel, geo)
    Ω = Triangulation(cutgeo, PHYSICAL)
    
    # ... FE setup and solve ...
    
    return (max_stress=σ_max, max_disp=u_max, mass=m)
end

# Use in optimization loop
for r in 0.05:0.01:0.2
    result = solve_bracket(r, 0.1, 1e6)
    println("r=$r, σ_max=$(result.max_stress/1e6) MPa")
end
```

### 1.7 Other Open-Source Options

| Library | Language | Immersed Support | Notes |
|---------|----------|------------------|-------|
| **FEniCSx** | Python/C++ | Experimental | Most mature FEM, but CutFEM needs work |
| **MFEM** | C++ | Manual | Very fast, DOE-backed, shifted boundary |
| **deal.II** | C++ | Partial | German powerhouse, some cut cell support |
| **scikit-fem** | Python | DIY | Simple, hackable, but build from scratch |

**Recommendation:** Start with Gridap.jl for prototyping. It has the best out-of-box immersed support.

---

## Part 2: Reduced Order Models (ROM/POD)

### 2.1 Concept

**Problem:** FEA is slow because we solve huge systems (100k+ DOFs).

**Insight:** For parametric variations, solutions live in a low-dimensional subspace.

**Solution:** Find that subspace, solve there instead.

```
Full Order Model:               Reduced Order Model:

K u = F                         K_r u_r = F_r
                                
K: 100,000 × 100,000           K_r: 20 × 20
u: 100,000 unknowns            u_r: 20 unknowns

Solve time: minutes             Solve time: milliseconds
```

### 2.2 How POD Works

#### Step 1: Collect Snapshots (Offline)

Run N full FEA simulations with different parameters:

```python
snapshots = []
for params in sample_parameter_space(n=100):
    result = nastran_solve(geometry(params), loads, bcs)
    snapshots.append(result.displacement_field)

U = np.array(snapshots)  # Shape: (n_samples, n_dofs)
```

#### Step 2: SVD to Find Basis (Offline)

```python
from scipy.linalg import svd

U_centered = U - U.mean(axis=0)
Phi, S, Vt = svd(U_centered, full_matrices=False)

# Keep top k modes
k = 20  # Usually 10-50 captures 99%+ variance
basis = Vt[:k, :]  # Shape: (k, n_dofs)
```

The singular values decay rapidly - most "information" is in first few modes:

```
Singular Value Spectrum (typical):

    S │
  100 │██
      │██
   50 │██░░
      │██░░░░
   10 │██░░░░░░░░░░░░░░░░░░░░░░
      └─────────────────────────
           1   5   10  15  20  ...
                Mode index
```

#### Step 3: Project and Solve (Online - FAST)

```python
def rom_solve(K, F, basis):
    """
    Solve in reduced space.
    
    Args:
        K: (n_dofs, n_dofs) stiffness matrix
        F: (n_dofs,) force vector
        basis: (k, n_dofs) POD modes
    
    Returns:
        u: (n_dofs,) approximate displacement
    """
    # Project to reduced space
    K_r = basis @ K @ basis.T  # (k, k)
    F_r = basis @ F            # (k,)
    
    # Solve tiny system
    u_r = np.linalg.solve(K_r, F_r)  # Instant!
    
    # Project back
    u = basis.T @ u_r  # (n_dofs,)
    
    return u
```

### 2.3 Complete ROM Implementation for Atomizer

```python
# atomizer/solvers/rom_solver.py

import numpy as np
from scipy.linalg import svd
from scipy.sparse import issparse
from typing import List, Dict, Callable, Tuple
import pickle

class ROMSolver:
    """
    Reduced Order Model solver for fast parametric optimization.
    
    Usage:
        # Training (offline, one-time)
        rom = ROMSolver()
        rom.train(parameter_samples, nastran_solver, n_modes=20)
        rom.save("bracket_rom.pkl")
        
        # Inference (online, fast)
        rom = ROMSolver.load("bracket_rom.pkl")
        u, error_est = rom.solve(new_params)
    """
    
    def __init__(self, n_modes: int = 20):
        self.n_modes = n_modes
        self.basis = None           # (n_modes, n_dofs)
        self.singular_values = None # (n_modes,)
        self.mean_solution = None   # (n_dofs,)
        self.parameter_samples = None
        self.K_func = None          # Function to build stiffness
        self.F_func = None          # Function to build force
        
    def train(self, 
              parameter_samples: List[Dict],
              solver_func: Callable,
              K_func: Callable = None,
              F_func: Callable = None,
              verbose: bool = True):
        """
        Train ROM from full-order snapshots.
        
        Args:
            parameter_samples: List of parameter dictionaries
            solver_func: Function(params) -> displacement field (n_dofs,)
            K_func: Function(params) -> stiffness matrix (for online solve)
            F_func: Function(params) -> force vector (for online solve)
        """
        self.parameter_samples = parameter_samples
        self.K_func = K_func
        self.F_func = F_func
        
        # Collect snapshots
        if verbose:
            print(f"Collecting {len(parameter_samples)} snapshots...")
        
        snapshots = []
        for i, params in enumerate(parameter_samples):
            u = solver_func(params)
            snapshots.append(u)
            if verbose and (i+1) % 10 == 0:
                print(f"  Snapshot {i+1}/{len(parameter_samples)}")
        
        U = np.array(snapshots)  # (n_samples, n_dofs)
        
        # Center the data
        self.mean_solution = U.mean(axis=0)
        U_centered = U - self.mean_solution
        
        # SVD
        if verbose:
            print("Computing SVD...")
        
        _, S, Vt = svd(U_centered, full_matrices=False)
        
        # Keep top modes
        self.basis = Vt[:self.n_modes, :]
        self.singular_values = S[:self.n_modes]
        
        # Report energy captured
        total_energy = (S**2).sum()
        captured_energy = (self.singular_values**2).sum()
        if verbose:
            print(f"Captured {captured_energy/total_energy*100:.2f}% of variance "
                  f"with {self.n_modes} modes")
        
        return self
    
    def solve(self, params: Dict) -> Tuple[np.ndarray, float]:
        """
        Fast solve using ROM.
        
        Args:
            params: Parameter dictionary
            
        Returns:
            u: Approximate displacement field
            error_estimate: Estimated relative error
        """
        if self.K_func is None or self.F_func is None:
            raise ValueError("K_func and F_func required for online solve")
        
        # Build matrices for this parameter set
        K = self.K_func(params)
        F = self.F_func(params)
        
        # Convert sparse to dense if needed (for small reduced system)
        if issparse(K):
            # Project while sparse
            K_r = self.basis @ K @ self.basis.T
        else:
            K_r = self.basis @ K @ self.basis.T
        
        F_r = self.basis @ F
        
        # Solve reduced system
        u_r = np.linalg.solve(K_r, F_r)
        
        # Project back
        u = self.basis.T @ u_r + self.mean_solution
        
        # Error estimate (residual-based)
        residual = F - K @ u
        error_estimate = np.linalg.norm(residual) / np.linalg.norm(F)
        
        return u, error_estimate
    
    def solve_interpolated(self, params: Dict) -> np.ndarray:
        """
        Even faster: interpolate in coefficient space.
        Requires parameter_samples to be stored with coefficients.
        """
        # Find nearest neighbors in parameter space
        # Interpolate their ROM coefficients
        # This avoids even building K and F!
        raise NotImplementedError("TODO: Implement RBF interpolation")
    
    def estimate_error(self, u_rom: np.ndarray, K: np.ndarray, F: np.ndarray) -> float:
        """Residual-based error estimate."""
        residual = F - K @ u_rom
        return np.linalg.norm(residual) / np.linalg.norm(F)
    
    def should_fallback_to_full(self, error_estimate: float, threshold: float = 0.05) -> bool:
        """Check if ROM error is too high, suggesting full solve needed."""
        return error_estimate > threshold
    
    def save(self, filepath: str):
        """Save trained ROM to file."""
        data = {
            'n_modes': self.n_modes,
            'basis': self.basis,
            'singular_values': self.singular_values,
            'mean_solution': self.mean_solution,
        }
        with open(filepath, 'wb') as f:
            pickle.dump(data, f)
    
    @classmethod
    def load(cls, filepath: str) -> 'ROMSolver':
        """Load trained ROM from file."""
        with open(filepath, 'rb') as f:
            data = pickle.load(f)
        
        rom = cls(n_modes=data['n_modes'])
        rom.basis = data['basis']
        rom.singular_values = data['singular_values']
        rom.mean_solution = data['mean_solution']
        
        return rom


class AdaptiveROMSolver(ROMSolver):
    """
    ROM with automatic fallback to full solver when error is high.
    """
    
    def __init__(self, n_modes: int = 20, error_threshold: float = 0.05):
        super().__init__(n_modes)
        self.error_threshold = error_threshold
        self.full_solver = None
        
    def set_full_solver(self, solver_func: Callable):
        """Set the full-order solver for fallback."""
        self.full_solver = solver_func
        
    def solve_adaptive(self, params: Dict) -> Tuple[np.ndarray, float, bool]:
        """
        Solve with automatic fallback.
        
        Returns:
            u: Displacement field
            error: Error estimate (or 0 if full solve)
            used_full: Whether full solver was used
        """
        # Try ROM first
        u_rom, error = self.solve(params)
        
        if error < self.error_threshold:
            return u_rom, error, False
        
        # Fallback to full
        if self.full_solver is not None:
            u_full = self.full_solver(params)
            return u_full, 0.0, True
        else:
            # Return ROM result with warning
            print(f"Warning: ROM error {error:.3f} exceeds threshold, "
                  f"but no full solver available")
            return u_rom, error, False
```

### 2.4 Integration with Atomizer Optimization Loop

```python
# atomizer/optimization/rom_optimization_runner.py

from atomizer.solvers.rom_solver import AdaptiveROMSolver
from atomizer.solvers.nastran_solver import NastranSolver
import optuna

class ROMOptimizationRunner:
    """
    Optimization using ROM for speed, Nastran for validation.
    """
    
    def __init__(self, problem_config: dict):
        self.config = problem_config
        self.rom = None
        self.nastran = NastranSolver()
        
    def train_rom(self, n_samples: int = 100, n_modes: int = 20):
        """Train ROM on problem-specific data."""
        
        # Sample parameter space
        samples = self._latin_hypercube_sample(n_samples)
        
        # Run Nastran for each sample
        def nastran_solve(params):
            result = self.nastran.solve(
                geometry=self._build_geometry(params),
                loads=self.config['loads'],
                bcs=self.config['bcs']
            )
            return result['displacement_field']
        
        # Train ROM
        self.rom = AdaptiveROMSolver(n_modes=n_modes)
        self.rom.train(samples, nastran_solve)
        self.rom.set_full_solver(nastran_solve)
        
        return self
    
    def optimize(self, n_trials: int = 500):
        """Run optimization using ROM."""
        
        def objective(trial):
            # Sample parameters
            params = {}
            for name, bounds in self.config['parameters'].items():
                params[name] = trial.suggest_float(name, bounds[0], bounds[1])
            
            # Fast ROM solve
            u, error, used_full = self.rom.solve_adaptive(params)
            
            # Compute objective (e.g., max stress)
            stress = self._compute_stress(u)
            max_stress = stress.max()
            
            # Log
            trial.set_user_attr('rom_error', error)
            trial.set_user_attr('used_full_solver', used_full)
            
            return max_stress
        
        # Run optimization
        study = optuna.create_study(direction='minimize')
        study.optimize(objective, n_trials=n_trials)
        
        # Validate best designs with full Nastran
        return self._validate_top_designs(study, n_top=10)
    
    def _validate_top_designs(self, study, n_top: int = 10):
        """Run full Nastran on top candidates."""
        
        top_trials = sorted(study.trials, key=lambda t: t.value)[:n_top]
        
        validated = []
        for trial in top_trials:
            result = self.nastran.solve(
                geometry=self._build_geometry(trial.params),
                loads=self.config['loads'],
                bcs=self.config['bcs']
            )
            validated.append({
                'params': trial.params,
                'rom_prediction': trial.value,
                'nastran_result': result['max_stress'],
                'error': abs(trial.value - result['max_stress']) / result['max_stress']
            })
        
        return validated
```

### 2.5 When ROM Works Well vs. Poorly

```
WORKS GREAT:                              WORKS POORLY:

✓ Parametric variations                   ✗ Topology changes
  (thickness, radii, positions)             (holes appear/disappear)

✓ Smooth parameter dependence             ✗ Discontinuous behavior
                                            (contact, buckling)

✓ Similar geometry family                 ✗ Wildly different geometries

✓ Linear problems                         ✗ Highly nonlinear

✓ Fixed mesh topology                     ✗ Mesh changes between solves
```

**For Atomizer's current parametric optimization: ROM is ideal.**

---

## Part 3: Comparison of Approaches

### 3.1 Speed Comparison

| Method | Typical Solve Time | Speedup vs Nastran | Training Cost |
|--------|-------------------|-------------------|---------------|
| **Nastran (baseline)** | 10 min | 1x | None |
| **Intact IMM** | 30 sec | 20x | None |
| **ROM/POD** | 10 ms | 60,000x | 100 Nastran runs |
| **Neural Surrogate (MLP)** | 10 ms | 60,000x | 50+ FEA runs |
| **L-BFGS on Surrogate** | <1 sec (20-50 evals) | 100-1000x vs TPE | Trained surrogate |

### 3.1.1 L-BFGS Gradient Optimization (IMPLEMENTED)

**Status:** ✓ Implemented in `optimization_engine/gradient_optimizer.py`

L-BFGS exploits the differentiability of trained MLP surrogates for ultra-fast local refinement:

```python
from optimization_engine.gradient_optimizer import run_lbfgs_polish

# Polish top candidates from surrogate search
results = run_lbfgs_polish(study_dir, n_starts=20, n_iterations=100)
```

**Performance:**
| Method | Evaluations to Converge | Time |
|--------|------------------------|------|
| TPE (surrogate) | 200-500 | ~30 min |
| CMA-ES (surrogate) | 100-300 | ~15 min |
| **L-BFGS (surrogate)** | **20-50** | **<1 sec** |

**When to use:** After Turbo mode identifies promising regions, before FEA validation.

**Full documentation:** `docs/protocols/system/SYS_14_NEURAL_ACCELERATION.md` (Section: L-BFGS Gradient Optimizer)

### 3.2 Capability Comparison

| Capability | Nastran | Intact | ROM | Neural |
|------------|---------|--------|-----|--------|
| Parametric opt | ✓ (slow) | ✓ | ✓✓ | ✓✓ |
| Topology opt | ✗ | ✓✓ | ✗ | ? |
| New geometry | ✓ | ✓ | ✗ (retrain) | ✗ (retrain) |
| Certification | ✓✓ | ✗ | ✗ | ✗ |
| Nonlinear | ✓ | ? | Limited | Limited |
| Dynamic | ✓ | ? | Limited | Limited |

### 3.3 Implementation Effort

| Approach | Effort | Risk | Payoff |
|----------|--------|------|--------|
| ROM/POD | 2-4 weeks | Low | Immediate 1000x speedup |
| Gridap prototype | 2-3 months | Medium | Topology opt capability |
| Intact license | $$ + 2-4 weeks | Low | Production meshfree |
| AtomizerField | Ongoing | Medium | Field predictions |

---

## Part 4: Recommended Atomizer Architecture

### 4.1 Three-Tier Fidelity Stack

```
┌─────────────────────────────────────────────────────────────────────────┐
│                      ATOMIZER SOLVER STACK                               │
├─────────────────────────────────────────────────────────────────────────┤
│                                                                          │
│  TIER 1: AtomizerField / ROM  (Screening)                               │
│  ──────────────────────────────────────────                              │
│  Speed:     ~10 ms per evaluation                                        │
│  Accuracy:  ~90-95%                                                      │
│  Use case:  Design space exploration, optimization inner loop            │
│  Output:    Promising candidates                                         │
│                                                                          │
│  TIER 2: Intact IMM / Gridap  (Refinement)  [FUTURE]                    │
│  ──────────────────────────────────────────────────                      │
│  Speed:     ~30 sec per evaluation                                       │
│  Accuracy:  ~98%                                                         │
│  Use case:  Topology optimization, complex geometry                      │
│  Output:    Optimized designs                                            │
│                                                                          │
│  TIER 3: NX Nastran  (Validation)                                       │
│  ────────────────────────────────                                        │
│  Speed:     ~10 min per evaluation                                       │
│  Accuracy:  100% (reference)                                             │
│  Use case:  Final validation, certification                              │
│  Output:    Certified results                                            │
│                                                                          │
└─────────────────────────────────────────────────────────────────────────┘

Workflow:

  Design Space ──▶ Tier 1 ──▶ Candidates ──▶ Tier 2 ──▶ Best ──▶ Tier 3
  (1M points)      (1 sec)    (1000 pts)    (8 hrs)    (10)    (2 hrs)
                                               │
                                               └── Skip if parametric only

**Current Implemented Workflow (Turbo + L-BFGS):**

```
┌─────────────────────────────────────────────────────────────────────────┐
│ 1. FEA Exploration    │  50-100 Nastran trials                          │
│    └─► Train MLP surrogate on FEA data                                  │
├─────────────────────────────────────────────────────────────────────────┤
│ 2. Turbo Mode         │  5000+ NN predictions (~10ms each)              │
│    └─► Find promising design regions                                    │
├─────────────────────────────────────────────────────────────────────────┤
│ 3. L-BFGS Polish      │  20 multi-start runs (<1 sec total)   [NEW]    │
│    └─► Precise local optima via gradient descent                        │
├─────────────────────────────────────────────────────────────────────────┤
│ 4. FEA Validation     │  Top 3-5 candidates verified with Nastran       │
│    └─► Final certified results                                          │
└─────────────────────────────────────────────────────────────────────────┘
```

### 4.2 Implementation Phases

```
PHASE 1: NOW (Q4 2024 - Q1 2025)
────────────────────────────────
[✓] Core optimization engine (Optuna/CMA-ES/NSGA-II)
[✓] NX/Nastran integration
[✓] Dashboard
[✓] MLP Surrogate + Turbo Mode
[✓] L-BFGS Gradient Optimizer  ← IMPLEMENTED 2025-12-28
[ ] ROM/POD implementation
[ ] AtomizerField training pipeline

PHASE 2: MEDIUM-TERM (Q2-Q3 2025)
────────────────────────────────────
[ ] Gridap.jl prototype
[ ] Evaluate Intact licensing
[ ] Topology optimization capability
[ ] Multi-fidelity optimization

PHASE 3: FUTURE (Q4 2025+)
──────────────────────────
[ ] Production meshfree solver (Intact or custom)
[ ] Full AtomizerField deployment
[ ] Cloud/HPC scaling
[ ] Certification workflow integration
```

---

## Part 5: Key Resources

### 5.1 Intact Solutions
- Website: https://intact-solutions.com
- PyIntact docs: https://intact-solutions.com/intact_docs/PyIntact.html
- Contact for pricing: info@intact-solutions.com

### 5.2 Gridap.jl
- Main repo: https://github.com/gridap/Gridap.jl
- Embedded: https://github.com/gridap/GridapEmbedded.jl
- Tutorials: https://gridap.github.io/Tutorials/stable/
- Paper: "Gridap: An extensible Finite Element toolbox in Julia" (2020)

### 5.3 ROM/POD References
- Benner et al., "Model Reduction and Approximation" (2017)
- Quarteroni et al., "Reduced Basis Methods for PDEs" (2016)
- pyMOR library: https://pymor.org (Python model order reduction)

### 5.4 Immersed FEM Theory
- Burman et al., "CutFEM: Discretizing geometry and PDEs" (2015)
- Schillinger & Ruess, "The Finite Cell Method: A Review" (2015)
- Main & Scovazzi, "The Shifted Boundary Method" (2018)

---

## Appendix A: Quick Reference - Immersed FEM Math

### Weak Form with Nitsche BCs

Standard elasticity:
```
Find u ∈ V such that:
a(u,v) = l(v)  ∀v ∈ V

a(u,v) = ∫_Ω ε(u):C:ε(v) dΩ
l(v) = ∫_Ω f·v dΩ + ∫_Γ_N t·v dΓ
```

With Nitsche (weak Dirichlet on Γ_D):
```
a(u,v) = ∫_Ω ε(u):C:ε(v) dΩ
       - ∫_Γ_D (σ(u)·n)·v dΓ       # Consistency
       - ∫_Γ_D (σ(v)·n)·u dΓ       # Symmetry  
       + ∫_Γ_D (β/h) u·v dΓ        # Penalty

l(v) = ∫_Ω f·v dΩ + ∫_Γ_N t·v dΓ
     - ∫_Γ_D (σ(v)·n)·g dΓ         # Consistency (RHS)
     + ∫_Γ_D (β/h) g·v dΓ          # Penalty (RHS)

where β ~ E (penalty parameter), h = cell size, g = prescribed displacement
```

### Ghost Penalty Stabilization

For small cut cells, add:
```
j(u,v) = Σ_F γ h² ∫_F [∂ⁿu/∂nⁿ]·[∂ⁿv/∂nⁿ] dF

where F = interior faces, [·] = jump across face, γ = stabilization parameter
```

---

## Appendix B: Julia/Gridap Cheat Sheet

```julia
# Installation
using Pkg
Pkg.add(["Gridap", "GridapEmbedded", "GridapGmsh"])

# Basic imports
using Gridap
using GridapEmbedded
using GridapGmsh  # For external meshes

# Geometry primitives
sphere(r)                    # Sphere of radius r
disk(r)                      # 2D disk
cylinder(r, h)               # Cylinder
box(corner1, corner2)        # Box

# Boolean operations
union(geo1, geo2)
intersect(geo1, geo2)  
setdiff(geo1, geo2)          # geo1 - geo2

# Custom level set
geo = AnalyticalGeometry(x -> my_sdf(x))

# Cutting
cutgeo = cut(bgmodel, geo)
Ω = Triangulation(cutgeo, PHYSICAL)
Γ = EmbeddedBoundary(cutgeo)

# FE spaces
reffe = ReferenceFE(lagrangian, Float64, 2)           # Scalar, quadratic
reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, 2)  # 3D vector

# Measures
dΩ = Measure(Ω, 4)           # Quadrature degree 4
dΓ = Measure(Γ, 4)

# Weak form
a(u,v) = ∫(∇(u)⋅∇(v)) * dΩ
l(v) = ∫(f*v) * dΩ

# Solve
op = AffineFEOperator(a, l, U, V)
uh = solve(op)

# Export
writevtk(Ω, "output", cellfields=["u"=>uh])
```

---

## Appendix C: ROM/POD Implementation Checklist

```
□ Data Collection
  □ Define parameter space bounds
  □ Latin hypercube sampling (or similar)
  □ Run Nastran for each sample
  □ Store displacement fields consistently (same mesh!)

□ Basis Construction
  □ Center snapshot matrix
  □ SVD decomposition
  □ Select number of modes (energy criterion: 99%+)
  □ Store basis and mean

□ Online Solver
  □ Stiffness matrix assembly for new params
  □ Force vector assembly for new params
  □ Project to reduced space
  □ Solve reduced system
  □ Project back to full space

□ Error Estimation
  □ Residual-based estimate
  □ Fallback threshold selection
  □ Adaptive solver wrapper

□ Integration
  □ Save/load trained ROM
  □ Plug into optimization loop
  □ Validation against full solver
```

---

*Document maintained by: Atomizer Development Team*
*Last updated: December 2024*